Obtención de datos

Se obtuvieron 227 registros de C. odorata a partir de las bases de datos del Inventario Nacional Forestal y Suelos (INFyS) 2004-2009 (^{CONAFOR, 2009}), del libro de “Árboles tropicales de México” (^{Pennington & Sarukhán, 2005}) y de la plataforma Global Biodiversity Information Facility (^{GBIF, 2017}) (Figura 1). Las bases de datos se depuraron dejando un registro por celda de 1 km² (30 segundos de arco), se revisó su correcta posición eliminando registros duplicados, mal georreferenciados, así como los ubicados en zonas urbanas y cuerpos de agua. Estos procedimientos se realizaron en la plataforma Niche ToolBox de la Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) (^{Osorio-Olvera, Vijay, Narayani, Soberón, & Falconi, 2016}) y en el programa ArcMap versión 10.3 (^{Environmental Scientific Research Institute [ESRI], 2014}). Posterior al proceso de depuración, 216 registros válidos se obtuvieron para la calibración de los modelos. Para una mejor evaluación de los datos, se generaron 50 modelos con diferentes criterios de modelación (Cuadro 1).

Variables bioclimáticas actuales y futuras

Se descargaron las 19 variables bioclimáticas (http://visualcrow.com/descarga-mexico.html) generadas para México del periodo 1910-2009 (^{Cuervo-Robayo et al., 2014}). Para la evaluación de la distribución potencial a futuro (2030), se utilizaron las variables bioclimáticas de los modelos de circulación general GFDL_CM3 y MIROC_ESM (http://www.ccafs-climate.org/data_spatial_downscaling/) proyectados a futuro cercano (2030), los cuales simulan trayectorias de forzamiento radiactivo, en las emisiones de gases de efecto invernadero, similares a las actuales (RCP 4.5 W·m^-2). Estas variables bioclimáticas tienen una resolución de un km² (30 segundos de arco); además, se incluyeron las variables pendiente (%), orientación (°) y altitud (m), obtenidas a partir del modelo digital de elevación con una resolución de 90 m en formato GeoTIFF. El tipo de vegetación se obtuvo de la serie VI “Uso de suelo y vegetación” del ^{INEGI (2016)} con una escala de 250 m. La edafología (tipo de suelo) de la serie I “Conjunto de datos vectoriales edafológico”, con una escala de 250 m, se obtuvo de la ^{CONABIO (1995)}. Estas últimas dos variables se descargaron en formato vectorial (shapefile). Todas las variables se estandarizaron y transformaron en formato ASCII a una escala de 1 km² (30 segundos de arco) con ayuda del programa ArcMap versión 10.3 (^{ESRI, 2014}).

Delimitación del área M

Las variables se ajustaron al tamaño del área de modelación espacial M (península de Yucatán), el cual se ha descrito como el espacio donde una especie está o se asume que podría estar, con base en el conocimiento biológico y de su capacidad de dispersión (^{Martínez-Méndez, Aguirre-Planter, Eguiarte, & Jaramillo-Correa, 2016}; ^{Soberón & Peterson, 2005}). A partir del área M, las variables ambientales se recortaron y ajustaron a un mismo tamaño de píxel de 1 km² (30 segundos de arco) con ayuda del programa ArcMap versión 10.3 (^{ESRI, 2014}).

Modelado de la distribución potencial actual y futura

Los modelos se generaron en el algoritmo MaxEnt versión 3.3.3k, el cual usa una técnica exploratoria e identifica sitios con valores similares de acuerdo con el principio de máxima entropía (^{Phillips, Anderson, & Schaphire, 2006}). Este algoritmo es uno de los más eficientes para modelar la distribución de especies utilizando solo datos de presencia (^{Elith et al., 2011}). Con el fin de aumentar la confiabilidad de los resultados de MaxEnt, 50 modelos de distribución se generaron y probaron bajo 10 criterios de modelación (Cuadro 2), utilizando diferentes reglas de aplicación de umbral y un número máximo de réplicas de 1 000, con un límite de convergencia de 0.00001 (^{Manzanilla et al., 2019}; ^{Phillips et al., 2006}). Esto permitió determinar cuál modelo se ajusta y explica mejor la distribución actual de la especie.

Los modelos se calibraron empleando 75 % de los datos de presencia de la especie en formato csv, seleccionados al azar para entrenar los modelos, y 25 % restante para la prueba de validación (^{Alba-Sánchez et al., 2010}; ^{Martínez-Méndez et al., 2016}; ^{Phillips et al., 2006}). La contribución relativa (%) de cada variable en los modelos se estimó mediante la prueba de jackknife (^{Shcheglovitova & Anderson, 2013}). La salida fue del tipo logística, el cual índica un índice de similitud (probabilístico) de condiciones idóneas para la especie con valores de 0 a 1, donde valores cercanos a 1 indican condiciones idóneas para el desarrollo y presencia de la especie (^{Coitiño, Montenegro, Fallabrino, González, & Hernández, 2013}; ^{Phillips et al., 2006}).

Para generar la modelación a futuro (2030), los parámetros del modelo de distribución actual con el mejor ajuste y desempeño estadístico se transfirieron al programa MaxEnt (^{Morrone & Escalante, 2016}).

Validación de los modelos

Los modelos se evaluaron por medio de la prueba del área bajo la curva (AUC), la cual se obtiene directamente de la técnica Receiver Operating Characteristic (ROC estándar) (^{Coitiño et al., 2013}; ^{Phillips et al., 2006}); no obstante, esta validación ha sido cuestionada por no considerar datos de ausencias verdaderas (^{Peterson, Papes, & Soberón, 2008}). Por tanto, fue necesario realizar un análisis de mayor robustez mediante una prueba de ROC parcial en el programa Tool for Partial-ROC versión 1.0 (^{Narayani, 2008}), la cual fue diseñada para contrarrestar las deficiencias del AUC (^{Peterson et al., 2008}).

Los análisis ROC parcial se realizaron con 50 % de los registros con una confiabilidad de 95 %, con 1 000 iteraciones mediante un remuestreo por bootstrap y fijando un error de omisión de 5 %. La prueba genera valores de 1 a 2; mientras más se aproxime a 2, el modelo será más confiable (^{Martínez-Méndez et al., 2016}; ^{Peterson et al., 2008}). Asimismo, se realizó una prueba de Z, donde el valor Z calculado de las proporciones AUC de ROC parcial debe ser mayor al valor Z de tablas (Z > 2.3 = 99 %), lo cual brindó un mayor soporte a los modelos (P < 0.01) (Martínez-Méndez et al., 2016; ^{Monterrubio-Rico et al., 2016}).

Áreas de distribución potencial

Los valores de la salida logística del modelo con el mejor ajuste estadístico se reclasificaron en tres categorías de calidad de hábitat con intervalos iguales (baja, media y alta), a través del módulo reclass del programa ArcMap 10.3 (^{ESRI, 2014}). El valor de la categoría de calidad de hábitat alta se utilizó para transformar los modelos de continuos a binarios (apto y no apto). Con base en la reclasificación, la superficie ocupada por las áreas de distribución potencial de C. odorata se estimó para los periodos actual y futuro.

Áreas de conservación y producción de semillas

Las áreas de conservación y de producción de semillas se seleccionaron con base en la identificación y delimitación de las áreas de mayor calidad de hábitat de los periodos actual y futuro. Un criterio adicional para determinar localidades con posibilidad de servir como refugio climático, para la conservación de la especie, fue la ubicación de las áreas naturales protegidas (ANP) (^{SEMARNAT & Comisión Nacional de Áreas Naturales Protegidas [CONANP], 2017}), cuyos límites se utilizaron y combinaron con las áreas de distribución determinadas para ambos periodos.

La selección de las áreas idóneas para la producción de semillas se basó en la evaluación de la calidad de estacion de los rodales (^{Huang, Titus, & Wiens, 1992}), los cuales reflejan la capacidad productiva de los sitios a través de la altura dominante de los árboles y de la relación alométrica diámetro-altura. Con base en lo propuesto por ^{Huang et al. (1992)} y ^{Yuancai y Parresol (2001)} se ha observado que los sitios con árboles de diámetros mayores suelen presentar dominancia en altura total, por lo que, acorde con lo reportado en la literatura, los sitios con mayores diámetros tendrían, teóricamente, mayor rendimiento en la producción de semillas. Las áreas idóneas para la producción de semillas se identificaron mediante los datos del diámetro máximo en formato raster del INFyS, a partir de los datos obtenidos de los conglomerados rectangulares de muestreo de 400 m², diseñados para selvas (^{CONAFOR, 2009}). Dicho archivo se transformó a formato vectorial en el programa ArcMap 10.3 (^{ESRI, 2014}); posteriormente, los polígonos de la distribución actual y futura se utilizaron para identificar las áreas donde la especie se encuentra. Finalmente, las áreas para la producción de semillas se eligieron con base en la reclasificación de los valores de la categoría diamétrica ≥50 cm; acorde con ^{Huang et al. (1992)} y ^{Yuancai y Parresol (2001)}, ese valor correspondería a un árbol de 20 m de altura, el cual se encontraría en etapa reproductiva.

Variables relevantes

La prueba jackknife indicó que las variables con mayor contribución en la generación de los modelos de distribución actual de C. odorata fueron la vegetación con 34.7 %, seguido de la precipitación del mes más húmedo con 14.6 %, edafología con 8.8 %, temperatura media del trimestre más frío con 8.6 % y pendiente con 7 %. En el análisis de la distribución potencial a futuro de los modelos GFDL_CM3 y MIROC_ESM, las variables ambientales más relevantes fueron la vegetación con 84 % y 83.7 %, temperatura media del trimestre más frío con 6.5 % y 9 %, y edafología con 2.8 % y 2.5 %, respectivamente (Figura 2). Estas variables son las que determinan el nicho ecológico actual y a futuro (2030) de la especie en la península de Yucatán.

Figura 2 Contribución de las variables ambientales en la distribución actual y futura (2030) de Cedrela odorata en la península de Yucatán. Veg = vegetación, Bio13 = precipitación del mes más húmedo, Edaf = edafología, Bio11 = temperatura media del trimestre más frío, Pend = pendiente. 

De acuerdo con la contribución de las variables ambientales, C. odorata requiere suelos del tipo Leptosol con precipitación del mes más húmedo de 216 mm y oscilación en la temperatura media del trimestre más frío de 22 a 24 °C. La especie tiene preferencia por los ecosistemas con vegetación secundaria arbórea y arbustiva de selvas medianas caducifolias, subcaducifolias, perennifolias y subperennifolias. Estos resultados concuerdan con los requerimientos ambientales de una especie forestal de zonas tropicales, los cuales han sido descritos en la literatura, lo que brinda un grado de confiabilidad de que la especie se encuentra en áreas congruentes a sitios con condiciones de hábitat similares a las descritas por ^{Rodríguez (2003)}, y ^{Pennington y Sarukhán (2005)}. Los resultados de las variables ambientales ya mencionadas son similares a los reportados por ^{Gómez et al. (2007)} y ^{Garza-López et al. (2016)}, quienes mencionan que las especies tropicales como cedro rojo y caoba suelen desarrollarse mejor durante las estaciones secas (no mayores a cuatro meses), en sitios con temperatura media anual de 24 a 28 °C, temperatura máxima del mes más cálido de 24 a 32 °C y temperatura mínima del mes más frío de 11 a 22 °C. ^{Andrade y Solís (2004)} señalan a C. odorata como una especie secundaria y primaria a la vez, debido a que suele ser pionera en zonas de vegetación secundaria en selvas medianas y altas, y también como elemento frecuente en el estrato superior de las selvas maduras. ^{Salazar et al. (2000)} mencionan que C. odorata requiere precipitaciones de 1 200 a 2 500 mm anuales y un gradiente altitudinal de 0 hasta 1 500 m para su óptimo desarrollo y crecimiento, condiciones que se localizan dentro de las zonas boscosas más densas hacia el sureste de la península de Yucatán. Las variables con mayor relevancia para el futuro (2030), en los dos modelos evaluados, fueron vegetación, edafología y temperatura media del trimestre más frío, las cuales se mantuvieron como las principales variables determinantes en la distribución de la especie. Lo anterior indica que los cambios o alteraciones en la vegetación asociada con C. odorata tendrían gran impacto en las áreas de distribución a futuro de la especie.

Ditribución potencial actual y futura

Los resultados indican una distribución potencial actual de C. odorata de 405 065 ha, equivalente a 2.9 % del área de estudio, donde dicha superficie reúne las condiciones ambientales más altas de calidad de hábitat, idóneas para la especie. La calidad de hábitat media representa una superficie de 19.2 % y la calidad de hábitat baja representa 77.9 %, la cual representa áreas marginales de casi nula idoneidad. Los resultados de esta clasificación permitieron generar un mapa de calidad de hábitat actual disponible para C. odorata en la península de Yucatán (Figura 3).

Figura 3 Clasificación de las áreas de calidad de hábitat de Cedrela odorata en la península de Yucatán. 

Por otra parte, los resultados de los modelos a futuro GFDL_CM3 y MIROC_ESM, proyectados a 2030, estiman un área potencial idónea de 277 037 y 226 138 ha, equivalentes a 1.6 y 2 % del área de estudio, respectivamente. Estos resultados representan una reducción importante entre 31 y 44.8 % del hábitat idóneo para la especie con respecto a su distribución potencial actual estimada (Figura 4).

Figura 4 Distribución potencial actual y futura al 2030 para Cedrela odorata en la península de Yucatán. 

Pese a que los escenarios de cambio climático predicen una reducción importante en el nicho ecológico de la especie, estos también indican que las poblaciones naturales de C. odorata encontrarían un refugio climático hacia las partes más densas de las selvas medias y altas, localizadas en el sureste del estado de Campeche, específicamente en la Reserva de la Biosfera de Calakmul. Los resultados sobre disminución de nicho ecológico concuerdan con las proyecciones de ^{Estrada-Contreras, Equihua, Laborde, Meyer, y Sánchez-Velásquez (2016)}, quienes modelaron la distribución actual y futura de C. odorata bajo escenarios de cambio climático en México, y reportan una disminución en sus áreas de distribución de 3.8 % bajo el escenario CGCM3 B1 para el año 2030 con respecto a su distribución actual. ^{Garza-López et al. (2016)} estimaron 60 % de reducción en el hábitat climático de S. macrophylla en la península de Yucatán para 2030. Los resultados de estos autores y los obtenidos en el presente estudio sugieren una vulnerabilidad media-alta por parte de las especies tropicales analizadas en el sureste mexicano. La combinación de altas temperaturas y la disminución en la precipitación serían la causa principal en la reducción de las condiciones idóneas de C. odorata en gran parte de la península; dicha combinación haría a la especie vulnerable a los incendios forestales y al ataque de plagas y enfermedades.

^{Villers-Ruiz et al. (2000)} señalan que los bosques tropicales del sureste de México se podrían convertir en bosques secos. A su vez, ^{Manzanilla-Quiñones y Aguirre-Calderón (2017)} predicen aumento en la temperatura promedio anual de 0.41 a 0.83 °C y disminución en la precipitación anual entre 35 a 71 mm en zonas tropicales del sureste mexicano para 2030, por lo que es muy probable que condiciones similares se presenten en la península de Yucatán para el mismo año.

Debido al calentamiento global, ya se están registrando modificaciones en la distribución de las especies, con cambios significativos en los taxones que habitan regiones montañosas, donde tienden a desplazarse altitudinalmente para poder subsistir (^{Lenoir, Gegout, Marquet, de Ruffray, & Brisse, 2008}; ^{Manzanilla-Quiñones, Aguirre-Calderón, Jiménez-Pérez, Treviño-Garza, & Yerena-Yamallel, 2019}; ^{Parolo & Rossi, 2008}; ^{Telwala, Brook, Manish, & Pandit, 2013}). Por tal razón, es prioridad realizar y analizar las simulaciones de los escenarios de cambio climático para determinar cuáles poblaciones están en riesgo de desaparecer por los efectos del cambio climático.

Áreas de conservación y de producción de semillas

Los modelos a futuro señalan que C. odorata conservará un área con mayor idoneidad de 88 603 a 64 342 ha, dentro de las reservas de la Biosfera Calakmul, Los Petenes y Ría Celestún (Figura 5). Estas reservas se sitúan principalmente en los municipios de Hopelchén, Champotón y Calkiní del estado de Campeche, y en el municipio de Halachó estado de Yucatán, donde las condiciones ambientales fungirían como refugios climáticos para la especie a futuro cercano (2030). Con respecto a la Reserva de la Biosfera Calakmul, los resultados son similares a lo reportado por ^{Garza-López et al. (2016}, ²⁰¹⁸⁾ para S. macrophylla y Lysiloma latisiliquum L. Benth. Estos resultados pueden ser útiles en la selección de áreas para la migración o recolonización asistida con fines de conservación y restauración de C. odorata en la península de Yucatán, como sugieren ^{Sáenz-Romero, Rehfeldt, Crookston, Duval, y Beaulieu (2009)} y la Secretaría del Convenio sobre la Diversidad Biológica (^{CDB, 2009}); además, el manejo adecuado de estas áreas beneficiará directamente a la protección de la especie.

Figura 5 Áreas de conservación de Cedrela odorata para el 2030 dentro de la poligonal de las áreas naturales protegidas en la península de Yucatán. 

Las áreas idóneas para la producción de semillas se localizaron principalmente en la Reserva de la Biósfera de Calakmul y Los Petenes, en el estado de Campeche, con diámetros máximos entre 50 a 80 cm y alturas de 10 a 30 m dentro de una superficie de 341 481 ha (Figura 6). En estas áreas hay árboles de cedro rojo con condiciones morfológicas y fenotípicas aptas para la producción y recolección de semillas y con los cuales es factible establecer actividades de conservación y restauración in situ. ^{Navarro-Cerrillo, Clavero, Vidaña, Quero, y Duque-Lazo (2016)} mencionan que definir áreas productoras de germoplasma ayudarían a la conservación y salvaguardarían el material genético de la especie de interés.

Figura 6 Áreas aptas para la producción de semillas de Cedrela odorata, en función de su diámetro máximo, en la Reserva de la Biósfera Calakmul y los Petenes en la península de Yucatán.