Cell size and wall dimensions drive distinct variability of earlywood and latewood density in Northern Hemisphere conifers

Northern Hemisphere network of tree-ring data

This study includes ring width and density data retrieved from the International Tree Ring Data Bank (ITRDB; https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). The included data sets were selected by screening for sample sites where most of the parameters MXD, LWD, MND, EWD, LWW, EWW and TRW (Fig. 1) were measured.

Moreover, the TRD parameter (not archived) was computed as:

(Eqn 1)

(t, the year of each tree-ring property.) In the following, we often use the terms ‘latewood density’ for LWD and MXD, ‘earlywood density’ for MND and EWD, ‘ring density’ for TRD and ‘ring width’ for EWW, LWW and TRW. The Northern Hemisphere network of tree-ring data (NH-network) largely consists of boreal conifers, targeted by Fritz Schweingruber and colleagues during the last three decades of the 20^th Century for temperature reconstruction purposes (e.g. Briffa et al., 2002b). See Fig. 2 for the geographic distribution of chronologies and genera, and Supporting Information Notes S1 and S2 and Figs S1 and S2 for further details on the NH-network.

Analyses in this study focused on interannual to decadal timescales. Therefore, longer timescale variance in the tree-ring data, such as the biological age trends, or lower frequency variability driven by internal or external factors (including climate) was removed through standardization (Fritts, 1976). Cubic smoothing splines with a 50% frequency response cutoff at 35 yr (Cook & Peters, 1981) were fitted to each individual tree-ring series. The splines were either divided into (ratios for ring widths) or subtracted from (residuals for densities) the tree-ring series to produce indices (Cook & Peters, 1997). The resulting indices were averaged arithmetically into tree-ring chronologies, that is, time-series with annually dated tree-ring information. The common signal, that is, a synchronous behavior in each chronology's tree-ring series through time, was assessed with the Rbar statistic (Wigley et al., 1984; Fig. S3). The interrelationship among NH-network parameters was investigated with pair-wise Pearson correlations for each site and species over the full length of each chronology.

Climate data and climate correlations

The climate correlations in the NH-network were calculated using the Climatic Research Unit Time-Series version 3.22 (CRU TS3.22) data set, which consists of 0.5° gridded monthly temperature and precipitation data spanning 1901–2013 (updated from Harris et al., 2014). The data set was detrended analogously to the tree-ring data to extract high-frequency variations, to match the limited frequency range of the detrended tree-ring data and to reduce probabilities of spurious correlations resulting from trends in the data sets. We used the matlab function seascorr, developed by Meko et al. (2011), to quantify the climate signals in the tree-ring data. Pearson correlation coefficients were computed between each parameter chronology (for each site and species) and the temperature data from the nearest grid point. The relationships were assessed using monthly data in a 19-month window spanning the period from April of the year before ring formation to October of the year of growth. Correlations were computed for the period of maximum overlap between each site chronology and the observational climate record. Although a common period would provide increased consistency, the maximum overlap was chosen to attain as robust correlations as possible, as the common period in our NH-network only spans 1941–1968. As there was often strong covariation between temperature and precipitation (e.g. higher warm-season temperatures associated with lower precipitation; Trenberth & Shea, 2005), we computed partial correlations between precipitation and tree-ring data to determine if the confounding covariation of climate parameters is important, but also to determine if moisture alone is limiting for growth.

Moreover, the broad spatial and phylogenetic behavior of the temperature response across the NH-network was analyzed using a cluster analysis, here the K-Means algorithm (MacQueen, 1967). The analysis was based on the correlation matrix quantifying relationships between monthly temperatures and the latewood density records. The clustering was based on the current year MXD parameter, as it displayed the most significant temperature response. The distance from each observation to the centroid of the cluster (group) was quantified using the Euclidean distance. The optimal number of clusters in the data set was estimated by means of the gap statistics (Tibshirani et al., 2001) and the silhouette width (Kaufman & Rousseeuw, 2009).

Tracheid anatomical data

The anatomical data used in this study included features of tracheid radial diameter, cell wall thickness (Fig. 3) and number of tracheids per ring from eight data sites; see open circles in Fig. 2 for the geographic distribution. These data are further described in Notes S4.

The density of conifer wood can largely be determined by two different wood anatomical parameters: CWA, the area of wall within the tracheid, and TA, the area within the outer dimensions of the tracheid (Eqn 2; Fig. 3b; Vaganov et al., 2006). CWA and TA profiles were derived from the tracheid diameters and cell wall thickness measurements using model equations established in Fonti et al. (2013; Notes S4), at a similar radial step size to that commonly used in the NH-network (10 μm). From the anatomical profiles, we calculated density profiles and extracted parameters such as MXD, in a similar manner to how MXD parameter is extracted from X-ray profiles (Fig. 1). The density profiles were calculated as follows:

