This study was carried out on two mountains in central Mexico: “La Malinche” (LM) and “Nevado de Toluca” (NT) ( Figure 1 ). LM is located between the states of Tlaxcala and Puebla and has a maximum elevation of 4,461 m (Acosta-Pérez and Kong, 1991). NT is located in Mexico State, between the municipalities of Toluca and Tenango del Valle, and has a maximum elevation of 4,680 m (Körner and Paulsen, 2004). Both mountains belong to the Trans-Mexican Volcanic Belt, however, there is a distance between the mountains of 181.6 km. On both mountains there are widely different conditions between the lower and upper extremes of the elevational gradient. The data from the model ClimateNA (Wang et al., 2016) indicate that the mean annual temperature (MAT) of LM is 9.7°C at the lower end of the elevation gradient ( Table 1 ) and 6.6°C at the upper end. The mean annual precipitation (MAP) at the lower end of LM it is 1208.2 mm and at the upper end is 1132.1 mm. On NT, the MAT is 9.9°C at the lower end and 5.5°C at the upper, and MAP is 1458.3 mm at the lower end and 1299.1 mm at the upper ( Figure 1 ). In addition to P. hartwegii, which is the species with the highest elevation distribution range (3500 - 4200), other tree species can be found on both mountains, such as P. montezumae Lamb. (3,000 - 3,200 m) and Abies religiosa Kunth Schltdl. & Cham. (2,800 - 3,400 m); as well as grasslands, alpine zacatonal, and highland paramo (4,000 - 4,400 m) (CONANP, 2013a; CONANP, 2013b).

On each mountain, sampling plots were selected at the lower and upper limits of the elevational distribution gradient of P. hartwegii. Three sites were sampled at each elevation on LM and two per elevation on NT ( Table 1 ). The sites were selected based on the criteria of being low competition stands of mature P. hartwegii trees, without evidence of disturbance due to logging, wildfire, or pests. At each site, we selected between 25 and 31 dominant adult trees with a straight trunk, over 40 cm in diameter at breast height. For each selected tree, we extracted a wood core from the bark to the center of the trunk using a 5-mm diameter increment borer (HAGLOF ^®, Suiza) at 1.30 m height above the ground. We sampled a total of 274 trees.

The collected increment cores were air-dried at room temperature and longitudinally sawn to obtain 1.57 mm thick transverse sections. Then, the resin was extracted with a bath in a water-pentane solution (2:1) for 48 hours. Each sample was exposed to a source of X-rays for 25 min in the wood densitometry laboratory of the Institut National de Recherche sur L’agriculture, L’alimentation et L’environnement, INRAe Val de Loire Center, Orleans, France. This was based on the procedure developed by Polge (1978) and described by Mothe et al. (1998). Then, microdensity profiles were obtained by scanning the X-ray images at a resolution of 4,000 dpi with the software WinDendro, Regent Instruments Inc. (Guay et al., 1992). This software converts the gray levels of the pixels detected in the X-rays image to density values, calibrated with a standard of known physical and optical density (Schweingruber, 1996).

The microdensity profiles obtained from the growth cores were imported into R (“R Development Core Team”, 2018). Then, we calculated the following traits related to the width and density of the growth rings: overall ring width (RW: mm), earlywood ring width (EWW: mm) latewood ring width (LWW: mm), overall mean ring density (RD: g cm^-3), earlywood ring density (EWD: g cm^-3), latewood ring density (LWD: g cm^-3), maximum density (MAXD: g cm^-3) and minimum density (MIND: g cm^-3). This was done using R functions designed to obtain these variables (Rozenberg et al., 2002).

Cambial age influences width of growth rings and in general the properties of wood, especially during the juvenile stage (the first 10 to 25 rings). Although it is a natural effect of tree growth, this effect can introduce bias in the subsequent statistical analyses. We therefore removed the effect of cambial age by adjusting the raw data using the regional curve standardization (RCS) procedure, in the four elevation plots, as described below ( Figure 2 ) (Esper et al., 2003).

Following Esper et al. (2003) and Rozenberg et al. (2020), we compared three methods to eliminate the effect of cambial age: raw data (no adjustment), residual RCS, and RCS ratio. Of these, the residual methods had the best fit to eliminate the effect of cambial age; we obtained one age-related curve for the lower end and another for the upper end for each of the mountains (LM and NT) ( Figure 2 ). In each of the populations, we calculated the expressed population signal (EPS) statistic, which serves to determine the statistical quality of the chronology. The EPS value has a range of 0 to 1.0, where a value of 0.85 is considered acceptable (Buras, 2017). In all of our study sites, the EPS for RW was 0.88 or higher, above the reference point used to evaluate the quality of the chronologies. The residual chronologies and the time-dependent RCS ratios were obtained from the means of individual series. For the above analyses, we used the dendrochronology library for R, dplR (Bunn, 2008).

