Área de estudio

El estudio se llevó a cabo en el ejido Emiliano Zapata, ubicado en Chignahuapan, Puebla, México (19°42’30” latitud norte, -98°2’30” longitud oeste). El clima es templado subhúmedo con lluvias en verano, caracterizado por una precipitación anual de 665 mm y temperatura media de 13.4 °C. El área se caracteriza por presentar bosques templados dominados por coníferas como Pinus patula Schl. et Cham, Pinus ayacahuite Ehren. ex Schltdl. y Abies religiosa (Kunth) Schltdl. et Chamen (^{Lazcano-Hernández, 2006}). Los suelos se clasifican como luvisoles, los cuales poseen una textura franca a franco arenosa, alta porosidad, un pH ligeramente ácido y un contenido de materia orgánica que oscila entre 6% y 7.2% (^{Galicia et al., 2015}).

Selección de unidades de muestreo

A partir de la revisión cartográfica del programa de manejo forestal maderable (PMFM) vigente del ejido (^{Unión de Ejidos de Producción, Explotación, Comercialización, Industrialización Agropecuaria y Forestal de la Sierra Norte de Puebla [Uepeciafsnp], s/f}), se diseñó un esquema de muestreo sistemático utilizando unidades circulares de 1000 m², separadas por una distancia equidistante de 100 m. En total, se seleccionaron 34 sitios en áreas de aprovechamiento forestal dominadas por Pinus patula (Fig. 1). Estos sitios, con características ecológicas homogéneas, se distribuyeron en 11 rodales con dos métodos de manejo (MDS y MMOBI) y uno de referencia (Referencia), el cual no ha sido intervenido en las últimas décadas (Tabla 1). En total se seleccionaron 16, 13 y 5 sitios circulares, para MMOBI, MDS y Referencia, respectivamente. Para MMOBI, se evaluó el tratamiento silvícola de cortas de selección, mientras que para MDS, se consideraron cortas de aclareo, liberación y regeneración. Dado que no todos los tratamientos se aplicaron en la misma anualidad, esta se subdividió en rodales.

La letra R indica los rodales descritos en la tabla 1.

Figura 1 Ubicación de los sitios circulares de 1000 m² para el muestreo dendrocronológico dentro del ejido Emiliano Zapata, Chignahuapan, Puebla. 

Muestreo y fechado dendrocronológico

En cada sitio de 1000 m², considerando todos los individuos registrados, se seleccionaron cuatro árboles representativos con base en su altura y diámetro. A la altura del diámetro normal (DN), se extrajeron dos núcleos de crecimiento por árbol utilizando un taladro Pressler de 5.1 mm de diámetro. Dado que en los sitios bajo manejo por MDS, la estructura fue homogénea; en los sitios bajo MMOBI y de referencia, la selección se enfocó en arboles maduros representativos de las condiciones del sitio, incluyendo al menos un individuo de categoría juvenil (DN mínimo de 20 cm) por sitio de dimensiones fijas.

El tamaño de muestra se integró por 136 árboles de P. patula, evitando aquellos enfermos, plagados o con problemas de crecimiento. Las muestras dendrocronológicas fueron secadas, montadas, y pulidas con diversos grados de lija en el Laboratorio Nacional de Dendrocronología del Centro Nacional de Investigación Disciplinaria en Relación Agua, Suelo, Planta, Atmósfera [Cenid Raspa], perteneciente al Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias [Inifap), ubicado en Gómez Palacio, Durango. La medición del grosor de los anillos de crecimiento (anillo total) se realizó con un sistema de medición Velmex (resolución 0.001 mm). La calidad del fechado fue evaluada con el software COFECHA, al correlacionar sucesivamente segmentos de 50 años en los crecimientos radiales con traslapes de 25 años, y una intercorrelación entre series mínima de 0.3281 (p < 0.01) (^{Holmes, 1983}).

Determinación del incremento en área basal e índice de ancho de anillo

