Rates of Stemwood Carbon Accumulation Are Linked to Hydroclimate Variability in Mexican Conifers

Abstract

Evapotranspiration demand has increased rapidly as temperatures have risen, affecting forest productivity. Consequently, carbon (C) uptake by forests is being modified; therefore, a more refined knowledge of the relationships between C capture and hydroclimate variability is required, particularly in drought-prone regions. In this study, we analyzed the relationships between climate and C capture as stemwood through the dendroecological analyses of radial growth in 15 conifer species distributed along a broad bioclimatic gradient in Mexico. The C content was calculated using densitometry data; correlations and mixed models were then used to determine the influence of climatic variables (precipitation, mean maximum and minimum temperatures, SPEI drought index) on tree growth. Each species showed specific responses to the climate with Taxodium mucronatum being the most responsive species. Both precipitation and maximum temperature best explained changes in C capture, with minimum temperature and SPEI playing secondary roles. The winter before the growth period was the most important season for C capture, particularly through positive responses to wet-cool conditions. However, the climatic influences of the current fall and summer were also notable. Seasonal climatic influences have implications for C uptake and forest productivity in the face of the severe droughts that repeatedly affect the study region.

1. Introduction

The frequency and intensity of drought periods have increased worldwide with more widespread and negative impacts on forests [1,2]. Drought affects tree respiration, photosynthesis, and growth rates negatively, therefore influencing forest productivity by decreasing the carbon (C) uptake rate [3]. Consequently, it has been a constant concern worldwide to improve our understanding of drought–growth relationships, particularly in the face of predicted climate warming and forecasted aridification [4]. Regarding radial growth, tree-ring variables such as width, density or blue intensity have been used to reconstruct C pools or fluxes and to assess climate impacts on wood production [2], but we still lack multi-species assessments in seasonally dry areas. Such comparisons would allow deciphering which species are more sensitive to more arid conditions and show lower or more variable long-term stemwood C accumulation rates as compared with more resilient species showing higher or more stable C long-term uptake rates.

Forest ecosystems are known to represent important C sinks that have been used in environmental services market policies to mitigate global warming (Red++). The scientific community has been concerned with providing scientific knowledge that highlights the sensitivity of trees to climatic variability. For example, [5] found that some tree species are vulnerable to rising maximum temperatures, while others appear to benefit from global warming and CO₂ fertilization [6]. Such divergent responses have been attributed to several processes related to biome type, site climate conditions, species functional traits, genetic load, etc. [7]. Moreover, current climate models indicate warmer conditions based on rising atmospheric CO₂ concentrations [8], which highlight the need to discern in a more refined manner the interaction between C uptake by trees and its relationships with hydroclimate variability based on the measurement of different tree-ring variables (e.g., width, wood density). Improved understanding of these mechanisms would allow the design and implementation of better management programs according to species’ climatic responses, as well as contribute to the understanding of C dynamics in different forest types and under contrasting climate conditions.

The existence of regions with diverse forests under varied climate conditions is a strategic opportunity for this purpose. Mexican forests are an example of diverse ecosystems, distributed in distinct and well-differentiated bioclimatic ranges [9]. In Mexico, there are several conifer species with high dendroecological potential: from arid to tropical environments, along gradients ranging from 0 up to more than 3000 m a.s.l. [10]. They represent varied soil, climatic, and physiographic features and can be used to reconstruct stemwood C uptake [11,12]. Many of these species and forests provide relevant ecological and commercial services [13]. Against such a perspective, these sites constitute ideal natural settings to test the influence of hydroclimate on the dynamics of C capture. Current methods for C estimation include the use of growth rings as natural archives that, through allometric relationships and using densitometric techniques, have proven their efficiency to model C fluxes and pools [14,15].

Although multiple efforts have been developed to estimate changes through time in C uptake, multispecies studies are limited [14,16]. The complexity is increased by the diverse sites tree species inhabit, as well as variations in the adaptive strategies of tree species [16]. Consequently, this requires a rigorous statistical process and the measurement of different tree-ring variables which allow inferring the dependences among C uptake, tree radial growth and hydroclimate. Models that include random effects, i.e., mixed models, are therefore an appropriate approach because they allow the determination of whether C fixation differs among trees subjected to similar hydroclimatic conditions. The distribution of residuals from this inference explains, to a large extent, the variation in C capture as related to regional hydroclimate. The development of successful climate change mitigation strategies will depend on this knowledge, given the forecasted scenarios of forest productivity loss in countries subjected to increasing drought stress such as Mexico [17].

