Location of sampled populations

Natural populations of P. hartwegii were sampled in December 2018 at the limits of its altitudinal distribution, 3 500 m (lower) and 4 000 m (upper), in the Cofre de Perote, Pico de Orizaba, and Monte Tláloc mountains in the states of Veracruz, Puebla, and Mexico, respectively (Figure 1). Environmental conditions in the three mountains are similar regarding climate, soils, topography, and vegetation cover (Table 1). The climate of the study areas is type C(w) semi-cold humid with summer rainfall (^{Garcia, 2004}). Soils correspond to humic Andosol, umbric Andosol and Leptosol units (^{IUSS Working Group WRB, 2015}).

Pure natural stands were located on each mountain where mature and dominant trees of the species were selected. The specific elevation of each stand varied slightly, due to the natural spatial distribution of vegetation; for the purposes of the study, elevations are referred to as 'lower' and 'upper' (Table 1). Two sampling stands were selected at each site with a minimum spacing of 300 m in a straight line.

Figure 1 Geographic location of the three mountains of the Transmexican Neovolcanic System with natural populations of Pinus hartwegii sampled at the altitudinal limits (≈3 500-4 000 m) of their natural distribution (information taken and modified from ^{Gómez-Guerrero et al., 2013}). 

Sample collection and analysis

We selected 25 mature, dominant, healthy, adult trees at each stand, with a vigorous crown, straight stem and without defects or physical damage (^{Villanueva-Diaz et al., 2007}). From each tree a 5 mm diameter wood core was collected using an increment borer, at 1.30 m height, up to the center of the stem and in a perpendicular orientation to the soil slope to avoid compression wood. Wood cores were stored in open polycarbonate boxes to allow moisture loss, and dried at room temperature in the laboratory.

Samples were sawn about 2 mm thick, and resin was extracted with a pentane bath (C₅H₁₂) for 48 h. The samples were dried again and exposed to X-rays with the procedure described by ^{Mothe, Duchanois, Zannier, and Leban (1998}). X-ray images (negatives) were scanned at 4 000 dpi to acquire microdensity profiles using WINDENDRO® (Windendro 2008e Regent instruments Canada Inc.).

Microdensity profiles were visually verified, compared with the sawn samples, and dated to identify the year of formation of each growth ring. According to the procedures described by ^{Rozenberg et al. (2012}), for each ring, the following variables were obtained: total ring width (TRW) (mm), of earlywood (EWW), and latewood (LWW); and total ring density (TRD), of earlywood (EWD), latewood (LWD), minimum (MID), and maximum (MAD) density (g∙cm^-3). Five profiles were excluded from the analysis because the ring boundaries were undefined due to an oblique grain angle, which slightly reduced the sample size at some sites (Table 1).

Analysis of densitometric data

Densitometric data were fitted with the residual method to remove the effect of cambial age, using independent models for each location (^{Rozenberg et al., 2020}). Rings close to the pith were excluded in the analysis to avoid variation due to juvenile effect (^{Cerano-Paredes et al., 2016}). Since the last ring (2018) had not been fully formed, it was also not included in the analysis. Based on these criteria and availability of instrumental climate data, it was determined that the most appropriate period to distinguish the effects of climatic factors was from 1960 to 2017.

With the average values of the growth ring components per tree, a two-way ANOVA was carried out with the MIXED procedure of SAS version 9.4 (SAS Institute, 2013), to estimate average values per location and evaluate the effect of elevation, mountain, and their interaction. Subsequently, the Pearson correlation coefficient was estimated separately for each location with the purpose of analyzing the interrelationships between annual radial growth and the components of the growth ring. The homogeneity of the structure of correlations between elevations and mountains was compared using the Chi-square test (^{Rathgeber, 2017}).

Growth ring traits

The average value of growth ring components corresponding to the period 1960-2017 was different between mountains and altitudinal limits (P ≤ 0.05), except for EWD (Table 2). Table 3 shows that narrower growth rings (TRW and EWW), but with higher LWW and density (TRD, EWD, LWD, MID and MAD) were found on Pico de Orizaba than on the other mountains. The average value of ring traits in Cofre de Perote and Monte Tláloc were similar, except EWD and MID which are 6.7 % and 4.1 % lower in Monte Tláloc.

Regarding elevation, Table 4 shows higher average radial growth with wider growth rings at the lower altitudinal limit, due to greater earlywood and latewood formation than at the upper limit. Average ring density was also higher at the lower limit, mainly due to higher LWD and MAD, as EWD was unchanged and MID was slightly higher at the upper limit. The effect of elevation was similar in the three mountains, except for TRW, EWW, LWW and MID in Pico de Orizaba where elevation did not affect them (Figure 2).

