Introducción
El germoplasma forestal utilizado en México para la producción de plántulas para abastecer los programas de reforestación nacional proviene de poblaciones naturales o de plantaciones forestales sin manejo, en las que no se considera la calidad fenotípica y genotípica de los individuos para su recolección (Muñoz et al., 2014). Además, parte de la planta producida es utilizada de forma indiscriminada para la reforestación de áreas en condiciones edafoclimáticas diferentes a las zonas de procedencia del germoplasma (Secretaría de Economía, 2016).
La selección errónea de este material y el uso inadecuado de las plantas producidas repercuten en la supervivencia, producción y rendimiento de las reforestaciones o plantaciones forestales comerciales (PFC) (Alba et al., 2005). Por otro lado, al introducir plantas forestales en áreas bajo condiciones ambientales completamente diferentes a las requeridas por cada especie pueden convertirse en vectores de plagas y enfermedades (Vanegas, 2016), alterar relaciones tróficas y provocar pérdida de biodiversidad (Fernández-Pérez et al., 2013).
La Comisión Nacional Forestal (Conafor) desde 2001 ha promovido la reforestación para la restauración ecológica de áreas degradadas (Vanegas, 2016). Un factor clave para su éxito es el manejo y producción de semillas, ya que, el germoplasma de buena calidad asegura una mayor producción de plantas con características que garantizan una alta supervivencia en campo (Conafor, 2014).
Como parte de un esfuerzo por regular el uso y la movilización indiscriminada de germoplasma forestal, a finales del año 2016 entró en vigor la Norma Mexicana NMX-AA-169-SCFI-2016 (Secretaría de Economía, 2016), que establece las especificaciones técnicas que se deben cumplir para obtener la certificación durante el proceso de establecimiento y manejo de las unidades productoras de germoplasma forestal.
Una forma de contribuir a la regulación del movimiento del germoplasma y aumentar el porcentaje de supervivencia de las reforestaciones es modelar la distribución potencial de las especies forestales (Perosa et al., 2014). Esta técnica permite identificar áreas con las condiciones bióticas (por ejemplo, tipos de vegetación) y abióticas, como temperatura, pendiente y humedad favorables para la permanencia de un taxon; también puede utilizarse para identificar zonas aptas para introducir taxa con un alto valor ecológico o económico (Morales, 2012).
En años recientes, han surgido diversas herramientas que facilitan el modelado de la distribución potencial de especies, tales como: GARP, Bioclim y MaxEnt. La selección del algoritmo está en función del conjunto y complejidad de los datos con los que se cuente (Conabio, 2017). Los modelos mecanicistas basados en el nicho ecológico, como MaxEnt, usan hechos (datos) que permiten predecir la distribución potencial mediante algunos métodos estadísticos (Kearney, 2006).
Maxent ha demostrado tener una buena capacidad de predicción basada, únicamente, en datos de presencia (Elith et al., 2006; Navarro-Cerrillo et al., 2011). El modelo se fundamenta en el principio estadístico de máxima entropía (cercana a la uniforme) que facilita hacer predicciones con el uso de información incompleta, lo que representa una ventaja, ya que, para la mayoría de los taxones se carece de datos referentes a ausencias verdaderas (Phillips et al., 2006).
El estado de Chiapas es deficitario en la producción de semilla forestal; entre los años 2008 y 2014 se produjeron 1 053.21 kg, distribuidos entre especies de Pinus sp. (377.24 kg), Cedrela odorata L. (150.6 kg), Tabebuia rosea (Bertol.) Bertero ex A. DC. (20.73 kg) y Chamaedorea elegans Mart. (484.64 kg). En contraste, anualmente la Conafor produce alrededor de 19 millones de plantas para llevar a cabo los trabajos de reforestación en la entidad. Estos datos evidencian la necesidad de disponer de germoplasma forestal en cantidad y calidad suficientes para cubrir los requerimientos estatales. Por lo anterior, el objetivo del presente estudio consistió en delimitar áreas potenciales para el establecimiento de Unidades Productoras de Germoplasma Forestal (UPGF) de P. oocarpa Schiede ex Schltdl. y P. pseudostrobus Lindl. en el estado de Chiapas, mediante el uso de algoritmos MaxEnt.
Materiales y Métodos
Ubicación de área de estudio
El estado de Chiapas se localiza al sureste de la república mexicana y cuenta con una superficie de 73 272.3 km2 (Inegi, 2016). Su complejo relieve se enmarca en siete regiones fisiográficas: Llanura Costera del Pacífico, Sierra Madre de Chiapas, Depresión Central, Altiplanicie Central, Montañas del Oriente, Montañas del Norte y Llanura Costera del Golfo. Debido a sus condiciones topográficas, Chiapas es una de los entidades con mayor diversidad biológica (Conabio, 2013).
Con base en la clasificación de World Reference Base for Soil Resources (WRB), en el estado existen 19 unidades de suelo (Inegi, 2006).
En Chiapas predomina el clima cálido húmedo sobre 39.4 % del territorio, seguido del cálido subhúmedo con 34.9 %, el semicálido húmedo con 14.2 % y en menor porcentaje el templado subhúmedo y el templado húmedo con el 7.0 y 3.2 %, respectivamente (Inegi, 2008). La temperatura media anual varía según la región, de 18 °C en los Altos de Chiapas a 28 °C en la Llanura Costera del Pacífico. La precipitación total anual fluctúa entre 900 a 4 000 mm y las altitudes desde 0 hasta 4 000 m (Conabio, 2013).
Datos de presencia
La presencia de ambas especies se obtuvo mediante la consulta de bases de datos y plataformas de ejemplares depositados en el Herbario Nacional de México (MEXU) y la Red Mundial de Información sobre Biodiversidad (REMIB) elaborada por la Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (Conabio). Las bases fueron depuradas y se eliminaron aquellos registros que se ubicaron dentro de un área buffer de 200 m a la redonda de las zonas urbanas definidas en la capa vectorial del Marco Geoestadístico Nacional.
Para P. oocarpa se obtuvieron 220 registros y para P. pseudostrobus 52 (Figura 1); en ambos casos, se seleccionó, de forma aleatoria, 28 % del total para validar el modelo (Ibarra et al., 2012). Se realizó una depuración de datos con base en los parámetros máximos y mínimos para cada variable. El reducido número de registros empleado para la modelación, en este caso se justifica por la ausencia de datos completos para las especies en cuestión; sin embargo, el algoritmo puede ser eficaz, incluso cuando el número de sitios en los que se ha documentado la presencia es bastante bajo (Costa et al., 2010). Para ambos taxa, los pixeles ocupados por más de un registro se eliminaron.
Variables edafoclimáticas
Se utilizaron 25 variables; 17 derivadas de los valores mensuales de temperatura y precipitación obtenidas de la plataforma WorldClim y ocho de tipo topográfico, climático y edafológico.
Las variables utilizadas en la modelación fueron seleccionadas con base en los requerimientos ambientales más significativos para ambas especies, también se consideró la disponibilidad de información geoespacial. Las variables fueron: temperatura (Eguiluz, 1982; Sáenz-Romero et al., 2006), precipitación (Fierros et al., 1999; Sáenz-Romero et al., 2006), altitud (Perry, 1991; Sáenz-Romero et al., 2006; Viveros-Viveros et al., 2007), tipo de suelo (Eguiluz, 1982; Fierros et al., 1999), pH (Eguiluz, 1982; Rueda et al., 2006), clima (Fierros et al., 1999; Rzedowski, 2006) y textura del suelo (Fierros et al., 1999), las cuales se obtuvieron de la Serie IV de Inegi (Ed.2), a una escala de 1:250 000 (Inegi, 2006). Para este estudio, además de las variables edafoclimáticas, se incluyó el Índice de Vegetación de Diferencia Normalizada (NDVI) como un indicador del porcentaje de cobertura vegetal y del vigor de la vegetación (Kulloli y Kumar, 2014) (Cuadro 1).
Clave | Variable ambiental |
---|---|
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) |
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) |
Alt | Altitud (msnm) |
Clim | Clima (tipo) |
Eda | Suelos (tipo) |
pH | pH (H+) |
NDVI | Índice de vegetación diferencial normalizado |
Pp | Precipitación media anual (mm) |
Tem | Temperatura (°C) |
Tex | Textura (tipo) |
Dado que las capas presentaron diferentes tamaños de pixel, se estableció uno estándar de 15 m a través del remuestreo de las capas con la herramienta resample de ArcGis 10.4® (ESRI, 2017), mediante el método CUBIC.
El mapa de altitud se elaboró a partir del Modelo Digital de Elevación de Inegi (Inegi, 2016) con resolución espacial de 15 m. La cartografía de temperatura media, con el método de gradiente altitudinal (Fries et al., 2012), con base en el registro histórico de 30 años de temperatura y altitud de 173 estaciones del Servicio Meteorológico Nacional (SMN, 2010). El mapa de precipitación se obtuvo mediante la interpolación de los datos de la precipitación media anual de 173 estaciones del SMN con el método de Kriging (Delaney, 1999).
La capa de clima se derivó de los datos vectoriales de Inegi (Inegi, 2008), los tipos de suelo, la textura y pH de los datos vectoriales de los perfiles de suelos y el conjunto de datos de erosión del suelo a escala 1: 250 000 (Inegi, 2006).
El NDVI resultó de procesar siete escenas Landsat 8 con diferentes fechas (Pat/Row: 20/48 (09-02-2015), 20/49 (09-02-2015), 21/48 (02-03-2016), 21/49 (02-03-2016), 21/50 (01-12-2015), 22/48 (25-01-2016) y 22/49 (09-01-2016), con el cual se generó un mosaico que permitió cubrir el estado de Chiapas. El módulo Atmosc de Idrisi Selva se utilizó para hacer la corrección atmosférica de las imágenes, con el propósito de disminuir los efectos de nubosidad y transformar los valores de números digitales a reflectancia (Chávez, 1996). Para el cálculo del NDVI se aplicó la Ecuación 1.
Donde:
NIR = Banda 5
R = Banda 4
El NDVI se tomó de acuerdo a lo descrito por Huete et al. (1999) y utilizado por Kulloli y Kumar (2014) y Ruiz-Huanca et al. (2005), en el que los valores cercanos a 0 indican áreas con baja cobertura vegetal; por el contrario, valores cercanos a 1 denotan alta cobertura y buen vigor.
La transformación de formatos de las capas vectoriales, los recortes y la homogenización del tamaño de la cartografía, así como las interpolaciones se realizaron siguiendo el método de Kriging en el software ArcGIS 10.4® (ESRI, 2017).
Para evitar un sobreajuste de los modelos por multicolinealidad entre variables (Dorman et al., 2013), se aplicó un análisis de correlación de Pearson. En la selección de variables para el modelado final se consideraron aquéllas que tuvieron un coeficiente de 0.80 y -0.80 (Fuentes et al., 2016) y p<0.0001. A partir de eliminar las variables correlacionadas, se hizo la prueba Jackknife para conocer las que aportaron mayor información al modelo (Phillips et al., 2006).
El modelado se hizo con el software MaxEnt 3.3.3. Para definir el modelo se ejecutaron seis algoritmos de los más utilizados y recomendados por diversos autores (Elith et al., 2006; Plasencia-Vázquez et al., 2014); para cada caso se llevaron a cabo 50 réplicas:
Función Cumulative con un umbral de convergencia de 1.0E-5, 1000 iteraciones y método de remuestreo crossvalidate.
Función Cumulative con un umbral de convergencia de 1.0E-5, 500 iteraciones y remuestreo crossvalidate.
La función logística con un umbral de convergencia de 1.0E-5, 1000 iteraciones y método de remuestreo bootstrap.
La función logística con un umbral de convergencia de 1.0E-5, 500 iteraciones y método de remuestreo bootstrap.
Función logística con un umbral de convergencia de 1.0E-5, 500 iteraciones y método de remuestreo crossvalidate.
Función Cloglog con un umbral de convergencia de 1.0E-5, 500 iteraciones y método de remuestreo crossvalidate.
Para obtener el mapa de distribución potencial se utilizó la regla de umbral Equal training sensitivity and specificity, debido a que fue la que mejor delimitó el área de distribución potencial; Plasencia-Vázquez et al. (2014) obtuvieron buenos resultados al modelar dos especies de Psitácidos (Amazona xantholora (Gray, 1859) y A. oratrix (Ridgway, 1887); y Liu et al. (2005) lo clasifican como uno de los mejores para establecer los umbrales de presencia o ausencia de las especies modeladas.
La evaluación del modelo se realizó con el valor del Área Bajo la Curva (UAC). Se usó la prueba Jackknife para conocer las variables que aportaron mayor información al modelo (Phillips et al., 2006).
Con base en los resultados del modelo se identificaron áreas con altas probabilidades de ubicar zonas adecuadas para el establecimiento de UPGF. Se hicieron recorridos de campo con el propósito de validar los mapas y delimitar rodales aptos para la producción de semilla. Al interior de esos rodales se llevó a cabo un inventario forestal, con la finalidad de seleccionar y marcar los árboles tipo 1 y tipo 2 (Secretaría de Economía, 2016); además, dichos sitios se describieron en función de las variables mostradas en el Cuadro 1.
Resultados y Discusión
Al realizar las pruebas de correlación de Pearson entre las 25 variables, se observaron problemas en 20 con coeficientes de 0.80 y -0.80 y p<0.0001, por lo que se eliminaron. Para el caso de P. pseudostrobus las variables consideradas en el modelo fueron Bio2, Alt, Eda, pH y textura; mientras que, para P. oocarpa, Bio2, Bio14, Alt, NDVI y pp.
De las seis pruebas para la salida del modelo, se optó por utilizar el modelo logístico con un umbral de convergencia de 1.0E-5, 500 iteraciones y método de remuestreo crossvalidate, debido a que fue el que delimitó con más exactitud las áreas de presencia de las especies. De acuerdo a Phillips y Dudík (2008), la salida de tipo logística realiza una transformación de la tasa de ocurrencia relativa mediante la cual MaxEnt puede estimar la probabilidad de presencia del taxon. Por otro lado, este modelo ha sido validado por diversos autores: Norris (2014) llegó a buenos resultados al utilizarlo para Tapirus terrestris (Linnaeus, 1758) e Ibarra et al. (2012) al modelar el nicho ecológico de Microcystis sp., entre otros.
Araujo y Guisan (2006) clasifican la precisión de los modelos en cinco categorías, en función de los valores de AUC: de 0.50-0.60, se cataloga como insuficiente; 0.60-0.70 como pobre; 0.70-0.80 como valores promedio; 0.80-0.90 como bueno y 0.90-1, como excelente. Al respecto, Phillips et al. (2006) señalan que modelos con predicciones perfectas alcanzan valores de 1; sin embargo, cuando se emplean únicamente datos de presencia es común obtener valores de AUC inferiores a 1.
Los modelos generados para cada especie mostraron valores de AUC aceptables, 0.882 para P. oocarpa y 0.947 para P. pseudostrobus, lo cual demuestra que pueden utilizarse para predecir la distribución y presencia de ambas especies, con un nivel de confiabilidad alto (Phillips et al., 2006). En otros estudios se han registrado valores de entre 0.7 y 0.9 (Wan et al., 2015), por lo que los autores concluyen que alcanzaron resultados precisos con el uso de MaxEnt.
Las áreas potenciales para las dos especies se concentran, principalmente, en la región Sierra Madre de Chiapas y Altiplanicie Central (Figuras 2 y 3). Para P. oocarpa se delimitó una superficie de 874 695 ha, con un alto potencial para el establecimiento de UPGF; esta especie de amplia distribución, se localiza en altitudes de 300 a 3 000 m, en suelos pobres y temperaturas entre 3 y 35 °C (Eguiluz, 1982; Perry, 1991). Para P. pseudostrobus el área fue de 478 493 ha, su intervalo de distribución altitudinal es más restringido, de 1 600 a 3 200 m y se desarrolla bajo temperaturas mínimas de -9 °C y máximas de 40 °C (Fierros, 1999). Sin embargo, la escasa superficie predicha por el modelo se debe, fundamentalmente, a que P. oocarpa prospera con mayor éxito en altitudes de 900 a 2 600 m, temperaturas mínimas de 12 °C y máximas de 22 °C, en suelos con profundidades de moderadas a profundas (Eguiluz, 1982); P. pseudostrobus es favorecida en altitudes de 1 500 a 2 400 m y con temperaturas de 14 a 20 °C (Fierros, 1999).
Respecto a las variables utilizadas (Cuadro 2), la altitud fue la que aportó mayor porcentaje al entrenamiento del modelo de ambos taxa, lo que concuerda con Schumann et al. (2016), quienes al modelar la distribución de 14 especies, la altitud resultó ser la variable más importante para 11 de ellas. En otras investigaciones con diferentes taxones del género Pinus se ha determinado una correlación positiva entre la altitud y diferentes variables morfológicas (Sáenz-Romero et al., 2006; Sáenz-Romero et al., 2012; Viveros et al., 2013); por lo que es de esperarse que dicha variable contribuya de manera importante con información para ambos modelos.
Variable | Contribución | |
---|---|---|
P. oocarpa (%) | P. pseudostrobus (%) | |
Alt | 84.5 | 97.3 |
Bio2 | 6.9 | 0 |
Bio14 | 6.2 | - |
Eda | - | 0.3 |
pH | - | 2.4 |
NDVI | 0.6 | - |
pp | 1.8 | - |
tex | - | 0 |
El pH del suelo influye en la disponibilidad de la mayoría de los nutrientes (Alcántar y Trejo, 2010); sin embargo, es una variable poco utilizada en el modelado de nicho ecológico. No obstante, para P. pseudostrobus, el pH se incorporó en la modelación, ya que, por lo general, su distribución está restringida a suelos ácidos, con valores entre 4.5 a 6.5 (Eguiluz, 1982; Rueda et al., 2006), lo cual coincide con las limitantes de distribución de otras coníferas (Pérez et al., 2014).
A pesar de que la altitud, la temperatura, la precipitación y el tipo de clima son variables con frecuencia consideradas en el modelado de la distribución potencial de especies (Schumann et al., 2016; Qin et al., 2017), en este trabajo solo fueron utilizadas la Bio2, Bio14 y la pp, ya que tuvieron un bajo peso en la construcción del modelo, debido a la alta correlación entre ellas (García, 2004; Fries et al., 2012).
La ruta particular seguida por MaxEnt para obtener la solución óptima del modelo, origina resultados distintos, ya que un algoritmo diferente podría generar la misma solución por medio de una ruta distinta, lo cual conduciría a valores de contribución porcentual distintos (Phillips et al., 2006). En este caso en particular, MaxEnt consideró más importante la altitud en comparación con el resto de variables usadas para la modelación.
Al comparar los resultados del Cuadro 2, se confirma lo expuesto en el párrafo anterior, ya que la altitud por si sola tiene un peso relevante en la presencia o ausencia de las dos especies (>80 %). Resultado afín a lo registrado por Cruz et al. (2014) en el entendido de que dicha variable tiene un alto porcentaje (>85 %) de participación en el modelo de tres especies forestales. Para P. oocarpa, la Bio2 y la Bio14 mostraron una participación importante en el modelo, lo que es similar para Catopheria chiapensis A. Gray ex Benth., Quercus martinezii H. Mull., Telanthopora grandifolia (Less.) H. Rob. & Brettell y Viburnum acutifolium G. Bentham.
El NDVI reveló ser una variable con poca aportación porcentual al modelo; dicho comportamiento es atribuible a que ese índice de vegetación refleja niveles de verdor de la vegetación, al respecto, los valores altos de NDVI para el área de estudio, se obtuvieron de pastizal, vegetación secundaria de bosque de pino-encino, selva mediana, bosques de pino-encino y selva alta perennifolia. Otra desventaja del NDVI es la saturación ante valores de índice de área foliar superiores a 2 (Huete et al., 1999; Ruiz-Huanca et al., 2005); por ello, es recomendable usarlo en zonas áridas y semiáridas en donde Schumann et al. (2016) consignaron buena respuesta al modelar la distribución de especies en dichos ecosistemas.
La textura del suelo fue la variable menos relevante, lo cual es atribuible al nivel de descripción de la información contenida en la cartografía digital utilizada, ya que clasifica la textura en partículas gruesa, media y fina. Fierros et al. (1999) argumentan que P. pseudostrobus se puede establecer en suelos con textura franco-arenosa, franco-arcillo-arenosa, arcillosa y arcillo-arenosa; mientras que, P. oocarpa en suelos con textura arenosa, migajón-arenosa, areno-arcillosas con buen drenaje (Eguiluz, 1982); lo anterior ratifica que ambas especies se distribuyen de forma indistinta en suelos con texturas de fina a gruesa.
Los resultados del modelo y los recorridos de campo permitieron establecer y registrar dos Unidades Productoras de Germoplasma Forestal: la UPGF Juznajab (N-07-019-JUZ-001/17) con las dos especies de interés, y la UPGF Coapilla (N-07-018-COA-002/17) con P. oocarpa (Figura 4).
En Juznajab, 2 % del área fue catalogado como no apto para P. oocarpa y 98 % se ubicó en la categoría de potencial medio; para P. pseudostrobus 100 % del área se catalogó con potencial alto, que concuerda con las condiciones observadas en campo, lo que se atribuye al estado actual de la vegetación, por ser un área bajo manejo forestal con una cobertura de copa promedio de 58 %. Asimismo, se confirmó que la precipitación presente en el sitio (1 350 mm) es inferior al óptimo requerido por la especie.
Conclusiones
Los resultados del modelo respaldan la ubicación de las zonas aptas para el establecimiento de UPGF, ya que disminuyen los tiempos y costos en la búsqueda de áreas con condiciones óptimas. Además, se sugiere que para generar modelos confiables a nivel de variedad es de vital importancia contar con un gran número de puntos de presencia que correspondan específicamente a la variedad de interés. De igual forma, se requiere caracterizar a detalle las condiciones edáficas y representar espacialmente la información a una mayor resolución espacial que los utilizados en este trabajo.