Objetivos

El objetivo principal del presente estudio fue modelar la distribución actual y futura de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal. Se plantearon los siguientes objetivos particulares; 1) modelar la distribución actual y futura (2041-2060) del bosque subalpino de Pinus hartwegii Lindl, 2) determinar las variables ambientales importantes en los periodos actual y futuro, 3) estimar la superficie que se vería reducida por el cambio climático y 4) identificar las áreas susceptibles al cambio climático.

Área de estudio

El Eje Neovolcánico Transversal se localiza entre las coordenadas geográficas 17° 30’ y 20° 25’ latitud norte y 96° 20’ y 105° 20’ longitud oeste y tiene una longitud de 930 km (^{Demant y Vincent, 1978}) (Fig. 1).

Figura 1 Área de estudio 

Datos de presencia

Se descargaron 247 registros (coordenadas GPS) disponibles de Pinus hartwegii Lindl, para México, de la aplicación Niche ToolBox (^{Osorio-Olvera, Vijay, Narayani, Soberón y Falconi, 2018}), de la Comisión Nacional para el Uso y Conocimiento de la Biodiversidad (Conabio) y se obtuvieron 96 registros GPS de la especie en campo.

Utilizando la aplicación de Niche Toolbox, se suprimieron las coordenadas duplicadas y con una distancia de separación menor a 1 km entre cada uno de los registros, dejando así, una sola coordenada por cada celda de 1 km². Este paso ayudó a evitar la sobreestimación de los modelos por los efectos de la autocorrelación espacial (^{Peterson y Nakazawa, 2008}; ^{Monterrubio-Rico et al., 2016}). Se dejaron solo los registros ubicados dentro del Eje Neovolcánico Transversal y se eliminaron aquellos ubicados en áreas urbanas y por debajo de los 3000 m snm, que es la altitud señalada por ^{Farjon et al. (1997)}, donde la especie forma bosques puros. Al final del proceso de depuración, se obtuvo una base de datos con un total de 65 registros para Pinus hartwegii Lindl; 60 de ellos fueron extraídos de la plataforma GBIF a través de la aplicación Niche Toolbox y cinco obtenidos en campo.

Obtención de información bioclimática actual y futura

Se descargaron las 19 variables bioclimáticas de Worldclim ver. 2.0 para el periodo presente (1970-2000) (^{Fick y Hijmans, 2017}) y futuro de los modelos de circulación global CRNMCM5 (Centro Nacional de Investigaciones Meteorológicas de Francia) y HADGEM2_ES (Centro Meteorológico Hadley del Reino Unido). Estos MCG fueron generados a partir de modelos regionales del proyecto de intercomparación de modelos acoplados fase 5 (CMIP5) con dos forzamientos radiactivos de 4.5 (trayectorias de CO₂ constantes) y 8.5 (trayectorias de CO₂ altas) para 2041-2060, tienen una resolución espacial de 1 km² y un formato GeoTiff (Tiff) (^{Worldclim, 2018}).

La Tabla 1 describe las 19 variables bioclimáticas descargadas para el presente y futuro.

Selección de variables

Debido a que las 19 variables bioclimáticas descargadas están correlacionadas entre sí, se procedió a eliminar la autocorrelación espacial mediante un análisis de multicolinearidad (^{Peterson y Nakazawa, 2008}; ^{Monterrubio-Rico et al., 2016}). La aplicación de dicho análisis ayudó a descartar aquellas variables con un coeficiente de correlación mayor o igual a 0.85, lo cual permitió maximizar la contribución de las variables en los modelos de distribución (^{Peterson y Nakazawa, 2008}; ^{Ortego, Gugger, Riordan y Sork, 2014}; ^{Martínez-Méndez, Aguirre-Planter, Eguiarte y Jaramillo-Correa, 2016}; ^{Monterrubio-Rico et al., 2016}). Se incluyeron también las variables pendiente (PEN; 1 km de resolución), suelo (SUE; escala 1; 250000), vegetación (VEG; serie V del Instituto Nacional de Estadística y Geografía [Inegi], escala 1; 250000) y altitud (DEM; modelo digital de elevación, resolución de 90 m) las cuales fueron descargadas en formatos vectorial y raster de las plataformas de Conabio e Inegi (^{Conabio, 2016}; ^{Conabio-Inifap, 1995}; ^{Inegi, 2013}; ^{Inegi, 2018}) y transformadas a archivos ASCII para integrarlas en la generación de los modelos de distribución.

