Abstract

The accelerating effects of climate change on Mexico’s temperate forests highlight the need for a precise and comprehensive assessment of ecosystem vulnerability. This study quantified how carbon density in aboveground biomass (cdAGB) in key forest species responds to bioclimatic gradients. Twelve taxa from the genera Pinus and Quercus, comprising six pine and six oak species, were assessed using data from the National Forest and Soil Inventory (2015-2020) together with 19 climatic variables. The analytical framework integrated Bayesian correlation analysis, Random Forest, and piecewise regression to identify inflection points or thresholds (ψ) characterizing the relationship between cdAGB and the bioclimatic predictors. Critical thresholds were identified in nine of the twelve species. Temperature-derived predictors exhibited the highest explanatory power. Davies test results showed that cdAGB in five of the studied species (three Pinus and two Quercus) exhibited significant bioclimatic thresholds characterized by abrupt declines within the central range of the predictor gradient (30th - 60th percentiles) Quercus showed a higher proportion of vulnerable sites (>60%) than Pinus (30.3% - 40.8%), and vulnerability probability followed complex, nonlinear spatial patterns structured mainly by latitude and longitude rather than elevation. Collectively, the results confirm a non-linear climate cdAGB relationship and underscore the importance of detecting ecological thresholds to refine projections of forest responses and optimize conservation strategies.

Keywords

Aboveground Biomass, Bioclimatic, Carbon Density, Thresholds, Vulnerability

Share and Cite:

1. Introduction

Worldwide, the total area under forest cover amounts to just over 4 billion hectares, representing approximately 31% of the Earth’s terrestrial surface. In this context, the Russian Federation, Brazil, Canada, the United States of America, and China possess the largest forest resources, collectively accounting for more than half of the global forest area (FAO, 2022). These ecosystems store an estimated 650 billion tonnes of carbon (≈650 Gt C), allocated as 45% in soils, 44% in biomass, and 11% in other carbon pools such as deadwood and leaf litter (FAO, 2010). Forests play a pivotal role in mitigating climate change by functioning as key regulators of the global climate system (Harris et al., 2021).

In Mexico, temperate forest ecosystems cover approximately 21% of the country’s land area and support over 7000 species, representing nearly 25% of the nation’s documented flora (Rzedowski, 1991). These ecosystems have undergone substantial biodiversity loss, with an estimated 25% of the original forest area having been converted to agricultural or livestock uses (Gómez & Arriaga, 2007). Moreover, Mexico harbors 43 of the 110 globally recognized pine species (Pinus spp.) and 161 of the 531 oak species (Quercus spp.) worldwide, exhibiting endemism levels of 55% in pines and 21% in oaks (Gómez & Arriaga, 2007; Sánchez, 2008).

Climate change is a global challenge driven by both natural and anthropogenic factors (Díaz, 2012). It manifests through species extinctions, pollution, natural disasters, and shifts in temperature and precipitation regimes (Barradas et al., 2011). According to Conde (2007), Mexico has lost up to 6.3 million hectares of forest in just over a century, ranking the country second in total forest loss across Latin America. Temperate forests worldwide are exhibiting complex responses to climate change, including variations in productivity associated with drought-induced stress resulting from rising temperatures (Pan et al., 2013).

According to the IPCC Sixth Assessment Report (IPCC, 2021), projections for the end of the century (2081-2100) indicate a best estimate of temperature increase ranging from 1.4˚C (SSP1-1.9) to 4.4˚C (SSP5-8.5) relative to the pre-industrial period (1850-1900). In addition, global land precipitation is projected to increase by between 0% and 13% compared with the 1995-2014 period, depending on the emission scenario. These conditions heighten the risk of extinction for vulnerable species (Watson et al., 2002; Wang et al., 2019). Temperature and precipitation are among the most influential climatic factors governing plant physiological processes such as photosynthesis, water and nutrient uptake, respiration, and leaf development (Wang et al., 2019; Ali et al., 2020; Bennett et al., 2020). Temperature regulates the rate of CO₂ assimilation and carbon losses through respiration, whereas precipitation influences stomatal conductance and water availability (Han et al., 2012).