(Eqn 2)

where ρ denotes density and thus ρ_CWA/TA represents the bivariate relative intra-annual wood density profile (Fig. S4a). We conducted similar parameter intercorrelation and climate correlation analyses as were performed for the NH-network and compared average broadscale features from the X-ray density data and the anatomical density data.

To identify the importance of either wall material or cell size for interannual variations in wood density, we applied a series of sensitivity analyses, a statistical method that aims to determine the relative influence of each input parameter on a modeled output parameter by alternately holding one input parameter constant. In this study, the modeled output parameter was anatomical wood density, and the input parameters were CWA and TA. Thus, the analysis involved the calculation of density profiles with Eqn 2, but also calculating interannual density profiles alternately using an average profile of either CWA or TA as follows:

(Eqn 3)

(Eqn 4)

where ρ_CWA denotes the univariate intra-annual wood density, where only CWA contributes to the interannual variation in density, and ρ_TA is the inverse case. Again, we extracted the density parameters from the profiles (Fig. S4a–c). Subsequently, bivariate and univariate densities were correlated to each other, but also to corresponding widths and to current-year temperatures. To describe the results of the sensitivity analysis, we henceforth usually adopt the terminology r[x, y], which denotes the correlation between the parameters x and y. Typically, x and y refer to widths, bivariate or univariate anatomical densities, or monthly temperature parameters.

Intercorrelation analyses of the NH-network

The intercorrelations between the ring width and density parameters in the NH-network revealed a strongly systematic pattern (Fig. 4).

We found a consistently negative association between ring width and earlywood density, a modest association between ring width and ring density, with both negative and positive correlations, and finally a predominantly positive correlation between latewood density and ring width (Fig. 4). We henceforth refer to this switch in the sign of correlations as the ‘EW–LW correlation-sign change’. The intercorrelation analysis also showed that the width parameters EWW and LWW were usually highly positively correlated to each other, but the earlywood density and latewood density parameters were usually only modestly correlated (Fig. S5).

We also observed a prominent gradient in correlation coefficients of many parameter pairs, from western North America to Eurasia. The gradient was most pronounced for r[ring width, latewood density] and strongly related to the uneven spatial sampling of Larix sp. and Picea sp. (Fig. 4b). The gradient feature is henceforth referred to as the ‘bimodal biogeographic correlation’, where Larix sp. exhibited a more pronounced ‘EW–LW correlation-sign change’ than Picea sp. (Fig. 4b). The ‘bimodal biogeographic correlation’ appeared for most parameter pairs (Figs 4b, S5) but was less pronounced for r[ring width, earlywood density] or r[EWW, TRW]. The pattern can partly be explained by referring to the fact that Picea sp. usually displayed lower within-chronology series intercorrelation (Rbar) (Fig. S3b). That is, Picea sp. ring-width series contained a lower degree of common variance – when noisier data are correlated with data with stronger signal (exhibited within density), the chance of attaining high correlation coefficients is reduced. It is, however, unlikely that the ‘bimodal biogeographic correlation’ was exclusively related to the lower Rbar for Picea sp. widths, for at least two reasons: the pattern was not discernible in the correlation between width and density in the earlywood, and these (negative) correlations were highly significant, pointing to a predictable, non-noise-driven process in Picea sp. ring widths; the pattern was maintained also when correlating only density parameters, for example r[latewood density, earlywood density] (Fig. S5). Thus, there probably also exists a systematic anatomical difference that can more comprehensibly explain the ‘bimodal biogeographic correlation’.

Climate correlation and grouping of the NH-network

The climate–growth relationships across the sites and species indicated that temperature is the most important climatic factor in the NH-network (Fig. 5), although weaker but consistent correlations with precipitation were also emerging in late summer for all density parameters at most latitudes (Fig. S6).

The temperature response of the ring-width parameters (Fig. 5) for sites above c. 55° latitude was generally positive in midsummer, and for sites around c. 40–55° N, the broadscale signature was a lagged negative correlation with previous years mid- or late summer temperatures. Latewood width showed reduced lagged correlations in comparison to earlywood width. The temperature correlations of the density parameters were usually stronger and more consistent over latitude. In the following, we explore the wood density temperature interaction in more detail.