The climate variables for each sampling location were obtained with the software ClimateNA model v6.30 (Wang et al., 2016), which includes monthly data as well as annual mean (MAT), maximum (T_MAX) and minimum (T_MIN) temperature and total annual precipitation (TAP) for the 1964-2016 period. We used the monthly temperature data to calculate the mean temperature of the growing season (T_GS; months of April through September), spring mean temperature (Tspring; months of April through June), and winter mean temperature (Twinter; months of January through March), The annual aridity index (AR_I) was calculated with the equation AR_I = TAP/(MAT+10) (Fniguire et al., 2014). To verify that the data from the ClimateNA software accurately reflected the conditions at the study locations, we calculated the correlation between the ClimateNA data, and the data recorded by local meteorological stations. The meteorological stations were “Amaxac de Guerrero” (No. 29042; 3,320 m of elevation) and “Acxotla del Monte” (No. 29161; at 2,443 m of elevation) for LM and “Tenango” (No. 15122; at 2,858 m) and “Nevado de Toluca” (No. 15062; at 4,283 m) meteorological stations for NT. In all cases, there were very strong correlations between MAT and TAP values (r > 0.97), except in the case of TAP with meteorological station No. 29161, which had a moderate correlation (r = 0.59).

All the statistical analyses were done using the adjusted data for the growth ring traits for the 1964–2016 period. To reduce the number of variables to analyze (RW, EWW, LWW, RD, EWD, LWD, MAXD, MIND), we estimate the Pearson correlation coefficients among the growth ring traits, using the mean values per sampling site per year (showed in annexes). Pairs of variables that were strongly correlated (r > 0.9) were consolidated by selecting only one representative variable; this resulted in the selection of five variables for subsequent analyses: RW, LWW, RD, EWD and MAXD.

The RN for each ring trait to each climate variable was calculated for each individual tree based on the linear model Y_i = β₁ + β₂X_i + ϵ_ij, where Y_i is the phenotypic trait (ring trait); X_i is the climate variable, β₁ is the intercept, β₂ is the slope (the PP value), and ϵ_ij is the model residuals. The model was fit using the lm function of R, from the R stats package (R Core Team, 2018). For each RN, we obtained the probability value (p) as well as the adjusted R² value as an indicator of the goodness of fit of the linear model. When p < 0.05 we considered that the slope of the RN was significantly different from zero, demonstrating that the ring trait showed PP, and the regression coefficient (β₂) was a quantitative estimate of the amount of PP of the ring trait to the climate variable for the corresponding tree. For the following statistical analyses, we only worked with the climatic variables: “T_GS” and “Twinter”, because they were the ones with which P. hartwegii showed greater PP through the traits of its growth rings. It should also be noted that a previous research (Carrillo-Arizmendi et al., 2022b) carried out correlations between these two climatic variables indicating that there is a relatively weak relationship (0.460 ≤ |r| ≤ 0.619), such that these two variables are not providing redundant information, so we consider it is appropriate to include both in order to compare the plastic response of trees to these two climate variables and determine which one might be more reliable across populations. The relative frequency of trees with significant RN (RN_S) for each growth ring trait-climate variable was calculated for each mountain/elevation combination (hereafter, “population”), to compare the overall “sensitivity” in PP of ring traits to different climate variables at the lower and upper end of the elevation gradient on both mountains. In addition, to identify possible relationships and/or compensatory effects on PP of growth ring traits, a correlation coefficient analysis and a principal component analysis (PCA) of the PP values for the ring traits were carried out. We also calculated the proportion of trees with positive and negative slope of RN_S for each ring trait. With these data, we determined whether there were differences in the proportion of trees with positive and negative RN_S between the two ends of the elevation gradient (lower and upper) and between mountains (LM and NT). The Glimmix procedure in SAS (SAS 9.4, 2019), with the link function for binary variables (proportions) was used for this purpose. A two-way ANOVA (SAS 9.4) was also done to evaluate the effect of elevation, mountain, and their interaction on the average slope value (β₂) of the positive and negative RN_S for each ring trait and climate variable to determine whether the populations differed in the magnitude of PP.