The aboveground live biomass carbon density (cdAGB) refers to the carbon stock per unit area contained in living trees, encompassing stems, branches, leaves, and seeds (Liu et al., 2014; Penman et al., 2018). It represents a key variable in the assessment of forest ecosystem services. Quantifying cdAGB enables the estimation of carbon storage capacity and, consequently, the evaluation of the role of forests in climate change mitigation. According to estimates reported by Rodríguez et al. (2016), the amount of carbon stored in aboveground live biomass across Mexico is equivalent to 1.69 gigatonnes of carbon (Gt C).

According to Conde (2007), Mexico’s forests rank among the ecosystems most vulnerable to climate change. Consequently, it is essential to examine the relationship between the aboveground live biomass carbon density (cdAGB) of temperate-forest species and climatic variables to determine their vulnerability under future climate-change scenarios. Vulnerability is defined as the degree to which a system is susceptible to, and has limited capacity to cope with, the adverse effects of climate change; it is determined by exposure, sensitivity, and adaptive capacity (IPCC, 2007). Vulnerability assessments enable the identification of species at greatest risk of population decline, range shifts, or local extinction (Foden et al., 2013).

The objective of this study was to determine the relationship between carbon density in aboveground biomass and bioclimatic variables in species of the genera Pinus and Quercus using piecewise regression, and to identify climatic vulnerability thresholds. We analyzed 12 species (six Pinus species and six Quercus species) employing data from the National Forest and Soils Inventory (2015-2020). This research contributes by identifying cdAGB thresholds that may be used to prioritize conservation measures and adaptive forest-management actions in response to climate change.

2. Materials and Methods

2.1. Description of the Study Area and Species

In Mexico, coniferous forests occur from 150 to 4000 meters above sea level (m a.s.l.), predominantly in mountainous regions where annual temperatures range from 6˚C to 28˚C and annual precipitation varies between 350 and more than 1000 mm (Rzedowski & Huerta, 2006). According to Galicia et al. (2016), these forests are found on seven principal soil types—Leptosol, Regosol, Luvisol, Phaeozem, Cambisol, Umbrisol, and Andosol—which have developed from igneous, sedimentary, and metamorphic parent materials.

For this study, twelve species were selected—six pines and six oaks: Pinus devoniana Lindl. (hereafter referred to as Pdev), Pinus douglasiana Martínez (Pdou), Pinus leiophylla Schiede ex Schltdl. & Cham. (Plei), Pinus oocarpa Schiede ex Schltdl. (Pooc), Pinus patula Schiede ex Schltdl. & Cham. (Ppat), Pinus pseudostrobus Lindl. (Ppse), Quercus crassifolia Bonpl. (Qcra), Quercus laurina Humb. & Bonpl. (Qlau), Quercus magnoliifolia Née (Qmag), Quercus radiata Trel. (Qrad), Quercus resinosa Liebm. (Qres), and Quercus rugosa Née (Qrug).

2.2. Data

The original data from the National Forest and Soil Inventory of Mexico (INFyS, 2015-2020) were preprocessed to produce a database structured by species and plot (ID). Each record contains geographic coordinates (longitude, latitude; datum WGS84) and the number of subplots sampled (sites). The variables considered were: n (number of individuals measured), DBH (mean diameter at breast height measured at 1.30 m; cm), HT (mean total height; m), BA (basal area; m² ha⁻¹), AGB (total aboveground biomass; Mg ha⁻¹), and cdAGB (carbon density in aboveground biomass; Mg C ha⁻¹).

The 19 bioclimatic variables were obtained from the WorldClim database, version 2.0 (Hijmans, 2024), at a spatial resolution of 30 arc-seconds for the period 1970-2000. Values for each of the 19 bioclimatic variables (Bios) at each site (ID) and for each species were extracted from the WorldClim raster layers using the terra package (Hijmans, 2024).

2.3. Bayesian Correlation and Climatic Predictors

To assess the correlation between cdAGB and each bioclimatic variable, and owing to the violation of the assumption of bivariate normality, the variables were rank-transformed and a Bayesian Spearman correlation coefficient (ρBayes) was estimated using the BayesFactor package (Morey & Rouder, 2024). This approach was chosen for its robustness to non-normal distributions and because Bayesian inference provides more flexible estimates and an explicit probabilistic framework for hypothesis comparison. In addition to the posterior median, we calculated credible intervals (CrI), the probability of direction (pd), the percentage in the region of practical equivalence (% in ROPE), and the Bayes factor (BF).