The cluster analysis divided the NH-network into five groups with characteristic climate response patterns. Fig. 6 shows the spatial repartition of the sites and species belonging to each cluster, with the average temperature response of the density parameters (see Fig. S7 for complete results).

Clusters 1, 2, 4 and 5 described distinct temperature signals, while cluster 3 grouped sites without a clear latewood density temperature correlation. Cluster 1 sites were predominantly located at high latitudes or altitudes in Scandinavia, European Russia and Canada, and was primarily composed of Picea sp., Larix sp. and Pinus sp. Its composition was not substantially different from the original fractionation of genera in the entire NH-network (Fig. S7a) and was therefore classified as geographically determined. Cluster 2 was composed of lower latitude sites, located mainly throughout the Alps, in which the genera Picea sp. and Abies sp. were overrepresented. Cluster 4 suggested a more or less genus-specific separation and is predominantly constituted by Picea sp. scattered across the North American continent. The majority of the sites in cluster 5 were found in Siberia, and contained predominantly Larix sp.

The temperature correlation for each cluster revealed striking complexities and similarities, among clusters but also among parameters (Fig. 6). The correlation structure was broadly characterized by a reliable positive correlation with current-year spring (March–May) temperature, with the temperature correlation being strongly dampened for all parameters at most sites (or even negative for the MND parameter) in midsummer (June–July), here termed the ‘midsummer decline’ (Figs 5, 6). Late summer was again characterized by positive correlations, particularly when gradually moving from the earlywood to the latewood density parameters. Furthermore, the earlywood densities showed systematic lagged correlations with previous-year late summer temperatures. The lagged correlations were increasingly dampened moving to ring density and were largely absent in the latewood density parameters. These patterns were observed, albeit to varying extents, in all of the five clusters (Fig. 6). However, the northernmost sites (Fig. 5), or cluster 5 (Fig. 6), did not indicate a ‘midsummer decline’ in the MXD parameter. Moreover, the southernmost earlywood density site chronologies (Fig. 5) displayed more prominent positive spring temperature signals starting already in February, accompanied by an extended period of negative precipitation signals from the previous October to the current June (Fig. S6).

Validation of the anatomical data set

Moving from the X-ray-derived density data to the cell anatomy-derived density data, the same broadscale features as described in Figs 4-6 were observed (Figs 7a, 8, upper panel).

To summarize, the ‘EW-LW correlation-sign change’ (Fig. 7a) was a prominent feature in the anatomical data; the ‘bimodal biogeographic correlation’ was also represented by the stronger ‘EW-LW correlation-sign change’ in Larix sp. (Fig. 7a); the temperature correlation was very similar to the NH-network results, even featuring a pronounced latewood density ‘midsummer decline’ for Picea sp. but not for Larix sp. (Fig. 8, upper panel). The only discrepancy between the cell anatomical data set and the tree-ring data set was the lack of a previous-year lagged temperature correlation for the earlywood anatomical densities (results not shown).

Sensitivity analyses of the influence of tracheid anatomy on wood density

We found that ring width measurements were consistently negatively correlated with the univariate tracheid area density (r[width, ρ_TA] < 0). Furthermore, we found positive coefficients for r[width, ρ_CWA] (Fig. 7a). Thus, both univariate versions of density lacked the ‘EW–LW correlation-sign change’. To attain an ‘EW–LW correlation-sign change’ in the observed direction (Fig. 4), we must then infer that TA dominated the variability in earlywood densities, whereas CWA was the primary control on the latewood density.

The bivariate density was positively correlated with the univariate CWA density (r[ρ_CWA/TA, ρ_CWA] > 0) throughout the entire tree ring and increased in strength moving from earlywood to latewood. By contrast, r[ρ_CWA/TA, ρ_TA] was positive only in the earlywood, and became negative in the latewood (Fig. 7b), and absolute values of correlation coefficients were rather unchanged, being higher than r[ρ_CWA/TA, ρ_CWA] only in the earlywood. This again suggested that the earlywood was dominated by fluctuations in TA, however it did not entirely exclude the role of CWA. The situation was more complex in the latewood. The fact that r[ρ_CWA/TA, ρ_TA] was negative in the latewood means that we counterintuitively observed a bivariate density increase when the cell size was increased. This finding suggests an even greater increase in CWA that compensates for fluctuations in TA. This is supported by the nearly constant lumen size in the latewood (Fig. S8). Fig. 3(c) conceptualizes these interannual anatomical differences between earlywood and latewood density.

