Mapping the Future: Climate-Induced Changes in Biomass Density Across Mexico’s Coniferous Forests

1. Introduction

Forests occupy approximately 31% of the Earth’s terrestrial surface, equivalent to around 4.06 billion hectares, and play a pivotal role in climate regulation by acting as major carbon sinks. It is estimated that global forests sequester approximately 662 gigatonnes (Gt) of carbon, distributed as follows: 44% in living biomass, 45% in soil organic matter, and 11% in dead wood and litter [1].

In Mexico, forests cover approximately 64.8 million hectares, accounting for about 33% of the national land area [1]. These forests encompass a wide range of types, including temperate, broad-leaved, mixed, and tropical forests [2], and are estimated to store around 1.69 Gt of carbon [3]. Changes in temperature and precipitation patterns driven by climate change have been shown to adversely affect forest growth, biomass accumulation, and carbon sequestration at a global scale [4,5].

The scientific foundation for studying climate change is unequivocal, and the corresponding findings are irrefutable. Between 2011 and 2020, the global average surface temperature rose by 1.09 °C compared to the 1850-1900 baseline, with more pronounced warming observed in the Northern Hemisphere, particularly over land (1.59 °C) compared to oceans (0.88 °C). Precipitation patterns are shifting, with increases observed at higher latitudes and decreases in the subtropics. Climate extremes are projected to become more frequent and intense, while the efficiency of carbon sinks is expected to diminish [6].

In Mexico, the average temperature increased by 0.31 °C per decade between 1971 and 2020, with the most pronounced effects occurring in the northern plateau during the summer months. According to projections under the SSP3-7.0 climate scenario, the average temperature is expected to rise by 0.82 °C between 2020 and 2039, and by 1.63 °C between 2040 and 2059. Moreover, certain northern states are likely to experience significant reductions in precipitation [7].

Precipitation and temperature are well-established as primary drivers of forest biomass productivity. Precipitation plays a pivotal role in plant physiology by regulating processes such as transpiration, nutrient uptake, stomatal conductance, and nutrient availability [8,9], in addition to shaping species-specific growth strategies. Variations in precipitation availability also impact water use efficiency [10,11], creating a direct relationship between precipitation levels and biomass productivity [2].

Likewise, temperature exerts significant control over plant growth by modulating fundamental physiological processes, including photosynthesis, respiration, and transpiration [4,12]. It also governs the rate of chemical reactions and CO₂ assimilation [13]. While rising temperatures may stimulate biomass accumulation in boreal forests, they tend to suppress growth in tropical forests [14]. Excessively high temperatures reduce growth rates, alter leaf pigmentation, impair root system development, and induce water stress, thereby disrupting normal growth patterns [15].

Plants from temperate climates exhibit a degree of cold tolerance, but extremely low temperatures negatively affect processes such as cell division, photosynthesis, and metabolic activity, leading to reduced productivity [16].

Climatic variations are driving forests to become increasingly dynamic systems [17], resulting in changes in tree species composition within ecosystems [11]. Simultaneous shifts in temperature and precipitation lead to either reductions or increases in ecosystem biomass, thereby altering global forest distribution patterns [18].

It has been well established that temperature exerts both positive and negative influences on aboveground forest biomass [4,19,20]. Furthermore, precipitation has been consistently shown to have a positive correlation with aboveground biomass [12,21]. These effects have been observed across various forest biomes, including temperate, tropical, and boreal ecosystems [14], as well as at the species level [8].

To estimate aboveground biomass in forest ecosystems, both dynamic variables (related to climate) and static variables (such as diameter at breast height, basal area, stand density, slope, aspect, and elevation, among others) have been jointly applied [8,22,23,24]. These mixed models, which incorporate intrinsic and extrinsic factors influencing aboveground biomass, have demonstrated moderate predictive capacity, with determination coefficients (R²) ranging from 0.19 to 0.86 [8,22,23]. In some cases, these models include nine or more explanatory variables [8].

While these models are useful in identifying correlations between biomass and environmental variables, they are limited in their ability to quantify the direct impacts of climate change on aboveground biomass [8,23,25]. Furthermore, forests are highly dynamic systems that undergo temporal changes due to natural disturbances, anthropogenic activities, and climatic fluctuations, which introduces additional complexity to the modeling process [24,26].