In this study, climatic sensitivity was defined as the magnitude and direction of the statistical response of cdAGB to bioclimatic gradients, with rank-based Bayesian correlations and segmented regression slopes serving as functional proxies. Vulnerability was defined functionally as a statistically significant reduction in cdAGB associated with adverse climatic conditions or bioclimatic thresholds, reflecting climate-driven constraints on carbon accumulation capacity rather than direct physiological failure, mortality, or demographic risk.

2.4. Response Variable

Predictor variables were selected using a two-step framework integrating multicollinearity control (VIF < 10) and Random Forest-derived importance metrics (Kuhn et al., 2008). The dataset was partitioned into training (80%) and testing (20%) subsets via stratified random sampling to preserve cdAGB distributions. Model performance and generalizability were evaluated using 5-fold cross-validation and independent testing, with predictive accuracy quantified by RMSE, R², and MAE. The relative importance of bioclimatic predictors for estimating cdAGB was determined by comparing model-internal importance with permutation-based importance derived from the testing set.

2.5. Breakpoint Analysis and Mapping

We applied a segmented regression model to statistically identify the presence of a breakpoint (ψ) in the response of aboveground biomass carbon density (cdAGB) to a previously derived bioclimatic predictor (Z). This model introduces a piecewise linear relationship by adding a nonlinear term to the Generalized Linear Model (GLM) linear predictor. The fundamental equation for the segmented relationship is as follows:

$\text{link}\left(E\left|{\text{cdAGB}}_i\right|\right)={\beta }_0+{\beta }_1{z}_i+{\beta }_2{\left(z_i-\psi \right)}_{+}+\text{Error}$ (1)

where:

E|cdAGB_i|: Is the expected value of Aboveground Biomass Carbon Density (cdAGB) for observation i.

Z: Is the independent bioclimatic variable (segmented variable).

ψ: Is the estimated change-point (or breakpoint), representing the bioclimatic threshold value.