The objective of this study is to assess changes in stemwood radial growth and wood density to reconstruct stemwood C accumulation rates. We compared data from 15 Mexican conifers growing under different climate conditions to identify the main climatic variables regulating C capture. We hypothesize that tree species will present different responses in their stemwood C uptake rates to hydroclimate variability, specifically to winter–spring water availability which mainly affects tree growth rates in the study area [4,5].

2. Materials and Methods

2.1. Study Area and Sampling Details

To encompass contrasting climatic conditions, 15 conifers located along bioclimatic gradients were sampled in Mexico, which is considered one of the most biodiverse regions of the planet [18] (Figure 1). The elevation at which these species were sampled ranged from 713 to 3500 m a.s.l. The dominant species in each stand occupied at least 70% of the basal area. From a random and centrally located tree, the 14 closest mature trees were selected to represent the stand, in contrast to the selective sampling of old trees growing in marginal conditions used to emphasize climate sensitivity (see [19]). Two cores were taken from each tree using a 5-mm Pressler increment borer at a height of 1.3 m above ground level, where the diameter at breast height was measured (DBH).

2.2. Tree-Ring Data

Once in the laboratory, wood samples were air dried, mounted, polished, and visually cross-dated following dendrochronological methods [20]. The visual cross-dating was checked using COFECHA software which calculates moving correlations between individual tree-ring width series and the site mean series [21]. Descriptive statistics related to growth data were calculated, specifically: the mean ring width and its standard deviation, the first-order autocorrelation (AR1), mean sensitivity (MS), and expressed population signal (EPS) (

Table 1. Descriptive statistics of the sampled species. Values are means ± SE.

Species (Code)	nc	Time Span	Tree-Ring Width (mm)	AR1	MS	EPS
Pinus caribaea var. hondurendis (Sénécl.) W.H. Barrett and Golfari (Pincar)	27	1982–2020	0.57 ± 0.05	0.31 ± 0.03	0.43 ± 0.01	0.88
Pinus oocarpa Schiede ex Schltdl. (Pinooc)	26	1986–2020	3.16 ± 0.20	0.47 ± 0.04	0.40 ± 0.01	0.89
Taxodium mucronatum Ten. (Taxmuc)	30	1987–2020	7.67 ± 0.25	0.63 ± 0.01	0.30 ± 0.01	0.98
Pinus jeffreyi Grev. and Balf. (Pinjef)	27	1975–2019	2.73 ± 0.16	0.51 ± 0.02	0.49 ± 0.01	0.95
Pinus cembroides subsp. lagunae (Rob. Pass.) D.K. Bailey (Pincem)	29	1970–2019	1.28 ± 0.07	0.36 ± 0.03	0.43 ± 0.01	0.86
Pinus engelmannii Carr. (Pineng)	30	1951–2019	2.13 ± 0.08	0.39 ± 0.03	0.45 ± 0.01	0.97
Picea martinezii T.F. Patterson (Picmar)	28	1862–2020	1.53 ± 0.10	0.70 ± 0.02	0.36 ± 0.01	0.92
Pinus maximartinezii Rzedowski (Pinmax))	29	1976–2020	2.41 ± 0.14	0.34 ± 0.04	0.48 ± 0.02	0.93
Pinus montezumae Lamb. (Pinmon)	29	1985–2019	2.01 ± 0.11	0.63 ± 0.04	0.42 ± 0.01	0.95
Pinus lumholtzii B.L. Rob and Fernald (Pinlum)	23	1917–2020	1.03 ± 0.09	0.55 ± 0.04	0.39 ± 0.01	0.93
Pinus patula Schltdl. et Cham. (Pinpat)	30	1952–2020	2.46 ± 0.09	0.68 ± 0.02	0.36 ± 0.01	0.96
Abies hickelii Flous and Gaussen. (Abihic)	29	1956–2020	2.84 ± 0.14	0.66 ± 0.04	0.34 ± 0.01	0.94
Juniperus deppeana Steud. (Jundep)	27	1986–2020	2.16 ± 0.19	0.39 ± 0.04	0.59 ± 0.02	0.95
Abies religiosa Kunth Schltdl. et Cham. (Abirel)	23	1960–2019	1.63 ± 0.13	0.35 ± 0.04	0.46 ± 0.01	0.94
Pinus hartwegii Lindl. (Pinhar)	20	1975–2019	1.64 ± 0.11	0.53 ± 0.04	0.35 ± 0.02	0.86