Delimitación del área M

El área de modelamiento espacial, M, es el espacio geográfico y ambiental donde se ha registrado la presencia de una especie y se delimita de acuerdo con el conocimiento biológico que se tiene de ella y con su capacidad de dispersión (^{Soberón y Peterson, 2005}; ^{Martínez-Méndez et al., 2016}). Aunado a este proceso, el uso de barreras biogeográficas en modelación espacial evita generar áreas de distribución sobreestimadas de la especie de interés (^{Martínez-Meyer, 2005}).

Para delimitar el área M, se emplearon las subprovincias fisiográficas en donde se ubicaban los registros de presencia de la especie, las cuales fueron Lagos y Volcanes de Anáhuac, Neovolcánica Tarasca y Volcanes de Colima (^{Inegi, 2001}). El área M sirvió para recortar las variables ambientales a un mismo tamaño.

Modelos de distribución de especies

MaxEnt es un algoritmo que funciona mediante el principio de máxima entropía, el programa calcula la probabilidad de áreas idóneas para una especie a partir de coordenadas puntuales de presencia, no requiere datos de ausencia y genera buenas predicciones, aun cuando las muestras son pequeñas (^{Phillips, Anderson y Schaphire, 2006}; ^{Pearson, Raxworthy, Nakamura y Peterson, 2007}).

Para generar los modelos de distribución de Pinus hartwegii Lindl se utilizaron los registros de presencia de la especie y las variables ambientales seleccionadas del análisis de multicolinearidad (BIO1, BIO2, BIO3, BIO4, BIO12, BIO15, BIO16, BIO18, PEN, SUE, VEG y DEM) en formato ASCII.

Estas variables se cargaron al programa MaxEnt Ver. 3.3.3, utilizando 75% de los registros de presencia seleccionados al azar para generar los modelos y 25% restante para validarlos (^{Phillips et al., 2006}). Se desactivaron las casillas Extrapolate y Do clamping con el fin de evitar sobreajustar los modelos (^{Elith et al., 2011}), la salida de los modelos fue del tipo logística.

Se generaron y probaron 10 modelos de distribución con el fin de determinar, por medio de las pruebas estadísticas del área bajo la curva (ABC), Roc parcial y Z, cuál modelo explica mejor la distribución actual de la especie. Dichos modelos se construyeron a partir de diversos criterios, los cuales se especifican en la Tabla 2; siendo M5 y M10, los modelos generados a partir de los parámetros de omisión del programa (^{Phillips et al., 2006}) pero con diferente replicado interno.

Validación de los modelos de distribución

Los modelos se evaluaron por medio de los valores del área bajo la curva (ABC) del análisis de las características operativas del receptor (ROC). Los valores de esta prueba oscilan de 0 a 1: valores entre 0.7 y 0.9 clasifican el desempeño de los modelos como bueno y valores mayores a 0.9 los califican como excelente (^{Peterson et al., 2011}; ^{Coitiño, Montenegro, Fallabrino, González y Hernández, 2013}); sin embargo, esta validación ha sido cuestionada por no considerar datos de ausencias verdaderas (^{Peterson, Papes y Soberón, 2008}).

Se utilizó el programa Tool for Partial-ROC versión 1.0 (^{Narayani, 2008}). Donde los dos subconjuntos de datos (presencia e idoneidad de hábitat) fueron contrastados con el modelo reclasificado en probabilidades de MaxEnt para obtener los valores puntuales de adecuación o idoneidad. También se extrajo el número de pixeles de cada valor de probabilidad del modelo reclasificado. Los archivos preparados en formato CSV se ingresaron al programa Tool for Partial-ROC ver. 1.0. Para el análisis, se usó una proporción de puntos para Bootraps de 50%, con un número de muestras con reemplazo de 1000 réplicas y un umbral de omisión de 95%. Se promedió el valor de radio de las 1000 réplicas de los modelos. Los valores de esta prueba oscilan de 1 a 2: un valor con promedio de radio de 1 representa un modelo igual al azar (^{Peterson et al., 2008}; ^{Martínez-Méndez et al., 2016}; ^{Garza-López et al., 2016}). Se realizó la prueba de hipótesis z.