β₁: Is the left slope (the initial rate of change before ψ.

β₂: Is the difference-in-slopes. The slope after the threshold is β₁ + β₂.

β₀: Is the intercept term.

(z_i − ψ)₊: Is the segmentation term, defined as $\left(z_i-\psi \right)\times \left(z_i>\psi \right)$ .

I(·): Is the indicator function (equal to 1 if z_i > ψ, and 0 otherwise).

The estimation of the breakpoint and model parameters was performed using the iterative algorithm implemented in the R package segmented (Muggeo, 2008) which translates the nonlinear problem into an approximate standard linear framework. To test whether the relationship is truly segmented versus a simple straight line (i.e., whether the breakpoint exists), the Davies Test was employed for the null hypothesis $H_0:{\beta }_2\left(\psi \right)=0$ .

If $H_0$ was rejected, the data for the bioclimatic variable (Z) were classified into two categories: (a) Vulnerable (V), representing sites where x < ψ, and (b) Non-vulnerable (N-V), representing sites where Z > ψ, based on the statistical significance (p < 0.05) of the slope in the segmented regression.

Breakpoint robustness and uncertainty (ψ) were quantified using a non-parametric bootstrap approach (1000 resamples) applied to the segmented regression model, with percentile-based confidence intervals providing an empirical, assumption-independent measure of breakpoint stability. Breakpoint-based classifications and spatial analyses were performed only when the Davies test indicated a statistically significant segmented relationship (p < 0.05); otherwise, relationships were treated as linear and all sites were conservatively classified as Non-vulnerable. When segmentation was supported, pre- and post-threshold slopes were estimated, and slope-change significance was assessed via the segmented-term t-statistic, enabling a quantitative evaluation of the magnitude and direction of climate sensitivity shifts across the breakpoint.

2.6. Spatially Controlled Climatic Analyses and Vulnerability Mapping of cdAGB

Partial correlation analyses were conducted to quantify associations between cdAGB and bioclimatic variables while controlling for latitude and longitude, thereby accounting for spatial structure. Altitude was excluded to avoid overcontrol of climate-driven signals, and Spearman’s rank correlation was employed to ensure robustness against non-normality and outliers.

Residual spatial dependence was evaluated using Moran’s I computed on segmented-regression residuals with a k-nearest neighbor weighting scheme (k = 5), where non-significant values indicated adequate control of spatial autocorrelation.

Subsequently, the sites were spatialized through a Vulnerability Map in a Geographic Information System (GIS) using the ggplot2 (Wickham et al., 2025) and sf (Pebesma, 2018) r libraries. The spatial vulnerability of cdAGB for each species was assessed as a binary outcome (Vulnerable vs. Non-vulnerable sites) using binomial generalized additive models (GAMs) with a logit link function, with the objective of estimating the spatial probability of cdAGB vulnerability. The models incorporated a two-dimensional spatial smoother (longitude-latitude) together with a one-dimensional elevation smoother, both estimated using restricted maximum likelihood (REML) and penalization to prevent overfitting. All statistical analyses were performed in R version 4.3.1 (R Core Team, 2023).

3. Results

Based on INFyS data, Quercus magnoliifolia exhibited the highest number of inventoried trees (n = 12,230) across 423 sites, Quercus radiata was the least frequent species, with 143 individuals recorded across 21 sites. Among pines, Pinus oocarpa recorded the largest number of trees (n = 6735, 419 sites), while Pinus devoniana recorded the lowest (n = 966, 96 sites; Table 1).

Pinus patula exhibited the highest cdAGB (49.1 Mg ha⁻¹), whereas Quercus radiata showed the lowest (2.7 Mg ha⁻¹). Mean cdAGB per genus was 23.75 Mg ha⁻¹ for Pinus and 13.21 Mg ha⁻¹ for Quercus.

Table 1. Dendrometric information for 12 temperate-forest species in Mexico used for the vulnerability analysis of aboveground carbon density with respect to bioclimatic variables.

Notes. n: number of trees; DBH: diameter at breast height (cm); HT: total height (m); cdAGB (Mg ha⁻¹): carbon density of aboveground live biomass (Mg ha⁻¹).

3.1. Bayesian Correlations of Carbon Density with Climate and Species

The response of cdAGB to bioclimatic variables was complex, exhibiting both positive and negative relationships and varying across species (Figure 1).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 1. Bayesian correlation between 19 bioclimatic variables and the carbon density of aboveground live biomass (cdAGB) for 12 temperate-forest species in Mexico. The vertical line separates temperature-related bioclimatic variables (left) from precipitation-related variables (right). Blue bars denote negative correlations and green bars denote positive correlations. The value above each bar indicates the percentage in the ROPE (Region of Practical Equivalence). Panels: (a) P. devoniana; (b) P. douglasiana; (c) P. leiophylla; (d) P. oocarpa; (e) P. patula; (f) P. pseudostrobus; (g) Q. crassifolia; (h) Q. laurina; (i) Q. magnoliifolia; (j) Q. radiata; (k) Q. resinosa; (l) Q. rugosa.

The three species showing the highest mean absolute Bayesian correlations (|ρ^Bayes|) between cdAGB and all bioclimatic variables were P. douglasiana (ρ = 0.24 ± 0.17, pd = 0.91, BF₁₀ > 10⁶), Q. resinosa (ρ = 0.23, pd = 0.90, BF₁₀ > 10⁵), and Q. rugosa (ρ = 0.17, pd = 0.99, BF₁₀ > 10⁸).

In contrast, the lowest absolute correlations were found in P. devoniana (ρ = 0.06, pd = 0.70, BF₁₀ = 0.34), P. leiophylla and P. oocarpa (ρ = 0.11, pd = 0.92, BF₁₀ > 10⁵), and Q. laurina (ρ = 0.13, pd = 0.91, BF₁₀ = 2.40).

3.2. Importance of Bioclimatic Variables

The analysis indicates that thermal variables predominate over precipitation variables. Temperature seasonality (bio4) was the most important predictor of cdAGB, exhibiting 100% relative importance in four of the 12 species (P. douglasiana, P. leiophylla, Q. crassifolia and Q. resinosa). Annual temperature range (bio7) was the second most important variable, showing 100% relative importance for three species (P. oocarpa, P. patula and Q. resinosa) (Figure 2).

By contrast, precipitation of the driest month (bio14) showed the lowest importance (0%) in 8 of the 12 species (four Pinus and four Quercus) (Figure 2).

Based on the average importance of the variables, the species showing the greatest climatic dependence were Q. magnoliifolia (67.2%), Q. rugosa (64.8%), P. oocarpa (63.3%), and P. pseudostrobus (61.1%), whereas the least sensitive species were Q. radiata (17.8%), P. patula (22.7%), and P. devoniana (36.4%).

Figure 2. Relative importance of 19 bioclimatic variables for predicting the carbon density of aboveground live biomass (cdAGB) for 12 temperate-forest species in Mexico. Each cell indicates the importance value estimated by a Random Forest model. From left to right: P. patula, P. devoniana, Q. radiata, Q. crassifolia, P. douglasiana, P. leiophylla, Q. resinosa, Q. laurina, P. oocarpa, Q. rugosa, P. pseudostrobus, and Q. magnoliifolia.

3.3. Vulnerability Thresholds in Pinus and Quercus Species

Statistically significant thresholds or breakpoint (ψ) were detected in 41.6% of the analyzed species, corresponding to those explicitly denoted by the † symbol (Table 2). Species-level segmented regression analyses indicated that five taxa (P. leiophylla, P. oocarpa, P. patula, Q. magnoliifolia, and Q. rugosa) exhibit statistically significant nonlinear responses of cdAGB across bioclimatic gradients, as evidenced by the Davies test (p ≤ 0.05).

Table 2. Species-level segmented regression analyses of carbon density in aboveground biomass as a function of bioclimatic variables, integrating collinearity, threshold detection, segment-specific slopes, and spatial autocorrelation.

VIF: variance inflation factor; Bio+: bioclimatic predictor with the highest relative importance; Davies-p: p-value from the Davies test used to detect a breakpoint; ψ: estimated bioclimatic threshold (breakpoint) of the segmented regression model; CI ψ: 95% confidence interval for ψ; U1.x − t: t statistic associated with the post-threshold change in slope; β₁ and β₂: slopes of the first and second segments of the model, respectively (Mg ha⁻¹); β₁ − t and β₂ − t: t statistics associated with β₁ and β₂; Corp-p: p-value of the partial correlation; IM: Moran’s I; p(IM): p-value of Moran’s I test; †: statistically significant breakpoint (Davies test, p < 0.05). BIO1, annual mean temperature (˚C); BIO2, mean diurnal temperature range (˚C); BIO4, temperature seasonality (SD ×100, unitless); BIO7, annual temperature range (˚C); BIO15, precipitation seasonality (coefficient of variation, unitless); and BIO19, precipitation of the coldest quarter (mm).

Multicollinearity was controlled using variance inflation factors (VIF), which remained below 10 across all species-specific models (maximum in Quercus laurina, VIF = 9.14; minimum in Quercus magnoliifolia, VIF = 1.49). The U1.x−t statistic approached and exceeded the critical threshold (|t| > 3) in species exhibiting a breakpoint, most notably Pinus patula (t = 3.96) and Quercus rugosa (t = 3.87), providing strong evidence of highly significant slope shifts following the threshold.

In Pinus oocarpa, cdAGB exhibited a pronounced asymmetric response across the thermal threshold identified by segmented regression (ψ = 13.28; p < 0.05). CdAGB remained stable below ψ but declined significantly beyond this threshold (β₂ = −3.209 Mg ha⁻¹ ˚C⁻¹), indicating increased vulnerability to thermal variability (Table 2; Figure 3(a)).

(a) (b)

(c) (d)

(e)

Figure 3. Segmented regression relationships between aboveground live biomass carbon density (cdAGB; Mg ha⁻¹) and the most influential bioclimatic predictor. Only species exhibiting statistically significant bioclimatic thresholds according to the Davies test (p ≤ 0.05) are shown. Points represent INFyS sampling sites. Red points indicate the range of the predictor in which significant changes in cdAGB are detected according to the segmented regression model (i.e., a threshold-driven response), whereas green points indicate ranges where no significant changes in cdAGB are detected. Solid black lines represent the fitted segmented regression, and shaded areas correspond to the 95% confidence intervals. Panels correspond to Pinus leiophylla (a) Pinus oocarpa (b), Pinus patula (c), Quercus magnoliifolia (d), and Quercus rugosa (e).

In Pinus patula, threshold effects manifested prior to ψ = 16.80, with a sharp pre-threshold decline in cdAGB (β₁ = −43.470 Mg ha⁻¹ ˚C⁻¹), followed by a comparatively stable response after ψ was exceeded (Table 2; Figure 3(b)).

In Quercus magnoliifolia and Q. rugosa, cdAGB responses were limited to the pre-threshold segment, with thresholds at ψ = 20.04 and ψ = 90.00, respectively Table 2; Figure 3(c), Figure 3(d). CdAGB increased in Q. magnoliifolia (β₁ = 2.206 Mg ha⁻¹˚C⁻¹), whereas Q. rugosa showed a weak decline (β₁ = −0.154 Mg ha⁻¹ mm⁻¹). The absence of post-threshold effects suggests stabilization of cdAGB beyond ψ.

Significant partial correlations (Corp-p ≤ 0.05) in three of the five species with identified breakpoints indicate an independent influence of the dominant bioclimatic variable on cdAGB. Low Moran’s I values (IM ≈ 0.001-0.212; Table 2) confirm the absence of residual spatial autocorrelation, supporting the robustness of the spatial model specification.

Segmented regression revealed breakpoints between the 30.3rd (P. patula) and 59.9th (Q. rugosa) percentiles of the predictor gradient, indicating that thresholds occur within the central portion of the data distribution and are not driven by outliers.

3.4. Spatial Distribution of Vulnerable and Non-Vulnerable Sites

According to the GAM-based models, the spatial component was highly significant for all species (p < 0.01) and was characterized by elevated estimated degrees of freedom (EDF > 1), indicating complex nonlinear spatial patterns in vulnerability probability across the geographic landscape (Table 3). By contrast, the effect

Table 3. Summary of species-specific binomial generalized additive models (GAMs) evaluating spatial and altitudinal effects on vulnerability status, reporting smooth-term penalization, inferential statistics, model discrimination (AUC), dispersion, and deviance explained.

EDF: Effective degrees of freedom of smooth terms (values > 1 indicate non-linear relationships, whereas values < 1 reflect penalized, near-linear effects); Statistic (χ²): Wald-type chi-square test statistic assessing the significance of smooth terms; p-value: Statistical significance of smooth effects; AUC: Area Under the Receiver Operating Characteristic (ROC) Curve, measuring the model’s discriminative performance; Dispersion (φ): Estimated dispersion parameter of the binomial GAM; values well below 1 indicate no evidence of overdispersion and reflect strong penalization and well-regularized model fit; Dev. Expl. (%): Proportion of deviance explained by the model relative to a null model. Only species exhibiting statistically significant bioclimatic thresholds according to the Davies test (p ≤ 0.05) are shown.

of elevation was weak or non-significant for most species (p > 0.05), with EDF values below 1 and only marginal levels of significance. A notable exception was Quercus magnoliifolia, for which elevation exerted a strong, significant, and nonlinear effect (p < 2 × 10⁻¹⁶).

The models exhibited excellent discriminative ability, with AUC values ranging from 0.94 to 0.99, and explained a substantial proportion of the deviance (63-93%), indicating a robust and reliable fit. Low dispersion parameter values (φ < 0.6) confirmed the absence of overdispersion and appropriate model regularization.

Figure 4. Spatial distribution of sites from the National Forest and Soil Inventory (2015-2020) for the genus Pinus and their vulnerability status. In this analysis, N-V represents non-vulnerable sites (green symbols), and V indicates vulnerable sites (red symbols). Circle size represents carbon density (Mg ha⁻¹), while the color gradient denotes the level of vulnerability.

Based on the GAM results, spatial patterns of vulnerable and non-vulnerable cdAGB sites relative to the dominant bioclimatic predictor are illustrated in Figure 4 (pine species) and Figure 5 (Quercus species). These patterns demonstrate a clear geographic segregation between vulnerable and non-vulnerable sites, underscoring the dominant role of spatial structure (latitude and longitude) over simple elevational gradients. In Pinus species, vulnerable sites predominantly occur at lower latitudes (Figure 4), whereas in Quercus species, vulnerability is mainly associated with lower longitudes (Figure 5).

The proportion of vulnerable sites was highest for Q. magnoliifolia (61.5%) and Q. rugosa (59.9%), and lower for P. patula (30.3%) and P. oocarpa (40.8%).

Figure 5. Spatial distribution of sites from the National Forest and Soil Inventory (2015-2020) for the genus Quercus and their vulnerability status. In this analysis, N-V represents non-vulnerable sites (green symbols), and V indicates vulnerable sites (red symbols). Circle size represents carbon density (Mg ha⁻¹), while the color gradient denotes the level of vulnerability.

4. Discussion

The aboveground carbon density (cdAGB) values obtained in this study, with maxima approaching 60 Mg ha⁻¹ for most species (Table 1), are consistent with previous reports for Mexico. Earlier, Cartus et al. (2014) generated cdAGB maps with values exceeding 50 Mg ha⁻¹. More recently, using data from the National Forest Inventory (2009-2014), Girón et al. (2024) reported that P. patula reaches accumulations of up to 53.59 Mg C ha⁻¹. These findings exceed values reported for other regions, including the evergreen needleleaf forests of northwestern Eurasia and western North America, where aboveground carbon accumulation potentials of approximately 21.45 ± 23.38 Mg C ha⁻¹—derived from biomass estimates of 42.9 ± 46.76 Mg ha⁻¹ (Chen et al., 2023)—have been documented.

In agreement with the literature, this study confirms that the correlation between aboveground carbon density (cdAGB) and bioclimatic variables is complex and inconsistent, depending on forest type and species (Stegen et al., 2011; Girón et al., 2024). At a global scale, various authors have reported a generally positive relationship between cdAGB and mean annual temperature (MAT, bio1) in coniferous forests (Ali et al., 2020; Chen et al., 2023; Liu et al., 2014; Reich et al., 2014). However, species-level results from this study reveal that this correlation can be either positive or negative (Figure 1), demonstrating a differentiated response. This variability is consistent with regional findings, such as those reported by Ma et al. (2023) and Zhang et al. (2023), who documented a negative relationship between MAT and aboveground biomass in coniferous forests in China.

Regarding precipitation, the literature suggests that its correlation with aboveground carbon density (cdAGB) is more consistent, being predominantly positive, although with values generally not exceeding 0.37 (Balima et al., 2021; García et al., 2024; Ma et al., 2023; Zhang et al., 2023), which aligns with the findings of the present study. Considering vulnerability as the degree of susceptibility of a system to adverse climate effects (IPCC, 2007), it can be inferred that a higher absolute correlation (|ρBayes|) between cdAGB and climatic variables indicates greater sensitivity and, consequently, higher vulnerability.

Under this framework, the most sensitive species were P. douglasiana (ρBayes = 0.25), Q. resinosa (ρBayes = 0.23), and Q. rugosa (ρBayes = 0.17). In contrast, the species exhibiting lower climate dependence and, therefore, lower vulnerability were P. devoniana (ρBayes = 0.06), P. leiophylla and P. oocarpa (ρBayes = 0.12), and Q. laurina (ρBayes = 0.13).

Frequency analysis of predictors identified by the Random Forest models showed a clear dominance of thermal bioclimatic variables, which comprised 75% of the predictors across the analyzed species (Table 2). Within this group, Bio2 (mean diurnal temperature range) was the most frequently selected predictor (33.3%), followed by Bio1 (annual mean temperature) (16.7%); collectively, these two variables accounted for 50% of all selected predictors. By contrast, precipitation-related variables (Bio12-Bio19) were markedly underrepresented, comprising only 25% (Table 2). Consistent with these results (Figure 1), previous studies have documented an inverse relationship between temperature-derived variables and aboveground biomass across forest types (Chave et al., 2014; Vieilledent et al., 2016).

The importance of bio4 in other models reaches up to 68% (Vieilledent et al., 2016). However, the significance of climatic predictor variables varies considerably across studies, highlighting the dependence on factors such as scale, forest type, and the species analyzed.

Notably, Bio2—identified as the most influential variable in this study—is rarely reported as a primary predictor. Other studies have identified bio9 (Bennett et al., 2020), bio17 (Stegen et al., 2011), bio12 (Chen et al., 2023), or bio1 (Wang & Ali, 2022) as the most influential variables for estimating aboveground biomass, thus evidencing the multifactorial nature of forest biomass distribution.

For the specific context of Mexico, Sandoval et al. (2024) reinforce this specificity, concluding that temperature variables are more robust predictors of aboveground carbon density (cdAGB) in coniferous forests than precipitation variables. Although in this study annual precipitation (bio12) was not the most important predictive variable (Figure 1), its positive correlation is consistent with reports from tropical forests and temperate forests in Durango (Balima et al., 2021; Ma et al., 2023; García et al., 2024). Likewise, the dual correlation (both positive and negative) observed with mean annual temperature (bio1) aligns with the variability reported in forests worldwide (Liu et al., 2014; Reich et al., 2014; Ma et al., 2023), highlighting the complexity of these interactions.

A central finding of this study is the non-linear nature of the relationship between aboveground carbon density (cdAGB) and bioclimatic variables (Table 2, Figures 3-5). This result has critical implications for modeling carbon reservoirs, suggesting that commonly used linear models may significantly under- or overestimate predictions under different climate scenarios. These findings are consistent with those reported by Vieilledent et al. (2016), who also identified climatic inflection points in the relationship between AGB, annual precipitation, and mean temperatures. More specifically, Guo et al. (2019) quantified a global forest threshold, reporting that the relationship is positive when mean annual precipitation is ≤1100 mm (with an increase of 7.1 Mg C ha⁻¹ per 100 mm), but becomes negative above this threshold (decreasing by 1.9 Mg C ha⁻¹ per 100 mm). This evidence underscores the need to adopt non-linear modeling approaches to accurately capture the vulnerability of forest biomass to climate change.

The vulnerability of forest ecosystems to climate change constitutes a global concern (Larjavaara et al., 2021), and the temperate regions of Mexico have been identified as particularly sensitive, with projections indicating risks of drastic reduction or even disappearance (Villers & Trejo, 1997; McKenney et al., 2007).

In this context, the present study advances understanding of this threat by identifying clear geographic patterns of cdAGB vulnerability driven by explicit spatial structure—latitude and longitude—but not by elevation (Table 3; Figure 4 and Figure 5). Vulnerable sites in Pinus species are predominantly associated with lower latitudes, whereas in Quercus species vulnerability is primarily linked to lower longitudes. This finding is particularly relevant because, although few studies have directly linked vulnerability to latitude at this scale, these same high latitudes were observed to exhibit the highest aboveground carbon density (cdAGB). This observation partially aligns with Lin et al. (2010), who reported an increase in biomass at mid-latitudes (30˚ to 75˚).

Finally, the latitudinal vulnerability identified is amplified when considering the observed and projected climatic trends for Mexico. Historically, the northern region of the country has experienced extreme warming (+2.7˚C between 1951 and 2017) and a drastic reduction in precipitation (up to −70 mm year⁻¹), a trend that contrasts with the southern region, which has been wetter (Cuervo-Robayo et al., 2020; Murray-Tortarolo, 2021). Climate projections exacerbate this scenario: by the end of the century (2070-2099), the northwest is expected to experience an additional warming of up to 3.5˚C and a precipitation decrease of 10 to 25% (Magaña et al., 2012).

This regional climatic divergence, which imposes severe water stress in the north, coincides with the location of high-biomass vulnerable Quercus sites. Therefore, the combination of the intrinsic sensitivity of these ecosystems (non-linear relationships) and the intensification of drought and extreme heat conditions confirms that the temperate forests of northern Mexico face a critical risk threatening their stability and their capacity as carbon sinks.

5. Conclusion

The response of aboveground carbon density (cdAGB) to bioclimatic variables is fundamentally non-linear. Thermal variables, particularly temperature seasonality (Bio2) and annual mean temperature (Bio1), were the principal predictors of cdAGB in the studied species, whereas precipitation-related variables exhibited substantially lower explanatory importance. The presence of critical thresholds in in five of the studied species demonstrates that carbon accumulation experiences an abrupt regime change at a climatic inflection point, most often characterized by a pronounced decline, with the response occurring either before or after the estimated threshold. At the genus level, Quercus exhibited a markedly higher proportion of vulnerable sites (>60%) than Pinus (30.3% - 40.8%), indicating that oak biomass is intrinsically more susceptible to climatic fluctuations. The probability of vulnerable and non-vulnerable site occurrence was primarily structured by latitude and longitude rather than elevation, with sites exhibiting pronounced spatial clustering. Collectively, the results reveal strong non-linear climate-cdAGB interactions, highlighting the critical role of ecological thresholds in refining forest dynamic models and guiding evidence-based conservation.

Acknowledgements

The authors gratefully acknowledge the Autonomous Agrarian University Antonio Narro for providing the necessary facilities, equipment, and personnel that made this research possible.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References