Abbreviations: nc, number of cores; AR1, first-order autocorrelation; MS, mean sensitivity; EPS, expressed population signal for the common timespan (1987–2019).

). The mean AR1 of ring widths measures year-to-year persistence in growth, while the mean MS of the ring-width accounts for relative changes in growth between consecutive years. The EPS is a measure of how coherent and replicated is a mean site ring-width series or chronology [22].

2.3. Series of Carbon Capture and Statistical Analysis

The two cores taken from each sampled tree were processed and measured in the tree-ring and wood anatomy laboratory at the University of São Paulo, Brazil. Initially, densitometry profiles were generated in an X-ray densitometry chamber (Faxitron X-Ray, Wheeling, IL, USA), with a measurement unit scale and a cellulose acetate wedge to calibrate wood density, according to established protocols [23,24]. Each radiographed sample was cross-dated to extract the annual density profile. Subsequently, the normal diameter was reconstructed from the total ring width (TRW) as an input for calculating the aboveground woody biomass (AWB) value for each species following [25] (Equation (1)). The AWB value allowed the annual C capture to be estimated (Equation (2)) considering the respective wood density (WD, in g cm^–3) and the fact that WD responds to seasonal changes in climate [26,27].

AWB = 0.0752 × DBH^2.4448 × 2.0331^WD

(1)

C = AWB × cf

(2)

where AWB is aboveground woody biomass, C is the annual carbon captured (in kg yr⁻¹), and cf is the carbon factor, obtained according to specific laboratory analyses [28].

