Wood Anatomical Traits Respond to Climate but More Individualistically as Compared to Radial Growth: Analyze Trees, Not Means

: Wood encodes environmental information that can be recovered through the study of tree-ring width and wood anatomical variables such as lumen area or cell-wall thickness. Anatomical variables often provide a stronger hydroclimate signal than tree-ring width, but they show a low tree-to-tree coherence. We investigate the sources of variation in tree-ring width, lumen area, and cell-wall thickness in three pine species inhabiting sites with contrasting climate conditions: Pinus lumholtzii in wet-summer northern Mexico, and Pinus halepensis and Pinus sylvestris in dry-summer north-eastern Spain. We quantiﬁed the amount of variance of these three variables explained by spring and summer water balance and how it varied among trees. Wood anatomical variables accounted for a larger inter-individual variability than tree-ring width data. Anatomical traits responded to hydroclimate more individualistically than tree-ring width. This individualistic response represents an important issue in long-term studies on wood anatomical characteristics. We emphasized the degree of variation among individuals of the same population, which has far-reaching implications for understanding tree species’ responses to climate change. Dendroclimatic and wood anatomical studies should focus on trees rather than on the mean population series.


Introduction
Ongoing global warming is expected to impair natural ecosystems by affecting forest growth and productivity-both directly and indirectly-through changes in temperature, rainfall amount and precipitation seasonal patterns [1]. Climatic stressors such as drought may also be exacerbated by climate warming and play an essential role in driving forest health [2]. Therefore, understanding how forests respond to those extreme climatic events is key to forecast how their carbon uptake capacity will change as climate warms and becomes more arid.
The study of growth rings of woody plants such as trees growing in seasonally dry forested biomes has allowed assessing how water shortage impacts wood production and forest productivity [3]. Tree-ring data are well-established environmental proxies that can be used to infer past growth conditions and glimpsing regional and local hydroclimate impacts on forests [4]. Therefore, the analysis of tree-ring width patterns has long been an excellent tool used by dendroclimatologists and dendroecologists to reconstruct hydroclimate variability and related forest dynamics [5]. For instance, tree-ring width data have allowed assessing drought impacts on radial growth and resilience of forests on a global scale [6,7]. Annually resolved tree-ring records, often coupled with remotely sensed indices, enabled the characterization of ontogenetic growth patterns [8] as well as the assessment of the combined negative effects of outbreaks and severe droughts on tree growth and mortality rates in drought-prone Mediterranean forests [9][10][11].
The analyses of xylem anatomical characteristics within tree rings such as lumen area and cell-wall thickness have long been thought to be an interesting tool for extracting annually to seasonally resolved environmental information from tree rings [12]. In recent decades, quantitative wood anatomy has been widely exploited in transdisciplinary global change research as a new class of metrics, opening an interesting prospect for dendroclimatological and dendroecological investigations including drought impacts on tree growth [13][14][15]. Compared to traditional tree-ring width data, the analysis of xylem features yields multiple advantages. First, anatomical variables are directly linked to the allometry and functioning of trees through processes such as tapering, hydraulic conductivity, and carbon allocation, which cannot be obtained from tree-ring width data [16]. Second, matching intra-ring cell position with the time of formation allows deciphering at finer temporal resolutions the relationships between growth and xylogenesis with climate at intra-annual resolution [17,18]. By matching seasonally or even monthly climate with information retrieved at the intra-annual scale from tree rings, we have expanded our understanding of the physiological mechanisms behind xylem vulnerability to drought [17,19] and of the impacts of climatic and biotic stressors on growth decline and tree mortality, particularly in areas subjected to seasonal water shortage [13,20].
Although wood anatomical features give detailed mechanistic insight into tree functioning and are more sensitive to specific intra-seasonal environmental factors than tree-ring proxies, they still present drawbacks, which limit its potential in dendroclimatology when compared to tree-ring variables. For instance, only lately have methodological advances such as improvements in computer speed and automated image-analysis systems enabled the processing and interpretation of large xylem cell datasets [21]. Such time-consuming procedures have impeded measurement strategies of wood anatomical properties, resulting in the limited replication of trees, rings, and cells per ring tested. However, a limitation probably even more relevant is the low between-tree correlation of wood anatomical features as compared with tree-ring width data [22], which could partly outweigh their hydroclimate signal and limit its application as climatic or ecological proxies. So, contrasting results emphasize relatively strong or weak common climatic signals found in wood anatomical variables which could depend on tree species, bioclimate, and number of sampled trees [22][23][24]. Such low common signal among conspecific and coexisting trees but strong responses to climate require additional investigations to validate the usefulness of wood anatomical parameters as environmental proxies.
The likely influences of both physiological processes and environmental factors on xylogenesis are plausible causes of the poor inter-tree correlation of wood anatomical series [25]. A particularly weak inter-series correlation for anatomical chronologies compared to ring-width series may also be related to a low year-to-year variability of xylem traits [22]. On the other hand, it was reported that a relatively weak inter-tree correlation or common signal in anatomical parameters can decouple them from intraannual climate variability, e.g., [15,26]. Nevertheless, it is widely accepted for numerous tree species and wood features that series of wood anatomical features can provide a strong hydroclimatic signal despite their low coherence among individuals. It remains a challenge to better quantify that coherence or similarity as a first step to perform attribution analyses and search for drivers of the high variability in wood anatomical responses among co-occurring trees. Here, we aimed to present statistical procedures designed to quantify: (i) the coherence among individuals of tree-ring and wood anatomical series, and (ii) the strength of the climatic signal contained in those series. We studied three pine species growing in contrasting bioclimates to test the hypothesis that the climatic signal is not related to inter-individual correlation in the tree variables. We expect a stronger tree-to-tree correlation in tree-ring width data as compared with wood anatomy data, but a higher response to climate of wood anatomy series. We compare the climate signal of growth and wood anatomy in Pinus lumholtzii L. Rob. and Fernald, Pinus halepensis Mill. and Pinus sylvestris L., three species which showed a marked sensitivity to seasonal drought in Mexico [27][28][29] and Spain [13,17,30]. We compare these three species because they encompass wide climatic and ecological gradients, i.e., from mountain mesic sites (P. sylvestris) to semi-arid (P. halepensis) to warm temperate conditions (P. lumholtzii). In addition, they withstand different dry seasons, with the typical Mediterranean summer drought affecting P. sylvestris and P. halepensis, whilst P. lumholtzii experiences summer monsoon rainfall and dry winter and spring. We aimed at examining the magnitude of intraspecific variability and the relationships between seasonal climate and tree-ring width and two major anatomical variables, transversal conduit diameter and cell-wall thickness. We focused on individual variability analyzed new data in the case of P. lumholtzii and re-analyzed previously taken data for P. sylvestris and P. halepensis coming from previous works which focused on mean species' series. Such combinations of contrasting climate conditions and species allowed us to address the following questions: (i) how variable are quantitative wood anatomical time series among coexisting trees of the same species as compared with tree-ring width data?, (ii) how is intraspecific variability structured between individuals and populations in the case of wood anatomy and tree-ring width data?, and (iii) do time series of anatomical traits contain a stronger individual hydroclimate signal than series of ring width?