El modelo de distribución con el mejor ajuste para Pinus hartwegii Lindl se eligió con base en el valor más alto del análisis Roc parcial, menor error estándar y valor significativo de z. Los valores de idoneidad del modelo con el mejor ajuste estadístico se reclasificaron a un modelo de calidad de hábitat de tres categorías con intervalos iguales: baja, moderada y alta. Se utilizó el valor de la categoría de calidad de hábitat alta (> 0.60); con ayuda del módulo reclass del programa Arcmap 10.3® (^{Environmental Systems Research Institute [ESRI], 2014}), el modelo se transformó de continuo a binario (apto-no apto), dicha transformación ayudó a calcular y delimitar la distribución actual del bosque subalpino puro de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal.

Modelación bajo escenarios de cambio climático

Para generar los modelos de distribución bajo escenarios de cambio climático, se transfirieron los parámetros del modelo de distribución actual con el mejor ajuste estadístico (^{Morrone y Escalante, 2016}) al programa MaxEnt Ver. 3.3.3. Se cargaron las variables ambientales en formato ASCII y los registros de la especie en formato CSV. Se utilizaron 75% de los registros de la especie para generar los modelos y 25% restante para validarlos (^{Phillips et al., 2006}), se desactivaron las opciones Extrapolate y Do clamping con el fin de evitar sobreajustar los modelos de distribución (^{Elith et al., 2011}).

Validación de los modelos bajo escenarios de cambio climático

Los modelos fueron evaluados por medio de los valores ABC del análisis ROC (^{Phillips et al., 2006}) y las pruebas de Roc parcial (^{Narayani, 2008}) y de Z (^{Monterrubio-Rico et al., 2016}). Se reclasificaron los valores de idoneidad en tres categorías; baja, moderada y alta y se utilizó el valor de la categoría de calidad de hábitat alta (> 0.60) como umbral para convertir los modelos continuos a binarios (apto-no apto).

Variables determinantes en la distribución actual y futura

Por medio de la prueba de Jackknife se determinó la contribución de las variables en la generación de los modelos de distribución (^{Phillips et al., 2006}). Dicha prueba permitió identificar los requerimientos ambientales adecuados para que Pinus hartwegii Lindl crezca, se desarrolle y subsista en la actualidad y a futuro en el Eje Neovolcánico Transversal.

Escenarios de cambio climático

De acuerdo con el forzamiento radiactivo empleado, los escenarios de cambio climático simulan la posible afectación de un conjunto de variables ambientales en la distribución de una especie. Las trayectorias se clasifican en conservadoras (RCP 4.5; trayectorias de CO₂ constantes) y extremas (RCP 8.5; trayectorias de CO₂ altas), las primeras predicen un aumento en temperatura promedio anual de 0.9 °C a 2 °C y las segundas un incremento de 1.4 °C a 2.6 °C en la temperatura promedio anual global para el 2050 (^{IPCC, 2013}; ^{Fernández-Eguiarte et al., 2015}). Los escenarios de cambio climático conservadores hacen referencia a una pérdida menor de superficie con las condiciones aptas para que la especie crezca y sobreviva, mientras que los escenarios extremos predicen pérdidas mayores en la distribución potencial de una especie.

Utilizando la información vectorial disponible de las áreas naturales protegidas federales (ANP) de la ^{Conanp (2017)} y sobreponiendo los archivos de las salidas binarias de los modelos actual-futuro con mejor desempeño estadístico en el programa Arcmap 10.3® (^{ESRI, 2014}), fue posible identificar las áreas susceptibles y estimar la distribución reducida por efectos del cambio climático en la poligonal de las ANP del área de estudio.

Modelación bajo escenarios de cambio climático

Los valores del ABC de la prueba de ROC fueron de 0.983 a 0.984 para el entrenamiento y de 0.960 a 0.970 para la validación, lo cual indica que el desempeño de los modelos fue excelente. Los resultados de las pruebas de Roc parcial y Z presentaron un desempeño estadístico excelente (Tabla 4).

Los resultados de las pruebas de Roc parcial, error estándar y Z indican un mejor ajuste en las trayectorias RCP para el modelo de circulación global CRNMCM5 en comparación con las trayectorias del modelo HADGEM2_ES.

Variables determinantes en la distribución actual y futura

Las variables ambientales relevantes actuales fueron BIO1, BIO18, BIO12, SUE y VEG, las cuales aportaron 84.9% de la variabilidad de los modelos (Fig. 3).