To evaluate the influence of the climate on the annual C capture rates of the study species, Pearson correlation coefficients between C and monthly climatic variables were calculated considering a 0.05 significance level and using a window from previous Septembers to the current September. Prior to correlation analyses, we checked that the analyzed variables followed normal distributions. The climatic variables considered were the monthly mean maximum (Tmax) and minimum (Tmin) temperatures, and total precipitation (PP). The climatic data were obtained from the CRU TS 4.06 database at a 0.5° spatial resolution. We also calculated the Standardized Precipitation-Evapotranspiration Index (SPEI) to assess changes in drought intensity and duration [29]; We extracted data from the SPEI webpage (https://spei.csic.es/, accessed on 10 May 2023) also at a 0.5° resolution with a temporal scale of 1 to 48 months and from January to December.

To summarize the relationships between climatic variables and C accumulation rates, a principal component analysis (PCA) was calculated on the covariance matrix containing 585 values corresponding to the 39 Pearson correlation coefficients (13 for each climatic variable) of the 15 study species. The number of principal components was selected by examining a scree plot. Tree species were grouped according to their loadings on the first three principal components.

Linear mixed-effects models were used to analyze the relationship between the seasonal environmental variables and the carbon storage, as follows: winter PP (from the previous December to February of the growing season), and spring PP (from March to May), summer PP (from June to September), and autumn PP (from October to November) of the growing season for the last three seasons. The general form of the mixed-effects models is written as follows (Equations (3) and (4)):

$$y_i=X_i\beta +Z_i{b}_i+e_i$$

(3)

$$\mathrm{with} e_i\sim N\left(0,{\sigma }^2{I}_{n_i}\right) \mathrm{and} b_i\sim N\left(0,Q\right)$$

(4)

where y_i is the n_i-dimensional vector of responses, β is the p-dimensional vector of fixed parameters, X_i and Z_i are the n_i × p and n_i × (q + 1)-dimensional design matrices from unknown variables, b_i is the (q + 1)-dimensional vector of Gaussian random effects, e_i is the n_i-dimensional vector of error terms, and Q is a positive semidefinite (q + 1) × (q + 1)-covariance matrix [30].

In model fitting, we used the carbon content (C) as a dependent variable. The fixed effects (predictors with interaction) of the model were the monthly climatic variables (Tmax, Tmin, PP) for each site or species. We selected the individual tree as the random effect and used a variance function to correct for heteroscedasticity. In this case, we used a power-type function to model the variance, expressed as follows (Equation (5)):

$$Var\left({ϵ}_{ij}\right)={\sigma }^2|TR{W}_{ij}{|}^{2\delta }$$

(5)

where σ² and δ are the scale and shape parameters, respectively, and TRW_ij is a predictor (tree ring width) of the variance function [31] that could be assigned to the variance given by the fitted values [32]. All parameters of the mixed-effects models were estimated via restricted maximum likelihood (REML) with an nlme R-package, using the function lme [33].

3. Results

3.1. Growth and Carbon Accumulation Rates

The species with highest and lowest growth rates were T. mucronatum (7.67 mm) and P. caribaea (0.57 mm), respectively (Table 1). The mean growth rate of all species was 2.35 mm. The species with highest and lowest AR1 values were P. martinezii (0.70) and P. caribaea (0.31), respectively. The species with highest and lowest MS values were J. deppeana (0.59) and T. mucronatum (0.30), respectively. The interval common to all species, i.e., with well-replicated tree-ring width series, was 1987–2019. All species’ chronologies showed EPS values above the 0.85 threshold during this common timespan, indicating that they showed coherent growth signals. We found notable differences between species in C accumulation over time (Figure A1). For example, T. mucronatum stood out above the rest in this regard with a mean value of 19.64 kg yr⁻¹, while P. cembroides presented the lowest value (1.58 kg yr⁻¹).

3.2. Responses of Carbon Accumulation Rate to Climate Variability

The P. oocarpa C accumulation series positively correlated with PP in the previous winter and current summer, while P. engelmannii was notable for its positive relationship (r = 0.42) with PP in September (Figure 2). Rainfall in the previous winter negatively affected P. montezumae (e.g., February, r = −0.46), P. caribaea, T. mucrunatum, and A. hickelii mainly, while P. montezumae and A. hickelii negatively responded to April PP (r = −0.36 to −0.40).

The Tmax in the previous winter positively correlated with C capture in P. cembroides, P. engelmannii, Picea martinezii, P. montezumae (e.g., February, r = 0.45), P. maximartinezi, P. lumholtzii, and J. deppeana (Figure 2). The same response was found for spring and summer Tmax for P. engelmannii, P. montezumae, P. patula, P. cembroides, P. martinezii (April, r = 0.38) and, particularly, for T. mucrunatum (e.g., August, r = 0.44). The summer Tmax negatively affected J. deppeana and P. jeffreyi, while P. engelmannii (September, r = −0.41) and P. montezumae were influenced by early-autumn Tmax.

The Tmin during the previous winter played a positive role in C capture in species such as A. hickelii, J. deppeana, P. jeffreyii, and P. montezumae (February, r = 0.40), as well as in spring-summer for P. engelmannii, P. jeffreyii, and P. martinezii (Figure 2). Low summer Tmin negatively affected J. deppeana, T. mucrunatum (e.g., July, r = −0.38) and P. jeffreyi.

Finally, the calculation of the PCA summarizing C-climate relationships led to three principal components accounting for 27.9%, 16.7% and 12.3% of the total variance. The species’ loadings most important for the first principal component corresponded to P. montezumae and P. maximartinezii (

Table 2. PCA loadings in the first (PC1), second (PC2) and third (PC3) principal components summarizing the relationships between climatic variables and C accumulation rates in the 15 study tree species. See species’ codes in Table 1.

Species	PC1	PC2	PC3
Pincar	−0.022	0.370	−0.053
Pinooc	0.022	−0.104	−0.104
Taxmuc	−0.010	0.788	−0.192
Pinjef	−0.124	0.083	0.552
Pincem	0.234	0.066	−0.431
Pineng	0.347	0.073	−0.087
Picmar	0.340	0.140	−0.149
Pinmax	0.401	−0.089	−0.108
Pinmon	0.604	−0.163	0.237
Pinlum	0.116	−0.028	0.053
Pinpat	0.222	−0.017	0.052
Abihic	0.095	0.296	0.337
Jundep	0.217	0.189	0.488
Abirel	0.199	−0.027	−0.053
Pinhar	−0.093	−0.188	−0.046

). In the case of the second principal component, the highest loadings corresponded to T. mucrunatum and P. caribaea, while in the third principal component the highest loadings were found for P. jeffreyii and J. deppeana, with P. cembroides showing the lowest loadings.

3.3. Seasonal Influences on C Accumulation

The behavior of the model residuals was improved by incorporating the variance model (Figure A2). In general, it was found that the fixed effects with interaction had different results for species, seasons, and climate variables, which indicates that C accumulation responses to climate differ between tree species (

Table A1. Estimates of the fixed parameters and covariance parameters of the mixed linear models for tree-ring carbon accumulation in the fifteen species studied. Only significant models (p < 0.05) are presented.

Season	Sign	Climate Variable	Species	Estimate	SE	t	p
Winter	Negative	Tmin	Pinus oocarpa	−1.5368	0.7296	−2.1062	0.0352
			Picea martinezii	−0.1308	0.0323	−4.0492	0.0001
			Pinus cembroides	−0.1228	0.0551	−2.2277	0.0259
			Pinus hartwegii	−0.3016	0.1161	−2.5969	0.0094
		Tmax	Abies hickelii	−0.6911	0.1883	−3.6698	0.0002
			Pinus hartwegii	−0.5182	0.091	−5.6959	0
			Pinus jeffreyi	−0.3339	0.1106	−3.0193	0.0025
			Abies religiosa	−0.18	0.0843	−2.135	0.0328
		PP	Taxodium mucronatum	−0.0343	0.0159	−2.1612	0.0307
			Abies hickelii	−0.0171	0.0035	−4.8889	0
			Pinus hartwegii	−0.0149	0.0037	−4.065	0
			Pinus oocarpa	−0.0058	0.0025	−2.3059	0.0211
			Pinus jeffreyi	−0.0032	0.0014	−2.3868	0.017
	Positive	Tmin	Abies religiosa	0.2212	0.1049	2.1075	0.0351
			Pinus jeffreyi	0.5046	0.1225	4.1211	0
			Pinus lumholtzii	0.7868	0.1446	5.443	0
			Taxodium mucronatum	1.0204	0.4261	2.3945	0.0167
			Abies hickelii	1.6024	0.2587	6.1943	0
		Tmax	Pinus lumholtzii	1.5527	0.108	14.3716	0
			Picea martinezii	0.1705	0.0296	5.7669	0
			Pinus cembroides	0.2214	0.061	3.632	0.0003
			Pinus montezumae	0.4707	0.2022	2.3283	0.0199
		PP	Pinus cembroides	0.0046	0.0016	2.8948	0.0038
			Picea martinezii	0.0058	0.0012	4.9885	0
			Pinus patula	0.0093	0.0034	2.7421	0.0061
			Pinus lumholtzii	0.0145	0.0049	2.9258	0.0034
Spring	Negative	Tmin	Pinus patula	−0.2487	0.1225	−2.0313	0.0422
			Abies hickelii	−1.1976	0.2997	−3.996	0.0001
		Tmax	Pinus lumholtzii	−1.1581	0.1263	−9.1675	0
			Taxodium mucronatum	−1.6198	0.5184	−3.1244	0.0018
			Pinus cembroides	−0.1741	0.0605	−2.8753	0.004
		PP	Taxodium mucronatum	−0.059	0.0163	−3.6169	0.0003
			Pinus lumholtzii	−0.0071	0.0033	−2.149	0.0316
			Pinus maximartinezii	−0.0051	0.0023	−2.2416	0.025
	Positive	Tmin	Pinus hartwegii	0.6474	0.1621	3.9931	0.0001
			Pinus maximartinezii	0.4826	0.1116	4.3254	0
			Pinus cembroides	0.2111	0.0734	2.8766	0.004
			Pinus jeffreyi	0.3072	0.1461	2.1024	0.0355
		Tmax	Abies religiosa	0.3252	0.1057	3.0769	0.0021
			Pinus engelmannii	0.2203	0.0767	2.8719	0.0041
			Picea martinezii	0.1044	0.0428	2.4421	0.0146
			Abies hickelii	0.5974	0.2391	2.4985	0.0125
		PP	Abies religiosa	0.0063	0.0025	2.5159	0.0119
			Pinus engelmannii	0.0063	0.003	2.0743	0.0381
			Pinus montezumae	0.0273	0.0099	2.7596	0.0058
Summer	Negative	Tmin	Taxodium mucronatum	−4.4819	0.7411	−6.0481	0
			Pinus montezumae	−1.791	0.4078	−4.3919	0
			Pinus lumholtzii	−0.6931	0.236	−2.9364	0.0033
			Pinus jeffreyi	−0.449	0.1686	−2.6628	0.0078
		PP	Pinus hartwegii	−0.0017	0.0007	−2.3339	0.0196
	Positive	Tmin	Pinus engelmannii	0.7208	0.135	5.3374	0
		Tmax	Pinus cembroides	0.1858	0.0763	2.4353	0.0149
			Pinus lumholtzii	0.3907	0.1519	2.5722	0.0101
			Pinus maximartinezii	0.5895	0.1069	5.5144	0
			Pinus jeffreyi	0.7097	0.1415	5.0162	0
			Taxodium mucronatum	2.784	0.5612	4.9607	0
			Pinus montezumae	0.8204	0.4037	2.0325	0.0421
		PP	Pinus oocarpa	0.0022	0.0009	2.4589	0.0139
			Abies religiosa	0.0027	0.0012	2.3231	0.0202
			Pinus engelmannii	0.004	0.0008	5.1523	0
			Pinus lumholtzii	0.0041	0.0011	3.6587	0.0003
			Taxodium mucronatum	0.0119	0.0052	2.2956	0.0217
Fall	Negative	Tmin	Picea martinezii	−0.1043	0.0435	−2.3945	0.0167
			Taxodium mucronatum	−2.0539	0.487	−4.2174	0
		Tmax	Pinus maximartinezii	−0.3713	0.0773	−4.8055	0
	Positive	Tmin	Pinus hartwegii	0.5715	0.1372	4.1665	0
		Tmax	Picea martinezii	0.0822	0.0389	2.1136	0.0346
			Pinus cembroides	0.1456	0.0738	1.9738	0.0484
			Abies hickelii	0.5509	0.2499	2.2049	0.0275
			Taxodium mucronatum	2.0133	0.4782	4.21	0
			Pinus hartwegii	0.8879	0.1062	8.3609	0
		PP	Pinus maximartinezii	−0.0034	0.001	−3.3869	0.0007
			Abies hickelii	−0.0029	0.0012	−2.4512	0.0142
			Picea martinezii	0.0012	0.0003	4.0476	0.0001
			Pinus cembroides	0.0013	0.0004	3.7593	0.0002
			Pinus patula	0.0037	0.001	3.5108	0.0004
			Pinus engelmannii	0.0051	0.0013	3.7902	0.0002
			Pinus lumholtzii	0.0057	0.0015	3.833	0.0001
			Taxodium mucronatum	0.0399	0.0063	6.3051	0

).

In spring, precipitation had a significantly positive effect on C capture in species such as P. montezumae, A. religiosa, and P. engelmannii, while for T. mucronatum, P. maximartinezii, and P. lumhotlzii, the effect of precipitation was negative. The C capture in P. lumholtzii, T. mucronatum, and P. cembroides subs. lagunae was negatively affected by Tmax, while for A. religiosa, P. engelmannii, A. hickelii, and P. martinezii, C capture showed a positive relationship with Tmax. The Tmin had a positive effect on C capture in P. maximartinezii, P. hartwegii, P. cembroides subs. lagunae, and P. jeffreyi, while A. hickelii and P. patula were negatively affected.

In summer, PP had a positive influence on C capture for P. maximartinezii, P. lumholtzii, P. oocarpa, A. Religiosa, and T. mucronatum, while only P. hartwegii presented a negative influence on the accumulation of C. Regarding Tmax, P. jeffreyi, P. maximartinezii, T. mucronatum, P. lumholtzii, P. cembroides subs. lagunae, and P. montezumae had a significantly positive influence on C accumulation, while in the case of Tmin T. mucronatum, P. montezumae, P. lumholtzii, and P. jeffreyi were negatively affected in terms of C storage, while P. engelmannii showed a positive influence of Tmin on C accumulation.

During autumn, PP positively affected T. mucronatum, P. lumholtzii, P. engelmannii, P. martinezii, P. cembroides subs. lagunae, and P. patula, while P. maximartinezii and A. hickelii were negatively affected. The Tmax had a positive effect on T. mucronatum, P. hartwegii, A. hickelii, P. martinezii, and P. cembroides subs. lagunae, and only P. maximartinezii was negatively affected. On the other hand, Tmin showed a positive effect on P. hartwegii; however, T. mucronatum and P. martinezii were negatively impacted.

In winter, PP presented a significantly negative influence on the accumulation of C in P. hartwegii, A. hickelii, P. jeffreyi, P. oocarpa, and T. mucronatum. However, for P. martinezii, P. lumholtzii, P. cembroides subs. lagunae, and P. patula, the effect of PP was positive. The Tmax positively affected C capture in P. lumholtzii, P. martinezii, P. cembroides subs. lagunae, and P. montezumae, while the species that presented negative influences were P. hartwegii, A: hickelii, P. jeffreyi, and A. religiosa. Finally, C capture was also negatively affected by Tmin in P. martinezii, P. hartwegii, P. cembroides subs. lagunae, and P. oocarpa, whereas P. jeffreyi, P. lumholtzii, A. hickelii, T. mucronatum, and A. religiosa presented positive values.

3.4. Carbon–SPEI Associations

Regarding associations between C accumulation rates and the SPEI drought index, we found positive relationships in some species (e.g., P. oocarpa, P. lumholtzii, P. hartwegii, A. religiosa) and negative relationships in other species (e.g., P. caribaea, T. mucronatum, P. martinezii, P. montezumae, A. hickelii, J. deppeana) (Figure A3). The two species with the highest and lowest C-SPEI correlations were P. hartwegii and P. montezumae at mid (6–9 months) to short (1–5 months) temporal scales, respectively. In both cases, these correlations were found from the previous winter to current spring.

4. Discussion

Our analyses improve the understanding of drought–growth relationships, in the context of forecasted climate warming, and how different impacts on radial growth and wood density would modify stemwood C accumulation rates under more arid conditions. The inter-annual ring width and wood density values were suitable proxies for modeling the temporal accumulation of C in Mexican conifers. The variability in C capture was influenced by the hydroclimate in this multi-species study, which is an unprecedented finding for Mexican forestry and tree-ring sciences.

Temporal rates of C capture among species differed, which is attributed to the adaptive discrepancies specific to their particular environments [16,34]. For example, the differences in C-accumulation rates between T. mucrunatum and P. cembroides could be directly related to the different soil moisture of the sites where they grow [35]. This highlights the importance of discerning the relevant role played by hydroclimatic variability as a driver of C accumulation, which is consistent with previous studies [36]. However, a greater radial growth does not necessarily imply C gain [37]. Environments with high moisture availability favor greater radial growth [38], whereas, in arid environments, water limitation reduced growth but may increase wood density and, potentially, C uptake [14,39].

The variable associations between monthly climate variables and C can be attributed to the adaptive conditions of each species, microsite conditions, plasticity, ontogenic trends, edaphic traits, topographic position, stand structure and forest type [36]. These findings are in line with those of previous studies such as [40], which argue that tree species develop specific response mechanisms to cope with hydroclimatic variability in drought-prone sites.

The combined influence of precipitation (PP) and mean maximum temperature (Tmax) in promoting C capture is evident. The effect of available soil water associated with high temperatures is reflected in the physiology of plants which show a more active cambium cell division and tracheid enlargement [41]. In spring and winter, the fitted models indicate positive influences on A. religiosa and P. engelmannii, whereas this occurs during summer in the case of P. maximartinezii, T. mucronatum, P. lumholtzii, P. martinezii, P. lumholtzii, and P. cembroides subs. lagunae. This concurs with other studies that reported seasonal humidity prior to the onset of growth as a driver of growth, perhaps by recharging soil moisture, while the summer climate influences the thickening of the tracheid cell walls being promoted by high temperatures [42]. In autumn, T. mucronatum, P. martinezii, and P. cembroides also benefitted from warm and wet conditions, which are related to the monsoon summer rains and delayed dormancy [43].

However, the humid and warm conditions also had a negative effect on C capture, as was the case in T. mucronatum and P. lumholtzii in spring, P. maximartinezii in autumn and P. hartwegii, A. hickelii, P. jeffreyi and P. oocarpa in winter. Previous studies have confirmed how hydroclimatic regimes have shifted and impacted tree growth [42]. Trees in neighboring sites have shown bimodality in growth [44], providing strong evidence of the extension of the growing season and changes in cambium phenology.

The minimum temperature (Tmin) plays a significant role in changes in wood density and tracheid maturation in latewood [45], which has implications for C accumulation rates. In contrast, Tmax reduces photosynthetic rates due to stress related to heatwaves and soil drying, which typically occur during drought events [46]. These dry conditions often reduce C uptake due to the decrease in available non-structural carbohydrates leading to carryover or legacy effects [47].

Regarding species behavior, T. mucronatum was very responsive throughout the year, which supports previous studies that report the species as highly sensitive to hydroclimate variability [48]. The summer and autumn PP and Tmax benefitted its C accumulation rate. Despite high summer evapotranspiration rates, this species seems to take advantage of its proximity to water bodies. In recent years, some tree species have been shown to present a bimodal growth pattern [49], which may also prolong their season of carbohydrate synthesis and wood formation (see [4]). Interestingly, spring PP and Tmax negatively affected C accumulation in this species, which can be attributed to cloudy conditions or increased evaporative demand [36].

The combined effect of rainfall and high temperatures at the beginning of the growing season may stress trees through the negative impact of elevated evapotranspiration rates [50,51]. The effects of increased temperatures with severe precipitation deficits have caused forest dieback and tree mortality due to high transpiration rates and hydraulic failure [1]. In this way, low biomass growth in drought-prone sites could reduce the competitive ability of some species to uptake C [14]. If this occurs, not only will net C gain be reduced, but carbon starvation may even occur and contribute to dieback. Interestingly, our findings are novel because drought-amplified stress may result in increased C pools [52]. In contrast, cold winter temperatures favor C accumulation by promoting the formation of more tracheids with thicker walls, whereas it negatively affects C accumulation by limiting tracheid expansion.

Three pines (P. maximartinezii, P. cembroides, and P. lumholtzii) were the species that showed significant relationships with seasonal climatic variation during three seasons. In all these seasons, the combination of PP and Tmax was the main driver of growth and its effect depended on the considered species. Overall, the magnitude of the statistical significance of the climatic influence was most notable during autumn, followed by summer (Figure 2). These results highlight the fact that the growing season may be extended after summer. According to [53], monsoon rains in conjunction with warm temperatures cause seasonal adjustments in the mechanisms of C synthesis and hydraulic functioning. Forecasted changes in the hydroclimate regimes will therefore have an impact on C fluxes, at least for these three species. Our results concur with those of previous studies that report the dependence of C production on thermal conditions and atmospheric demand in summer and autumn [36].

The capture of C in P. martinezii, P. engelmannii, A. religiosa, and P. jeffreyi was sensitive to summer and autumn hydroclimate conditions (combined influence of PP and Tmax, followed by Tmin). According to several authors, minimum temperatures regulate C allocation in tree organs (e.g., the stem, branches, shoots, leaves, and roots), which could explain this pattern [54].

The C accumulation in species such as J. deppeana and P. caribaea was less influenced by the hydroclimate (see Table A1). Under management scenarios, these species seem to be less affected by hydroclimate variability, although their resilience remains to be evaluated.

The mixed models improved the variance of climatic factors, but we must acknowledge that other variables should be incorporated into future studies, including genetic load, stand structure, site topography, soil features, disturbances, and xylogenesis data. A spatial analysis would refine the analyses of hydroclimate influences on C accumulation and allow identifying geographical patterns at the level of ecoregions and biomes [55]. Other limitations of our study are that we reconstructed the radial growth and aboveground woody biomass using general allometric equations, but we lacked detailed data on changes in belowground biomass and species-specific allometric relationships developed for the study sites. Our database nevertheless represents a solid input for adapting new algorithms relating changes in tree biomass and climate (see [16,56,57]) and, if necessary, replicating the analyses in other bioclimatic zones. Another important aspect would be to explore the use of other models such as generalized additive models, which would allow a better understanding of nonlinear influences of the hydroclimate on C capture.

5. Conclusions

There are differences in the C accumulation rate of the 15 conifer species studied across Mexico. Our findings reveal the critical role of the climate in the response of C accumulation. The combination of seasonal PP and Tmax is the main driver of C accumulation, suggesting a positive role in C capture of an extended growing season. The season to which most species respond is winter; however, the magnitude of the climatic influence of PP and Tmax on C capture is greatest during autumn, followed by summer. This means that the hydroclimate changes in autumn, and possibly late summer, are modulating stemwood C accumulation through changes in wood density and tracheid cell-wall thickness. Thus, droughts play a prominent role in the dynamics of C fluxes in the studied species. The modeling strategy allowed us to identify the seasonal climatic factors that affect the dynamics of C capture. This knowledge is crucial in the face of more frequent and severe climate extremes, such as hot droughts, to develop new strategies to improve forest resilience and C uptake.