Study Sites and Tree Species
We selected three pine species subjected to seasonal drought and inhabiting Spain and Mexico ( Figure 1). First, Aleppo pine (Pinus halepensis Mill.) was sampled from a mixed pine-juniper Mediterranean forest situated in the Vedado de Peñaflor, Aragón, northeastern Spain (41 • 47 N, 0 • 43 W, 560 m a.s.l.), in the semi-arid Middle Ebro Basin. The tree cover is 55%. Gypsum and marls comprise the parent rock material, while soils are Regosol with loamy texture of the subsurface horizons. The topography of the terrain consists of small hills and plateaus at or below 5 degrees. Second, Scots pine (Pinus sylvestris L.) wood samples were collected in a mixed pine-oak-juniper forest located in Corbalán, Aragón, northeastern Spain (40 • 26 N, 0 • 58 W, 1303 m a.s.l). The tree cover is 58%, the terrain slope is 25 degrees, and soils are Cambisol with loamy texture. The P. halepensis and P. sylvestris sites were selected because they were severely impacted by the 2012 drought [31]. Third, stem wood samples were collected from sad or Lumholtz's weeping pine (Pinus lumholtzii B.L. Rob. and Fernald) trees growing in a pine-oak forest located in Sierra Madre Occidental, Chihuahua state, northern México (27 • 04 N, 107 • 08 W, 1390 m a.s.l.). The site vegetation is formed by mixed pine and oak-pine forests with junipers and other shrubs [32] with 60% crown cover. Terrain slope is 1.7 degrees, and soils are mainly Leptosols and Regosols of volcanic origin.