The research conducted by [27] indicated that, among ten conifer species, precipitation demonstrated a positive correlation with aboveground biomass density (0 ≤ ρ ≤ 0.20), whereas temperature exhibited a negative correlation (−0.20 ≤ ρ ≤ 0). Importantly, the magnitude and direction of the correlation between climatic variables and aboveground biomass in forest ecosystems are highly contingent upon several factors, such as spatial scale [28], forest type [4], and species-specific traits [12,29].

In order to deepen our understanding of how climate change specifically influences aboveground biomass, it is imperative to develop new modeling frameworks. These models should integrate dynamic climatic data to directly assess the effects of climate variability and enable future projections—an objective that cannot be fully achieved by models that rely solely on static stand or site-level variables.

Building upon the aforementioned considerations, this study's primary objective was to develop bioclimatic and dynamic models specifically designed to evaluate the impacts of climate change on the geospatial distribution of aboveground live biomass carbon density in Mexico's pine forests. The analysis incorporated four distinct climate scenarios based on Representative Concentration Pathways (RCPs) 2.6, 4.5, 6.0, and 8.5 W/m², projected for the years 2050 and 2070. These models are intended to provide highly accurate information that can enhance forest management practices and contribute to climate change mitigation efforts. It is hypothesized that the future climate projections will exert a significant influence on the aboveground biomass carbon density in the pine forests of Mexico.

2. Materials and Methods

2.1. Study Area Description

Conifers in Mexico are found from sea level up to over 4,000 m in altitude, with the highest diversity concentrated in the Sierra Madre Occidental (SMOc) and the Sierra Madre Oriental (SMO), and Mexico hosts 49 species of pine, representing 40% of the approximately 120 pine species worldwide [30]. They span approximately from 32 °N north latitude to 19 °N, and from 105 °W west longitude to 98 °W, on soils derived from igneous and metamorphic rocks [31]. These ecosystems are characterized by a temperate to cold climate, with annual precipitation ranging from 350 to 1200 mm and average temperatures ranging from 6 to 28 °C [32], with rainy summers and dry, cold winters, dominated by coniferous trees such as pines (Pinus spp.), firs (Abies spp.), cedars (Cupressus spp.), and junipers (Juniperus spp.) [33].

2.2. Data Acquisition and Cleaning Process

A total of 22 parameters including type, topographic, species richness, biomass, and carbon (aboveground and belowground), health indices, climate, geographic, and dasometric factors were obtained from [2]. From these, we used data on carbon density in aboveground live biomass (henceforth referred to as _cdAGB), geographic coordinates (lat, long), and ID of each cluster, exclusively from the parameter type='coniferous forests'. The _cdAGB is the carbon stock per unit area of live biomass contained in living trees, including stems, branches, leaves, and seeds [34]. These data were calculated from the National Forest and Soil Inventory (INFyS) of Mexico (2009-2012) [2]. Additionally, we utilized 19 bioclimatic variables from WorldClim version 2.0, covering the period from 1970 to 2000 (Appendix A), in raster format with a resolution of 30 arcseconds, published in January 2020 [35]. These variables have been previously used in similar studies [8,9,12,14,28,36].

Using the 'raster' library in R [37], we employed the geographic coordinates of each cluster to extract values for the 19 bioclimatic variables (Bios). To identify erroneous or atypical data in _cdAGB, we applied a Principal Component Analysis (PCA) to the centered and scaled matrix of _cdAGB and 19 Bios, using the 'FactoMineR' v.2.9 library [38]. We considered outliers as those falling outside the 95% confidence ellipse of the PCA, which were subsequently excluded from the dataset to ensure the coherence and reliability of the results.

2.3. Predictive Modeling of _cdAGB

Due to the significant climatic variability recorded in the pine forests of the country [31], and aiming to enhance predictions of _cdAGB, a bioclimatic stratification of Mexico's pine forest areas was conducted. We utilized 19 bioclimatic variables in raster format, employing the "GeoStratR" library [39], to create strata (geographical space) in the same format. The clusters from INFyS located within each generated stratum were separated for independent analysis.