Posterior al fechado, las series de crecimiento se están-darizaron utilizando una curva tipo Spline para calcular índices de ancho de anillo (IAA), utilizando el paquete dplR en R (^{Bunn, 2008}). El proceso de estandarización es necesario para eliminar efectos de edad, geométricos, competencia, entre otros conocidos como ruido, que no son de interés para el investigador. Para determinar las tasas de crecimiento radial, se calculó el incremento en área basal (IAB). El IAB es indicador temporal de la productividad forestal (cm² año^-1), calculado a nivel árbol, de acuerdo con la siguiente formula:

donde:

R = radio del fuste (cm)

n = año de formación del anillo

Finalmente, los valores de IAB fueron agrupados a escala de árbol, sitio, rodal, tratamiento silvícola o método de manejo, de acuerdo con el análisis estadístico realizado.

Información climática y su relación con el crecimiento radial

Se utilizó información del Sistema de Información Hidrológica proporcionado por la Comisión Nacional del Agua (Conagua) para obtener datos diarios del área de estudio en el periodo de 1973 a 2023. Los datos se obtuvieron de cinco estaciones en operación con condiciones climáticas similares, ubicadas a una distancia máxima de 28 km de los sitios de muestreo. Las estaciones utilizadas fueron 29032 (Tlaxco, 96.7% datos completos), 29052 (El Rosario, 95.7% datos completos), 29165 (El Epazote, 86.5% datos completos), 21021 (Tetela de Ocampo, 97.5% datos completos) y 21140 (Chignahuapan, 78.6% datos completos). Los datos diarios recopilados incluyeron temperatura (máxima, mínima y promedio), precipitación (milímetro por día) y evaporación (milímetro por día) de cada estación. Para garantizar la calidad de las series de datos meteorológicos, se empleó el paquete Climatol (^{Guijarro, 2019}) en R. Este paquete permite la homogenización y el análisis de las series temporales climáticas, además de la corrección de errores e interpolación de datos faltantes (^{Addou et al., 2022}).

Para evaluar la influencia de las variables climáticas registradas diariamente en el crecimiento radial de P. patula, se utilizó el paquete dendrotools de R (^{Jevšenak y Levanič, 2018}). Este paquete permite calcular una matriz de correlaciones considerando una ventana móvil de datos diarios. A diferencia de otros métodos, el paquete está mejor adaptado al crecimiento radial continuo (datos diarios vs datos mensuales), lo que resulta en una selección más precisa de las relaciones entre clima y crecimiento. Para el ancho de la ventana, se establecieron 21 días y 240 días como límites inferior y superior, respectivamente, lo que significa que el tamaño mínimo de análisis de correlación fue inferior a un mes, mientras que el máximo fue equivalente a 8 meses. Las correlaciones comenzaron desde el día del año (DOY) 240 del año anterior (t^-1) (finales de agosto) hasta el DOY 300 del año actual (t) (finales de octubre). Se utilizó la función de agregación (suma) para la precipitación y evaporación; mientras que para la temperatura se utilizó el promedio. Esto quiere decir que, por ejemplo, para la ventana mínima, se consideró precipitación acumulada de 21 días, mientras que para temperatura se utilizó el promedio. Para mejorar la precisión del análisis y generar intervalos de confianza en el análisis de correlación, el análisis se replicó por método de manejo y tratamiento silvícola considerando un bootstrap de 100 veces y un α = 0.05.

Diferencia en el crecimiento radial de P. patula entre métodos y tratamientos silvícolas