Figura 3 Variables importantes en los modelos actuales 

Las variables ambientales relevantes para 2041-2060, de acuerdo con el modelo CRNMCM5 para las trayectorias RCP 4.5 y 8.5, fueron BIO1, SUE, BIO16, PEN y VEG, con una contribución de 94.2% y 92.6%, respectivamente (Fig. 4).

Figura 4 Variables importantes del modelo CRNMCM5 a futuro 

Por otra parte, las variables relevantes para 2041-2060 para las trayectorias RCP 4.5 y 8.5 del modelo HADGEM2_ES fueron DEM, BIO1, SUE, VEG y BIO12, con una contribución de 88.1% y 88.5%, respectivamente (Fig. 5).

Figura 5 Variables importantes del modelo HADGEM2_ES a futuro 

Los resultados de los modelos de distribución señalan a la BIO1 y DEM como las variables determinantes en la distribución actual y futura del bosque subalpino de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal.

Escenarios de cambio climático

Los resultados de la trayectoria RCP 4.5 prevén un aumento en temperatura promedio anual entre 0.7 °C (CRNMCM5) y 1.5 °C (HADGEM2_ES), dichos aumentos podrían provocar reducciones entre 11.7% (CRNMCM5) y 18.3% (HADGEM2_ES) en la distribución actual del bosque subalpino de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal para 2041-2060.

Los resultados de la trayectoria RCP 8.5 predicen un aumento en temperatura promedio anual entre 1.1 °C (CRNMCM5) y 2.1 °C (HADGEM2_ES); ante tales aumentos en temperatura promedio anual, la distribución actual del bosque subalpino de Pinus hartwegii Lindl se vería reducida en 18.3% (CRNMCM5) hasta 26.9% (HADGEM2_ES) para 2041-2060.

En la Tabla 5 se describe la superficie actual y la estimada para el bosque subalpino de Pinus hartwegii Lindl de los dos modelos de circulación global probados para 2041-2060.

Debido a la gran escala del área de estudio, esta se dividió en cuatro zonas, las cuales abarcan la distribución subalpina actual y futura de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal. La zona 1 abarca las áreas naturales protegidas Volcán Nevado de Colima, la zona 2 el Pico de Tancítaro, la zona 3 comprende el Nevado de Toluca e Izta-Popo y la zona 4 está conformada por La Malinche, El Cofre de Perote y El Pico de Orizaba. La clasificación del área de estudio por zonas permitió identificar las áreas naturales protegidas más susceptibles a los efectos del cambio climático.

La Tabla 6 muestra las superficies actuales y futuras estimadas por zonas para el bosque subalpino de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal.

En la Tabla 6 se puede observar que los modelos bajo escenarios de cambio climático predicen que las zonas 3 y 4 serían las más afectadas por los efectos del cambio climático a futuro.

La Figura 6 muestra la distribución de las condiciones ambientales idóneas actual y futura de los escenarios de cambio climático conservador (CRNMCM5 RCP 4.5) para el bosque subalpino de Pinus hartwegii Lindl, dentro de la poligonal de las áreas naturales protegidas del Eje Neovolcánico Transversal.

Figura 6 Distribución actual y futura del escenario conservador (CRNMCM5 RCP 4.5) para el bosque subalpino de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal.  

De acuerdo con los resultados del escenario climático conservador (MCG CRNMCM5 RCP 4.5) para 2041-2060, los sitios más susceptibles al cambio climático serían el Nevado de Toluca, con pérdidas de hasta 53.1% (2853.5 ha), el Pico de Orizaba con una reducción de hasta 33.5% (3292.1 ha) y La Malinche con una disminución de hasta 19.6% (54.3 ha) en la distribución actual de Pinus hartwegii Lindl. Sitios como el Volcán Nevado de Colima (+1.9%, 61 ha) y el Cofre de Perote (+3.5%, 160.8 ha) presentarían un aumento en la superficie idónea, mientras que el Pico de Tancítaro presentaría una ligera reducción (0.7%, 27.2 ha). Estos resultados indican que estos sitios no serían tan afectados por el cambio climático; sin embargo, dichas modificaciones sí cambiarían sus distribuciones para 2041-2060. Por otra parte, el Izta-Popo presentaría un aumento de hasta 37.4% (361.5 ha) para 2041-2060.

Variables determinantes para Pinus hartwegii Lindl