Climate Data
For all study sites, the climate data were obtained from the TerraClimate dataset (Available online: http://www.climatologylab.org/terraclimate.html, accessed on 2 May 2022), covering the period 1958-2020. TerraClimate is a high-spatial-resolution (~4 km) global dataset of monthly temperature, precipitation, and climatic water balance (P-PET, difference between precipitation and potential evapotranspiration), which is a combination of data from WorldClim and CRU TS4.0 climate datasets [33]. Selected study sites have a sharp difference in climate regimes. Note that all sites are subjected to severe seasonal drought and have negative annual water balances. The climate at the northeastern Spanish sites is typical Mediterranean climate characterized by dry summers and mild, wet winters. In particular, based on climate data for the period 1970-2020 ( Figure  2), the climate at the site at Vedado de Peñaflor (Pinus halepensis) is continental Mediterranean with an average yearly rainfall of 400 mm (semi-arid conditions), the lowest rainfall recorded in July. Annual mean temperature was 14.9 °C with maximum and minimum mean monthly temperatures of 31 °C (July) and 2.3 °C (January), respectively. The estimated annual water balance was −505 mm with water deficit occurring from May to September. At the Corbalán site (Pinus sylvestris), the climate is continental Mediterranean (Figure 2) with a mean annual temperature of 10.3 °C, with frosts occurring from December to January, and maximum and minimum mean monthly temperatures of 25.7 °C (July) and −1.4 °C (January), respectively. Total precipitation was 481 mm, with water deficit occurring from late June to end August. The estimated annual water balance was −353 mm. The climate in the Pinus lumholtzii study area is sharply contrasted with the European sites: climate is temperate with dry spring and hot but wet summers because of the North American Monsoon influence. The mean annual temperature was 15.2 °C, and the annual precipitation sum was 781 mm, with 46% of this precipitation falling from July to August, and the annual water balance was −400 mm.

Climate Data
For all study sites, the climate data were obtained from the TerraClimate dataset (Available online: http://www.climatologylab.org/terraclimate.html, accessed on 2 May 2022), covering the period 1958-2020. TerraClimate is a high-spatial-resolution (~4 km) global dataset of monthly temperature, precipitation, and climatic water balance (P-PET, difference between precipitation and potential evapotranspiration), which is a combination of data from WorldClim and CRU TS4.0 climate datasets [33]. Selected study sites have a sharp difference in climate regimes. Note that all sites are subjected to severe seasonal drought and have negative annual water balances. The climate at the north-eastern Spanish sites is typical Mediterranean climate characterized by dry summers and mild, wet winters. In particular, based on climate data for the period 1970-2020 (Figure 2), the climate at the site at Vedado de Peñaflor (Pinus halepensis) is continental Mediterranean with an average yearly rainfall of 400 mm (semi-arid conditions), the lowest rainfall recorded in July. Annual mean temperature was 14.9 • C with maximum and minimum mean monthly temperatures of 31 • C (July) and 2.3 • C (January), respectively. The estimated annual water balance was −505 mm with water deficit occurring from May to September. At the Corbalán site (Pinus sylvestris), the climate is continental Mediterranean ( Figure 2) with a mean annual temperature of 10.3 • C, with frosts occurring from December to January, and maximum and minimum mean monthly temperatures of 25.7 • C (July) and −1.4 • C (January), respectively. Total precipitation was 481 mm, with water deficit occurring from late June to end August. The estimated annual water balance was −353 mm. The climate in the Pinus lumholtzii study area is sharply contrasted with the European sites: climate is temperate with dry spring and hot but wet summers because of the North American Monsoon influence. The mean annual temperature was 15.2 • C, and the annual precipitation sum was 781 mm, with 46% of this precipitation falling from July to August, and the annual water balance was −400 mm.

Tree-Ring Width and Wood Anatomical Data
Original plant material of P. lumholtzii was collected in Mexico [27] and Spain for P. halepensis and P. sylvestris [13,17,30]. Briefly, to build site-specific chronology, stems of dominant trees were cored at breast height, i.e., at approximately 1.3 m above ground using either 5 or 10 mm Pressler increment borer to take 2 or 1 cores, respectively, per tree for tree-ring and anatomical analyses. The 5 mm radial wood cores were fist checked to exclude reaction wood, eccentricities, or wounding and then air-dried, glued onto wooden supports, and sanded with sandpapers of increasing grain until ring boundaries were conspicuous. Then, they were visually cross-dated, and their tree-ring widths (TRW) were measured under a binocular scope using a Lintab-TSAP measuring device (Rinntech, Heidelberg, Germany). The visual cross-dating was validated using the COFECHA software which calculates moving correlations between individual tree-ring width series and the site mean series [34]. for Pinus lumholtzii (f). In the in the climate diagrams, monthly precipitation below 100 mm is scaled 2:1 with mean monthly temperature (vertically hatched) and above 100 mm 15:1. Turquoise bars below the x-axis show the months where frosts can occur (when absolute monthly minimums are equal or lower than 0 °C).

Tree-Ring Width and Wood Anatomical Data
Original plant material of P. lumholtzii was collected in Mexico [27] and Spain for P. halepensis and P. sylvestris [13,17,30]. Briefly, to build site-specific chronology, stems of dominant trees were cored at breast height, i.e., at approximately 1.3 m above ground using either 5 or 10 mm Pressler increment borer to take 2 or 1 cores, respectively, per tree for tree-ring and anatomical analyses. The 5 mm radial wood cores were fist checked to exclude reaction wood, eccentricities, or wounding and then air-dried, glued onto wooden supports, and sanded with sandpapers of increasing grain until ring boundaries were conspicuous. Then, they were visually cross-dated, and their tree-ring widths (TRW) were measured under a binocular scope using a Lintab-TSAP measuring device (Rinntech, Heidelberg, Germany). The visual cross-dating was validated using the COFECHA software which calculates moving correlations between individual tree-ring width series and the site mean series [34].
For the scope of this study, we selected five mature trees for P. halepensis and P. sylvestris sampled during the 2011-2012 and 2013-2014 winters, respectively, and four mature trees for P. lumholtzii sampled in the summer 2013. Selected trees were well crossdated and covered similar time spans. We divided the cores into 3-5 cm long blocks to facilitate further processing. Then, we used a rotary microtome (Leica RM 2255; Leica Microsystems, Germany) to obtain transversal cross-sections (15-20 μm thick) of P. halepensis and P. sylvestris, while for P. lumholtzii samples, we used a sledge core microtome [35]. Wood sections were stained by mixing safranin (1%) and astra blue (0.5%) solutions. They were then fixed and permanently mounted onto glass microscope slides using a synthetic resin (Eukitt™; Merck, Darmstadt, Germany). Digital images were captured at 40-100× magnification using a digital camera attached to a light microscope (Olympus BH2). We created panoramas by stitching together multiple overlapping for Pinus lumholtzii (f). In the in the climate diagrams, monthly precipitation below 100 mm is scaled 2:1 with mean monthly temperature (vertically hatched) and above 100 mm 15:1. Turquoise bars below the x-axis show the months where frosts can occur (when absolute monthly minimums are equal or lower than 0 • C).
For the scope of this study, we selected five mature trees for P. halepensis and P. sylvestris sampled during the 2011-2012 and 2013-2014 winters, respectively, and four mature trees for P. lumholtzii sampled in the summer 2013. Selected trees were well cross-dated and covered similar time spans. We divided the cores into 3-5 cm long blocks to facilitate further processing. Then, we used a rotary microtome (Leica RM 2255; Leica Microsystems, Wetzlar, Germany) to obtain transversal cross-sections (15-20 µm thick) of P. halepensis and P. sylvestris, while for P. lumholtzii samples, we used a sledge core microtome [35]. Wood sections were stained by mixing safranin (1%) and astra blue (0.5%) solutions. They were then fixed and permanently mounted onto glass microscope slides using a synthetic resin (Eukitt™; Merck, Darmstadt, Germany). Digital images were captured at 40-100× magnification using a digital camera attached to a light microscope (Olympus BH2). We created panoramas by stitching together multiple overlapping images using PtGui free software (PTGui. Available online for download: https://ptgui.com/, accessed on 14 February 2022) ( Figure S1).
Tree-ring width and wood anatomy data were analyzed for the periods 1970-2012, 1950-2010, and 1970-2013 in the case of P. halepensis, P. sylvestris and P. lumholtzii, respectively (see Table 1). Images of P. halepensis and P. sylvestris wood sections were processed using the ROXAS software [21], while images of P. lumholtzii wood sections were processed using the recently developed AutoCellRow software [36]. These analyses provided measurements of transversal conduit area and double cell-wall thickness of tracheids performed on radial direction within the dated annual ring. For each annual ring, we then calculated the mean conduit area (MCA) and cell-wall thickness (CWT). Table 1. Statistical characteristics of raw tree-ring width and xylem trait chronologies (Rbar and mean EPS statistics were calculated from the detrended chronology). SD is the standard deviation; AR1, first-order autocorrelation coefficient calculated from raw data; Rbar is the average inter-series correlation between all series from different trees; EPS is the expressed population signal; MS is the mean sensitivity of the raw time series; variables: TRW, tree-ring width (in mm); MCA, mean conduit area (in µm 2 ); CWT, cell-wall thickness (in µm).

Statistical Analyses
Several descriptive statistics widely used in tree-ring sciences [4] were calculated on raw data series, including mean, standard deviation, first-order autocorrelation (AR1), the mean sensitivity (MS), mean series intercorrelation (Rbar), and the expressed population signal (EPS) MS measures relative changes in a variable among consecutive rings. Rbar gives the measure of the common signal between time series [37]. The EPS, which evaluates how well-replicated is the mean series or chronology as compared with a theoretical, wellreplicated chronology [38]. Dendro statistics were calculated with the dplR R package [39].
We used linear mixed models (LMMs; [40]) to explore the species-specific variation in the considered traits as well as the effects of seasonal water balance on this variation. This modeling approach is an excellent tool to represent the hierarchical sampling design and has found wide application in ecological studies in recent years as an excellent tool to decompose the trait variation within grouped data (e.g., [41,42]).
We used a double-step approach: first, we built intercept only linear mixed effect models for each of the trait characteristics defined above, i.e., TRW, MCA, and CWT, which was entered as response variables to explore the variability attributed to the clustering species, individual trees, and inter-annual temporal variation. Such a design was implemented into the modeling approach that contained only tree species, tree ID, and cambial age as random effects, controlling variation in the intercept of the resulting models. To well represent the study design within the model structure, we nested tree ID within the tree species clustering variable. In all fitted models, the response variables were log-transformed (log(var+1)) to achieve normality assumptions.
By including these grouping factors as random effects in the models, we investigated the intraclass correlation coefficient (ICC) accounted by each factor, which is meaningful to understand how much of the overall variation in the response variable is explained simply by the clustering groups. The ICC is calculated by dividing the random effect variance, σ 2 i , by the total variance, i.e., the sum of the variance attributable to the random effect and the residual unexplained variance, σ 2 ε derived from the models using the VarCorr function in R. Being σ 2 Species , σ 2 Tree , σ 2 Age the variance accounted by the species, tree ID and cambial age random variable, respectively, the ICC for the Species group takes the form: and so on for all the clustering variables defined above. The ICC can be interpreted as the proportion of the variance explained by the grouping structure in the population, and it is generally between 0 and 1.0, where higher values reflect greater between-group variability.
Aside from the intercept-only effect model, we created additional models for each species separately by including seasonal water balance (P-PET) as fixed effects while removing the variable species from the random effects to test the role of climate variation within and between the grouping characteristics defined above and how this variation influenced xylem traits. In detail, for each tree species, LMMs were fitted to assess the effect of spring (March-April-May) P-PET MAM on TRW and MCA, and summer (June-July-August) P-PET JJA on CWT. The selection of that seasonal climatic predictors was based on the results of prior studies on the same forests showing a significant effect of spring hydroclimate conditions on radial growth and wood anatomy and a positive relationship between summer climate and CWT [13,17,27,29,30]. Individual trees and cambial age were considered random effects to account for repeated measures within a site. Before analysis, the response variable was log-transformed (log(trw+1)) to achieve normality assumptions; then, all fixed terms were centered and scaled to improve parameter estimates and allow direct comparisons of the regression coefficients. For all selected species, the period from 1970 to the most recent year was considered for the analysis.
We investigated the percentage of variance explained by each factor by quantifying the proportion of variance explained by fixed and random terms we calculated a full partitioning of variance into three components: variance attributable to climate predictors via fixed slopes, variance attributable to random terms, and unexplained variance attributable to residuals. Variances were estimated using restricted maximum likelihood (REML). From these models, marginal and conditional R 2 values [43] were derived to examine the variation explained by fixed and fixed plus random factors, respectively, using the r.squaredGLMM function in the MuMIn package [44]. A residual diagnosis was performed to check the validity of the model assumptions (normality and homoscedasticity of residuals). We calculated the variance components on these random-effects models using the r2mlm function from the r2mlm [45,46] R package.
All statistical analyses were carried out with R v4.0.4 [47], and parameter estimation was carried out by using the lme4 library [48,49], which provides unbiased variances for small sample sizes using an unstructured error covariance, and the lmerTest package to obtain p-values.

Statistics and Variability of Radial Growth and Wood Anatomy
The widest and narrowest rings were formed by P. halepensis and P. sylvestris, respectively, while the largest MCA and CWT values were found in P. lumholtzii (Table 1 and Figure 3). TRW presented the highest AR1 values, whereas CWT showed the lowest MS values. The MS was always higher for TRW than for wood-anatomical variables, and it reached a maximum value of 0.44 in P. sylvestris. The standard deviation (SD) and MS of MCA were the highest in P. sylvestris, which presented an elevated year-to-year variability in lumen area (Figure 3), with sharp drops corresponding to dry periods. For instance, there was a strong MCA reduction (ca. −50%) in response to the severe 2005 drought (Figure 3). The high SD value found for TRW in P. halepensis is probably due to the different ages of selected trees. Within each species, Rbar values were always higher for TRW and varied from 0.39 in P. lumholtzii to 0.69 in P. sylvestris, but in P. halepensis, the TRW Rbar (0.45) was lower than the MCA Rbar (0.61). In contrast, the MCA Rbar was very low (0.04) in P. lumholtzii. In the case of CWT, the Rbar ranged from 0.11 (P. lumholtzii) to 0.23 (P. sylvestris) ( Table 1). Most EPS were below the 0.85 threshold suggested by Wigley et al. [38] for well-replicated series, except MCA in P. halepensis and TRW in P. sylvestris.
for TRW and varied from 0.39 in P. lumholtzii to 0.69 in P. sylvestris, but in P. halepensis, the TRW Rbar (0.45) was lower than the MCA Rbar (0.61). In contrast, the MCA Rbar was very low (0.04) in P. lumholtzii. In the case of CWT, the Rbar ranged from 0.11 (P. lumholtzii) to 0.23 (P. sylvestris) ( Table 1). Most EPS were below the 0.85 threshold suggested by Wigley et al. [38] for well-replicated series, except MCA in P. halepensis and TRW in P. sylvestris.

Trait Variance Decomposition between Species and Individuals
Inspecting intercept only models we found, in terms of the proportion of the variance explained by the grouping structure aka ICC, that the between-species variability played a significant role in the overall variation in TRW (47%) and MCA (37%), and that variance in CWT (75%) was almost completely explainable by this parameter (Figure 4). This means

Trait Variance Decomposition between Species and Individuals
Inspecting intercept only models we found, in terms of the proportion of the variance explained by the grouping structure aka ICC, that the between-species variability played a significant role in the overall variation in TRW (47%) and MCA (37%), and that variance in CWT (75%) was almost completely explainable by this parameter (Figure 4). This means that species-specific differences play an important role in the cell wall thickness variability. Furthermore, the percentage of variance explained by variation among tree individuals in all traits was relatively smaller for TRW and CWT compared to MCA, indicating a much higher coherence among trees in these variables. that species-specific differences play an important role in the cell wall thickness variability. Furthermore, the percentage of variance explained by variation among tree individuals in all traits was relatively smaller for TRW and CWT compared to MCA, indicating a much higher coherence among trees in these variables. . Species, the clustering variable including P. halepensis, P. sylvestris, and P. lumholtzii; Tree, individual trees; Age, cambial age. ICC was defined as σ 2 i/(σ 2 i+σ 2 ε), using mixed model random variance terms, and represents the consistency of behavior across multiple measurements, where σ 2 i is the random effect variance term and σ 2 ε is the error variance term from the LMMs. Unexplained variance attributable to residuals is also reported. . Species, the clustering variable including P. halepensis, P. sylvestris, and P. lumholtzii; Tree, individual trees; Age, cambial age. ICC was defined as σ 2 i /(σ 2 i +σ 2 ε ), using mixed model random variance terms, and represents the consistency of behavior across multiple measurements, where σ 2 i is the random effect variance term and σ 2 ε is the error variance term from the LMMs. Unexplained variance attributable to residuals is also reported.

Effects of Seasonal Climate on Measured Traits and Variance Decomposition
Through fixed-effect mixed modeling fitted for each species separately, we were able to attribute an amount of variance in TRW, MCA, and CWT to seasonal water balance variation (fixed variable). Model results are presented in the Supplementary Materials (Tables S1-S3), and model-predicted relationships between tree-ring width and wood anatomical variables and P-PET by species are shown in Figure 5. Overall, variations in TRW were well-explained by P-PET MAM but with a low amount of explained variance, with some species-specific peculiarities, while distinct outcomes were found for wood anatomy models. variation (fixed variable). Model results are presented in the Supplementary Materials (Tables S1-S3), and model-predicted relationships between tree-ring width and wood anatomical variables and P-PET by species are shown in Figure 5. Overall, variations in TRW were well-explained by P-PETMAM but with a low amount of explained variance, with some species-specific peculiarities, while distinct outcomes were found for wood anatomy models.  For P. halepensis, a significant positive effect of P-PET MAM was detected only for TRW, while no significant effects were detected for MCA and CWT. However, for TRW, the proportion of variance explained by the fixed factor alone (i.e., marginal R 2 ) was low (0.03), as well as the proportion of variance explained by both the fixed and random factors, i.e., the conditional R 2 , which was 0.59 (Table S1). The sources of variation varied notably across measured traits and pine species (Figure 6) as well. For example, 7% of the variance in TRW and MCA of P. halepensis was due to differences between individuals. In contrast, this variance corresponded to 13% in the case of CWT, indicating a lower coherence among trees in this variable. Hight inter-annual temporal age variability was associated with TRW and MCA, while it was negligible for CWT, leaving more than 78% of trait variance unexplained. The variability between P. sylvestris individuals represented a large part of the variance (excluding unexplained residual variance) in MCA, accounting for most one-third of the total variance (38%), much higher as compared to TRW (9%) and CWT (17%) models. For this species, we found a positive effect of P-PET MAM on both TRW and MCA and P-PET JJA on CWT. All models shared a lower marginal R 2 compared to the conditional R 2 , i.e., they fall below 0.08 (Table S2).
variance in TRW and MCA of P. halepensis was due to differences between individuals. In contrast, this variance corresponded to 13% in the case of CWT, indicating a lower coherence among trees in this variable. Hight inter-annual temporal age variability was associated with TRW and MCA, while it was negligible for CWT, leaving more than 78% of trait variance unexplained. The variability between P. sylvestris individuals represented a large part of the variance (excluding unexplained residual variance) in MCA, accounting for most one-third of the total variance (38%), much higher as compared to TRW (9%) and CWT (17%) models. For this species, we found a positive effect of P-PETMAM on both TRW and MCA and P-PETJJA on CWT. All models shared a lower marginal R 2 compared to the conditional R 2 , i.e., they fall below 0.08 (Table S2). Figure 6. Variance in tree-ring width (TRW) and wood anatomical variables (MCA, mean conduit area; CWT, cell-wall thickness) explained by variation in seasonal water balance (PPET) between individual trees, and age for P. halepensis, P. sylvestris, and P. lumholtzii. Tree represents individual trees, and Age is the individual cambial age. Unexplained variance attributable to residuals is also reported.
For P. lumholtzii, a significant positive effect of P-PETMAM was detected only for TRW, while no significant regression coefficients were detected for MCA and CWT. For TRW, the marginal and conditional R 2 values were 0.03 and 0.52, respectively (Table S3). Without taking into account the unexplained residual variance, for TRW, MCA, and CWT, we identified differences among tree individuals, as the strongest contributor to this variation, accounting for 24%, 28%, and 31% of total traits variability, respectively. Interannual temporal age variability was more influential on TRW (24%) compared to the other traits.
When we included the seasonal water balance, i.e., from spring to winter in the mixed models, we observed a general increase in the percentage of explained variance by the fixed factors in both the P. sylvestris and P. halepensis for the variables TRW and MCA and a weak increase in explained variance for all traits in P. lumholtzii ( Figure S3). Figure 6. Variance in tree-ring width (TRW) and wood anatomical variables (MCA, mean conduit area; CWT, cell-wall thickness) explained by variation in seasonal water balance (PPET) between individual trees, and age for P. halepensis, P. sylvestris, and P. lumholtzii. Tree represents individual trees, and Age is the individual cambial age. Unexplained variance attributable to residuals is also reported.
For P. lumholtzii, a significant positive effect of P-PET MAM was detected only for TRW, while no significant regression coefficients were detected for MCA and CWT. For TRW, the marginal and conditional R 2 values were 0.03 and 0.52, respectively (Table S3). Without taking into account the unexplained residual variance, for TRW, MCA, and CWT, we identified differences among tree individuals, as the strongest contributor to this variation, accounting for 24%, 28%, and 31% of total traits variability, respectively. Inter-annual temporal age variability was more influential on TRW (24%) compared to the other traits.
When we included the seasonal water balance, i.e., from spring to winter in the mixed models, we observed a general increase in the percentage of explained variance by the fixed factors in both the P. sylvestris and P. halepensis for the variables TRW and MCA and a weak increase in explained variance for all traits in P. lumholtzii ( Figure S3).

Discussion
Interpreting the climatic effects on Mediterranean woody species by mean series of tree rings and xylem functional traits is of paramount importance to assess the potential response of tree species to changing climate. However, compared to tree ring data, the functional interpretation of climate-xylem trait relationships still presents many drawbacks, which limit its potential in dendroclimatology when compared to tree-ring variables. In this study, we found that wood anatomical time series account for a large intraspecific and intraindividual variability as compared with tree-ring width data, and anatomical traits respond to hydroclimate, in terms of shape and strength, more individualistically than tree-ring width data.
The patterns in the variance components of the key tree-ring and wood-anatomy traits analyzed in this study are rich in information in many respects, and some of the findings from these results are unexpected. There are three main elements that we would like to emphasize: (1) explained variance of traits is not fairly distributed across all investigated species and traits, (2) the total amount of within-species variation depends on the con-sidered trait, and (3) the proposed models explained only a low percentage of the total observed variance.

The Percentage of Variance Explained for Each Trait Differs between Tree Species
Our results reflect the patchwork of the relationship between climate and growth, or wood anatomy and variation explained by models differing between species and traits well delineated by prior works. Overall, we expected a significant positive relationship between spring water balance (P-PET MAM ) and radial growth (TRW) for all droughtsensitive investigated species, although we found a low percentage of explained variance by the models due to individual variability in the responses. The negative impact of limited water availability on these species has been exacerbated due to the increase in growingseason temperatures and evapotranspiration rates over the past decades, showing greater sensitivity to drought in dry regions of Spain and Mexico, e.g., P. sylvestris populations in Camarero et al. [50] and Bose et al. [51]. Significant wood anatomical responses to spring-summer water balance controlling earlywood lumen area, and therefore hydraulic conductivity, and CWT in P. sylvestris has been also highlighted in previous works, e.g., [13]. The significant relationships previously reported between climate (water balance) and xylem traits in P. halepensis and P. lumholtzii based on mean series [17,[28][29][30] were not observed considering the full spectrum of individual variability. However, there are some significant differences between the three investigated species and the two study sites that are worth emphasizing as we compared very different seasonally dry biomes subjected to wet (Mexico) or dry (Spain) summer conditions with different impacts on xylem development [17,30]. In fact, including the effect of the water balance in the fitted models, we remarked a significant increase in the variance explained in tree-ring width and lumen area of P. sylvestris and P. halepensis ( Figure S3), suggesting coherent drought signals in the trait patterns of the species, typical of Mediterranean dry biomes, e.g., [50]. It is likely instead that the low explained variance by fixed climate in all traits of P. lumholtzii is attributable to the weak seasonal climate signal of considered variables in our analysis. The growth rate of P. lumholtzii is indeed shown to be markedly related to the water balance of the prior year wet seasons [27].
A non-uniform distribution of the climate response among species suggests the role of non-common ecological processes in determining trait variability among individuals at each study site. For instance, the limiting effect of reduced water availability on growth and tracheid expansion in one species in its environment is likely to not have the same implications on another species growing in another environment due to different climatic or soil conditions. This is because, beyond the possible morpho-anatomical adjustments faced differentially by the species (see below discussion on the phenotypic plasticity), it is probable that the limiting effect of an environmental variable on wood formation can be amplified or smoothed by the joint effect of other site-specific-abiotic or biotic-variables such as soil water retention capacity, tree size and allometry (e.g., sapwood/leaf area ratio, rooting depth) or tree-to-tree competition, among others.
These differences would partly explain the common reduced variance explained by the models for the three pine species considered. This would mean that, on average, all components of traits variation ought to be considered in any study to disentangle the species' behavior in response to climate, although in practice, this is often not feasible. It is thus critical to decide which components of variation should be taken into account considering the ecological question addressed. For instance, among the source of variability of wood traits, the cambial age or the tree height [52] may play significant roles in trait variation for different tree species ( Figure 6, and variance decomposition in Figure S2 with age as fixed variable). Additionally, unfiltered data might include high first-order autocorrelation, indicating that radial growth and xylem traits throughout were strongly influenced by conditions in the preceding year, which can partly outweigh seasonal climatic signals.

The Percentage of Variance Explained Differs between Traits and Individuals
The percentage of explained variance differs between traits within a single species as found for P. halepensis and P. lumholtzii, in which cases tree-ring width was more informative than wood-anatomical variables (MCA, CWT). This is consistent with many studies that found that earlywood vessel cross-sectional lumen area and ring width in conifer and broadleaved species contained unique climate information [15,22], but that this variable was not very sensitive, or was deemed too weak when compared to tree-ring width features [13,24,53,54]. The lack of consensus and uniform protocol for selecting and defining the most "climatically relevant" vessels [24], which may contain different climatic information depending on vessel ontogenesis and location within the ring [22], is likely to constitute a limitation. For instance, intra-annual density fluctuations (IADFs) in P. halepensis caused by resumed cambial activity within the same growth season represent a valuable ecological proxy of late summer to early autumn water availability, although their drivers are not fully understood [17]. Alternatively, factors such as tree size or ontogeny, stand structure and dynamics may mask the climatic signal contained by wood anatomical series. Nevertheless, wood anatomical data provide unique ecological information on how trees modulate their potential hydraulic conductivity and the carbon investment in cell-wall thickening in response to hydroclimate variability [14,24,54], particularly at intra-annual resolution [17].
The lower common signal often found in wood-anatomical series may appear a small award in comparison to the time-consuming job of creating well-replicated lumen-area or cell-wall thickness series. However, this aspect cannot be generalized as species such as P. sylvestris displayed in this study a discrete climatic signal both in the growth ring and in the wood-anatomical variables, therefore containing valuable climatic information related to water balance. The high variability of vessel cross-sectional area found in P. sylvestris indicates that wood-anatomical features constitute a valuable source of information to assess how plastic growth and phenology responses are in response to local hydroclimate variability.

There Is Great Variability between Individuals of the Same Population for Wood Anatomical Traits
High levels of adjustment or plasticity in the considered traits have been described in the past, and, indeed, they are often controlled for in comparisons among species, e.g., [22]. However, this study clearly shows that within-species variation is a significant, not minor, component of the overall variance in these traits. The finding that the amount of within-species trait variation is greater than interspecific variation indicated a large overlap between species due to the high individual variability.
Within-individual variability accounted for an overwhelming proportion of variance in the response of wood anatomical traits to water balance ( Figure 6). Such high individual variability among individuals of the same population in the case of wood anatomy represents an important issue per se in long-term studies on the anatomical characteristics of wood, as it would also indirectly affect variability in the species-specific response to climate. This partially agrees with results from studies evaluating responses of xylem traits to climate suggesting that most of the observed variation in those xylem traits was related to within-individual responses to year-to-year environmental change, rather than to average site conditions, and that these responses were mainly driven by adjustments of tracheid size and cell wall thickness [25]. This endogenous regulation may be due to intrinsic features such as tree height, branch length, or ontogeny affecting tapering [52,55]. Indeed, the year-to-year xylem hydraulic functionality is preserved by the combination (cell-to-wall area, length, etc.) of vessel traits developed in different growth rings, which, in turn, can be compensated for at the organ level (i.e., whole branch or trunk). With this compensation, the xylem characteristics might not respond to environmental changes in the same way to yearly environmental fluctuation as leaf attributes do, demonstrating significant within-individual resilience [56]. Such a source of within-individual variability that does not depend on the sampling site but exclusively on ontogenesis or tree features (e.g., height) is carefully evaluated in dendroecological studies and may partially be controlled by the use of standardized sampling procedures [57]. A caveat to our results is the low sample size, which may limit the robustness of our conclusions. Further research could test our ideas by analyzing more trees per population and more tree populations experiencing contrasting climate conditions.

Conclusions
Comparing the climate signal of growth and wood anatomy of three species encompassing wide climatic and ecological gradients, our findings commonly highlight the need of determining the impact of within-individual vs. environmental variation in order to fully understand how the environment drives species' xylem and growth adjustments. Research has consistently highlighted that small changes in xylem characteristics can have a major influence on hydraulic functionality and consistent dependence on climatic variability even with extremely low inter-tree correlations. However, if these xylem features differ more within than across tree populations, it begs the question of whether such mean population series are useless or, at the best case, biased. If a successfully constructed tree-ring chronology is an expression of the common environmental drivers of the growth of a population, even if the environmental predictors are often difficult to identify; however, this may not be true for the dendroanatomical time series. To this, the analytical criteria by which the average series of anatomical features are often analyzed should be not the same as the classic dendrochronological series, as the high variability between individuals of the same species casts doubt on some of the dominant assumptions of dendrochronological studies, including the replication and limiting factor principles. Although the amount of individual "noise" can be minimized by sampling many trees, it should be assessed if the extra effort would be enough to reduce the amount of intra-individual variability. Furthermore, the presence of a limiting factor common to the investigated tree species that would affect the amount of wood produced by a tree in any single year could have no causal implication on the characteristics of the wood anatomical traits of the species. This study emphasizes that the degree of variance around a species' mean characteristic must be accounted for since it has profound implications for understanding within-species variability and xylem adjustments to climate. We argue that future climate-focused analyses on growth and wood anatomy should explicitly consider individual variability (trees) rather than emphasizing mean population series (chronologies).

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/f13060956/s1, Figure S1: Wood cross-section of Pinus halepensis (A), Pinus sylvestris (B), and Pinus lumholtzii (C), Figure S2: Sources of variation within species for model relationships between seasonal water balance, tree age (years) and raw tree-ring and xylem anatomical traits. Relative variance decomposition for: tree-ring width (MRW), mean conduit area (MCA), and cell-wall thickness (CWT). In the legend, sigma 2 is the variance attributable to residuals; mean variation is the variance attributable to the random ID, i.e., random intercept variation; fixed is the explained variance by fixed variable, i.e., seasonal P-PET and tree age., Figure S3. Variance in tree-ring width (TRW) and wood anatomical variables (MCA, mean conduit area; CWT, cell-wall thickness) explained by variation in spring, summer, autumn, and winter seasonal water balance (PPET) between individual trees, and age for P. halepensis, P. sylvestris, and P. lumholtzii. Tree represents individual trees, and Age is the individual cambial age. Unexplained variance attributable to residuals is also reported, Table S1: Multilevel model results for Pinus halepensis.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.