For the selection of predictors (Bios) of _cdAGB, within each stratum, the 'machine learning' (ML) technique employing 10-fold cross-validation was utilized, employing backward selection. A tuning grid ranging from 1 to 7 predictors (nvmax = 1:7) was considered [40]. From each subset of predictors identified by ML, the statistical significance of regression coefficients (p < 0.05) and the variance inflation factor (VIF) were evaluated to mitigate multicollinearity effects.

With the subset of predictors that met these criteria, the model was trained using the Generalized Linear Models (GLM) method with the 'identity' link function. For this purpose, 80% of the data (randomly obtained and by quantiles) from the previously cleaned matrix was used. This approach (GLM) is suitable for continuous and positive response variables [2,9] in this type of study.

The model fitting was performed using the 'caret' package in R [40]. Hypothesis tests ($\propto =0.05$) were conducted on the regression coefficients ($H_0:{\beta }_i=0vs{H}_1:{\beta }_i\ne 0,\dots ,{\beta }_{ij}$) of the final model. We utilized the "lmg" metric to calculate the importance of each bioclimatic variable (as a percentage of model explained variance), which provides a decomposition of the model's explained variance into non-negative contributions [41]. With the remaining 20% of the data, the model was validated using various techniques: Leave-One-Out Cross-Validation (LOOCV), Cross-Validation (CV; k=10), Repeated Cross-Validation (RCV, k=10, rep=10), and Bootstrap (reps=100). The evaluation included calculations of Root Mean Square Error (RMSE), Pseudo Coefficient of Determination (Rsquared), and Mean Absolute Error (MAE) to assess the model's performance. These validations were conducted using the 'caret' R package [40].

2.4. Current and Future Prediction and Rate of Change of _cdAGB

Using the generated bioclimatic models, current predictions of _cdAGB were made for each stratum, employing the corresponding bioclimatic predictors. The raster library of the R software [37] was utilized with the argument `type = "response"` to generate raster maps of _cdAGB predictions.

To predict future _cdAGB, raster layers from three different General Circulation Models (GCMs) were obtained: MIROC-6, GISS-E2, and CMCC-ESM2. From these models, only the bioclimatic predictors of _cdAGB for each stratum were selected. The GCMs are derived from future climate projections of the Coupled Model Intercomparison Project Phase 5 (CMIP5), at a resolution of 2.5 degrees. Four Representative Concentration Pathways (RCPs) developed by the Intergovernmental Panel on Climate Change [42] were used, characterized by their projected Total Radiative Forcing (TRF) for the year 2100, ranging from 2.6, 4.5, 6.0, to 8.5 W/m², considering the years 2050 and 2070. These RCP scenarios are used to project potential future trajectories of climate change based on different emission levels and mitigation efforts.

Through raster algebra, the 3 raster layers of each GCM, RCP, and year were averaged. Using this average and the corresponding bioclimatic model, future prediction of _cdAGB was conducted for each stratum. Similarly, the R raster library [37] was employed with the argument `type = "response"` to generate raster maps of future _cdAGB predictions. To identify changes in future _cdAGB relative to the current scenario, the following expression was used: _cdAGB (future) - _cdAGB (current). The results are as follows: 0 if pixel values are equal in both scenarios; positive if _cdAGB is higher in the future; and negative if _cdAGB is higher in the current scenario.

To visualize the vulnerable areas of Mexico's coniferous forests to changes in _cdAGB due to climate change, a 40 km radius buffer was created around each site from the National Forest and Soil Inventory. This delineated the forested areas potentially affected. Subsequently, the number of pixels for each _cdAGB change category or rate within these buffer polygons was quantified. The resulting change rate was then presented in map form for each RCP and year.