Las variables relevantes para ambos periodos analizados fueron BIO1, DEM, BIO18, BIO12, SUE y VEG. Una oscilación de la temperatura promedio anual de 6.8 °C a 13.1 °C es importante en el crecimiento y desarrollo de la especie, estos resultados son ligeramente menores a los obtenidos por ^{Gutiérrez y Trejo (2014)}, quienes mencionan que la temperatura promedio anual para el crecimiento de Pinus hartwegii Lindl es de 9.9 °C a 13 °C. También resultan ser menores a los obtenidos por ^{Martínez-Méndez et al. (2016)} para Abies religiosa (Kunth) Schltdl & Cham (12 °C a 16.4 °C), especie que coexiste con Pinus hartwegii Lindl en las partes altas del Eje Neovolcánico Transversal.

^{Martínez-Méndez et al. (2016)} mencionan que la DEM es un factor importante en la distribución de cinco especies de abetos mexicanos (Abies religiosa, Abies concolor, Abies guatemalensis, Abies hickelii y Abies durangensis var. coahuilensis), al igual que para estas especies, la altitud es un factor importante para la distribución de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal.

Para la BIO18, los mismos autores del párrafo anterior encontraron que esta variable es relevante en la distribución de cinco especies de abetos mexicanos (Abies concolor, Abies guatemalensis, Abies durangensis var. durangensis, Abies durangensis var. coahuilensis y Abies vejarii). ^{Sáenz-Romero et al. (2012)} mencionan que la precipitación de los meses de verano es importante en la distribución actual de Abies religiosa en el Eje Neovolcánico Transversal. Por lo que la variable (BIO18) es determinante en la distribución actual de especies de la familia Pinaceae que se localizan en las partes altas de México. De acuerdo con los resultados obtenidos en este estudio, los valores adecuados de la precipitación durante los meses de verano para el crecimiento de Pinus hartwegii Lindl son de 316 mm a 575 mm.

La precipitación total anual idónea para el desarrollo de Pinus hartwegii Lindl es de 1304 mm; de acuerdo con los resultados de los modelos a futuro, se prevé una disminución entre 2.8% y 3.8% en la precipitación promedio anual para 2041-2060. Estos resultados difieren de los estimados por ^{Sosa-Rodríguez (2015)} y ^{Sáenz-Romero et al. (2010)}, quienes mencionan que la precipitación promedio anual en México disminuiría hasta 20% para 2020 y 9% para 2060. ^{Manzanilla y Aguirre (2017)} estimaron una disminución de 2.6% a 5.3% en la precipitación promedio anual de la reserva de la biosfera Selva el Ocote, Chiapas para 2030; resultados similares a los obtenidos en el presente estudio.

Para la variable SUE, trabajos como el ^{Cruz-Cárdenas et al. (2016)} mencionan que el pH del suelo es importante en la distribución de seis especies de pináceas en el estado de Michoacán, pero no sustancial para Pinus hartwegii Lindl, aunque para este trabajo, el suelo fue ligeramente importante. Estudios como el de ^{Galicia et al. (2016)} sugieren que las pináceas del centro de México prefieren suelos del tipo Andosol, los cuales son ligeramente ácidos. En este trabajo, las zonas 2, 3 y 4 presentaron suelos del tipo Andosol.

En cuanto a la VEG, de acuerdo con ^{Inegi (2013)} y con lo obtenido en este trabajo, el bosque subalpino de Pinus hartwegii Lindl del Eje Neovolcánico Transversal coincide con los tipos de vegetación pradera alta de montaña, bosque de pino y bosque de oyamel. Trabajos como los de ^{Estrada-Contreras et al. (2015)} e ^{Islas-Báez et al. (2015)} mencionan que los tipos de vegetación son importantes en la distribución potencial actual y futura de una especie.

Bajo el contexto de protección legal de especies de la familia Pinaceae en México, estudios como el de ^{Aguirre y Duivenvoorden (2010)} mencionan que 70% de los 56 taxa del género Pinus modelados en México tienen solo 10% de su distribución natural resguardaba dentro de la poligonal de alguna área natural protegida federal.