Para evaluar diferencias entre métodos de manejo para el IAB, se construyó un modelo mixto, considerando una estructura de tres niveles. El primer nivel considera las mediciones repetidas en el tiempo por cada árbol de acuerdo con el método de manejo (e.g. IAB), seguida de una agrupación espacial (e.g. árboles anidados a sitio circular, árboles anidados a rodal y árboles anidados en sitio circular dentro de un rodal) y finalmente covariables a escala de sitio (área basal, densidad arbórea, elevación, pendiente, clima). Para evaluar la multicolinealidad de los factores, se evaluó el factor de inflación de la varianza (VIF < 5). Dado que IAB no se distribuyó de manera normal, los datos de IAB se transformaron a logaritmo base natural con el paquete bestNormalize en R (^{Peterson, 2021}). Para modelar la autocorrelación entre mediciones repetidas en el tiempo, se evaluó una estructura AR (covarianza) con diferentes estructuras (potencia, constante o combinación de ambas) (^{Everitt, 2005}). Así mismo, dado que los valores de IAB no suelen comportarse de manera lineal, se incluyó una función Spline de orden 2 para modelar su tendencia (Ecuación 2). El mejor modelo fue seleccionado consi-derado el menor criterio de información de Akaike (AIC) y criterio de información bayesiano (BIC). El análisis estadístico fue realizado en R 4.3.2 (^{R Core Team, 2023}) utilizando el paquete nlme (^{Pinheiro et al., 2018}).

donde:

log(IAB) _ijk =incremento en área basal transformado logarítmicamente para el árbol 𝑘, en el sitio 𝑗, en el tiempo 𝑖

𝛽₁, … , 𝛽₇ = coeficientes asociados a las variables independientes

bs(𝐴ñ𝑜 ) = función Spline de orden 2 para modelar el comportamiento temporal

𝑏_0𝑗= efecto aleatorio del sitio 𝑗 (intercepto aleatorio)

𝑏 _𝑘j ∙𝐴ñ𝑜= pendiente aleatoria del año dentro del sitio 𝑗 y el árbol 𝑘

𝜀 _𝑖jk = error aleatorio asociado a cada observación, modelado con una estructura de correlación AR(1) (dependencia temporal entre los años)

Edad del arbolado

Un total de 131 árboles de P. patula fueron datados correctamente, lo que representa 96% del total de las muestras obtenidas (136). La edad promedio del arbolado fue de

36.35 años, sin embargo, la edad máxima fue de 108 años, perteneciente a un árbol ubicado en un sitio de referencia (1916-2023) (Fig. 2). A nivel método, el arbolado más joven se encontró en MDS ([27.2 ± 2.19] años) seguido por MMOBI ([42.5 ± 1.0] años). Para los tratamientos silvícolas dentro de MDS, las cortas de regeneración tuvieron el arbolado más joven [(19.4 ± 1.08] años), mientras que el aclareo, el más longevo ([38.4 ± 4.74] años). A escala general, la intercorrelación entre todas las series fue de 0.48, superior al umbral mínimo requerido en estudios dendrocronológicos para un fechado correcto (^{Speer, 2010}).

Figura 2 Diagrama de edades en el arbolado de Pinus patula bajo manejo forestal en el ejido Emiliano Zapata, Chignahuapan, Puebla. 

^{Ortega-Rodríguez (2014)}, en Perú, señala valores de intercorrelación entre series de 0.69 para la misma especie, lo cual es superior a lo encontrado en este estudio. Esta diferencia posiblemente se deba al hecho de ser plan-taciones con una menor edad (25 años), lo que asegura una alta sincronía entre patrones de crecimiento. Es importante destacar que los valores de intercorrelación indican que, aunque este estudio no se centró en la reconstrucción climática, los árboles presentaron una señal común asociada al clima. Este hallazgo no es frecuente en bosques manejados, donde la densidad arbórea, la dinámica forestal y, en general, las condiciones del sitio suelen enmascarar su asociación con el clima (anillos de crecimiento complacientes vs sensitivos) (^{Fritts, 1976}). En este caso, P. patula, más allá de ser una especie de rápido crecimiento con interés comercial, muestra una clara sensibilidad al clima.

El potencial para utilizar P. patula en estudios dendrocronológicos ha sido documentado por ^{Quezada-Rivera (2013)} y ^{Belay (2016)}, mientras que ^{Rodríguez-Ortíz et al. (2011)} y ^{Ortega-Rodríguez (2014)} utilizaron información de sus anillos de crecimiento para evaluar efectos de aclareos en el crecimiento radial, lo que resalta la importancia de la dendrocronología en la evaluación del incremento radial de la especie.

