Área de estudio

La Faja Volcánica Transmexicana se localiza en las coordenadas 17° 30’ y 20° 25’ LN, y 96° 20’ y 105° 20’ LO (^{Ferrusquía-Villafranca, 1998}) (Figura 1); presenta altitudes de 256 a 5 650 m (^{Instituto Nacional de Estadística y Geografía [INEGI], 2018}). La topografía de la zona presenta relieve accidentado y con mucha pendiente en las montañas (^{Ferrusquía-Villafranca, 1998}). El clima dominante corresponde al tipo templado (húmedo y subhúmedo) con temperatura media anual de 12 a 18 °C y precipitación anual de 600 a 1 500 mm (^{García, 1998}).

Figura 1 Ubicación geográfica de la Faja Volcánica Transmexicana. 

Descripción de la especie

Abies religiosa es un árbol monoico, perennifolio, con altura de hasta 60 m y diámetro máximo de 1.80 m. El árbol tiene corteza grisácea, áspera y agrietada con grosor de 18 a 25 mm; ramas en forma de cruz; hojas simples y alternas de 15 a 35 mm de longitud por 1.5 mm de ancho; y conos solitarios de 8 a 16 cm de largo y 4 a 6 cm de ancho (^{Protectora de Bosques [PROBOSQUE], 2007}; ^{Rzedowski, 2006}).

En México, la especie se localiza en Morelos, Estado de México, Hidalgo, Puebla, Michoacán, Jalisco, Colima, Guerrero, Tlaxcala, Veracruz y Ciudad de México (^{Martínez-Méndez et al., 2016}) (Figura 2). En la Faja Volcánica Transmexicana, A. religiosa se distribuye en un gradiente altitudinal de 2 400 a 3 500 m y es la especie dominante del tipo de vegetación conocida en México como bosque de oyamel (^{Madrigal, 1967}; ^{Rzedowski, 2006}).

Figura 2 Distribución natural de Abies religiosa en México. 

Registros de Abies religiosa

Se descargaron 1 042 registros de A. religiosa de la plataforma Niche Toolbox (^{Osorio-Olvera, Vijay, Narayani, Soberón, & Falconi, 2017}). Con el objetivo de evitar la autocorrelación espacial y el sobreajuste de los modelos (^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}), se suprimieron las coordenadas duplicadas, mal georreferenciadas y cercanas entre sí (distancia >1 km), dejando un registro por cada celda de 1 km². Al final del proceso de depuración, se obtuvieron 341 registros de A. religiosa distribuidos en el área de estudio, con los cuales se generaron los modelos de distribución potencial para los periodos actual e histórico.

Información bioclimática actual e histórica

En el Cuadro 1 se muestran las 19 variables bioclimáticas disponibles en WorldClim versión 2.0 para el periodo actual (1970-2000) (^{Fick & Hijmans, 2017}) y de los modelos de circulación global (MCG) CNRMCM5 (Centro Nacional de Investigaciones Meteorológicas de Francia) y MIROC_ESM (Instituto de Investigaciones Oceánicas y Atmosféricas, Instituto Nacional de Estudios Ambientales y Agencia Japonesa para la Ciencia y Tecnología Marina y Terrestre) para el Holoceno medio. Estos últimos fueron generados a partir de Modelos Regionales del Proyecto CMIP5 (^{Coupled Model Intercomparison Project Phase 5, 2013}) del Panel Intergubernamental de Cambio Climático (IPCC). Todas las variables tuvieron resolución espacial de 1 km² y presentaron formato Geotiff (TIF). Además, se incluyeron las variables tipo de suelo (SUE; escala 1: 250 000) y altitud (DEM; modelo digital de elevación, resolución de 90 m), las cuales se descargaron en formato vectorial y ráster desde la plataforma del INEGI (^{INEGI, 2014}, ²⁰¹⁸). Al final, todas las variables se homogeneizaron a formato ASCII con escala de 1 km².

Selección de variables