The temperature correlation of the univariate tracheid area density (r[temperature, ρ_TA]) was almost identical to r[temperature, ρ_CWA/TA] in the earlywood parameters (cf the negative June correlation (Fig. 8) not observed for r[temperature, ρ_CWA]). By contrast, in the latewood, r[temperature, ρ_CWA] was almost identical to r[temperature, ρ_CWA/TA]. However, r[temperature, ρ_TA] was an almost perfect mirror image of r[temperature, ρ_CWA/TA].

Overall, these results indicated that earlywood density was predominantly controlled by fluctuations in tracheid size, whereas the fluctuations in latewood density were more difficult to disentangle because of the high degree of collinearity between TA and CWA (Fig. 7b).

Tracheid anatomy explains tree-ring data relationships

We have shown that TA is the main anatomical feature determining earlywood density. TA had a negative influence on density (Fig. 3c), but was positively correlated with the number of cells per ring (Fig. S9), which is the main determinant of ring width (Vaganov et al., 2006). Therefore, when correlating earlywood density and ring width it is very likely to show a negative relationship, as observed in this study (Fig. 4). The importance of TA in determining density was reduced in latewood compared with earlywood, but also, instead of having a negative influence on density as in the earlywood, TA had a positive association with density in the latewood (cf conceptual model of earlywood and latewood density in Fig. 3). As the contributions from both CWA and TA to latewood density appear to be positive, and both CWA and TA are positively correlated with the number of tracheids (Fig. S9), positive correlations between width and density are very likely to be obtained, as observed in this study (Fig. 4). Thus, the ‘EW–LW correlation-sign change’ is probably driven by the fact that the earlywood and latewood components maintain this fundamental shift in relative importance of cell anatomical features.

Inferring structure and functionality using tree-ring data

Further capitalizing on the combined use of anatomical data and tree-ring data, our results are also relevant because they provide evidence for a causal link between a widely used ecological and environmental proxy (density) (Chave et al., 2009; Stoffel et al., 2015; Wilson et al., 2016) and its fundamental structural–functional origin (anatomy).

The consistency of our discoveries (> 90% of sites had significant negative correlations between ring width and earlywood density (α < 0.05)) implies an almost universal mode of coupled environmental response of earlywood width and earlywood density. The hydraulic efficiency in conifers increases linearly with the number of conduits or cells (c. earlywood width), and increases exponentially with the lumen area of the cells (Sperry et al., 1994). Because earlywood density is mainly controlled by changes in tracheid size (TA and lumen area are interchangeable in the earlywood, i.e. almost perfectly correlated; Fig. S9), and changes in tracheid size have an exponential impact on conductivity, density and width should be interesting proxies for interannual variations in conductivity because of the currently superior replication and geographical coverage of tree-ring data. This prospect has previously also been suggested by Dalla-Salda et al. (2011), who conducted their study on 27 juvenile plantation Pseudotsuga menziezii trees.

Furthermore, there was also a near universal mode of coupled environmental response of latewood width and latewood density (> 80% of sites showed significant positive correlations between ring width and latewood density; α < 0.1). The mechanical function of the xylem is promoted by increased deposition of wall material (Chave et al., 2009) and an increased number of latewood cells. By maintaining the shift in the relative importance of tracheid features moving from the earlywood to the latewood, the same environmental conditions can consecutively promote the conductivity of the earlywood and the mechanical function of the latewood. The density of the latewood may thus compensate for any impairment in the mechanical function of the earlywood, complying with universal plant functional trade-offs (Baas et al., 2004).

Environmental responses in the functionality of either earlywood or latewood anatomy have, to our knowledge, only been established at the local scale (e.g. Bryukhanova & Fonti, 2013; Pritzkow et al., 2014; Castagneri et al., 2015). In this study, we indirectly demonstrate, for the first time, that conifers may increase their hydraulic efficiency in the earlywood, as well as their mechanical strength in the latewood, as a response to the same environmental conditions. Interestingly, these modes of responses are largely independent of phylogeny or prevailing climate regimes in the Northern Hemisphere.

Anatomical basis and sequential formation explain differences in temperature responses

Tracheid size, the major determinant of earlywood density, is a function of turgor pressure and hormonal control through cell wall relaxation occurring during the phase of cell enlargement (Cosgrove, 1993). For trees in cold-limited environments, conductivity, expressed by lumen area, is increased with temperature (Petit et al., 2011). The apparent paradox that temperature could drive tracheid enlargement can perhaps be explained by the fact that water availability generally is not a limiting factor in cold environments. This was indirectly observed in this study, where responses to precipitation were modest and responses to temperature were pronounced (Figs 5, S6). Because the temperature-sensitive cell size has a negative influence on earlywood density, the major response is negative in midsummer, coinciding with the period when most of the cells are formed and mature (e.g. Seo et al., 2008; Cuny et al., 2014). As TA does not exclusively determine earlywood density, there is probably also some influence of CWA. The positive temperature response preceding the negative temperature response, also found as a lagged correlation in the previous year, is presumably best explained as a wall deposition signal.

