Introducción
Dentro de la familia Meliaceae se reportan 50 géneros y alrededor de 1000 especies con una amplia distribución pantropical (Patiño, 1997; Pennington y Sarukhán, 2005). De las especies que destacan en esta familia están Cedrela odorata L. y Swietenia macrophylla King, por su alto valor comercial (Mendizábal-Hernández et al., 2009).
Particularmente Cedrela odorata tiene un rango altitudinal que varía entre 0 y 1200 m y tolera una amplia gama de condiciones climáticas, distribuyéndose en el continente americano desde el norte de México hasta el norte de Argentina (Gómez et al., 2007). En México se distribuye en el Pacifico, desde el sur de Sinaloa hasta Guerrero y la Depresión central y Costa de Chiapas, mientras que en el Golfo desde el sur de Tamaulipas hasta la Península de Yucatán (Patiño, 1997; Romo-Lozano et al., 2017).
Sus áreas de distribución y las poblaciones han disminuido considerablemente (Patiño, 1997; Rodríguez et al., 2003) por la extracción intensiva de su madera, en el país y de forma general en Mesoamérica, durante más de 200 años (Cavers et al., 2003). Dicha situación ha llevado a incluir a Cedrela odorata en el listado de la Convención sobre el Comercio Internacional de Especies Amenazadas de fauna y flora silvestres (CITES, 2007) y, en México, en la Norma Oficial Mexicana NOM-059-SEMARNAT-2010 sobre protección ambiental y especies nativas de flora y fauna como especie sujeta a protección especial (SEMARNAT, 2010).
La sobreexplotación de los recursos y la reconversión de la cubierta vegetal, además del cambio climático, reconocido como el principal reto que enfrenta la humanidad (Rodríguez y Mance, 2009), han ocasionado la disminución del hábitat de las especies (Cushman, 2006; Uribe, 2015), y el deterioro ambiental y la pérdida de biodiversidad (Challenger y Dirzo, 2009; Pérez y Ferreira, 2016). Para la reducción de estos problemas se deben tomar acciones de restauración, conservación y manejo de las especies (ONU, 2016) como lo son la inclusión nacional (SEMARNAT, 2010) e internacional (CITES, 2007) en algún estatus de protección a las especies más susceptibles para su aprovechamiento sostenible o su conservación futura. Para ello, los Sistemas de Información Geográfica (SIG) (Barrios, 2006; Gañan et al., 2015) y la modelación de nichos ecológicos con cualquiera de los métodos existentes (Phillips et al., 2006; Pliscoff y Fuentes-Castillo, 2011) son herramientas de gran utilidad.
En la actualidad y con el fin de contribuir al conocimiento de las especies y su interacción con el ecosistema, se han realizado trabajos de modelación de nichos ecológicos en combinación con los SIG para diferentes especies y con diversos enfoques. Por ejemplo, 1) definir áreas de distribución potencial para su conservación y restauración como es el caso de las palmas Thrinax radiata Lodd. ex Schult. & Schult. f. y Cryosophila argentea Bartlett en la Península de Yucatán (Hernández-Ramos et al., 2015) o la delimitación de nichos ecológicos para las cactáceas Ferocactus histrix (DC.) G.E. Linds., Mammillaria bombycina Quehl y M. perezdelarosae Bravo & Scheinvar en Aguascalientes, México (Meza-Rangel et al., 2014), 2) delimitar áreas de mayor probabilidad de presencia de Dendroctonus mexicanus Hopkins en bosques de coníferas de la Meseta Purhépecha en Michoacán, México, con la finalidad de emprender acciones de prevención y control de plagas (Martínez-Rincón et al., 2016), 3) delimitar áreas con mayor potencial de éxito en la reforestación y establecimiento de plantaciones forestales para Pinus sylvestris L., P. nigra J.F. Arnold, P. pinaster Aiton, P. halepensis Mill., Quercus ilex L. subsp. ilex, Crataegus monogyna Jacq. y Acer opalus Mill. subsp. granatense (Boiss.) Font Quer & Rothm. en la Provincia de Granada, España (Navarro-Cerrillo et al., 2016) y 4) identificar áreas de conservación y refugio ante escenarios de variación climática en Swietenia macrophylla en el sureste de México (Garza-López et al., 2016) o zonas futuras de cultivo con escenarios de cambio climático para Coffea arabica L. en Nicaragua (Läderach et al., 2012), entre otros. Sin embargo, en los bosques tropicales, los cuales son los más antiguos, diversos y ecológicamente complejos (Whitmore, 1993; Meli, 2003), hace falta mucha información ecológica que contribuya a entender las interacciones climáticas para el desarrollo de las especies que crecen en estos ecosistemas. Los bosques tropicales juegan un papel indispensable para enfrentar y reducir la degradación ambiental, así como para la generación de conocimiento que contribuya a la conservación, manejo y aprovechamiento de las especies incluidas en CITES (CITES, 2007) y la NOM-059-SEMARNAT-2010 (SEMARNAT, 2010). Sin embargo, son los ecosistemas más susceptibles de perder biodiversidad por los cambios en las condiciones ambientales (Pérez y Ferreira, 2016).
Bajo ese escenario, y considerando la importancia económica de la especie, se planteó el objetivo de determinar mediante modelos de simulación de nicho ecológico la distribución histórica y actual de Cedrela odorata en México, de igual forma modelar dos escenarios futuros de distribución con la finalidad de poder tener aplicación en los programas de conservación, manejo o aprovechamiento de la especie, así como en el establecimiento de plantaciones forestales comerciales (PFC) o delimitación de áreas productoras de germoplasma.
Materiales y Métodos
Área de estudio
El área de estudio se delimitó de acuerdo con el área de distribución natural en México para Cedrela odorata que reporta Patiño (1997), y la distribución mencionada del bosque tropical y bosque mesófilo considerado por Romo-Lozano et al. (2017) por ser una especie característica de estos ecosistemas.
Base de datos
Se obtuvo una base de 1747 coordenadas geográficas de Cedrela odorata construida a partir de la información contenida en 1) Inventario Nacional Forestal y de Suelos de 2004-2009 (INFyS 2004-2009) realizado por la Comisión Nacional Forestal (CONAFOR, 2016), 2) Red Mundial de Información sobre Biodiversidad (REMIB) elaborada por la Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) (CONABIO, 2015), 3) Base de datos de Missouri Botanical Garden (TROPICOS, 2018), 4) Herbario Nacional de México (MEXU) (UNIBIO, 2011) y 5) recorridos de campo en los años 2016 y 2017 en el área de distribución natural reportada. La base de datos fue depurada siguiendo algunos de los procedimientos que realizaron Monterrubio-Rico et al. (2016) como son la eliminación de puntos duplicados, referencias ubicadas en poblados, vías de comunicación principales, plantaciones forestales de la especie y lugares en los que las condiciones altitudinales o climáticas hicieran imposible su presencia, y se verificaron en campo algunas coordenadas particulares.
Modelación de distribución histórica y actual
Para la modelación se emplearon 19 variables climato-ambientales (BioClim) del periodo 1950 a 2000, obtenidas de Worldclim (WorldClim, 2016) a una resolución de 1 km2 por pixel (Hijmans et al., 2005), y las capas de altitud, orientación, pendiente y tipo de suelos de la Serie IV de INEGI (INEGI, 2014) a una escala homologada de 1: 1,000,000 obtenidas del Sistema Nacional de Información sobre Biodiversidad de CONABIO (Cuadro 1).
Clave | Variable ambiental |
---|---|
Bio1 | Temperatura promedio anual (°C) |
Bio2 | Oscilación diurna de la temperatura (°C) |
Bio3 | Isotermalidad (°C) |
Bio4 | Estacionalidad de la temperatura (desviación estándar *100) (°C) |
Bio5 | Temperatura máxima promedio del periodo más cálido (°C) |
Bio6 | Temperatura mínima del mes más frío (°C) |
Bio7 | Oscilación anual de la temperatura (°C) |
Bio8 | Temperatura media del mes más húmedo (°C) |
Bio9 | Temperatura media del mes más seco (°C) |
Bio10 | Temperatura media del cuatrimestre más cálido (°C) |
Bio11 | Temperatura media del cuatrimestre más frio (°C) |
Bio12 | Precipitación anual (mm) |
Bio13 | Precipitación del periodo más húmedo (mm) |
Bio14 | Precipitación del periodo más seco (mm) |
Bio15 | Estacionalidad de la precipitación (Coeficiente de variación, CV) |
Bio16 | Precipitación del trimestre más húmedo (mm) |
Bio17 | Precipitación del trimestre más seco (mm) |
Bio18 | Precipitación del cuatrimestre más cálido (mm) |
Bio19 | Precipitación del cuatrimestre más frío (mm) |
Pendiente | Pendiente (%) |
Altitud | Elevación (msnm) |
Orientación | Orientación (°) |
Tipo de suelo | Suelos (tipo) |
Los registros históricos mensuales se manipularon a través de las capas en formato raster de las 19 variables BioClim en ArcMap 10.4 (ESRI, 2017), realizando un recorte con una máscara de México definido por las capas de la Serie IV de INEGI. Para evitar un sobre ajuste de los modelos por colinealidad como lo mencionan Dormann et al. (2013) o que se tenga un error de tipo II (cada presencia no cuenta con un grado de libertad) de acuerdo a Miller et al. (2007), se realizó un análisis de correlación entre las 23 variables en ArcMap 10.4 con la herramienta Spatial Analysist Tools>Multivariate>Band Collection Statistics, utilizando las fórmulas siguientes para determinar las matrices de covarianza (Cov) (1) y correlación (Corr) (2) (ESRI, 2017):
y
Donde Z es el valor por celda, i y j son las capas raster, µ es la media de la capa, N representa el número de celdas, k es la denominación de cada celda en particular y δiδj es la desviación típica de las variables i y j.
La selección de las variables se realizó considerando solo las que tuvieran un coeficiente menor a 0.80 y -0.80 (Rissler y Apodaca, 2007; Fuentes et al., 2016) y que fueran las más significativas a nivel biológico para esta especie (Romo et al., 2013). Posteriormente se transformaron a un formato ASCII para su manipulación en MaxEnt (Steven et al., 2017).
La modelación de nicho ecológico se realizó en el programa MaxEnt versión 3.4.1 (Steven et al., 2017), el cual está fundamentado en el concepto de máxima entropía, donde el resultado de distribución es la respuesta de la especie a las variables medio ambientales de su entorno (Phillips et al., 2006). Este programa es considerado un algoritmo con un buen desempeño (Elith et al., 2011; Miranda et al., 2016), debido a que en la distribución proyectada utiliza la restricción de que cada valor esperado de las variables ambientales debe concordar con el promedio de su valor empírico, lo cual hace eficientes las predicciones a partir de información incompleta (Monterrubio-Rico et al., 2016).
Se utilizó 75% de la muestra para el entrenamiento (Martínez-Méndez et al., 2016) a través de 100 réplicas, 500 iteraciones (Loaiza y Roque, 2016) y un límite de convergencia de 0.00001 (Garza-López et al., 2016). Se aplicó la regla de umbral de máxima prueba de sensibilidad (plus) más especificidad y no la mínima presencia de entrenamiento que viene por default (Loaiza y Roque, 2016). También los tres tipos de modelos clásicos: Logistic (Porcentaje de la distribución de MaxEnt, cuyo valor crudo es al menos r), Cumulative (Registro crudo de r) y Raw (Modelo exponencial de MaxEnt). De igual manera se incluyó la corrección al modelo Logistic propuesta por Phillips et al. (2017) Cloglog, para visualizar el efecto que tiene cada uno de los diferentes modelos en la representación del nicho ecológico de la especie (Ortíz-Yusty et al., 2014). La validación se realizó empleando 25% de la muestra (Cartaya et al., 2016; Miranda et al., 2016) y la técnica Crossvalidate.
Los modelos generados fueron sometidos a las pruebas de la curva de respuesta con el análisis de omisión o comisión, de sensibilidad ROC (Receiver Operating Characteristic), el área bajo la curva (AUC por sus siglas en ingles) y la prueba de Jackknife para verificar la confiabilidad y el efecto de cada variable dentro del modelo (Scheldeman y van Zonneveld, 2010).
La validación de significancia estadística de predicción se realizó mediante la prueba binominal de omisión para todas las corridas utilizando los 11 test binomiales: valor fijo acumulado 1, valor fijo acumulado 5, valor fijo acumulado 10, el mínimo de puntos de presencia de entrenamiento, el percentil 10 de los puntos de presencia de entrenamiento, prueba de igualdad de la sensibilidad y especificidad de entrenamiento, valor máximo de la sensibilidad más la especificidad de entrenamiento, prueba de igualdad de la sensibilidad y especificidad de prueba, valor máximo de la sensibilidad más la especificidad de prueba, balance de los puntos de omisión de entrenamiento, área predicha y valor umbral; comparación de la entropía de las distribuciones originales y consideradas a un umbral determinado (Phillips et al., 2006). Para considerar admisible la modelación se requiere que todos sean significativos con un valor de p<0.01 (Romo et al., 2013).
La selección del mejor modelo se hizo de acuerdo con Peterson et al. (2008), donde los modelos que tengan valores en AUC cercanos a 1 y probabilidad 0 son los más consistentes. Se tomó en cuenta lo propuesto por Merow et al. (2013) y Ortíz-Yusty et al. (2014), quienes mencionan que la elección de los resultados se deberá hacer a la luz de los conocimientos biológicos de cada especie. En apego de la distribución reportada para Cedrela odorata, se tomó como referencia lo publicado por Pennington y Sarukhán (2005), representando los resultados bajo un modelo de consenso (Araújo y New, 2007; Marmion et al., 2009).
Modelación de distribución futura
La modelación futura a 2050 (promedio de 2041 a 2060) se realizó empleando las capas climáticas con una resolución espacial de 30 segundos de grado, información cartográfica resultado del Programa de Diagnóstico del Modelo Climático e Intercomparación (PCMDI por sus siglas en inglés) en su quinta fase para el modelo del sistema climático global versión 4 (CCSM4 por sus siglas en ingles), las cuales ya fueron corregidas por WorldClim versión 1.4 como clima de referencia actual (WorldClim, 2016).
Se emplearon las variables de pendiente (%), altitud (m), orientación (°) y suelos (tipo) con dos variantes de modelación. La primera modelación se efectuó usando las variantes de precipitación total mensual (pr) con variables climato-ambientales, y la segunda a través de variables bioclimáticas (bi) en el escenario futuro más conservador denominado RCP2.6 (Representative Concentration Pathways). En este escenario se plantea que el aumento de la radiación solar forzada será de más de 2.6 W/m2, que se eleve la concentración de CO2 de entre 350 y 400 ppm, y de 445 a 490 en el CO2 equivalente, con un incremento en la temperatura a nivel global de entre 2 y 2.4 °C (Weyant et al., 2009). La elección de las variables a utilizar fue mediante el análisis de covarianza y correlación (ESRI, 2017), y la regresión para obtener la distribución final se realizó con las combinaciones del mejor modelo resultante de la modelación histórica.
Elaboración de mapas de distribución histórica, actual y futura
Con el modelo que más se acercó a la realidad se elaboraron los mapas de distribución histórica, actual y futura reclasificando los resultados de 0 a 100% en ArcMap 10.4 (Hernández-Ruíz et al., 2016), tomando en consideración el valor máximo del AUC.
El mapa de distribución histórica se realizó con todas las probabilidades resultantes, mientras que el mapa de distribución actual solo con la probabilidad mayor de 50% (York et al., 2011). Los mapas de distribución futura se prepararon considerando todas las probabilidades para los dos escenarios.
Resultados
Base de datos y selección de variables
Se obtuvieron, posterior a la depuración, 1200 pares de coordenadas de Cedrela odorata, de las cuales se emplearon 900 datos para la modelación y 300 para la validación dentro del mismo programa. Al realizar las pruebas de covarianza y correlación entre las 23 variables, se observaron problemas de correlación debido a que los valores de las capas Bio1, Bio4, Bio10, Bio11, Bio13, Bio16 y Bio17 fueron superiores a 0.80 y -0.80 por lo cual se eliminaron, quedando 16 variables para modelar la distribución histórica y actual de la especie en México, sin tener problemas de sub o sobre ajuste en el modelo empleado.
Modelación y mapas de distribución histórica y actual
Al seleccionar las variables (16) se realizó la modelación con los cuatro tipos de regresiones de MaxEnt: Logistic, Cumulative, Raw y Cloglog, en las cuales no se observaron diferencias en los estadísticos de ajuste, pero sí en las distribuciones proyectadas en los mapas. Para los cuatro casos, las tasas de omisión en los registros de presencia de entrenamiento y prueba, y el área predicha en función del umbral acumulativo, son cercanas; además, ambas presentan una tendencia hacia la línea recta (Fig. 1A).
Los resultados de la aplicación de la técnica Receiver Operating Characteristic (ROC) indicaron que los cuatro modelos son adecuados, ya que el valor del cálculo del área bajo la curva (AUC) para los datos de entrenamiento (línea roja) y para los datos de prueba (línea azul) es de 0.926 y 0.921, respectivamente. Estos resultados están orientados a la parte superior izquierda por lo cual, y de acuerdo a la literatura, mientras ambas curvas se sitúen lo más juntas posible, mejor es el ajuste del modelo, lo que indicaría que no habría ningún error de omisión. De igual forma al emplear una clasificación aleatoria por el programa se obtuvo el valor máximo (0.5) (Fig. 1B).
Al validar la significancia estadística de predicción mediante la prueba binominal de omisión para todas las corridas utilizando los 11 test binomiales, todos ellos son significativos p<0.01, por lo cual se acepta el modelo empleado (Cuadro 2).
Umbral acumulativo | Umbral logístico | Descripción | Área predicha fraccional | Tasa de omisión de entrenamiento | Tasa de omisión de prueba | Valor de p |
---|---|---|---|---|---|---|
1 | 0.041 | Valor acumulado fijo 1 | 0.305 | 0.009 | 0.013 | <0.0001 |
5 | 0.181 | Valor acumulado fijo 5 | 0.218 | 0.035 | 0.045 | <0.0001 |
10 | 0.314 | Valor acumulado fijo 10 | 0.181 | 0.076 | 0.084 | <0.0001 |
0.57 | 0.023 | Presencia mínima de formación | 0.334 | 0 | 0.006 | <0.0001 |
12.4 | 0.356 | 10 percentil de presencia en el entrenamiento | 0.168 | 0.1 | 0.11 | <0.0001 |
18.321 | 0.45 | Sensibilidad y especificidad de entrenamiento iguales | 0.143 | 0.143 | 0.149 | <0.0001 |
5.714 | 0.209 | Sensibilidad máxima de entrenamiento más especificidad | 0.211 | 0.037 | 0.045 | <0.0001 |
17.286 | 0.436 | Sensibilidad y especificidad de prueba iguales | 0.147 | 0.134 | 0.149 | <0.0001 |
5.966 | 0.217 | Máxima sensibilidad de prueba más especificidad | 0.209 | 0.043 | 0.045 | <0.0001 |
0.57 | 0.023 | Equilibrar la omisión del entrenamiento, el área prevista y el valor umbral | 0.334 | 0 | 0.006 | <0.0001 |
5.152 | 0.185 | Equivale a la entropía de las distribuciones umbral y original | 0.217 | 0.035 | 0.045 | <0.0001 |
Al analizar la prueba de Jackknife, las variables con mayor porcentaje de contribución para modelar la distribución de Cedrela odorata en México fueron la precipitación anual (mm) (Bio12: 38.3%), la precipitación del periodo más seco (mm) (Bio14: 32.2%), la oscilación diurna de la temperatura (°C) (Bio2: 9.4%) y la altitud (m) (8.2%), mientras que las demás variables incluidas están por debajo de 2.3%. En el caso del factor de importancia de las variables con mayor peso, vuelven a ser altitud (27.2%), precipitación anual (Bio12: 16.1%) y la oscilación diurna de la temperatura (Bio2: 13.3%). Sin embargo, la oscilación anual de la temperatura (°C) y la precipitación del cuatrimestre más cálido (mm) tienen un factor mayor de 10% (Bio7 y Bio18) (Cuadro 3), situación que señala que las interacciones entre las variables son indispensables para la presencia o ausencia de la especie.
Variable | Porcentaje de contribución (%) | Factor de importancia (%) | Variable | Porcentaje de contribución (%) | Factor de importancia (%) |
---|---|---|---|---|---|
Bio12 | 38.3 | 16.1 | Tipo de suelo | 1.4 | 1 |
Bio14 | 32.2 | 0.7 | Bio5 | 1.2 | 0.6 |
Bio2 | 9.4 | 13.3 | Bio3 | 1 | 3.5 |
Altitud | 8.2 | 27.2 | Bio7 | 0.9 | 12 |
Bio8 | 2.2 | 5.1 | Bio15 | 0.8 | 3.7 |
Bio18 | 1.8 | 10.2 | Orientación | 0.6 | 1.9 |
Bio9 | 1.8 | 4.4 | Pendiente | 0.1 | 0.3 |
En la representación gráfica de los modelos de MaxEnt para Cedrela odorata los colores cálidos indican una alta probabilidad de presencia disminuyendo paulatinamente conforme el color cambia a azul. Los resultados muestran que el empleo de la regresión Cloglog es la que más se apega a la distribución potencial para la especie (Fig. 2A), caso contrario al emplear la modelación de tipo Raw, ya que la proyección del área potencial es muy restrictiva y compacta (Fig. 2B), seguido de la regresión Cumulative (Fig. 2C) y la regresión Logistic (Fig. 2D), siendo estas un punto intermedio para la modelación de la distribución histórica de Cedrela odorata en México.
Al definir que la mejor modelación es la obtenida por el modelo Cloglog se construyó el mapa de distribución actual para Cedrela odorata en México al emplear solo la probabilidad >0.50 (Fig. 3).
Modelación y mapas de distribución futura
Al realizar los análisis de covarianza y correlación en las capas de información futura se tiene que para la modelación con las variantes de precipitación (pr) solo nueve variables fueron las adecuadas para modelar la distribución futura, mientras que al emplear las variantes climáticas se modeló con 14 variables. Ambos escenarios fueron realizados con la regresión Cloglog, la técnica Crossvalidate y el umbral de máxima prueba de sensibilidad más (plus) especificidad. En los dos escenarios planteados la tasa de omisión, calculada para los registros de presencia de entrenamiento y de prueba, y el área predicha en función del umbral acumulativo son cercanas; además, en ambas su tendencia es hacia la línea recta (Figs. 4A, B). Los resultados de la prueba ROC indicaron ajustes adecuados, ya que el valor de AUC para los datos de entrenamiento (línea roja) y para los datos de prueba (línea azul) es de 0.924 y 0.926 para pr, respectivamente (Fig. 4C), y de 0.927 y 0.921 para bi (Fig. 4D). De igual forma al emplear una clasificación aleatoria en los casos se obtuvo el valor máximo (0.5).
Las pruebas binominales de omisión en las dos distribuciones son significativas p<0.01. Para RCP2.6 pr las variables que mayor aporte presentaron a la modelación futura de C. odorata en México (Fig. 5A) fueron temperatura media del cuatrimestre más cálido (°C) (58.2%), altitud (m) (12.1%), temperatura mínima del mes más frío (°C) (11.3%) y estacionalidad de la temperatura (desviación estándar *100) (°C) (10.4%). De igual forma la variable oscilación anual de la temperatura (°C) en el factor de importancia contribuyó en 26.7%. En el caso del modelo RCP2.6 bi las variables precipitación del periodo más seco (mm), precipitación anual (mm) y altitud (m) fueron las de mayor aporte a la modelación futura bajo esta variante con 36.8%, 30.7% y 8.2% respectivamente, mientras las variables de oscilación anual de la temperatura (°C) (24.5%) y oscilación diurna de la temperatura (16.3%) (°C) son las dos más altas en el factor de importancia (Fig. 5B).
Discusión
Al realizar la selección de las variables que no están correlacionadas entre sí, mejora la precisión de los resultados evitando el sobre o sub ajuste de la distribución histórica y futura de la especie, tal y como lo mencionan Miller et al. (2007) y Dormann et al. (2013) en sus estudios de análisis de colinealidad y correlación de las variables espaciales, y de lo aplicado por Obregón et al. (2014) al modelar la distribución de las especies de mariposa Pseudophilotes abencerragus Pierret y P. panoptes Hübner en España, así como Elith et al. (2006) al emplear la prueba de correlación de Spearman para mejorar los métodos de predicción de la distribución de las especies a partir de datos de presencia.
Los resultados muestran una buena eficiencia de los modelos aplicados para predecir la distribución histórica de C. odorata, siendo superior la regresión Cloglog en la modelación por apegarse con mayor fidelidad a la realidad, por lo cual fue empleada en la modelación futura. La superioridad de esta técnica de regresión se debe a que es una corrección de la regresión Logistic que deriva de la interpretación de MaxEnt versión 3.4.0 como un proceso de Poisson inhomogéneo (IPP), dándole una justificación teórica más fuerte a la transformación (Phillips et al., 2017).
Al comparar la tasa de omisión y el área predicha en función del umbral acumulativo con la línea de referencia (negra), son muy cercanas tanto en la modelación histórica como futura; por lo cual se define con eficiencia el umbral acumulativo para C. odorata (Figs. 1, 4). Estos resultados son similares a los reportados por Álvarez et al. (2013) al definir la distribución potencial y la probabilidad de presencia de cinco especies: Pinus sylvestris, P. nigra, P. pinaster, P. halepensis y P. pinea L. en Andalucía, España, empleando cinco variables orográficas, 17 climáticas y 15 edáficas, y por Hernández-Ruíz et al. (2016) al modelar con 22 variables la distribución potencial y características geográficas de poblaciones silvestres de Vanilla planifolia Andrews (Orchidaceae) en Oaxaca, México. Sin embargo, difiere con lo obtenido por Loaiza y Roque (2016) al modelar con 14 variables la distribución potencial de Armatocereus brevispinus Madsen (Cactaceae) en la región sur en Ecuador. Esta diferencia se asume que se debe al tipo de variables utilizadas en este último estudio, ya que emplea variables de nubosidad anual, frecuencia de heladas y presión de vapor anual que en el presente trabajo no se integraron.
El comportamiento de la curva ROC fue adecuado para la modelación histórica y para las dos variantes del escenario futuro, esto confirma que la modelación realizada con MaxEnt tuvo una buena capacidad predictiva, ya que es muy semejante la curva resultante al aplicar los datos de entrenamiento y el test aplicado por el programa. Al evaluar el área bajo la curva ROC (AUC) se registró un valor cercano a 1 en la clasificación, utilizando todos los datos de entrenamiento, y de 0.5 al emplear una clasificación aleatoria por el programa que es el valor máximo (Álvarez et al., 2013). Lo anterior concuerda con lo expuesto por Araújo et al. (2005), Araújo y Guisan (2006) y lo aplicado por Miranda et al. (2016), donde un valor mayor en AUC de 0.9 es satisfactorio en este tipo de modelaciones. De igual forma cumple con lo expuesto por Newbold et al. (2009) al mencionar que un modelo con valor mayor de 0.7 corresponde a un modelo de elevada precisión o alta discriminación.
Las modelaciones con MaxEnt mostraron significancia estadística en todos los test binomiales (p<0.01), lo cual señala la confiabilidad estadística de los resultados a 99% tal como lo indican Romo et al. (2013) al evaluar 76 test binomiales para modelar la distribución de las especies de mariposa: Boloria dia L., B. eunomia Esper, B. selene Denis & Schiffermüller, B. euphrosyne L., B. pales Denis & Schiffermüller y B. napaea Hoffmannsegg en la Península Ibérica, en el sudoeste de Europa, y Loaiza y Roque (2016) al evaluar solo los primeros tres test binomiales al realizar la distribución potencial de la cactácea Armatocereus brevispinus en la región sur del Ecuador.
La técnica de regresión que representa mayor fiabilidad, sin subestimar la superficie geográfica de distribución de C. odorata, fue de tipo Cloglog con el método de validación Crossvalidate, empleando el umbral de máxima prueba de sensibilidad más (plus) especificidad. Nuestros resultados difieren con lo aplicado por Cartaya et al. (2016) y Figueroa et al. (2016), donde la combinación de la regresión de tipo Cumulative y la técnica de Bootstrap fue la adecuada al modelar el nicho ecológico del roedor Cuniculus paca Brisson en Ecuador, y del oso Tremarctos ornatus Cuvier en Perú, respectivamente. Sin embargo, son semejantes con lo expuesto por Loaiza y Roque (2016) al modelar la distribución potencial bajo el umbral de mínima presencia de entrenamiento de Armatocereus brevispinus en Perú, donde mencionan que el empleo de la regresión tipo Logistic es superior por su fácil conceptualización e interpretación de los resultados por su escala binaria. Estas diferencias se asumen por el comportamiento biológico y los requerimientos mínimos de cada especie para su proliferación y distribución.
Esta combinación de procedimiento en la modelación realizada para C. odorata en México (Cloglog+Umbral de máxima prueba de sensibilidad más (plus) especificidad+ Crossvalidate) difiere con lo reportado por Merow et al. (2013) donde se menciona que para un análisis de idoneidad de hábitat y de influencia de las variables climáticas en la distribución potencial de las especies se deberá utilizar el formato Raw por ser un formato de salida crudo. Sin embargo, en este estudio se observó en el mapa de distribución potencial resultante con Raw una subestimación del área potencial de la especie en cuestión, indicando probabilidades de presencia muy restringidas a áreas específicas. Caso contrario es el empleo del umbral de máxima prueba de sensibilidad más (plus) especificidad, el cual ha demostrado ser uno de los métodos más robustos para generar este tipo de mapas de distribución a partir de mapas de probabilidades continuas (Liu et al., 2005; Ortíz-Yusty et al., 2014).
El área de distribución histórica modelada para C. odorata concuerda con lo reportado por Gómez et al. (2007) para el estado de Hidalgo, al modelar su distribución y algunos escenarios climáticos en esta región. Al comparar los resultados de este estudio con la distribución de C. odorata reportada por Pennington y Sarukhán (2005), para la selva mediana subperennifolia, muestran una disminución de su área de distribución. Por otro lado, si hacemos una comparación visual de los reportes anteriores, particularmente lo expuesto por Patiño (1997) y Rodríguez et al. (2003), donde se menciona que las poblaciones vegetales y su área de distribución han disminuido considerablemente, los resultados de este estudio podrían ser tomados como base para la actualización de la distribución histórica en México de esta especie así como para la implementación de planes de conservación del nicho ecológico de C. odorata y/o especies asociadas, tal como lo realizaron Hernández-Ruíz et al. (2016) para Vanilla planifolia en el estado de Oaxaca. De igual forma, el área de distribución potencial obtenida para C. odorata es muy similar a lo reportado por Garza-López et al. (2016) para Swietenia macrophylla en la región sur de México, al modelar la distribución contemporánea y futura de esta especie.
En el caso de la distribución futura bajo dos variantes (pr y bi) del escenario climático RCP2.6 a 2050, se observa que en ambos casos disminuye el área de distribución general de la especie. Particularmente en el escenario con la variante de precipitación (pr) en la región sur de México la disminución del área se ve marcada en el sur de Veracruz y Chiapas, Tabasco y la Península de Yucatán. El área de distribución de la especie se compacta en el norte de Chiapas y en la parte centro de Veracruz, mientras que en la región centro del país se mantienen las condiciones favorables de crecimiento. Al modelar la distribución futura empleando variantes climáticas (bi), la tendencia de disminución general del área de distribución y las áreas de compactación se mantienen, aunque en menor grado; sin embargo, este escenario señala un área de refugio muy marcada en la parte sur de Campeche y Quintana Roo.
Estos resultados pueden ser la pauta para realizar planes y programas de manejo y conservación a futuro de la especie, tal y como lo proponen Obregón et al. (2014) al estudiar la ecología, biología y distribución de Pseudophilotes abencerragus y P. panoptes empleando un escenario futuro, Ortíz-Yusty et al. (2014) al modelar la posible fluctuación de la distribución potencial de la tortuga Podocnemis lewyana Duméril en escenarios de cambio climático en Colombia, y Pliscoff y Fuentes-Castillo (2011) al realizar una revisión de las nuevas herramientas y enfoques empleados al modelar la distribución de las especies y ecosistemas. También podrán servir como base en estudios posteriores de esta especie en su migración asistida como lo han documentan Sáenz-Romero et al. (2009) y la Secretaría del Convenio sobre la diversidad Biológica (2009) en otras especies.
Cabe mencionar que la presencia y/o ausencia de cualquier especie está sujeta a limitaciones históricas, ya sea de tipo ambiental u orográfico, que determinan su dispersión o distribución geográfica, y que son el resultado de un proceso evolutivo completo (Maciel-Mata et al., 2015). Por lo anterior, una especie no siempre se encontrará en todas las áreas potenciales de distribución como lo mencionan Soberón y Peterson (2005) al explicar el diagrama de BAM, y Broennimann et al. (2006) al describir que este tipo resultados que representa el nicho ecológico idóneo y la amplitud de distribución de las especies son de tipo probabilísticos.
La modelación de la distribución histórica, actual y futura, así como sus áreas potenciales de presencia y/o proliferación, contribuyen a entender la ecología de distribución de cualquier especie. Particularmente en el caso de especies raras o con algún estado de conservación, tal y como lo mencionan Leal-Nares et al. (2012) al generar un modelo espacial basado en el conocimiento ecológico para Pinus martinezii E. Larsen, García-Aranda et al. (2012) al modelar la distribución actual y potencial para Taxus globosa Schltdl., y Cartaya et al. (2016) al identificar la distribución geográfica potencial del hábitat de Cuniculus paca, pues se da la pauta para poder planear actividades o realizar planes de conservación, restauración o manejo acordes con sus condiciones ambientales específicas (Meza-Rangel et al., 2014; Hernández-Ramos et al., 2015). También se pueden definir áreas productoras de germoplasma y/o zonas de establecimiento de poblaciones artificiales con fines de aprovechamiento o conservación (Navarro-Cerrillo et al., 2016), o bien realizar planes futuros de acuerdo con modelos de escenarios de distribución de las especies (Läderach et al., 2012; Garza-López et al., 2016).
Cabe mencionar que aunque los resultados obtenidos en el estudio fueron adecuados estadísticamente (Araújo et al., 2005; Liu et al., 2005; Araújo y Guisan, 2006; Miller et al., 2007; Newbold et al., 2009; Álvarez et al., 2013; Dormann et al., 2013; Merow et al., 2013; Romo et al., 2013; Ortíz-Yusty et al., 2014; Cartaya et al., 2016; Figueroa et al., 2016; Loaiza y Roque, 2016; Phillips et al., 2017), y concuerdan con otras aproximaciones para esta especie y el hábitat donde se desarrolla (Patiño, 1997; Rodríguez et al., 2003; Pennington y Sarukhán, 2005; Gómez et al., 2007; Garza-López et al., 2016), se deberá tomar en cuenta que las capas empleadas tienen un retraso de 17 años, por lo que la actualización constante de esta proyección, o bien el crear capas más actuales, mejoraría los resultados de este trabajo.
Conclusiones
Los resultados de la modelación de nicho ecológico con MaxEnt sugieren que es posible usar los mapas históricos para actualizar y determinar la distribución de Cedrela odorata en México. A partir de las variables analizadas con la prueba de Jackknife se mostró que la precipitación anual, la precipitación del periodo más seco, la oscilación diurna de la temperatura y la altitud (m), contribuyeron con más información para modelar la distribución de esta especie.
De los cuatro modelos usados, la regresión Cloglog es el más apegado con la distribución natural de la especie. La distribución histórica obtenida a partir de este modelo con respecto a la distribución actual reportada mostró que C. odorata tiene valores muy altos de probabilidad de distribuirse en la parte sur de la Península de Yucatán, norte y sur de Chiapas, y la llanura costera del Golfo en el estado de Veracruz. Esto sugiere que en la actualidad existe una reducción de la distribución de la especie en el país. La modelación futura bajo las dos variantes del escenario de cambio climático señala que las áreas de distribución de C. odorata seguirán disminuyendo con el paso de los años, por lo cual es indispensable planear acciones que mitiguen los efectos del cambio climático global en los bosques tropicales donde se desarrolla.