La autocorrelación espacial se eliminó mediante un análisis de multicolinearidad (^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}), donde se descartaron aquellas variables con coeficiente de correlación ?#8805;0.85 (^{Manzanilla et al., 2019}; ^{Monterrubio-Rico et al., 2016}), para maximizar la contribución de las variables en los modelos de distribución (^{Martínez-Méndez et al., 2016}; ^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}).

Delimitación del área M

El área de modelamiento M es el espacio geográfico ambiental donde una especie se ha reportado o donde se supone que podría estar de acuerdo con el conocimiento biológico disponible (^{Martínez-Méndez et al., 2016}; ^{Soberón & Peterson, 2005}). Para delimitar el área M se utilizaron las subprovincias fisiográficas donde se localizaban los registros de A. religiosa: Llanuras y Sierras de Querétaro e Hidalgo, Chapala, Mil Cumbres, Lagos y Volcanes de Anáhuac, Neovolcánica Tarasca y Volcanes de Colima (^{INEGI, 2001}). El área M sirvió para recortar las variables a un mismo tamaño y con ello evitar la generación de áreas de distribución sobrestimadas (^{Monterrubio-Rico et al., 2016}).

Generación de los modelos de distribución

Los modelos de distribución se generaron con los registros de presencia de A. religiosa en formato CSV y con las variables BIO1, BIO2, BIO3, BIO4, BIO12, BIO14, BIO15, BIO18, BIO19 en formato ASCII (Cuadro 1), tipo de suelo (SUE) y altitud (DEM).

Las variables fueron introducidas al programa de máxima entropía (MaxEnt), donde se utilizó 75 % de los registros de la especie para realizar la prueba de entrenamiento y el 25 % restante para la prueba de validación de los modelos (^{Martínez-Méndez et al., 2016}; ^{Monterrubio-Rico et al., 2016}; ^{Phillips, Anderson, & Schaphire, 2006}). Se aplicó un replicado interno por validación cruzada y se desactivaron las casillas Extrapolate y Do clamping para evitar el sobreajuste de los modelos (^{Elith et al., 2011}). Las salidas de los modelos fueron de tipo logística, la cual representa un índice ambiental de probabilidades con valores de 0 a 1. Los valores cercanos a 0 señalan condiciones ambientales no aptas, mientras que los valores cercanos a 1 sugieren condiciones ambientales excelentes para el crecimiento y desarrollo de la especie (^{Coitiño, Montenegro, Fallabrino, González, & Hernández, 2013}; ^{Phillips et al., 2006}).

Se generaron y probaron cinco modelos de distribución con el fin de determinar por medio de las pruebas ROC (análisis de las características de funcionamiento del receptor del área bajo la curva [ABC]), Roc Parcial y Z, cuál modelo ajustaba con la distribución actual de la especie (Cuadro 2). Al modelo con mejor desempeño estadístico se sobrepuso el shapefile áreas naturales protegidas federales decretadas para la región (^{Comisión Nacional de Áreas Naturales Protegidas [CONANP], 2017}). La superficie actual de oyamel dentro de la poligonal de las áreas naturales protegidas se estimó con el programa ArcGIS 10.3 (^{Environmental Systems Research Institute [ESRI], 2014}).

Los parámetros del modelo de distribución actual con el mejor ajuste estadístico (^{Morrone & Escalante, 2016}) se transfirieron al programa de máxima entropía (MaxEnt), para generar los modelos de distribución en el pasado.

Validación de los modelos

Los modelos se evaluaron por medio de los valores del ABC de ROC, donde los valores entre 0.7 a 0.9 se consideran buenos y los mayores de 0.9 excelentes (^{Coitiño et al., 2013}; ^{Peterson et al., 2011}). No obstante, este tipo de validación ha sido cuestionado por no considerar ausencias verdaderas (^{Lobo, Jiménez, & Real, 2008}; ^{Peterson, Papes, & Soberón, 2008}), por lo cual fue necesario efectuar un análisis Roc Parcial en el programa Tool for Partial Roc versión 1.0 (^{Narayani, 2008}), para contrarrestar las deficiencias del ABC (^{Peterson et al., 2008}).