In the latewood, where both TA and CWA were found to positively influence density, we did not observe, and should not expect, negative responses to temperature. Studying the positive temperature correlations for the density parameters (Figs 5, 6), it is clear that their timing did not fully match the timing of ring/cell formation (e.g. Seo et al., 2008; Cuny et al., 2014). Earlywood densities showed consistently positive responses to temperature in the previous year. Latewood densities had pronounced responses to early spring temperatures. The cell wall contribution to wood density, regardless of whether in earlywood or latewood, must therefore be realized by capitalizing on resources obtained over an extended period, including the time before wood formation. This is supported by Kagawa et al. (2006), who reported that earlywood formation strongly relies on photosynthates stored from the previous year in Larix gmelinii, and by Kuptz et al. (2011), who concluded that older pools of photosynthates are also used for secondary stem growth in Picea abies. Vaganov et al. (2009) also found a carry-over effect of the δ¹³C (the ratio of the stable isotopes ¹³C and ¹²C) of latewood in the previous year on that of earlywood in the current year in P. abies and Pinus sp. Although the use of resources obtained over an extended period is compatible with latewood results from the anatomical data set, the previous-year lagged response in the earlywood was not detectable (Fig. 8). This absence of a lagged correlation in the anatomical data does not, however, rule out a lagged use of resources from the previous year to build earlywood tissue, but it dilutes the support for a lagged resource-use.

To summarize, the differentiated relative importance of anatomical features causes major differences between the temperature responses of earlywood and latewood density. Moreover, these differences may also be amplified by a possible dependence on previous year's reserves for earlywood density.

Similarities in wood density temperature responses are still elusive

In addition to distinct differences, we found conspicuous similarities in the temperature signals of earlywood and latewood densities. These were represented by synchronous positive early season responses followed by the ‘midsummer decline’ in the latewood, which was negative for the earlywood (Fig. 6). In the earlywood, we conclude that this temporal evolution of correlation coefficients is created by the opposing influences of the positive cell wall signal and the negative cell size signal. We do not rule out the possibility that such a mechanism is also present in the latewood. However, this possibility was not directly supported by our analyses, because ρ_CWA also displayed a ‘midsummer decline’ (Fig. 8). Recurring midsummer temperature-induced drought, during which high temperatures become a stress factor, is also questionable, because density was slightly negatively affected by increased precipitation just after midsummer (Fig. S6), and thus difficult to reconcile with midsummer temperature stress. The role of resource allocation is also inconclusive. A depletion of midsummer assimilates by concurrent cell division and cell enlargement may leave only resources accumulated in early and late summer for latewood formation. This could explain reduced correlations with midsummer temperatures. Cuny et al. (2015) found that the cell division and enlargement processes require astonishingly small amounts of carbon, compared with the wall deposition process, and a depletion when assimilation is at its peak (Bourdeau, 1959) is therefore unlikely.

Concluding remarks

We found a transition in the relative importance – from cell size towards cell wall dimensions – that drive distinctly different year-to-year variability of earlywood and latewood density, respectively, in Northern Hemisphere conifers. This transition explains the relationships between radial stem growth and wood density; earlywood densities were negatively correlated with ring widths, while latewood densities were positively correlated with ring widths. It also explains the complex temperature responses of wood density (notably varying intraseasonally in strength or even sign). These findings substantially improve our understanding of interannual variations in conifer growth and are vital for the interpretation of existing and future temperature reconstructions relying on the MXD parameter (e.g. Briffa et al., 2002b; Wilson et al., 2016).

Furthermore, because we could so easily integrate the anatomical data patterns with the phylogenetically and geographically diverse tree-ring data patterns, and identify a mechanism explaining why we observe these patterns, we conversely infer that the anatomy of earlywood and latewood tissues varies according to some simple universal rules. The year-to-year variability of earlywood tissue is dominated by variation in cell size, which supports hydraulic−functional responses to climate. By contrast, the year-to-year variability of latewood tissue is more dependent on cell wall variations, which permit mechanical–functional responses to climate. By combining the detailed anatomical data with the large-scale tree-ring data, we illustrate how Northern Hemisphere conifers organize investments in earlywood and latewood tissue, to continuously guarantee an adequately balanced and sustainable functioning of the xylem.