Long-Term Wood Micro-Density Variation in Alpine Forests at Central M é xico and Their Spatial Links with Remotely Sensed Information

: Ongoing climate variability strongly a ﬀ ects high-elevation forests, inﬂuencing the wood formation process (e.g., xylogenesis). Furthermore, spatio-temporal studies to establish links of wood properties and tree performance are needed. Using linear mixed-e ﬀ ects models, empirical cumulative distribution functions, and spatial analysis, we explore time trends and space connections of wood density of Pinus hartwegii Lindl. to remotely sensed variables (Moderate Resolution Imaging Spectro-radiometer MODIS-derived) in two high-elevation forests in M é xico, Tl á loc (TLA) and Jocotitl á n (JOC) Mountains. Results indicated that elevation and cambial age e ﬀ ects are important factors explaining wood density variation. Minimum earlywood—MID, average—AVE, and maximum latewood density—MXD were statistically similar between mountains ( p > 0.05), but TLA showed a signiﬁcant increase in MID over time with higher values after 1950. Wood density values and spatial correlations were site-dependent with TLA exhibiting the highest correlations between MXD and the Normalized Di ﬀ erence Vegetation Index (NDVI) of the spring season (r = 0.59, p < 0.05). Overall, correlations to remotely sensed information were positive with MXD, negative for MID and divergent for AVE. Historical temperature deﬁnes MID along the elevation gradient, while MXD was related to soil moisture only at low-elevation sites where soils are deeper. We found that two high-elevation forests, 115 km away from each other, with similar climate, soil, and vegetation, behaved di ﬀ erently regarding their xylogenesis, indicating the potential of using the link between wood micro-density and remotely sensed information to understand forest response to climate change e ﬀ ects. J.J.V.-H.; Resources, and W.R.H.; P.R.; Visualization, W.R.H.; Writing—original draft, A.C.-D. and A.G.-G.; Writing—review and editing, A.G.-G., W.R.H. J.J.V.-H.


Introduction
Wood micro-density (hereinafter named wood density) is linked to the water and nutrient transport processes, plant mechanical support and in the long-term is related to the carbon-storing capacity of trees [1,2]. Furthermore, the retrospective analysis of tree-ring density is useful to capture the historical plant responses to climate variation [3] since the memory of carbon gain is imprinted in the tree-ring