Growth ring traits were more “sensitive” to the climate variables T_GS and Twinter; that is, they had the overall largest number of significant RN with these two climate variables ( Table 2 ) and a high percentage of trees with significant RN for more than one ring trait ( Figure 3 ). Overall, the percentage of trees with significant RN for a growth ring trait was around 50% for both T_GS and Twinter across elevations and mountains, except for the RN related with Twinter at low elevation in LM mountain, were it was only 44.3% ( Table 2 ). None of the other climate variables had similar percentages of significant RN across all populations. For instance, there was a high percentage of significant RN with mean annual temperature (MAT) in three populations but not at high elevation in LM. T_MAX was also important in three populations, but not at low elevation in LM and TAP was important only at the upper end in NT mountain. Given that T_GS and Twinter had the highest proportion of trees with significant RN, the remaining analyses were done only for these two climate variables. Looking at individual growth ring traits, the proportion of trees with significant RN to T_GS and Twinter differed only for RW, between mountains at the lower end, with a lower percentage on LM mountain ( Table 2 ). Despite these differences among populations, the overall proportion of trees with significant RN in response to T_GS and Twinter was high on both mountains; over 70% of trees had significant RN for at least one growth ring trait, and over 40% of trees in all populations showed significant RN for three or more growth ring traits ( Figure 3 ).

The correlation matrix and the first two PC obtained from the PCA showed that for both climatic variables (T_GS and Twinter), the PP of growth ring traits was separated in three distinctive groups: ring and latewood width in one group, ring and earlywood density in other, and maximum density in the third group ( Figure 4 ). In the correlation matrix of PP values associated with both T_GS and Twinter, RW and LWW were strongly correlated, as well as RD and EWD.In general, the PP of growth ring traits in response to the two climate variables (T_GS and Twinter) showed a similar structure. Thus, the PP of each ring variable in response to T_GS was strongly correlated (r > 0.87) with the PP of the same variables in response to Twinter. The exception was LWW, where the correlation was slightly weaker (r = 0.76).

A wide variation in PP was found for growth ring traits (RW, RD, and MAXD) across populations in both T_GS and Twinter-related NR ( Figure 5 ). Even though significant RN ranged from negative to positive slopes for all traits, negative trends were more common for RW and positive trends for RD and MAXD, except for MAXD at low elevation in NT, where negative trends were predominant. In addition, there were significant elevation and mountain effects on both the proportion of trees with negative and positive RN and the mean absolute PP value for a given growth ring trait ( Figure 5 ). The proportion of trees with negative trends in their RN generally decreased from the lower to the upper end for all traits in both mountains, although the effect of elevation was more pronounced for RW and MAXD at the NT mountain for both climate variables ( Figure 5 ). Overall, mean absolute PP values were higher on LM than on NT for both positive (in RD) and negative (in RW and RD) RN ( Figure 5 ). Also, at the lower end of the elevation gradient mean absolute PP values were higher than at the upper end for both positive (in RD) and negative (RW) RN. The effect of elevation on the absolute PP value of negative RN for RW was larger at the LM mountain ( Figure 5 ). Mean absolute PP values for MAXD showed a similar trend, but differences between mountains and elevations were not significant.

Overall, radial growth of P. hartwegii trees was more sensitive to temperature than to precipitation. This coincides with a previous report by Soto-Carrasco (2021) who analyzed the PP of P. hartwegii on three other mountains in Mexico (Mount Tláloc, Pico de Orizaba, and Cofre de Perote). In that study, precipitation was the climate variable with the fewest significant RN for growth ring traits. Liu et al. (2015) mentioned that high-mountain trees generally have low sensitivity to precipitation, although they did not propose any explanation for this finding. Hendrickson et al. (2004) and Storkey (2004) mention that, under high mountain conditions, temperature is the main factor influencing most physiological processes and phenotypic changes in plants; precipitation seems to be a limiting factor only for some processes associated with xylogenesis, particularly during extreme drought events (Pompa-García and Venegas-González, 2016; Camarero and Gutiérrez, 2017; Correa-Díaz et al., 2020; Pompa-García et al., 2021).

In our study, the differences in tree sensitivity to climatic variables could be also due to differences in the degree of micro-spatial variation in temperature versus precipitation at the study locations. Ambient temperature is more homogeneous within a single location, associated mainly with site elevation, such that the value recorded at the locality level better reflects the growing conditions of all trees on it, while rainfall (and soil moisture) is very heterogeneous over short distances (Gómez-Plaza et al., 2000; Wundram et al., 2010; Yang et al., 2018). Therefore, total annual precipitation values recorded at the site level do not necessarily reflect the moisture conditions for each tree, potentially making it more difficult to detect significant RN for them. It is worth mentioning that, although our analyses do not consider the possibility of analyzing multi-annual or lagged effects, it would be interesting to do so in future research and thus learn how individual trees within a population differ in the reaction norms of ring growth traits in response to interannual fluctuation of climatic variables, especially precipitation.