Figure 2 Average value of width (a) and density (b) components of growth rings from 1960-2017 at the lower (≈3 500 m) and upper (≈4 000 m) limits of the altitudinal distribution of Pinus hartwegii in the Cofre de Perote (P), Pico de Orizaba (O) and Monte Tláloc (T) mountains. TRW = total ring width, EWW = earlywood width, LWW = latewood width, TRD = ring density, EWD = earlywood density, LWD = latewood density, MID = minimum ring density and MAD = maximum ring density 

Correlation structure between growth ring traits

Table 5 shows the correlation matrix between ring width and ring density traits; the Chi-square test indicated that the correlation structure was not homogeneous at the altitudinal distribution limits of P. hartwegii. The strong positive correlation of TRW with EWW (r ≥ 0.95) was similar at both altitudinal limits, but the correlation of these with LWW was higher at the lower limit.

Correlations of TRD with EWD (r ≥ 0.82), EWD with MID (r ≥ 0.92) and LWD with MAD (r ≥ 0.91) were also homogeneous at the two altitudinal limits. At the upper limit, correlation of TRD with MID (r ≥ 0.72) was higher and with LWD (r ≥ 0.20) and MAD (r ≥ 0.07) was heterogeneous (Table 5). Correlations of width variables (TRW, EWW and LWW) with density variables (TRD, EWD, LWD, MID and MAD) showed a different structure at the two limits. At the lower limit, TRW, EWW and LWW had a weak, positive or negative correlation (-0.35 ≤ r ≤ 0.35) with TRD, EWD and MID, while, at the upper end, correlations with TRD and EWD were negative (-0.82 ≤ r ≤ -0.17); furthermore, correlations of LWD and MAD with TRW, EWW and LWW at the lower limit of Pico de Orizaba and with LWW at the upper limit of Monte Tláloc were not significant (P > 0.05). At the upper limit, correlations between growth and density were negative for earlywood components (EWD and MID), and positive for latewood components (LWD and MAD).

The structure of correlations was also heterogeneous between mountains; Cofre de Perote and Monte Tlaloc contributed more to this heterogeneity at the lower limit, and Pico de Orizaba and Monte Tlaloc at the upper limit. A more homogeneous structure of correlations between altitudinal limits was found in Pico de Orizaba (Table 5).

Differences in growth ring traits

The effect of elevation on radial growth and growth ring components of P. hartwegii was significant, with rings 50 % wider and 5 % denser at the lower limit (≈3 500 m) of the altitudinal gradient from 1960-2017. ^{Morgado-González et al. (2019}) describe similar results for this species, with smaller average ring width and density with increasing site elevation. For Larix decidua Mill., ^{Rozenberg et al. (2020}) also reported higher average ring width and density at the lower limit of the altitudinal gradient. However, in Fagus silvatica Lipsky a compensatory effect of elevation was found, reflecting higher average growth and lower wood density at the lower limit of the gradient (^{Topaloğlu, Ay, Altun, & Serdar, 2016}).

Dissimilarities in radial growth between altitudinal limits may be due to differences in temperature and growing season. Data derived from the Climate NA model (^{Wang et al., 2016}) show that, in the study mountains, at the lower altitudinal limit of P. hartwegii forests, mean temperature is 3 °C higher and the frost-free period in the year is 64 days longer than at the upper limit. Although high mountain species are adapted to low temperatures, mean annual temperature (6 °C) at the upper limit may not be enough to reach the optimal level of photosynthesis, which only occurs at times of maximum irradiance (^{Körner, 2012}; ^{Sanfuentes, Sierra-Almeida, & Cavieres, 2012}). A suboptimal temperature, together with a shorter favorable growth period, is reflected in narrower growth rings and lower average density, as pointed out by ^{Cerrato et al. (2019}) in Pinus cembra L. The results show that there are better growth conditions at the lower limit of the gradient, which allow trees to modulate the xylem function more effectively, favoring hydraulic conductivity, growth, structural support, and carbon accumulation (Cerrato et al., 2019). In contrast, at the upper limit, the response emphasizes the hydraulic function, similar to that reported by ^{Camarero, Rozas, and Olano (2014}) in Juniperus thurifera L.

Studies in other conifers show that reduced soil moisture (due to lower precipitation, water freezing, or lower retention capacity) induces latewood formation, with smaller diameter tracheids and greater cell wall thickness, to maintain the xylem hydraulic function (^{Camarero et al., 2014}; ^{Moyes, Castanha, Germino, & Kueppers, 2013}).