Por otra parte, ^{Martínez-Méndez et al. (2016)} hallaron que nueve de los diez taxa del género Abies registrados en México presentan su área de distribución fuera la poligonal de las áreas naturales protegidas. En el presente trabajo, el bosque subalpino actual de Pinus hartwegii Lindl se encuentra protegido en 86.5% por las áreas naturales protegidas que se encuentran en el Eje Neovolcánico Transversal. Sin embargo, dichas áreas no están exentas de los efectos del cambio climático. Los modelos a futuro prevén que la temperatura promedio anual en el bosque de Pinus hartwegii Lindl aumente entre 0.7 °C y 2.1 °C por lo que, de acuerdo con el ^{IPCC (2014)}, aumentos de 1 °C a 2 °C representarían un daño moderado a la biodiversidad y arriba de 3.5 °C sería perjudicial para muchas especies.

^{Garza-López et al. (2016)}, utilizando trayectorias RCP de 6.0 (promedio de 18 Modelos de Circulación Global), estimaron una reducción de 60% en el hábitat climático de la caoba (Swietenia macrophylla) en la Península de Yucatán para 2030. En este trabajo, se prevén reducciones entre 11.7% (escenario RCP 4.5 MCG CRNMCM5) y 26.9% (escenario RCP 8.5 MCG HADGEM2_ES) en la distribución actual de Pinus hartwegii Lindl. Siendo el Nevado de Toluca y el Pico de Orizaba los sitios más susceptibles al cambio climático de acuerdo con los modelos proyectados a futuro.

^{Villers y Trejo (1998)} y ^{Sáenz-Romero et al. (2012)} mencionan que los bosques secos aumentarían su distribución (7.4%), mientras que los bosques templados y fríos reducirían sus nichos ecológicos entre 65% y 87.6%, respectivamente; además, sitios como la reserva Mariposa Monarca presentarían condiciones más cálidas (2.3 °C) para 2060 (^{Sáenz-Romero et al., 2010}), condiciones climáticas similares a las que se llegarían a presentar, de acuerdo con los resultados obtenidos en este trabajo (0.7 °C a 2.1 °C) para 2041-2060 para la misma región geográfica. Sáenz-Romero et al., 2010 mencionan para que una especie arbórea pueda subsistir a los cambios en el clima deberá migrar altitudinalmente 100 m por cada 0.5 °C de aumento. Los resultados del presente trabajo sugieren un aumento en temperatura promedio anual entre 0.7 °C y 2.1 °C para 2041-2060; antes tales cambios, Pinus hartwegii Lindl estaría obligado a desplazarse entre 100 m y 500 m arriba, por lo que los bosques puros de la especie tendrían que reubicarse entre 3100 m y 3500 m snm y no en 3000 m snm, que señala actualmente ^{Farjon et al. (1997)}.

^{Sáenz-Romero et al. (2012)} indican que el aumento en temperatura media del mes más cálido reduciría en 87.6% la cubierta forestal de Abies religiosa para 2060 en el Eje Neovolcánico Transversal; por otra parte, el aumento en temperatura promedio anual para 2050 reduciría hasta 26.9% de la superficie de Pinus hartwegii Lindl para esa misma región geográfica.

^{Estrada-Contreras et al. (2015)} hallaron un aumento de nicho de 3.7% para Veracruz para 2050, resultado similar a lo obtenido en este estudio para el Cofre de Perote (+3.5%) para 2041-2060. ^{Cruz-Cárdenas et al. (2016)} prevén una reducción de 38.3% en la superficie de Pinus hartwegii Lindl en Michoacán para el horizonte 2015-2039, resultado que es menor a lo obtenido en el presente estudio para el Nevado de Toluca (53.1%). ^{Islas-Báez et al. (2015)} mencionan que la variable más importante en la distribución de Abies religiosa para 2030 y 2050 en la Reserva de la Biosfera Mariposa Monarca es la precipitación anual, variable que difiere con las variables determinantes (BIO1 y DEM) obtenidas de los modelos generados para este estudio.

De acuerdo con lo analizado en los trabajos anteriores y el presente estudio, los escenarios de cambio climático predicen panoramas no muy alentadores para el bosque subalpino de Pinus hartwegii Lindl en el Eje Neovolcánico Transversal. Ante la posible pérdida en la superficie boscosa subalpina, una de las consecuencias inmediatas en dicho ecosistema sería la reducción y función de los servicios ecosistémicos como la captura de carbono, recarga de acuíferos y regulación del clima local, lo cual afectaría directamente a las poblaciones rurales aledañas a los ecosistemas de alta montaña del centro de México.