On the other hand, the sensitivity of ring growth traits to the inter-annual climate fluctuations also shows differences between mountains. There was a higher overall proportion of significant RN on NT than on LM ( Table 2 ). The difference between mountains could be due to differences in the range or variation coefficient of climate variables, differences in age of trees, or to other site-level factors affecting the response of radial growth to fluctuation of climate variables. Although the variation coefficient of TAP was similar between the two mountains (14%), the average value of TAP was higher on NT and there was a larger difference between the lower and higher elevations (1458.3 mm and 1299.1 mm, respectively, compared to 1132.1 mm at the lower end and 1208.2 mm at the upper end of the elevation gradient on LM). On NT there was also a higher variation coefficient of MAT than on LM, which could have facilitated the detection of the PP response of trees in this mountain. In addition, trees sampled at NT were older, on average, than at LM ( Table 1 ), which could suggest that sensitivity of radial growth to climate fluctuation increases with age in P. hartwegii trees. Similar findings have been reported in other studies, particularly for high mountain species (Voelker, 2011; Dominik et al., 2019). Besides, the proportion of trees with significant RN was similar at both ends of the altitudinal range, despite the broader inter-annual climate fluctuation at high elevation. At the upper end of the gradient, there was a higher variation coefficient in MAT and TAP, than at the lower end. However, environmental conditions at the upper end of the gradient are far more restrictive for vegetative growth (sub-optimal temperature and lower depth and moisture retention capacity of soil), limiting the response of trees to fluctuations in temperature and precipitation (Holtmeier and Broll, 2005; Silva et al., 2016). Despite the results obtained, we suggest that for further research, the variation of conditions within a year or the global variation between years should be considered, to corroborate that phenotypic plasticity does not refer only to the response to average climatic conditions from one year to another.

The similar correlation matrix and separation of growth ring traits in three distinctive groups, based on the PP response to T_GS and Twinter, shows a strong relationship between the temperature variables at different periods of the year and their effect on the growth ring traits. In other words, the PP values for any growth ring trait in response to T_GS and with Twinter were strongly correlated. However, the separation of ring width traits from ring density traits, and these from maximum density indicates a different plastic response, globally, of these traits to the temperature variables. The positive loadings of PC2 to RD and EWD, and negative loadings to RW and LWW show contrasting RN in these traits to increasing temperature; a positive response for the first group and a negative one for the second. MAXD also showed, globally, a negative response to increasing temperature, more pronounced in the RN to Twinter. This contrasting response in RN for width vs. density seems to indicate a compensatory effect in resource allocation at the global (population) level in response to increasing temperature, reducing cell division and size, and increasing cell wall deposition (Rozenberg et al., 2020). A reduction in the rates of cell division and/or the diameter of the cells (less ring width) and an increase in the thickness of the cell wall (higher wood density) with increasing temperature, could allow an increase in the functionality and hydraulic security of the xylem in response to increased water stress (Rozenberg et al., 2020). Several authors have pointed out that during episodes of drought or increased temperature, evapotranspiration increases, and with it, the water demand of the trees, increasing the vulnerability of the xylem to hydraulic failures due to cavitation and embolism (Bréda et al., 2006; Meinzer et al., 2010; McDowell et al., 2011).

On the other hand, the correlation analysis also indicated a low and generally non-significant correlation between the PP values for width and density traits. Thus, despite the potential compensatory effect between growth ring traits in response to temperature described at the population level, at the tree level this compensatory effect between growth ring traits was not shown. This is corroborated by the broad within-population variation in PP values found in all growth ring traits in the P. hartwegii populations sampled. Thus, trees differed in their PP response; they seem to have a predominant response in one or few growth ring traits, and for those with a simultaneous significant RN for several growth ring traits, the PP response was relatively independent for width and density traits. Soto-Carrasco (2021), found similar results, arguing that trees exposed to similar variations in temperature can respond differently in the PP of their radial growth traits. Similarly, Pélabon et al. (2013), indicate that high intra-population variation is common and reflects variation in the specific physiological processes associated with the phenotypic traits involved. Escobar-Sandoval et al. (2021), mention that variation in the direction of the slope of the RN among trees within the same population is possible because different trees are positioned at different points along the response curve to the climate variable, i.e., they inherently have different RN. In addition, it has been reported that the RN of the majority of phenotypic traits analyzed in plants over wide temperature intervals are generally expected to be curvilinear (Arnold et al., 2019). As such, the detection of individual trees whose significant RN have signs opposite to the general trend could suggest differences among the trees in their position on the temperature curve. However, the results could also be due to the effect of other micro-environmental factors intrinsic to each population (e.g., micro-spatial variation in temperature or soil moisture), as has been described by other authors (Diamond and Kingsolver, 2012; Tammaru and Teder, 2012; Smoliński et al., 2020). Natural tree populations are characterized by high phenotypic variation in adaptive traits, due to their inherently high genetic variation and high pollen dispersal capacity, which enables them to cope with environmental changes. If phenotypic plasticity has any adaptive role, it would also be expected to be highly variable; thus, genetic variation in climatic sensitivity for radial tree growth is also a possible cause of the variation in PP within populations at both mountains (Rolland et al., 1999; Lebourgeois et al., 2010).