Se siguieron las recomendaciones de ^{Peterson et al. (2008)} usando 1 000 réplicas por bootstrap entre los archivos de presencia e idoneidad de hábitat y estableciendo un error de omisión de 5 %. La prueba genera valores de 1 a 2, donde un valor con promedio de radio 1.0 equivale a un modelo al azar (^{Garza-López et al., 2016}; ^{Lobo et al., 2008}; ^{Peterson et al., 2008}). Para determinar si los modelos eran válidos estadísticamente se hizo una prueba de Z entre las proporciones del ABC de Roc Parcial (^{Martínez-Méndez et al., 2016}). El mejor modelo se eligió con base en el valor más alto de las pruebas de Roc Parcial, desviación estándar menor y un valor confiable de Z (P < 0.01).

Los valores de las salidas logísticas de los modelos se reclasificaron en tres categorías de calidad del hábitat con intervalos iguales (baja, moderada y alta) con la herramienta reclass del programa ArcMap 10.3 (^{ESRI, 2014}). Los valores de la categoría de calidad de hábitat alta se emplearon para transformar los modelos continuos a binarios (apto-no apto). Con base en la reclasificación, se estimó la superficie de la distribución potencial de A. religiosa en la Faja Volcánica Transmexicana para ambos periodos analizados.

Las proyecciones de los modelos con el mejor desempeño estadístico para cada periodo se visualizaron y plasmaron en mapas de distribución potencial con el programa ArcMap 10.3 (^{ESRI, 2014}). Por último, las variables más importantes de cada periodo analizado se determinaron con la prueba de Jackknife, la cual permitió cuantificar la contribución de las variables ambientales en los modelos de distribución (^{Phillips et al., 2006}).

Distribución potencial histórica y actual

Dada la diversificación reciente del género Abies en el mundo (^{Xiang et al., 2015}) es muy probable que las condiciones climáticas de hace 6 000 a 12 000 años fuesen más cálidas (+2 °C) en zonas templadas y frías de Norteamérica (^{Caballero et al., 2010}; ^{Svensson et al., 2008}). Esto posiblemente delimitó la distribución de los bosques de coníferas hacia las zonas templadas y frías de la Faja Volcánica Transmexicana (^{Caballero et al., 2010}; ^{Madrigal, 1967}; ^{Martínez-Méndez et al., 2016}; ^{Rzedowski, 2006}).

De acuerdo con la prueba de Jackknife, las variables ambientales limitantes con mayor importancia en la distribución de A. religiosa en la Faja Volcánica Transmexicana fueron altitud (DEM con 56.5 %), precipitación anual acumulada (BIO12 con 18.4 %) y precipitación del trimestre más cálido (BIO18 con 8.3 %). Estas tres variables explicaron 83.2 % de la variabilidad del modelo. Estos resultados son similares a los encontrados por ^{Farjon y Filer (2013)}, ^{Madrigal (1967)} y ^{Rzedowski (2006)}, quienes indican que la humedad anual constante y lluvias abundantes en verano son esenciales para la distribución de este abeto. En el presente estudio se observó que el gradiente altitudinal idóneo para el crecimiento y desarrollo de A. religiosa fue de 2 129 hasta 3 687 m. Estos resultados son similares a los que se reportan en la literatura, ya que trabajos previos mencionan que la distribución de A. religiosa oscila de 2 400 hasta 3 660 m (^{Calderón de Rzedowski & Rzedowski, 2001}; ^{Madrigal, 1967}; ^{Nieto de Pascual-Pola, 1995}, ^{Rzedowski, 2006}). Cabe resaltar que, de acuerdo con los resultados de este estudio, la distribución del bosque de oyamel se puede encontrar 271 m por debajo de la altitud mínima y 27 m por encima de la altitud máxima reportada en la literatura. En cuanto a la precipitación anual acumulada, los modelos actuales estimaron un promedio de 1 125 mm. Este valor fue ligeramente mayor que el reportado por ^{Madrigal (1967)} y similar a lo mencionado por ^{Rzedowski (2006)} (1 000 mm y >1 000 mm, respectivamente). ^{Martínez-Méndez et al. (2016)} y ^{Sáenz-Romero et al. (2012)} mencionan que la precipitación de los meses de verano es determinante en la distribución de la especie, lo que concuerda con lo obtenido en el presente estudio; sin embargo, los estudios dendrocronológicos realizados para la especie, dentro de la zona de estudio, indican que la precipitación de invierno-primavera es de alta relevancia para el crecimiento de A. religiosa (^{Cerano-Paredes et al., 2014}; ^{Huante, Rincón, & Swetnam, 1991}). Lo anterior difiere con lo obtenido en este estudio, pues los modelos evaluados no consideraron a la precipitación de invierno-primavera como un factor clave en la distribución del oyamel. No obstante, dicha variable podría ser de mayor relevancia que la precipitación en los meses de verano, durante la etapa de crecimiento radial de la especie.