In addition, uncertainty in future predictions of _cdAGB was quantified using the standard error: EE(θ) = √[Var(θ). This involved repeated predictions of _cdAGB from MIROC-6 + GISS-E2 + CMCC-ESM2, for each RCP and year, which were averaged using the raster library (Hijmans, 2010), and generating a raster map of EE(θ).

2.5. Bioclimatic Predictors: Current Analysis and Future Projections

In order to determine the influence of future bioclimatic predictors on the changes in _cdAGB of Mexican forests, the non-parametric Wilcoxon signed-rank test was employed. This test was used to evaluate whether there are significant differences between the medians of two samples: the current bioclimatic predictor and the future bioclimatic predictor (e.g., current Bio1 vs. Bio1 2050, RCP 4.0). A significance level of 95% was used.

For this purpose, the current and future values (average of 3 GCM, RCP, and year) of the predictors for each stratum (Bios) were extracted using the geographic coordinates of each cluster and employing the raster library [37]. These analyses were represented in violin plots, with temperature-derived Bios shown in orange and precipitation-derived Bios in blue, with statistical significance indicated by the * symbol at the top.

All statistical analyses, figures, and geographic processing were performed using R software [43].

3. Results

3.1. Distribution Carbon Density Of Aboveground Live Biomass

The stratification process of the coniferous forests in Mexico resulted in three strata: Stratum I, in the northwest of the country (Figure 1a; n = 60), Stratum II, in the central-southern part of the country (Figure 1b; n = 450), and Stratum III, located in the Sierra Madre Occidental, SMO (Figure 1c; n = 463).

3.2. Models For Predicting Carbon Density in the Aboveground Live Biomass

The bioclimatic variables selected by the algorithms were statistically significant (p < 0.05) for predicting _cdAGB. The predictors (Bios) of _cdAGB were representative of both precipitation variables (3 out of 6) and temperature variables (3 out of 6). However, at the stratum level, temperature variables are more important (Imp., Table 1) than precipitation variables. Bio 5 (Max Temperature of Warmest Month) was the best predictor of _cdAGB in 2 strata (I and II). No model exhibited multicollinearity (VIF < 1.16), which ensures that _cdAGB predictions are not overestimated.

The equations for estimating _cdAGB in each stratum are as follows:

Stratum I: _cdAGB (Mg C ha^-1) = 64.8332 + -0.1897*Bio 5 (°C x10) + 0.0586*Bio 18 (mm)

(1)

Stratum II: _cdAGB (Mg C ha^-1) = 87.6362 + -0.2572*Bio 5 (°C x10) + 0.0190*Bio 12 (mm)

(2)

Stratum III: _cdAGB (Mg C ha^-1) = 26.6379 + -0.0959*Bio 10 (°C x10) + 0.0933*Bio 13 (mm)

(3)

The mean of _cdAGB does not follow a pattern among strata; that is, it is not lower in stratum I and higher in stratum III. The mean _cdAGB in stratum II represents, on average, up to 1.6 times more than in strata I and III; and it ranges from 23.14 Mg C ha^-1 (stratum I) to 42.57 Mg C ha^-1 (stratum II), with maximum values reaching almost 180 Mg C ha^-1 (Table 2), but it also has the highest variability (CV > 60%). None of the variables showed normality (p < 0.0001) according to the Shapiro-Wilk and Anderson-Darling tests (Table 2).

3.3. Validation of Predictive Models for _cdAGB

According to the pseudo R² derived from the validation procedure, bioclimatic variables explain an average of 19% of the _cdAGB (Table 3). The LOOCV method calculates the lowest R² (0.10), while the other methods yield an average R² of 0.22. According to the validation, when using the models to predict aboveground biomass carbon density with independent data from the training set, the estimation errors (RMSE) could range from 13.18 Mg C ha⁻¹ (stratum III) to 42.42 Mg C ha⁻¹ (stratum I), with the best estimates observed in stratum III.

3.4. Current And Future Prediction of Carbon Density of Aboveground Live Biomass

Using the bioclimatic models generated here, at the INFyS sites, predictions range from 7.48 to 34.79 Mg C ha^-1; from 16.68 to 71.29 Mg C ha^-1; and from 10.98 to 41.91 Mg C ha^-1 for strata I, II, and III, respectively (Figure 2a-Figure 2c). The highest _cdAGB is observed/predicted in stratum II (Figure 2b) in the central region of the country. The highest differences between observed _cdAGB – estimated _cdAGB are -19 and +30; -52 and +124; -29 and +56 Mg C ha^-1, respectively for strata I, II and III.

According to future predictions of _cdAGB up to 2050 and 2070 (Figure 3, Figure 4 and Figure 5) using the bioclimatic model for each stratum (equations 1-3), under any climate scenario (RCP), only reductions in _cdAGB are expected. Forests in stratum II (Figure 4a - Figure 4h, situated in the central and southern regions of the country) would be the most affected, with decreases ranging from -5 to -10 Mg ha^-1, potentially up to -20 Mg C ha^-1 by 2070, whereas those in stratum III (Figure 5a – Figure 5h, SMO) would experience lesser impacts ranging from 0 to -5 Mg C ha^-1. The most significant losses of _cdAGB are expected under RCP 8.5 and are anticipated to be greater by 2070.

Table 4. Number of pixels contained within the 40 km radius buffer around each sampling site of INFyS to differentiate changes in carbon density of aboveground live biomass in the conifer forests of Mexico under future climatic scenarios. Each pixel represents 0.98 x 0.98 km on each side.

Stratum	Changes in cdAGB	RCP 26		RCP 45		RCP 60		RCP 85
		2050	2070	2050	2070	2050	2070	2050	2070
I	(-20 – -15 | )	0	0	0	0	0	0	0	0
	(-15 – -10 | )	0	0	0	0	0	35	0	107
	(-10 – -5 | )	0	240	2361	2280	906	2885	3011	3180
	(-5 – 0 | )	3287	3047	926	1007	2381	367	276	0
	(0 – 5 | )	0	0	0	0	0	0	0	0
II	(-20 – -15 | )	0	0	0	0	0	0	0	1148
	(-15 – -10 | )	0	30	6	1526	0	3928	4796	22466
	(-10 – -5 | )	17750	19473	23821	22317	23571	19912	19047	229
	(-5 – 0 | )	6093	4340	16	0	272	3	0	0
	(0 – 5 | )	0	0	0	0	0	0	0	0
III	(-20 – -15 | )	0	0	0	0	0	0	0	0
	(-15 – -10 | )	0	0	0	0	0	0	0	0
	(-10 – -5 | )	0	6	0	0	0	0	1290	3281
	(-5 – 0 | )	9680	9819	9825	9825	9825	9822	8535	6544
	(0 – 5 | )	145	0	0	0	0	3	0	0

Where: Observed carbon density of aboveground live biomass (Mg C ha^-1); RCP = Representative Concentration Pathways. The color of the box in column two corresponds to the color of the pixel in Figure 3, Figure 4 and Figure 5.

Table 4. Number of pixels contained within the 40 km radius buffer around each sampling site of INFyS to differentiate changes in carbon density of aboveground live biomass in the conifer forests of Mexico under future climatic scenarios. Each pixel represents 0.98 x 0.98 km on each side.

Stratum	Changes in cdAGB	RCP 26		RCP 45		RCP 60		RCP 85
		2050	2070	2050	2070	2050	2070	2050	2070
I	(-20 – -15 | )	0	0	0	0	0	0	0	0
	(-15 – -10 | )	0	0	0	0	0	35	0	107
	(-10 – -5 | )	0	240	2361	2280	906	2885	3011	3180
	(-5 – 0 | )	3287	3047	926	1007	2381	367	276	0
	(0 – 5 | )	0	0	0	0	0	0	0	0
II	(-20 – -15 | )	0	0	0	0	0	0	0	1148
	(-15 – -10 | )	0	30	6	1526	0	3928	4796	22466
	(-10 – -5 | )	17750	19473	23821	22317	23571	19912	19047	229
	(-5 – 0 | )	6093	4340	16	0	272	3	0	0
	(0 – 5 | )	0	0	0	0	0	0	0	0
III	(-20 – -15 | )	0	0	0	0	0	0	0	0
	(-15 – -10 | )	0	0	0	0	0	0	0	0
	(-10 – -5 | )	0	6	0	0	0	0	1290	3281
	(-5 – 0 | )	9680	9819	9825	9825	9825	9822	8535	6544
	(0 – 5 | )	145	0	0	0	0	3	0	0

Where: Observed carbon density of aboveground live biomass (Mg C ha^-1); RCP = Representative Concentration Pathways. The color of the box in column two corresponds to the color of the pixel in Figure 3, Figure 4 and Figure 5.

3.5. Evaluation of Climate Projections on Variables Predicting _cdAGB

The Wilcoxon test provided sufficient evidence (p<0.0001) to reject the null hypothesis (H0) of equal medians between the current Bio and the same Bio under a future climate scenario. In all cases, the median of the estimated temperature Bios for 2050 and 2070, projected by three General Circulation Models (GCMs), in any climate scenario (RCP; left, Figure 6), will be higher than the current median (Table 2). Generally, the temperature is expected to increase by 2.66 °C (2050) and 3.36 °C (2070) compared to the present (Figure 6a, Figure 6c, and Figure 6e). The coniferous areas in stratum II are predicted to be vulnerable due to a likely increase of +2.7 °C by 2050, whereas the less vulnerable areas will be in stratum III (+2.5 °C).

Generally, the test demonstrated that in the coniferous forests of Mexico, precipitation levels in 2050 and 2070 could statistically be the same (p < 0.05) as those currently recorded in stratum I (Figure 6b). In stratum II, precipitation is projected to decrease (p < 0.0001) by 5.3% (2050) and 6.4% (2070) (Figure 6d), while in stratum III, it could either increase or decrease (Figure 6f). When analyzed individually (Current vs Bio-RCP-Year), significant differences exist only for precipitation variables (right, Figure 6), and not for temperature variables (left, Figure 6).

The uncertainty (represented by the standard error) estimated for _cdAGB shown here (Figure 7), for the most critical scenario (RCP 85 | 2070), indicates that the standard error ranges from 2.5 Mg C ha^-1 (stratum III) to 6 Mg C ha^-1 (stratum II), well below the RMSE (Table 3) estimated with independent data. This uncertainty is similar in less critical climatic scenarios (RCP 26 | 2050; Figures not shown).

4. Discussion

The selection of _cdAGB predictors through automated algorithms (stepwise, Machine Learning, neural networks, etc.) is generally efficient. These algorithms have been tested in natural forests [44,45], tropical rainforests [46], and even in plantations [19]. Although the variables selected by these algorithms are statistically significant, p<0.05 [36], the models generated here exhibited multicollinearity effects (VIF > 10), consequently, these models had to be evaluated to mitigate this effect and enhance the accuracy of the predictions.

Some authors [47] have used linear procedures (lm) for predicting NPP (Net Primary Production) in grasslands, with favorable results. We tested this technique; however, it was demonstrated that in this type of study, it is difficult to meet all the assumptions of a regression model (${\epsilon }_i~N\left(0,{\sigma }^2\right)$, which is why the 'glm' procedure was ultimately employed, as has been done in such studies [2,9]. In addition to these procedures, Random Forest [8,48,49] and Bayesian models [20] have also been used for predicting AGB (Above-Ground Biomass) from temperature and precipitation variables.

4.1. Predictor Variables of Aboveground Biomass

In this study, the algorithms selected both temperature variables (Bio5 and Bio10) and precipitation variables (Bio12, Bio13, and Bio18) as the best predictors of _cdAGB. Temperature variables were the most important in the model (contributing 8 to 10% to the pseudo R²) for strata I and II (Table 1, Figure 1a and Figure 1b), but not for stratum III. According to some authors [21,47], on a global scale, precipitation and temperature variables are the best predictors of AGB. Specifically, annual mean temperature (Bio1) and temperature of the warmest quarter (Bio10) are climatic variables associated with biomass distribution on broad scales. Metrics demonstrate that temperatures appear to be more important than precipitation variables [28].

It has been demonstrated that bioclimatic variables of temperature (Bio1 and Bio5) and precipitation (Bio12) are associated with the accumulation of AGB in forests [20], in boreal ecosystems [19], temperate seasonal ecosystems [18], tropical rainforests [50], and tropical seasonal ecosystems [4,5]. As can be noted, Bio5 and Bio12 (Appendix A) are good predictors of AGB across different types of ecosystems. For instance, a study conducted in the SMO (stratum II of this same study) indicates that average temperature (Bio1) is the most important predictor for AGB in temperate forests [24].

In Australian forests, encompassing 15 types of forests predominantly dominated by eucalypts [8], it has been demonstrated that climatic variables are better predictors of AGB than soil variables, with Bio 9 (Mean Temperature of Driest Quarter) being the most important variable. Regardless of the metrics, as in our study, temperature variables appear to be the most important for predicting AGB (Table 1).

However, the relationship between bioclimatic variables and AGB is entirely dependent on the type of ecosystem, the species, and the region of the world. For this reason, other authors emphasize that mean annual precipitation (MAP) has a relatively greater importance (0.19%) than mean annual temperature (MAT) (0.05%) for predicting AGB in Larix plantations in northern and northeastern China. Furthermore, its importance is also dependent on the model structure [19].

Continuing with the previous narrative, the relationship between precipitation and AGB can be complex, as different responses can be observed depending on the type of forest and climatic conditions. This variable significantly influences the accumulation of different components of AGB (branches, stems, roots, and needles) in conifer plantations. Therefore, it is important to consider it in predictive models and in the evaluation of the climate-forest relationship [12]. However, climatic variables such as temperature and precipitation together can improve AGB estimates [21]. In tropical forests [51], 13 predictor variables of AGB were used, including geographical, topographical, hydrological, soil, and even species-related variables (coverage), etc. They found that the relative influence of MAP on AGB is 37.6%, making it the most important, while MAT has a relative influence of less than 1%.

4.2. Correlation Between Bioclimatic Variables and Carbon Density

The correlation of _cdAGB with temperature variables (Bio 5 and Bio 10) in temperate forests of Mexico is negative but positive with precipitation variables (Table 1). Similar findings were reported by [29], who observed a negative correlation (-0.46 > r < -0.63; -0.43 > r < -0.60) between the stem biomass of Larix gmelinii and Betula platyphylla with MAT. Likewise, other studies demonstrate that MAT is negatively correlated with AGB, either at the species level [12] or at the stand level [36].

Similar to this study, in temperate and dry forests, a positive correlation between _cdAGB and MAP has been observed [26]. This same relationship is found in other studies; for instance, [29] discovered that MAP positively correlates with the stem biomass of Larix gmelinii (0.84 > r < 0.92) and Betula platyphylla (0.76 > r < 0.88). Additionally, [12] observed that MAP positively and significantly (p < 0.05) correlates with the AGB of Pinus koraiensis Siebold & Zucc., Larix olgensis A. Henry, and Pinus sylvestris var. mongolica Litv.

Furthermore, the study conducted by [28] reveals that on a global scale (> 6200 forests and 61 countries), MAT is positively correlated with foliage biomass, although the geographical patterns of correlation are inconsistent. [14] found that while temperature correlates positively in boreal forests, the opposite occurs in tropical forests. Studies conducted by [4] in boreal forests demonstrate that MAT can explain (R²) between 26% and 45% of carbon density and is positively correlated; whereas MAP can explain between 28% and 67% of carbon density (both aboveground and belowground); however, the correlation between these two variables is positive when precipitation is between 0 to 1000 mm, and negative when it exceeds 1000 mm.

In contrast to this study, in boreal and temperate forests, a positive relationship has been found between _cdAGB and MAT [18], but a negative one in humid regions [26].

4.3. Predictive Capacity of Bioclimatic Models

Model validation is crucial for assessing the predictive capacity of a model based on new data. In reality, within this context (AGB-climatic predictors relationship), few studies undertake this process [19,47,48,52]. As known, when generating a model using the 'glm' procedure, statistics such as R², RMSE, MAE, etc., are not computed. Upon validating our models, we were able to calculate these metrics and evaluate the predictive capacity of each model. We observed that bioclimatic variables can explain up to 22% of _cdAGB in these types of forests (Table 3), a reasonably significant value at the eco-region level [4], but not at the species level [12,29], where variables can explain up to 84% of AGB at this scale.

In general, the 'cross-validation' technique has been the most commonly used to validate models for estimating AGB [19,36,52,53]; in this study, we used it along with three other techniques (Table 3). It is noteworthy that, when dealing with climatic predictors, these can explain around 20% (Table 3), whereas predicting AGB from predictors derived from vegetation indices, e.g., normalized difference vegetation index (NDVI), spectral information [52], or satellite optical images and unmanned aerial vehicles [53], can explain (R2) up to 80% of AGB. The inclusion of different predictor variables such as MAP, MAT, clay content, pH, dryness index, and stand age can explain up to 44.4% of _cdAGB in temperate forests [4].

However, exclusive stand predictors (diameter at breast height, age, stand density) can explain up to 98% of the ABG variance [19]. The RMSE metric is entirely dependent on the units and scale of the dependent variables, which is why such distant differences are observed across studies.

4.4. Current and Future Projection of _cdAGB

Under any climate scenario (RCP - year), our models predict losses in _cdAGB ranging from 5 Mg C ha^-1 (2050) to 20 Mg C ha^-1 by 2070 (Figure 3, Figure 4 and Figure 5) in Mexico's coniferous forests. Globally, positive changes in Total Carbon Density (TCD) are projected for temperate forests, averaging 2.23 Mg C ha⁻¹ (RCP 26, 45, and 85; 2050) and 1.99 Mg C ha⁻¹ by 2070. Specifically for Mexico, these authors [4] predict changes of ±20 Mg C ha⁻¹. It's important to note that, for the country, only two plots were considered in [4], whereas our study was conducted at a finer scale, generating models from n=48 to n=370 (Table 2).

A study conducted by [54] under various simulation scenarios, altering temperature and precipitation, in Yunnan Province, China, suggests that the combined effects of these variables are more complex than anticipated. They can result in both gains and losses (as observed in this study) in carbon sequestration across different forest types, attributed to decreased precipitation and increased temperature. Our study shows a decrease in _cdAGB, primarily due to a temperature increase of 1 to 3°C and a precipitation decrease of approximately 10%. Nearly three decades ago [55], the vulnerability of Mexico's forests to climate change was assessed, revealing that under these conditions (+2°C and -10% precipitation), the area of humid and dry temperate forests would significantly decrease (to less than half of their current size).

In a study conducted in the Brazilian Atlantic Forest (AF) [56], it was found that in 34.7% of the existing forest fragments, AGB could increase, while in 2.6% it could decrease by the year 2100. Models predict an 8.5% increase in total carbon stock; additionally, 76.9% of AF would be suitable for a potential increase in AGB by 2100 under an RCP 4.5, solely due to climate change. This contrasts with findings here, possibly due to the type of forests and geographical location, but is similar to what was found by [36] in subtropical evergreen forests in China, showing a decrease in AGB by 2050-2070, varying according to different climate scenarios (RCP 2.6 > RCP 4.5 > RCP 6.0 > RCP 8.5).

5. Conclusions

The relationships between _cdAGB and climate are more complex than they appear. Bioclimatic variables, especially those related to temperature, have proven to be significant predictors of _cdAGB in the developed models and can explain up to 19% of _cdAGB. Bio 5 (Max Temperature of Warmest Month) emerged as the most robust predictor. Predictions using bioclimatic models suggest that all strata of the country's coniferous forests, under any climate scenario (RCP - Year), will experience declines in _cdAGB (in the absence of anthropogenic activities) by 2050 and 2070, particularly under the RCP 8.5 scenario. It is expected that forests in stratum II will be the most affected, with significant reductions that could reach up to -20 Mg C ha^-1 by the year 2070. The temperature projected by the GCMs for 2050 and 2070 will be significantly higher than current levels, with increases of up to 3.55 °C in stratum I by 2070. Precipitation could remain unchanged in stratum I, decrease (-10%) in stratum II, and vary (±10%) in stratum III. These future climate projections are expected to lead to a redistribution of the _cdAGB in Mexican forests. The conclusions suggest an urgent need for adaptation and mitigation strategies to conserve biodiversity and ecosystem services provided by conifer forests in Mexico. Forest management should consider not only biodiversity conservation but also the carbon storage capacity in these ecosystems.