Study Sites
Two high-elevation mountain areas across the Trans-Mexican volcanic belt were studied, namely Tláloc Mountain (TLA) (19.39 (Figure 1a). The dominant tree species is Pinus hartwegii Lindl, forming pure, uneven-aged open stands [43]. Given that the summit of each mountain is different (Table 1), we focused on selecting sites formed from pure stands of P. hartwegii concentrating on those with the most contrasting elevations. We sampled P. hartwegii trees at two elevation positions (3500-3900 m asl for TLA and 3700-3800 m asl for JOC) with two contrasting aspects (Northwest and Southwest). Mean annual temperature is similar at both sites (14.8 ± 0.6 and 14.2 ± 0.4 for TLA and JOC, respectively). Precipitation is slightly higher at TLA, mainly in the period July-August (Figure 1b,c) [44]. The climate at both sites is classified as temperate subhumid (C(w 2 )) according to Köppen climate classification as modified by García [45] for México conditions. Soils at both sites are Andisols, developed from volcanic material resulting in high organic matter content epipedons with a high water-holding capacity [46].  (Figure 1a). The dominant tree species is Pinus hartwegii Lindl, forming pure, uneven-aged open stands [43]. Given that the summit of each mountain is different (Table 1), we focused on selecting sites formed from pure stands of P. hartwegii concentrating on those with the most contrasting elevations. We sampled P. hartwegii trees at two elevation positions (3500-3900 m asl for TLA and 3700-3800 m asl for JOC) with two contrasting aspects (Northwest and Southwest). Mean annual temperature is similar at both sites (14.8 ± 0.6 and 14.2 ± 0.4 for TLA and JOC, respectively). Precipitation is slightly higher at TLA, mainly in the period July-August (Figure 1b,c) [44]. The climate at both sites is classified as temperate subhumid (C(w 2 )) according to Köppen climate classification as modified by García [45] for México conditions. Soils at both sites are Andisols, developed from volcanic material resulting in high organic matter content epipedons with a high water-holding capacity [46].  [44]; and (c) General view of Mexican high-elevation forest at TLA.

Sample Collection and Density Measurements
At both mountain sites (TLA and JOC), we collected 80 increment core samples at breast-height, using a 5-mm diameter Pressler borer. The average number of trees sampled per site condition (a combination of elevation × aspect) was 10 trees, sampling only dominant and healthy mature trees. Once the increment cores were oven-dried, they were sawn using a twin-blade (1.5 mm thick cross-sections) and then, resin-extracted with in a pentane bath. To obtain X-ray data, each core was exposed about 25 min to X-ray densitometry [47,48]. Then, the X-ray microfilms were scanned (4000 dpi) and analyzed using WinDendro (Regent Instruments Inc., Quebec, Canada). WinDendro software converts the gray levels of the X-ray microfilms to density values by comparing them to a standard of known physical and optical density through a calibration curve [49]. We discarded some trees in the process because of blurred information in some micro rings (≈20% of total) but not less than six trees represented a site and condition, thus, a total of 63 trees were used for the analysis. Finally, we calculated MID, AVE, and MXD annually for each tree. All samples were prepared and processed at INRA Val de Loire Centre, Orléans France.

Mountain and Topographic Gradient Effect on Density Profiles
The topographic effect in tree-ring density variables (MID, AVE, and MXD) was tested by linear mixed-effects models for repeated measured data using R ('nlme' package) [23,50]. An autoregressive process of order 1 was included for the correlational structure of repeated measurements (AR1). First, we tested significant differences for the mean wood density profiles between mountains, using trees and cambial age as random factors. Since we found no significant differences between mountain profiles, we decided to combine them into a single database to test the elevation, aspect, and cambial age effects at a regional level. Here, we used a nested effect of trees within mountains as random factors, blocking any effect by the mountain source. Random intercepts and slopes were set to measure the changes of every tree. Finally, statistical differences among topographic conditions were tested by pairwise comparison using the post hoc Bonferroni's adjustment using emmeans library in R [51].

Unbiased Trends of Density Profiles and Temporal Stability
To obtain unbiased MID, AVE, and MXD profiles by the cambial age effect, we adjusted the density measurements by site (elevation × aspect) with a linearly non-parametric method similar to the residuals RCS method [52]. The use of the subtraction operator instead of the ratio allows producing time series in the original density scale (g cm −3 ) and a straightforward interpretation. Afterward, temporal trends were tested by linear mixed-effects models as described above, using the slope and significance as indicators of change over time.
Empirical cumulative distribution functions (CDF) were calculated for adjusted MID, AVE, and MXD time series to probe the stability over time using a breaking point in 1950. We chose this year since it matched with a notable forest growth reduction (basal area increment) observable at TLA mountain [34]. Then, we used the nonparametric Kolmogorov-Smirnov (KS) test to prove that the distribution function of the period <1950 was equal to the more recent period (>1950); this method is particularly useful when data does not come from a normal distribution as in this case [21].

Climatic Influence on Tree-Ring Wood Density Variables
We correlated monthly temperature (maximum and minimum), precipitation and soil moisture with the adjusted density profiles by elevation level (e.g., 3500, 3700, 3800, 3900 m asl) using the bootstrapped Pearson's correlation from "treeclim" in R [53]. All climatic variables, excluding soil moisture, were obtained using a downscaled time series (according to 250-m digital elevation model) generated by software package ClimateNA v550 from 1902 to 2016 [44,54], based on the Climate Research Unit database (CRU ts 4.01) ( Table 1). Soil moisture was retrieved using the data from the project Climate Change Initiative (CCI) from the European Space Agency (ESA). We used the combined

Linking Remotely Sensed Variables from MODIS to Wood Density
We compared the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) 16-day composite time series (250-m spatial resolution) [55], and Leaf Area Index (LAI) 8-day composite time series (500-m spatial resolution) [56,57] with the adjusted tree-ring density values from 2000 to 2016. Preprocessing analysis included pixel quality analysis (avoiding snow-and cloud-contaminated pixels), harmonic interpolation (when pixel gaps were presented) and maximum value composite for LAI data (16-day composite) [58,59].
A Pearson product-moment correlation coefficient was used to explore the link between temporal and spatial variables as described by Correa-Díaz, et al. [39] where the yearly tree-ring derived measurements (e.g., MXD) are spatially correlated to remotely sensed maps from November of the previous year until September of the current year. Once the correlation maps were calculated, not statistically significant pixels were masked out according to their p-values (p < 0.05). This procedure allows for detecting the highest spatio-temporal association between yearly and bi-weekly data.

