Área de estudio y especies estudiadas

México se ubica en el hemisferio norte, posee una superficie de 1 964 375 km² (^{INEGI, 2009}), se distribuye a ambos lados del Trópico de Cáncer, con elevaciones oscilan entre 1500 y 5000 msnm; predomina el clima cálido subhúmedo ((Am(f)) con temperaturas desde 2 a 30 °C, y precipitación anual de 0 a 4 000 mm (^{García, 1998}).

Las especies seleccionadas (símbolo - número de registros) fueron (Cuadro 1): Pinus arizonica Engelm. (Par - 541), P. ayacahuite Ehren. (Pay - 738), P. cembroides Zucc. (Pce - 1,244), P. devoniana Lindl. (Pde - 239), P. douglasiana Martínez (Pdo - 299), P. durangensis Martínez (Pdu - 1,217), P. engelmannii Carr. (Pen - 990), P. greggii Engelm. (Pgr - 33), P. hartwegii Lindl. (Pha - 64), P. herrerae Martínez (Phe - 338), P. lawsonii Roezl. (Pla - 61), P. leiophylla Schltdl. et Cham. (Ple - 906), P. lumholtzii Robins & Fern. (Plu - 812), P. maximinoi H. E. Moore (Pma - 76), P. montezumae Lamb. (Pmo - 127), P. oocarpa Schiede (Poo - 1,028), P. patula Schl. et Cham. (Ppa - 130), P. pringlei Shaw. (Ppr - 76), P. pseudostrobus Lindl. (Pps - 473) y P. teocote Schltdl. et Cham. (Pte - 830). La distribución altitudinal de muchas de estas especies varía entre 50 a 3000 msnm (^{Eguiluz, 1982}), incluso algunas pueden alcanzar el límite superior de la vegetación arbórea de hasta 4,000 msnm (^{Yeaton, 1982}). En los pinares la temperatura media anual fluctúa entre 6 y 28 °C y la precipitación desde 350 a más de 1000 mm año^-1 (García, 1988).

Base de datos: climáticas, topográficas y biológicas

Del sitio web WorldClim se obtuvieron 19 cubiertas geográficas de variables climáticas en formato raster de 1 km² de resolución (Cuadro 2), información derivada de temperatura y precipitación (^{Hijmans et al., 2005}) del periodo 1950 al 2000. La altitud, exposición y pendiente se calcularon en ArcMap 10.2 a partir del Continuo de Elevaciones Mexicano 3.0 (CEM 3.0) de INEGI. Un total de 10 222 registros biológicos de presencia (coordenadas geográficas) de las especies se obtuvieron de la base de datos del Inventario Nacional Forestal y de Suelos (INFyS 2009 - 2014) proporcionados por la Comisión Nacional Forestal (CONAFOR) y antes de su análisis se verificaron para que no existieran registros repetidos y que su distribución coincidiera sobre las áreas boscosas.

Modelación y validación de la distribución de especies

Para generar los modelos de distribución potencial se utilizó MaxEnt versión 3.3.3k (Maximum Entropy Species Distribution Modeling) (^{Phillips et al., 2006}), el cual solo requiere datos de la presencia de la especie. Para generar el modelo de cada especie se utilizaron todas las capas Bioclim y las topográficas, solicitando al software 500 interaciones y salida de tipo logística, este formato permite interpretarse como la estimación de probabilidad (de 0 a 1) de especies condicionada a variables ambientales (^{Phillips y Dudík, 2008}), el 50 % del total de registros se usó para entrenamiento y el resto para prueba. Los modelos fueron evaluados acorde con pruebas curva de respuesta (análisis de omisión/comisión) y sensibilidad Curva Operada por el Receptor (ROC) - Área bajo la curva (AUC) (^{Elith et al., 2006}; ^{Aguirre y Duivenvoorden, 2010}) y pruebas Jackknife para medir el efecto de cada variable en el modelo (^{Hijmans et al., 2005}; ^{Ramos-Dorantes et al., 2017}). Los modelos obtenidos en MaxEnt se reclasificaron en ArcMap 10.2, para obtener superficies en los diferentes niveles de probabilidad y con ello obtener el perfil bioclimático de cada especie. Un análisis de componentes principales se realizó usando la matriz estandarizada de contribución de cada variable bioclimática y topográfica para conocer su comportamiento en función a variables Bioclim.

Evaluación de los modelos de predicción

La prueba estadística área bajo la curva ROC muestra que los modelos de predicción son sensibles para aquellas especies con menos de 340 registros (Figuras 1A, B, C, E, F, G, H, K, L y O), al observar que las curvas de entrenamiento y de prueba se encuentran muy dispersas del umbral de predicción, es decir, sobrestimando o subestimando las proyecciones. Para especies con más de 470 registros el ajuste fue muy bueno, puesto que en el análisis de omisión/comisión (Figuras 1D, I, J, M, N, Ñ, P, Q, R y S) la tasa de omisión es cercana a la omisión estimada (línea diagonal), tanto de la muestra de entrenamiento como la de prueba. La evaluación de la predicción de los modelos mediante la técnica ROC ha sido aceptada como un método estándar y proporciona una medida simple del desempeño del modelo (^{Benito de Pando y Peñas de Giles, 2007}; ^{Aguirre y Duivenvoorden, 2010}).

Figura 1 Análisis del ajuste de modelos de distribución potencial para 20 especies de pino en México, mediante la técnica ROC. 

El análisis de resultados muestra que tanto el valor de AUC como el valor relativo de probabilidad (VRP) disminuyen significativamente (p≤0.001) conforme aumenta el número de registros de especies (Figura 2A y 2B). ^{Stockwell y Peterson (2002)} demostraron que la modelación de cualquier especie depende de la cantidad y calidad de los datos; de esta última, la precisión geográfica de su ubicación, incluso sugieren que para que la predicción sea robusta, se utilice un mínimo de 50 registros. Aunque se ha realizado modelación de especies con menos de 6 registros (^{Aguirre y Duivenvoorden, 2010}), ésta se ha reportado como no significativa. En nuestro estudio los registros variaron desde 33 en P. greggii hasta 1,244 en P. cembroides (Figura 1), y es claro que la modelación de la primera especie pudo no ser robusta por tener pocos registros, y el resto de las especies presenta más de 60 datos. De aquí se deriva que el ajuste del modelo aumenta con menor número de registros, pero subestimado y sobreestimado las predicciones (Figura 2A). ^{Elith et al. (2006)} demostraron que la estimación de parámetros en MaxEnt es crítica cuando el tamaño de muestra es pequeño, y viceversa.

Figura 2 Relación entre valores del área bajo la curva, AUC (A), valor relativo de probabilidad (B) y número de registros. Tendencia de valores de AUC con respecto al número de registros de cada especie (C). 

El promedio de los AUC de todas las especies fue 0.9589, con P. pringlei (0.982), P. patula (0.980) y P. maximinoi (0.971) como las especies con mayor ajuste, contrario a lo observado con P. oocarpa (0.922), P. leiophylla (0.942) y P. teocote (0.945) (Figura 2C), incluso estos últimos muy aceptables ya que ^{Phillips y Dudík (2008)} establecen que cuanto más próximo a 1 esté el valor de AUC mejor es el modelo. Además, este promedio supera el valor (0.75) que sugieren ^{Elith et al. (2006)} y ^{Aguirre y Duivenvoorden (2010)}, quienes aseguran que valores de AUC mayores a 0.75 indican que las predicciones basadas en solo datos de presencia son lo suficientemente precisas para establecer planes de manejo.

La modelación de la distribución potencial de coníferas de nuestro estudio mediante MaxEnt se puede considerar confiable. ^{Aguirre y Duivenvoorden, (2010)} registraron un promedio de AUC de 0.9238, menor al promedio de nuestro estudio (0.9589) pero variando de 0.9994 a 0.7660. En P. pringlei, P. patula, P. maximinoi, P. teocote, P. leiophylla y P. oocarpa estudiadas por estos autores, la modelación de todas ellas se mejoró en nuestro estudio con AUC desde 0.928 hasta 0.967. Con 391 registros de P. herrerae, ^{Ávila et al. (2014)} reportaron un valor de AUC de 0.9730, similar a nuestro estudio (Figura 1). Estas variaciones en ajuste del modelo y de zonas de predicción se deben básicamente al número de registros usados y su ubicación geográfica. En 11 especies modeladas por ^{Ramos-Dorantes et al. (2017)} y consideradas en nuestro estudio, el AUC promedio fue 0.9660, levemente superior al de nuestro estudio (0.9589). En el 54.5 % de ellas, nuestra investigación mejoró la modelación (P. ayacahuite, P. greggii, P. hartwegii, P. patula, P. pseudostrobus y P. teocote). El ajuste de los modelos no siempre es bueno, y en Pinus oocarpa Schiede ex Schltdl. el valor de AUC fue solo 0.548 y 0.445 en Abies guatemalensis Rehd. (^{Cruz-Cárdenas et al., 2016}), pero usando otras variables bioclimáticas; un valor cercano a 0.5 es realmente un valor al azar (^{Felicísimo et al., 2012}).

Variables ambientales que modelan la distribución de las especies

Al considerar las primeras tres variables en el modelo de cada especie, las que más contribuyeron a predecir su distribución son altitud, Bio10 y Bio1, en 10, 9 y 6 especies, respectivamente (Cuadro 3). Sin embargo, una sola variable (Bio4) puede explicar el 57.4 % de la distribución potencial de P. lawsonii. El perfil bioclimático es diferente entre especies (Cuadro 3) y depende de su distribución geográfica en el país. Al respecto, ^{Eguiluz (1982)} asoció la distribución de 77 taxas del género Pinus en México, y mostró que varias de ellas dependen de temperatura máxima (11 especies donde se registran hasta 45 °C, entre ellas P. pinceana, P. pseudostrobus var. estevezii Mart., P. greggii y P. oocarpa), mínima (hasta -23 °C, donde crecen P. ayacahuite, var. brachyptera, P. lumholtzii, P. engelmannii y Pinus arizonica), precipitación (P. radiata var. binata Lem y P. greggii en 100 y 2900 mm de precipitación anual), pero también de altitud (P. muricata D. Don y P. hartwegii Lind. crecen entre 120 hasta 4000 msnm, respectivamente).

El análisis de compontes principales muestra la existencia de dos grupos claramente definidos: un grupo compacto de especies (Pinus cembroides, P. durangensis, P. arizonica, P. lumholtzii, P. leiophylla, P. herrerae, P. leiophylla y P. ayacahuite) distribuidas especialmente a Bio1, Bio11, Bio18 y Bio19; y un segundo grupo destacando P. devoniana, P. douglasiana y P. oocarpa, dependiendo de variables Bioclim derivadas de precipitación, Bio13, Bio15 y Bio16 (Figura 3A). Pinus greggii demostró una respuesta diferente las otrasespecies, y su distribución depende de precipitación del mes más seco (Bio14) y exposición. Al considerar una contribución mayor al 10 %, la altitud aparece en la estructura del modelo en más de 55 % de las especies, le siguen Bio10 (50 %) y Bio1 (35 %); las variables que no explican más de ese valor son Bio6, Bio9, Bio15, Bio17, Bio19 y exposición (Figura 3B).

Figura 3 Análisis de componentes principales de la contribución de variables bioclimáticas y topográficas a la modelación de 20 especies de coníferas en México (A) y frecuencia en los modelos (B). La simbología se explica en Materiales y Métodos. 

Similar a nuestro estudio y derivado de la prueba de Jacknife, ^{Téllez et al. (2005)} reportan que las variables que explican la distribución de P. arizonica, P. devoniana, P. durangensis y P. pseudostrobus son Bio1 y Bio12. En 12 especies de coníferas evaluadas por ^{Ramos-Dorantes et al. (2017)}, Bio1 fue la variable más importante, y 11 de estas especies se incluyeron en este mismo estudio respondiendo de esta misma variable bioclimática. En Pinus herrerae Martínez, la temperatura máxima y altitud (1985 - 2227 msnm) determinan su distribución (^{Ávila et al., 2014}), y concuerda con nuestro estudio de que la exposición no explica la distribución de especies de coníferas. Una investigación realizada por ^{Cruz-Cárdenas et al. (2016)} a escala regional en ocho especies incluidas en nuestro estudio, contrasta con lo reportado aquí: Bio6 es la variable más importante para P. douglasiana (92.3 %), P. lawsonii (76.7 %), P. pseudostrobus (92.4 %) y P. teocote (92.3 %); la concordancia existe solo con P. montezumae (48.9 %) al observar que la altitud determina su distribución.

Superficie y mapas de distribución potencial

Las tres especies con mayor superficie potencial estimada por MaxEnt acorde a una probabilidad mayor o igual que 0.70 son: P. montezumae, P. devoniana y P. pseudostrobus con 14,744.8, 14,436.1 y 11,594.8 km² (Figuras 4 y 5O, L y Ñ),, mientras que las de menor superficie fueron P. durangensis, P. lumholtzii y P. engelmannii con 0.1, 275.9 y 819.1 km² (Figuras 4 y 5N, J y M). Oaxaca (+) y Durango (*) son los estados que registran superficie idónea para nueve especies (Figura 4), Chihuahua para siete, Guerrero para seis, Michoacán y Puebla para cinco (Figura 4), siete estados solo registran condiciones ideales para una especie, pero diferente. El 50 % de las especies registraron un valor de probabilidad mayor a 0.9 (Figuras 5A, B, C, E, F, G, H, L, O y S), y solo P. durangensis no posee superficie dentro de este rango (Figura 5M).

Figura 4 Superficie estimada por especie en un umbral de probabilidad mayor a 0.70 y los tres estados del país con mejores condiciones de idoneidad de ocurrencia de cada especie. 

Figura 5 Mapas de distribución potencial para veinte especies de Pinus en México. 

La comparación de las superficies potenciales de especies de coníferas en México, entre estudios, resultaría ambigua porque los estudios pueden estar representados con diferente escala de información y análisis espacial (^{Ávila et al., 2014}; ^{Cruz-Cárdenas et al., 2016}; ^{Ramos-Dorantes et al., 2017}); además, ^{Téllez et al. (2005)} y ^{Aguirre y Duivenvoorden (2010)} proporcionan superficies a nivel de especie y taxa sin indicar número de registros, y presentan superficies en diferentes rangos de probabilidad. De estos últimos autores, el perfil bioclimático y distribución geográfica es similar al encontrado aquí en P. arizonica, P. devoniana, P. durangensis y P. pseudostrobus.^{Stockwell y Peterson (2002)} coinciden que el número de registros es relevante en la calidad del modelo de predicción y por ende en la estimación de superficies la cual es importante, ya que ^{Cruz-Cárdenas et al. (2016)} indican que especies como P. leiophylla y P. teocote bajo diferentes escenarios de cambio climático reducirán su superficie. La delineación de la distribución geográfica de las especies estudiadas aquí, incluye perfectamente las cinco zonas geográficas reportadas por ^{Sánchez-González (2008)}, y son útiles en la evaluación de la distribución de las especies por efectos del cambio climático (^{Phillips y Dudík, 2008}).