Por otro lado, las variables más importantes en la predicción de distribución histórica (modelo MIROC_ESM) fueron temperatura media anual (BIO1 con 65.9 %), precipitación anual acumulada (BIO12 con 18.6 %) y oscilación diurna en temperatura (BIO2 con 4.6 %); estas tres variables contribuyeron con 89.1 % en la variación del modelo. Para el modelo CNRMCM5, las variables de mayor importancia fueron la temperatura media anual (BIO1 con 62.4 %), la precipitación anual acumulada (BIO12 con 16.2 %) y la oscilación diurna en temperatura (BIO2 con 6.5 %), las cuales contribuyeron con 85.1 % del modelo. Los resultados del presente estudio concuerdan con los de ^{Madrigal (1967)} y ^{Rzedowski (2006)}, quienes mencionan que la oscilación diurna en la temperatura (11 a 16 °C) es determinante para el crecimiento de A. religiosa. Por lo tanto, se corroboró que esta variable ha sido un factor limitante desde hace 6 000 años (10.6 a 17.2 °C) en la distribución del oyamel en la Faja Volcánica Transmexicana. Con respecto a la temperatura media anual, los valores entre 4.6 y 15.7 °C son congruentes con lo reportado por ^{Rzedowski (2006)}, quien menciona que fluctúan entre 7 y 15 °C. Por otra parte, según las condiciones climáticas estimadas por los modelos de circulación global, la precipitación anual durante el Holoceno medio fue de 80 a 224 mm mayor que la actual y la temperatura promedio anual fue 1 °C más fría y no más cálida como lo mencionan ^{Svensson et al. (2008}) y ^{Caballero et al. (2010)}.

Debido a que las superficies estimadas y las condiciones ambientales descritas en la literatura y en este estudio presentaron similitudes, se puede suponer que la distribución potencial de A. religiosa ha permanecido estable a través del tiempo (según la teoría del conservadurismo de nicho ecológico; ^{Peterson, 2011b}). Así, en el presente estudio, las evidencias encontradas por los modelos actual e histórico parecen indicar que la distribución de la especie en la Faja Volcánica Transmexicana no ha cambiado en 6 000 años. Esto probaría (con la limitante de que no se cuentan con registros palinológicos) que el nicho ecológico de A. religiosa se ha mantenido durante los periodos evaluados. Cabe resaltar que, aunado a que la distribución actual e histórica no cambió significativamente, se logró identificar que sitios como el Nevado de Toluca y la Ciudad de México presentaron superficies mayores del oyamel en el pasado: +2 658.4 ha y +20 129.7 ha, respectivamente (Cuadro 4; Figura 5). Por lo tanto, se puede inferir que las condiciones climáticas de hace 6 000 años fueron más frías y húmedas en las zonas montañosas templadas y frías de la Faja Volcánica Transmexicana, lo que provocó que la distribución del oyamel se mantuviera ahí hasta el presente (^{Madrigal, 1967}; ^{Martínez-Méndez et al., 2016}; ^{Ramírez-Barahona & Eguiarte, 2013}; ^{Rzedowski, 2006}).

Figura 5 Distribución potencial (actual e histórica) del oyamel (Abies religiosa) en el Nevado de Toluca (A) y la Ciudad de México (B).