Tree-Ring Density Statistics and the Cambial Age Effect
Minimum earlywood density (MID) ranged from 0.16 to 0.72 g cm −3 , with a mean of 0.

Linking Remotely Sensed Variables from MODIS to Wood Density
We compared the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) 16-day composite time series (250-m spatial resolution) [55], and Leaf Area Index (LAI) 8-day composite time series (500-m spatial resolution) [56,57] with the adjusted tree-ring density values from 2000 to 2016. Preprocessing analysis included pixel quality analysis (avoiding snow-and cloud-contaminated pixels), harmonic interpolation (when pixel gaps were presented) and maximum value composite for LAI data (16-day composite) [58,59].
A Pearson product-moment correlation coefficient was used to explore the link between temporal and spatial variables as described by Correa-Díaz, et al. [39] where the yearly tree-ring derived measurements (e.g., MXD) are spatially correlated to remotely sensed maps from November of the previous year until September of the current year. Once the correlation maps were calculated, not statistically significant pixels were masked out according to their p-values (p < 0.05). This procedure allows for detecting the highest spatio-temporal association between yearly and bi-weekly data.

Tree-Ring Density Statistics and the Cambial Age Effect
Minimum earlywood density (MID) ranged from 0.16 to 0.72 g cm −3 , with a mean of 0.  Raw tree-ring widths were negatively associated to MID and AVE at TLA (r = −0.50 and r = −0.33, p < 0.05) and JOC (r = −0.45 and r = −0.39, p < 0.05), and positively related to MXD only at TLA (r = 0.29, p < 0.05) ( Figure S1). Contrary to the expected plateau in wood density for mature trees, we found a constantly increasing trend in MID and AVE at TLA and an evident reduction in MXD at JOC (Figure 3a-c). According to the cambial age, at both mountains, a peak in MXD occurs about the (r = 0.29, p < 0.05) ( Figure S1). Contrary to the expected plateau in wood density for mature trees, we found a constantly increasing trend in MID and AVE at TLA and an evident reduction in MXD at JOC (Figure 3a-c). According to the cambial age, at both mountains, a peak in MXD occurs about the age of 40 y (in the juvenile stage); however, TLA showed a second peak of MXD around to the age of 150 y ( Figure 3c).
Forests 2020, 11, x FOR PEER REVIEW 6 of 18 age of 40 y (in the juvenile stage); however, TLA showed a second peak of MXD around to the age of 150 y (Figure 3c).

Elevation and Aspect Variations for Tree-Ring Wood Density
We found that all tree-ring wood density variables showed an elevation (p < 0.05) and cambial age effect (p < 0.01) with no differences between mountains (p = 0.81). This means that these factors affected wood density in the same way regardless of geographical location. The aspect was not statistically significant for AVE (p = 0.07), and MXD (p = 0.24), however, it was relevant for MID (p = 0.029). The interaction elevation × aspect was marginally significant for MID and MXD (p < 0.08) ( Table 2). Table 2. Results of the analysis of variance on linear mixed-effects models fitted to the tree-ring wood variables of Pinus hartwegii. In the analysis, trees and mountain sites were tested as random effects. DF, degrees of freedom. Significance levels † p < 0.1, * p < 0.05, and ** p < 0.001.