Efecto del clima en el crecimiento

Todas las variables climáticas evaluadas resultaron asociadas con en el crecimiento radial de P. patula. De manera general, la precipitación fue la variable con mayor influencia positiva (r = 0.73) mientras que la temperatura mínima fue la de menor influencia (r = 0.37). Esto concuerda con ^{van Zonneveld et al. (2009)} y ^{Ortega-Rodríguez (2014)}, quienes señalan que la precipitación y la temperatura máximas del mes más cálido son las variables climáticas que mejor definen el hábitat natural y crecimiento de P. patula. A nivel del método de manejo, MMOBI fue, en su mayoría, el método donde se identificó la mayor señal climática, seguido por el de referencia, y MDS el de menor señal (Tabla 2). Este resultado podría explicarse por las bajas intensidades de corta realizadas en estos rodales, así como las condiciones topográficas, caracterizadas por mayor elevación y suelos más someros y pedregosos observados en campo, que favorecen la formación de anillos más sensitivos (^{Fritts, 1976}). Esto apoya parcialmente la hipótesis de este trabajo, ya que el arbolado con mayor intensidad de manejo forestal mostró los valores de correlación más baja con el clima.

En MDS, se estimó una densidad arbórea promedio de 502 ha^-1, mientras que en MMOBI fue de 443 ha^-1, pero el comportamiento del volumen maderable fue diferente, de tal manera que en MDS, se registraron 248 m3 ha^-1 con un área basal de 25 m² ha^-1, en tanto que, en MMOBI el volumen promedio fue 328 m3 ha^-1 con 29 m² ha^-1. Entonces, en MDS, el volumen promedio por árbol fue de 0.49 m3 y, en MMOBI, de 0.74 m3, es decir, en este, en términos de parámetros forestales agregados, la ocupación del sitio resultó mayor que en el primero. Aun así, conforme a la guía de densidad de P. patula, propuesta por ^{Santiago-García et al. (2013)}, en ambos casos, los árboles se encuentran en el espacio de crecimiento libre sin mortalidad. Reyes-Cortés et al. (2020) sugieren que, en plantaciones de P. patula, el espacio disponible para el crecimiento y la competencia entre individuos afectan tanto el crecimiento radial como su relación con factores ambientales.

Para precipitación, el valor más alto de correlación se encontró en MMOBI (r = 0.73), representando el periodo de lluvias acumulado de enero 18 a mayo 28, mientras que, en MDS, el valor fue menor (r = 0.55), de noviembre 18 del año previo a mayo 20. Esto quiere decir que, la precipitación acumulada del invierno previo hasta la primavera (inicio de estación de lluvias) es fundamental para el crecimiento de P. patula en la zona de estudio (^{Quezada-Rivera, 2013}). Este periodo es importante no solo para P. patula sino también para otras especies forestales estudiadas en el norte de México (^{Cerano-Paredes et al., 2009}; ^{Stahle et al., 2016}). Sin embargo, el hecho de incorporar datos diarios en este estudio permitió determinar con mayor precisión el periodo de alta influencia en el crecimiento radial. Por ejemplo, ^{Quezada-Rivera (2013)} en la zona de Tlaxcala-Puebla determinó la precipitación estacional de enero a junio como la más relevante; mientras que en este estudio fue de enero 18 a mayo 28 (Tabla 2). Notablemente, este periodo representa poco menos de 20% de la precipitación total, lo que implica que, en el área de estudio, más que la precipitación de verano, donde se distribuye la mayor parte de la precipitación, la humedad al inicio de la estación de crecimiento es esencial para un mayor incremento radial anual para P. patula.