Differences in radial growth and density among mountains are also related to differences under environmental conditions. Narrower rings were found at Pico de Orizaba, but with a higher proportion of latewood and higher average density and component densities, including minimum and maximum ring density (Table 3), which could be related to an earlier transition to latewood formation and a shorter growth period, as suggested by the estimated climatic data. This mountain had a mean annual temperature (and during the frost-free period) 1 to 2 °C lower than the others; in addition, the frost-free period was 56 days shorter, starting 44 days later and ending 12 days earlier than in the other mountains, which may explain the higher relative proportion of latewood and higher ring density. The lower temperature at the top of Pico de Orizaba could explain the higher MAD (at the end of the growing season), to reduce cavitation risks at that time of the year, when the greatest temperature drop occurs (^{Arias, Bucci, Scholz, & Goldstein, 2015}; ^{Scholz, Bucci, Arias, Meinzer, & Goldstein, 2012}).

Interrelationships between radial growth and ring components

In most cool-temperate conifers it is common to find a strong positive relationship between ring width and earlywood and latewood components, as reported by ^{Rozenberg et al. (2020}) for L. decidua (r ≥ 0.7). ^{Bjorklund et al. (2017}) also indicate positive relationship between average ring density and that of its components for 27 conifer species (r ≥ 0.5), in addition to a weak negative correlation between ring width and ring density and its components. Negative correlations between wood density and ring width are mentioned in P. hartwegii by ^{Morgado-González et al. (2019}); however, ^{Correa-Díaz et al. (2020}) found that maximum ring density had a positive, although weak (r = 0.3), correlation with ring width at Monte Tláloc.

In general, the pattern of correlations between growth ring traits in P. hartwegii is similar to that described for other conifers (^{Björklund et al., 2017}; ^{Dalla-Salda, Martinez-Meier, Cochard, & Rozenberg, 2011}); however, in this species, the pattern is modified with elevation. At the higher limit, the correlation of TRW with LWW was weaker and the negative correlations of TRW, EWW and LWW with TRD, EWD and MID were stronger (Table 5). This effect is consistent with differences under climate conditions at the gradient limits; low temperatures and short length of the favorable growing season at the upper limit cause latewood to have less relative importance in growth ring width, as has been shown in other species (^{Correa-Díaz et al., 2020}; ^{Düthorn et al., 2016}).

On the other hand, the lower photosynthetic capacity at the upper limit of the gradient, associated with lower average temperature, generates greater competition for resources, especially during earlywood formation when shoot growth also occurs (^{Antonucci et al., 2017}; ^{Li et al., 2017}). As environmental conditions become harsher, processes associated with energy capture, such as photosynthesis, are reduced to a greater degree (^{Pompa-García & Venegas-González, 2016}). In the upper part, lower temperature and lower soil water availability limit photosynthetic activity as noted by ^{Pompa-García, Camarero, Colangelo, and Gallardo-Salazar (2021)} in Pinus leiophylla Schl. et Cham. This could explain the stronger negative correlation of TRW and EWW with TRD and EWD at the upper limit of the gradient. This could also explain the difference in magnitude of the correlations of TRW and EWW with LWD and MAD at the upper limit of Pico de Orizaba compared to Cofre de Perote and Monte Tláloc. These correlations were stronger on Pico de Orizaba and Monte Tláloc, where the upper part is cooler. When conditions improve in this more unfavorable zone, earlywood density decreases and latewood density increases; that is, the tree can distribute its resources more efficiently between earlywood and latewood of the ring, to favor conduction functions with larger cells in earlywood, and xylem support functions, with thicker-walled cells in latewood. Under more limiting conditions of Pico de Orizaba, there are fewer possibilities than in Cofre de Perote and Monte Tláloc to distribute resources efficiently.

The lower radial growth of P. hartwegii at the upper limit, linked to lower wood density, reflects the effect of limiting factors and response mechanisms of the species to adjust to environmental conditions during wood formation, maintaining xylem functionality (^{Park & Spiecker, 2005}; ^{Rozenberg et al., 2020}). These differences at the altitudinal gradient limits of the species help to understand the differential response of its populations to increased temperature and environmental fluctuation associated with climate change. However, more detailed studies regarding phenology of cambial activity are required to understand the process of xylogenesis and to identify possible compensatory mechanisms among the components of the growth ring. These mechanisms would be required to maintain the functions of mechanical support and xylem hydraulic conductivity (^{Gindl, Grabner, & Wimmer, 2000}) in a scenario of greater climate fluctuation occurring at the altitudinal gradient limits of natural distribution of the species.