Variable
Fixed Significantly higher MID and AVE were seen at 3700 m asl, regardless of aspect. Thus, statistical differences (using Bonferroni's adjustment for multiple comparisons) for MID were calculated between the 3700-NW site (elevation-aspect) to 3500-NW, 3800-SW, and 3900-NW sites (p < 0.05) (Figure 4a). For AVE, differences were between the 3700-SW site with 3800-SW, and 3900-NW (p <

Elevation and Aspect Variations for Tree-Ring Wood Density
We found that all tree-ring wood density variables showed an elevation (p < 0.05) and cambial age effect (p < 0.01) with no differences between mountains (p = 0.81). This means that these factors affected wood density in the same way regardless of geographical location. The aspect was not statistically significant for AVE (p = 0.07), and MXD (p = 0.24), however, it was relevant for MID (p = 0.029). The interaction elevation × aspect was marginally significant for MID and MXD (p < 0.08) ( Table 2). Table 2. Results of the analysis of variance on linear mixed-effects models fitted to the tree-ring wood variables of Pinus hartwegii. In the analysis, trees and mountain sites were tested as random effects. DF, degrees of freedom. Significance levels † p < 0.1, * p < 0.05, and ** p < 0.001. Significantly higher MID and AVE were seen at 3700 m asl, regardless of aspect. Thus, statistical differences (using Bonferroni's adjustment for multiple comparisons) for MID were calculated between the 3700-NW site (elevation-aspect) to 3500-NW, 3800-SW, and 3900-NW sites (p < 0.05) (Figure 4a). For AVE, differences were between the 3700-SW site with 3800-SW, and 3900-NW (p < 0.05), and marginally with 3500-NW and 3900-SW sites (p < 0.10) (Figure 4b). Finally, differences for MXD were found between the 3700-SW site with 3800-NW, 3900-NW, 3800-SW and 3900-SW sites (p < 0.05) (Figure 4c).

Trends and Stability in Adjusted Tree-Ring Wood Density Profiles
MID statistically increased over time (p < 0.001) (calendar year effect) at TLA mountain ( Figure  5a). The average increment in MID for TLA is equivalent to 33 kg m −3 decade −1 , this is approximately an increase of 50 kg m −3 throughout the life of the trees. However, at the tree level, we found increases up to 100 kg m −3 . The rate of increase in MID was significantly higher (p = 0.001) at high elevation sites (3900-NW and 3900-SW). In contrast, MID, AVE, and MXD did not change over time at any mountain (p > 0.10) (Figure 5b,c). The difference between MXD and MID was similar at both mountains, but this amplitude increased over time at TLA (p = 0.03), possibly caused by a rising MID, denoting higher reaction response at this site.

Trends and Stability in Adjusted Tree-Ring Wood Density Profiles
MID statistically increased over time (p < 0.001) (calendar year effect) at TLA mountain (Figure 5a). The average increment in MID for TLA is equivalent to 33 kg m −3 decade −1 , this is approximately an increase of 50 kg m −3 throughout the life of the trees. However, at the tree level, we found increases up to 100 kg m −3 . The rate of increase in MID was significantly higher (p = 0.001) at high elevation sites (3900-NW and 3900-SW). In contrast, MID, AVE, and MXD did not change over time at any mountain (p > 0.10) (Figure 5b,c). The difference between MXD and MID was similar at both mountains, but this amplitude increased over time at TLA (p = 0.03), possibly caused by a rising MID, denoting higher reaction response at this site.

Trends and Stability in Adjusted Tree-Ring Wood Density Profiles
MID statistically increased over time (p < 0.001) (calendar year effect) at TLA mountain ( Figure  5a). The average increment in MID for TLA is equivalent to 33 kg m −3 decade −1 , this is approximately an increase of 50 kg m −3 throughout the life of the trees. However, at the tree level, we found increases up to 100 kg m −3 . The rate of increase in MID was significantly higher (p = 0.001) at high elevation sites (3900-NW and 3900-SW). In contrast, MID, AVE, and MXD did not change over time at any mountain (p > 0.10) (Figure 5b,c). The difference between MXD and MID was similar at both mountains, but this amplitude increased over time at TLA (p = 0.03), possibly caused by a rising MID, denoting higher reaction response at this site.  The empirical cumulative distribution functions (CDF) and the Kolmogorov-Smirnov test confirmed that all tree-ring wood density distributions were significantly different between the two periods analyzed, <1950 and >1950 (Table 3). Remarkably, we found the same pattern between mountains for MID and AVE, showing a shift to denser rings after the second half of the 20th century. However, for MXD divergent conclusions arose, whereas a larger CDF after 1950 (p < 0.001) was found in TLA (less dense rings), and a smaller CDF was found for JOC (denser rings > 1950) (p < 0.001) ( Figure S2).

Responses of Wood Density to Climate
MID responded positively to the temperature (maximum and minimum), mainly at low-elevation sites. Remarkably both mountains followed the same pattern, with the maximum correlations in April for minimum temperature (r = 0.55, p < 0.01) and during the summer season for maximum temperature (r = 0.54, p < 0.01) (Figure 6a). In contrast, precipitation and soil moisture were not significantly related to MID except for few months at 3500 and 3700 m asl (Figure 6a). A divergent response was seen for maximum temperature correlated to MXD since the lowest elevation sites (i.e., 3500 m asl) showed negative correlations in April (r = −0.40, p < 0.01) and positive correlation for the sites with the highest elevation in JOC (r = 0.22, p < 0.01). In the same way, we found a negative (3500 m asl), and positive (3800 m asl) correlations with minimum temperature. Finally, MXD showed a negative relationship with soil moisture during the summer at the lowest elevation (3500 m asl) (Figure 6b).

The Connection between Remotely Sensed Information and Wood Density Traits
Significant correlations were detected between remotely sensed information and tree-ring wood density variables; however, the best spatial associations (referred to as a larger surface with high correlations) were at TLA for MXD ( Figure S3). Figure 7 summarizes the Pearson correlation with the total area for both mountains considering an α = 0.05. Overall, positive correlations were found for MXD, negative for MID and an opposite pattern for AVE.

The Connection between Remotely Sensed Information and Wood Density Traits
Significant correlations were detected between remotely sensed information and tree-ring wood density variables; however, the best spatial associations (referred to as a larger surface with high correlations) were at TLA for MXD ( Figure S3). Figure 7 summarizes the Pearson correlation with the total area for both mountains considering an α = 0.05. Overall, positive correlations were found for MXD, negative for MID and an opposite pattern for AVE. maximum association in early March (day of the year 065, r = 0.59, and 56% total area), and with LAI in late August (doy 241, r = 0.58, 37% total area (Figure 7). MID was negatively correlated to NDVI from previous winter to spring with a maximum association in middle-January at TLA (doy 017, r = −0.60, 44% total area), and with EVI from previous-December (previous doy 337, r = −0.54, 35% total area). At JOC, AVE was negatively correlated to EVI in early-April at JOC (doy 097, r = −0.55, 55% total area). Spatially, the northernmost part of TLA Mountain was more suitable to establish significant links between tree-ring wood density variables and remotely sensed data (Figure 8a). Therefore, MXD seems to be more reliable to explain spring photosynthetic activity with drops in 2003, 2011 and 2013.  In TLA, MXD significantly correlated to the NDVI from previous winter to spring with a maximum association in early March (day of the year 065, r = 0.59, and 56% total area), and with LAI in late August (doy 241, r = 0.58, 37% total area (Figure 7). MID was negatively correlated to NDVI from previous winter to spring with a maximum association in middle-January at TLA (doy 017, r = −0.60, 44% total area), and with EVI from previous-December (previous doy 337, r = −0.54, 35% total area). At JOC, AVE was negatively correlated to EVI in early-April at JOC (doy 097, r = −0.55, 55% total area). Spatially, the northernmost part of TLA Mountain was more suitable to establish significant links between tree-ring wood density variables and remotely sensed data (Figure 8a). Therefore, MXD seems to be more reliable to explain spring photosynthetic activity with drops in 2003, 2011 and 2013.

Hypothesis Testing Summary
The temporal wood density analysis provides new insights into how high-elevation forests modulate wood formation processes (i.e., xylogenesis) to cope with climate variation and assure survival. We showed that both mountains (TLA and JOC) had statistically the same mean wood density, but the elevation and cambial age are important factors. The notable increase in MID over time only at TLA rejects our hypothesis of the same response at different locations. Stronger evidence on the adjustment of wood density in TLA at high elevations suggests a higher response in this location. Probably the fact that JOC is an isolated volcano has also isolated this tree community and resulting in it being a less responsive forest. While the cumulative distribution functions (CDF) detected significantly higher wood density after 1950 (MID and AVE), progressive changes over time evaluated with linear mixed models were partially demonstrated. Unexpectedly, considering the well-known concept that MXD is highly correlated to summer temperatures, we found a poor association, and it was MID that demonstrated a higher correlation. MID was relevant because lower wood density relates to wider cells which occurs in the growing season. Finally, although we compared three different spatial variables that were highly correlated among them (NDVI, LAI, and EVI), the NDVI was the most important related to MXD during the spring season at TLA but no concurrent at JOC.

Hypothesis Testing Summary
The temporal wood density analysis provides new insights into how high-elevation forests modulate wood formation processes (i.e., xylogenesis) to cope with climate variation and assure survival. We showed that both mountains (TLA and JOC) had statistically the same mean wood density, but the elevation and cambial age are important factors. The notable increase in MID over time only at TLA rejects our hypothesis of the same response at different locations. Stronger evidence on the adjustment of wood density in TLA at high elevations suggests a higher response in this location. Probably the fact that JOC is an isolated volcano has also isolated this tree community and resulting in it being a less responsive forest. While the cumulative distribution functions (CDF) detected significantly higher wood density after 1950 (MID and AVE), progressive changes over time evaluated with linear mixed models were partially demonstrated. Unexpectedly, considering the well-known concept that MXD is highly correlated to summer temperatures, we found a poor association, and it was MID that demonstrated a higher correlation. MID was relevant because lower wood density relates to wider cells which occurs in the growing season. Finally, although we compared three different spatial variables that were highly correlated among them (NDVI, LAI, and EVI), the NDVI was the most important related to MXD during the spring season at TLA but no concurrent at JOC.

Tree-Ring Wood Density Patterns across Mountains, Elevations, and Aspects
Similar to other studies [63,64], we found that temporal variations in wood density are strongly influenced by cambial age and once the juvenile stage is finished there can be a declining or irregular trend (Figure 3). In fact, for TLA a clear increasing trend was observed for MID and AVE even with ages older than 100 years, which could suggest that trees growing at TLA have been adapting their wood morphological outcomes to climate variability more strongly than at JOC. The negative correlation found between raw tree-ring width and density (MID and AVE) ( Figure S1) suggest that any increase in radial growth results in lower wood density until it reaches a plateau. Lower density and wider cell lumen improve hydraulic conductivity but with a higher risk for cavitation [11]. On the other hand, the positive relationship with MXD for TLA could suggest the densification of woody tissues (higher carbon investment) within wider rings but only until a threshold is reached ( Figure S1) [3].
Elevation was an important factor explaining P. hartwegii wood density variation whereas aspect was a poor predictor except for MID (Table 2). Morgado-González, et al. [25] found similar results at TLA; however, they found that elevation × aspect interaction was statistically significant for AVE, while in this study only MID and MX were marginally significant (p < 0.10). Possible discrepancies could be attributed to the acquisition method since Morgado-González, et al. [25] tested X-ray computed tomography rather than X-ray densitometer that is demonstrated to be more accurate [65]. Many studies that focused on understanding wood density response with elevation have not found a clear pattern across regions. For example, an increase [66], a decrease [67] or even no change with elevation has been reported [68]. Thus, the increase in wood density at intermediate elevations (e.g., 3700 m asl) with a decrease at extreme elevations (e.g., 3500 and 3900 m asl) ( Figure 4) is a pattern previously reported by Topaloglu, et al. [69] for Fagus orientalis Lipsky.

The Temporal Stability of Tree-Ring Wood Density at Two High-Elevation Forests
Conversely to other studies showing a common decrease in average or latewood density during the second half of the twentieth century in Northern ecosystems [22,23,70], here we found a rising MID at TLA (using linear mixed models). Nevertheless, the temporal trends were dependent on the statistical method used, thus, the CDF analysis showed a shift to denser rings after 1950 (Table 3) for MID and AVE regardless of Mountain site. Although not all long-term trends in Figure 5 are statistically significant with linear mixed models (p < 0.05); they matched with slope direction in Kolmogorov-Smirnov test results. For example, less dense MXD values in TLA after 1950 agree with a negative slope observed in Figure 5c and denser MXD values in JOC agreed with a positive slope (Figure 5c). Andreu-Hayles, et al. [21] reported the same wood density adjustment with a shift to denser rings after 1950 for Picea glauca in Alaska.
Despite the statistical method used, the results suggest that, most notably at TLA, P. hartwegii modulate wood formation process according to climate variation (Figure 5a) as influenced by topographic position. We argue that climate sensitive tree-ring series are also indicative of the regulation of wood formation processes. In this regard, TLA exhibited a higher series intercorrelation and mean sensitivity than JOC (0.486 vs. 0.338 and 0.25 vs. 0.18, respectively), however, further information is needed to support this affirmation. Probably the fact that JOC is an isolated volcano expose this mountain to abrupt daily changes in temperatures, influencing cambial activity and leading lower correlation between wood density and climate and remote sensed variables.
The MID increase and the subsequent reduction of wood density range (e.g., the difference between MXD and MID) in TLA, could mean that trees are reducing water transport, despite high-elevation sites having better water status conditions than low-elevation sites [71]. Nevertheless, it is interesting to point out that not all trees followed the same trend which could denote either high phenotypic plasticity or a high genetic variation of P. hartwegii [72] (Figure 5).

Climate-Wood Density Responses
The maximum latewood density (MXD) is a well-known proxy for summer temperatures in northern forests [7,9,10,73,74]. Overall, MXD increases due to higher temperatures allow for production and deposition of more cell wall components due to longer growing seasons [7]. However, in Mexican temperate forests, an inverse relationship has been reported, Pompa-García and Venegas-González [26] found a negative correlation with maximum temperature during the spring season. In our study, we did not find a noticeable influence of temperature but rather soil moisture was important during the late growing season (Figure 6b). MXD is related to cell differentiation processes (thickening or lignification) in latewood tracheids at the end of the summer season [75]. The negative response found between MXD and soil moisture suggests that the growth of P. hartwegii is controlled by soil moisture at the lowest elevation (3500 m asl) which drives the thickening of cell walls and the production of total tracheid lumens [11,76]. The soil water content affects lumen size by influencing the turgor pressure during the tracheid's differentiation process [77,78] On the other hand, the positive correlations between MID and temperature (maximum and minimum) suggest that the higher values in MID are accelerated by high temperatures with an exception at the highest elevation (3900 m asl). In this regard, Camarero, et al. [11] found similar results for Juniperus thurifera with positively (temperature) and negatively (rainfall) correlations with MID during the May-June period.

Mountain Dependent Associations between Wood Density and MODIS Variables
During the last decades, a major effort has been made to evaluate the dynamics of wood formation (xylogenesis) at a tracheid level and assessing the rates of differentiation (cell enlarging and wall thickening), however, these studies are spatially restricted [79,80]. Thus, to assess xylogenesis at larger spatio-temporal scales, studies have explored relations between remotely sensed variables and the carbon allocation process. For example, several studies have correlated tree-ring width records with NDVI across different spatial and temporal scales [34][35][36]38]. However, less frequently, tree-ring density variables have also been explored [21,40,41]. From these studies, we can conclude that photosynthetic activity and carbon allocation (e.g., MXD) are broadly linked although they deal with different aspects of plant physiology.
Documented links between MXD and NDVI at northern latitudes showed that summer NDVI (May and June) is important to explain xylogenesis as related to the vigor canopy and growing season effect on the cell differentiation process [21,40]. However, we did not find the May-June period relevant to wood density, instead, we found that the previous winter until spring period showed a strong positive correlation (highest correlation in doy 065-March) with NDVI at TLA but not at JOC (Figure 7). This relationship implies that higher NDVI values during the first months of the year are translated into denser woody tissues at the end of the growing season but not all sites are suitable to establish robust links. It is well known that there is a lag time between processes, first with photosynthesis and then with the carbon allocation process. Thus, several intra and inter-specific factors could modify the first signal captured by remote sensed information such as edaphoclimatic conditions or the age of trees, that can alter the priority of carbon allocation, leaf area and light use efficiency [81]. For example, we saw a poor association at JOC for almost all remotely sensed variables, thus, although both mountains are strongly temperature-limited, JOC is probably restricted by another factor. The negative correlation observed between NDVI and MID (highest correlation in doy 017-January) indicates that high photosynthetic activity (NDVI) at this time is matched with low MID. Perhaps because carbon allocation at the start of the growing season is primarily directed towards shoot growth (leaf area) and the start of cambial activity, wood cell division and enlargement (i.e., wider xylem cells) have priority over cell wall thickening. Mäkinen, et al. [82] reported a decrease in wood density when well-watered and optimal nutrients conditions are supplied.
The NDVI is an appropriate indicator of canopy vigor and thus the health status of plants [83]. The fact that the MXD time series at TLA mirrored the NDVI activity but not at JOC, reinforces the idea that more studies are needed to improve our understanding of the relationship between the canopy and cambial activity at local scales and the potential impact of climate change on forests tree growth and productivity.

Conclusions
We demonstrated that two high-elevation forests closely located on different mountain sites and dominated by the same tree species (Pinus hartwegii Lindl.) could have statistically similar wood density profiles, but different wood density trajectories suggesting different behaviors in the wood formation process. Tree-ring wood adjustments coupled with a decrease in productivity (previously reported) can provide early-warning signals of stressed forests as water transport efficiency can be affected as a result of climate variation. Compared to MXD, we reported higher correlations between MID and temperature. A set of significant correlations between MXD and the Normalized Difference Vegetation Index-NDVI indicates the potential use of the link between wood micro-density and remotely sensed information to understand forest response climate change effects but not at all site conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/4/452/s1, Figure S1. Scatter plot between tree-ring width and tree-ring wood density by Mountain. The solid line represents a generalized additive model fitted by Mountain. Figure S2. Empirical cumulative distribution function (CDF) for adjusted tree-ring wood density variables by mountain. Figure S3. Spatial correlation (p < 0.05) between maximum latewood density-MXD, average density-AVE and minimum latewood density-MID and normalized difference vegetation index-NDVI at Tlaloc Mountain. "DOY" represents the day of the year for the current year while the word "doy" in lowercase represents the day for the previous year. Funding: This research received funding from The National Council for Research and Technology (CONACyT) and the ECOS NORD Program for French-Mexican collaboration, under the project number 0276777 "Impacto del cambio climático en la adaptación y plasticidad fenotipica del crecimiento en árboles forestales"-M16A01 "Impact du changement climatique sur la plasticité phénotypique de la microdensité du xylème" and the J. G. Boswell Endowed Chair in Soil Science.