En general, la temperatura máxima tuvo una influencia negativa, es decir temperaturas más altas están asociadas a un menor crecimiento, especialmente en MMOBI (r = - 0.59) durante los meses de marzo a mayo, periodo donde se presentan los valores más cálidos en el año. La temperatura mínima presentó resultados mixtos, ya que mientras para MDS y MMOBI se encontró un efecto negativo, r = -0.46 y r = -0.34; para Referencia, fue positivo (r = 0.47) (Tabla 2). Sin embargo, esto depende de la estación del año, ya que, por ejemplo, el efecto positivo está asociado con el inicio de la primavera, mientras que el efecto negativo es en verano, donde temperaturas más frescas están asociados a mayores crecimientos. La temperatura promedio, en todos los casos, mostró una correlación negativa, resultando la regeneración como el tratamiento con mayor impacto, cuyo valor fue de -0.55, en el mes de febrero.

Finalmente, la demanda evaporativa ejerce un importante control negativo sobre el crecimiento de P. patula (r = -0.67), especialmente desde finales de marzo hasta junio, que suelen ser los meses más calurosos y con menor humedad en el ambiente. Por ejemplo, para MMOBI, la demanda evaporativa del periodo de marzo 31 a junio 30 disminuyó notablemente el crecimiento radial de la especie. Resulta interesante que dentro de los tratamientos silvícolas del MDS; el arbolado localizado en las cortas de regeneración tuvo, para todas las variables climáticas, la señal correlativa más intensa que el determinado en las cortas de liberación, o aclareos; de esta manera, para el caso de la precipitación, la correlación de 0.66, resultó ser la segunda más importante del comportamiento general (Tabla 2).

Incremento en área basal y su modelación estadística

En promedio, el IAB más alto se encontró en arbolado bajo MMOBI ([15.9 ± 0.44] cm² año^-1), seguido del correspondiente para MDS ([13.0 ± 0.50] cm² año^-1) y el valor más bajo de IAB se observó para el sitio de referencia ([11.0 ± 0.80] cm² año^-1). En el análisis entre tratamientos, en MDS, el arbolado perteneciente a las cortas de Regeneración tuvo la mayor productividad ([18.8 ± 1.12] cm² año^-1), mientras que Liberación la menor ([9.6 ± 0.64] cm² año^-1). Las principales diferencias estadísticas entre rodales (p < 0.05) fueron entre el grupo de Liberación y Aclareo de MDS (Rodal 2 al 3) y los de selección de MMOBI (Rodal 9 al 11); así como los rodales 6 y 7 (Fig. 3). Notablemente, existió una alta sincronía en picos máximos y mínimos, por ejemplo, durante el año 2014, la mayor parte de los sitios mostró valores máximos de productividad, posiblemente asociado al tratamiento aplicado y las condiciones ambientales atípicas. La precipitación en ese año (910 mm) fue sustancialmente mayor en comparación con años anteriores (585 mm), lo que resultó en un incremento significativo en la disponibilidad de agua para los árboles. Los datos de precipitaciones registradas en las cinco estaciones evaluadas refuerzan esta hipótesis; entonces, la combinación de condiciones climáticas favorables de un año previo con las intervenciones de manejo forestal, parecen haber creado un ambiente propicio para el incremento sustancial de la productividad de los árboles.

Área gris representa el error estándar entre árboles de un mismo rodal. Línea punteada es el promedio de IAB por rodal.