The general trend of higher frequency of negative RN for ring width traits and positive RN for ring density traits was most noted at the upper end of the elevation gradient on both mountains. This partially coincides with the results of Soto-Carrasco (2021), who mentions that on the mountain Pico de Orizaba there was a higher proportion of positive RN for ring density traits at higher elevation. However, in her study, the percentage of positive and negative RN for the ring width traits did not differ between elevations. We also found several differences in the average PP values between the ends of the elevation gradient sampled. The highest values of PP for both positive and negative RN were found at the lower elevation on both mountains, for most growth ring traits in response to both T_GS and Twinter. These results coincide with the report by Escobar-Sandoval et al. (2021), where the absolute value of the mean slope of positive and negative RN were generally higher at the lower end of the gradient of the three mountains analyzed in that study, which was attributed to the differences in the mean temperature between the two extremes of the elevation gradient. In our study, the average value of T_GS during the study period was 2.0°C higher at the lower than at the upper end of the elevation gradient on LM and 4.3°C higher on NT. The difference in Twinter was even more marked, with 5.3°C and 4.5°C higher temperatures at the lower than at the higher elevation on LM and NT, respectively. A change in the slope of the RN (i.e., a change in the phenotypic plasticity) associated with a change in the range of T_GS and Twinter values implies that the RN associated with temperature are not linear (Arnold et al., 2019). This is the case in many biological systems (Camarero and Gutiérrez, 2004; Körner, 2012). It is common, for example, that as temperature increases, the rate of change (slope of the regression line) of metabolic processes is modified (Arnold et al., 2019; Escobar-Sandoval et al., 2021). It is therefore possible that the higher phenotypic plasticity observed at the lower end of the gradient is due to the average position of the population with respect to a non-linear RN of these traits to T_GS and Twinter, as has been pointed out in other studies with different growth traits in several plant species (Valladares et al., 2006; Arnold et al., 2019).

The elevational differences in PP were generally more pronounced on LM than on NT, associated in most cases with the higher average PP value at the lower end of LM ( Figure 5 ). A steeper slope of the negative RN for ring width and the positive RN for ring density indicates that trees on LM are adjusting more rapidly their radial growth traits to the temperature increase. Gilmore (1968), mentions that the factors that favor a reduction in ring width generally increase their density, mainly in the second half of the radial growth period (Büntgen et al., 2010; Bouriaud et al., 2015). The differences between mountains, therefore, suggest different needs and/or responses of trees, as a function of the micro-environmental conditions on them. Even though temperature at the lower end was similar on both mountains, LM has lower annual precipitation, particularly during the summer and early fall ( Figure 1 ), so trees might be experiencing higher water stress. This differential response is also shown in the wide intra-population variation in the sign and magnitude of the RN slopes found in all growth ring traits, with individuals that differed from the general population trends, discussed before. Although it is difficult to identify the main factors involved in this PP variation, the importance of PP variation as a buffering mechanism of the impact of climate change at the population level is evident, as has been mentioned by Escobar-Sandoval et al. (2021).

Therefore, P. hartwegii trees seem to have the capacity to show plastic responses with different growth ring traits under climate change-driven environmental variation. Populations at the lower extreme of the elevation gradient showed higher proportion of trees with negative trends in their RN to temperature and higher absolute values of PP for growth ring traits, so they seem to be experiencing higher impacts from global warming and climate change. However, the wide intrapopulation variation detected in PP of these traits offers positive perspectives on the short-term response capacity of these populations, as PP may serve as a mechanism to buffer the impact of interannual fluctuation and the gradual increase in temperature in the medium term. However, it is fundamental to obtain additional information with respect to the physiological mechanisms associated with this variation in PP, its possible genetic origin (Zacharias et al., 2022) and the adaptive implications (DeWitt and Scheiner, 2004; Ghalambor et al., 2007; Chown et al., 2008).