Figura 3 Series temporales de incremento en área basal (IAB) dividido por método de manejo (Método de Desarrollo Silvícola- MDS, Método Mexicano de Ordenación de Bosques Irregulares - MMOBI y referencia (Testigo). 

Dentro de la agrupación espacial evaluada en la modelación estadística del IAB (e.g., anidación), se encontró que los mejores modelos correspondieron a árboles anidados dentro de sitio circular y árboles anidados dentro rodal. Sin embargo, el modelo de árboles anidados dentro de sitio mostró un menor AIC en comparación del resto, y, por tanto, fue seleccionado para el modelo final (efectos aleatorios). De igual manera, con base en este estadístico de ajuste (AIC), la estructura de varianza que mejor reflejó la correlación temporal para los datos fue del tipo potencia.

El tiempo (años), el método de manejo, las características intrínsecas del rodal (área basal y densidad arbórea) y las climáticas (precipitación y evaporación) fueron las variables más relevantes para explicar el comportamiento de IAB (efectos fijos). El método de manejo fue un factor importante con diferencias estadísticas entre todos los métodos de manejo (p < 0.001, Tabla 3). El área basal tuvo un efecto positivo en el crecimiento radial (0.354, p < 0.001), mientras que la densidad arbórea, negativo (-0.001, p = 0.012). Por tanto, aquellos rodales con un gran número de árboles están asociados con valores pequeños de IAB. La remoción de arbolado derivado del manejo forestal suele estimular el crecimiento de la masa forestal residual a escala de rodal; sin embargo, a escala individual, algunos árboles pueden mostrar una respuesta contraria, posiblemente asociada a la ruptura de interacción con otros individuos. Esta información no podría obtenerse sin un análisis dendrocronológico detallado (^{Trujillo-Martínez et al., 2020}).

Notablemente, años con alta precipitación influyen positivamente en el crecimiento radial de P. patula, mientras en años secos y de alta demanda evaporativa ocurre un reducido crecimiento (Tabla 3). También es importante mencionar que, en su modelación, los efectos fijos explican alrededor de 50% de varianza, mientras que los aleatorios una proporción del 25%; es decir, existe una variabilidad en el crecimiento de los árboles que es difícil de capturar, al menos con las variables seleccionadas en este estudio.

De la figura 4a se destaca que el arbolado bajo MMOBI presentó una asíntota de máximo crecimiento alrededor de 2010 (alrededor de 30 cm² año^-1), mientras que para MDS, este punto aún no se ha alcanzado; ya que los árboles tanto en regeneración como en liberación aún están en desarrollo. La asíntota observada en MMOBI es semejante a la registrado por autores como ^{Rodríguez-Ortíz et al. (2011)} y ^{Santiago-García et al. (2016)}, quienes señalan que los árboles de P. patula llegan a un punto de inflexión en su crecimiento entre los 20 años y los 30 años. Así mismo, el arbolado de referencia, si bien ha mantenido su ritmo de crecimiento, está notablemente por debajo de lo observado con aquellos que presentan un manejo silvícola. La alta capacidad de crecimiento de P. patula ha sido documentada con anterioridad; por ejemplo, ^{Figueroa-Navarro et al. (2019)} destacan la rápida recuperación de biomasa en rodales con P. patula manejados con MDS, lo que denota su alta capacidad de resiliencia.

Figura 4 (a) Incremento en área basal predicho con el modelo seleccionado en la tabla 3 para los diferentes métodos de manejo. 

A diferencia de otros estudios donde se modela el crecimiento radial con curvas como Chapman-Richards o Weibull, el uso de modelos mixtos permite incorporar la autocorrelación temporal y espacial de los datos, la variabilidad inherente a cada árbol analizado y entre grupos de árboles, así como el efecto de covariables de interés (Tabla 3 y Fig. 4). Las variables consideradas dentro del modelo final presentaron un VIF < 4, por lo que no existió problema de multicolinealidad. Así mismo, los interceptos a escala de sitio, pendientes de los árboles y los residuales del modelo tendieron a la normalidad, lo que asegura una correcta interpretación del modelo final (Fig. 4b, c y d).

La importancia de este tipo de estudios en el manejo forestal, incorporando información dendrocronológica, es la de tener indicadores confiables acerca del efecto de diferentes tratamientos silvícolas en conjunto con el efecto del clima sobre los bosques, considerando a la par una perspectiva de cambio climático y resiliencia (^{Gómez-Mendoza y Arriaga, 2007}; ^{Correa-Díaz et al., 2023}). Por ejemplo, en los últimos años se han detectado, con mayor frecuencia, evidencia de procesos de declinación forestal, asociados a mortalidad masiva o ataque de plagas y presencia de enfermedades, donde P. patula no es la excepción (^{Sáenz-Romero et al., 2020}).