Heartwood Relationship with Stem Diameter in Pinus canariensis Plantations of Gran Canaria, Spain

: The development of stem heartwood and the factors that control it play an important role in tree physiology, thereby impacting demographic and ecological processes of woody species. We investigated the relationships of stem heartwood with site-and tree-level variables in Pinus canariensis plantations. A total of 30 plots were sampled in the island of Gran Canaria, Spain, over a large elevation range (995–1875 m) and on terrain with different slopes (4%–70%) and exposures. The 15 pines closest to each plot center were measured and cored to quantify growth rates and the size of heartwood, also known as “tea”. We used generalized linear mixed models (GLMMs) to account for both ﬁxed and random effects while evaluating the best predictors of heartwood presence. Stem diameter was the variable most correlated with heartwood radius, and allowing for a random slope and intercept of this relationship accounted for spatially related variability. Furthermore, the GLMM model became more effective when the relationship between stem diameter and heartwood was modeled using the presence/absence of the “tea” rather than its measured size. Other site-and tree-level variables either were not statistically signiﬁcant or improved the model relatively little. Because stem heartwood affects both wood quality and the amount of carbon that trees can store, our ﬁndings have implications for forest management and carbon-conscious policies.


Introduction
Xylogenesis, i.e., wood formation, is a complex process that involves cambium division, growth, elongation and thickening of cell walls, and in the case of tracheids and vessels, programmed cell death. The part of the xylem that contains living parenchyma cells with storage substances is called sapwood [1]. As the tree grows, the xylem cells die, giving rise to the formation of duramen or heartwood [2]. The sapwood performs important functions in various processes, such as storage of carbohydrates in the parenchyma cells, transport of water from the roots to the crown, and response to wounds, while the heartwood does not play a visible physiological role [3][4][5]. Detection of the proportion of sapwood in total stem biomass is necessary for large-scale assessments of respiration and transpiration of woody plants [6].
The heartwood determines tree resistance to physical and biological stressors, ranging from wind to insects, as well as multiple features of the wood (color, density, chemistry, etc.) [7,8]. Trees regulate xylem formation and sapwood to heartwood conversion to maintain efficient water transport as the tree height increases [9,10]. While heartwood formation is an actively regulated stage in tree growth [11], the role of environmental factors is still being debated, especially in connection with internal mechanisms, such as age and loss of vitality of the parenchyma cells [12], toxic effects of polyphenols [13], high levels The climate of the Canary Islands is characterized by cool trade winds and cold surface waters [42], which create moderately temperate conditions during the warm season on land, despite the low latitude of the islands (~28° N). In winter, the Azores anticyclone changes location, and the frequency of trade wind days decreases to about half of the total [43]. Large-scale climatic processes are translated into precipitation and temperature regimes on land that are heavily influenced by topography, especially elevation and exposure [43]. Overall, precipitation increases and air temperature decreases at higher elevations and on slopes directly exposed to north-easterly trade winds [44]. Monthly average total precipitation and mean air temperature for a representative station ( Figure S3) demonstrate limited seasonal differences in temperature, with August being the warmest month (mean of 21.1 °C) and January-February the coolest ones (mean of 12.4 °C).
Wet conditions are confined to the winter season (December-February), with monthly total precipitation being highest in December (38 mm) and lowest in June-August (6 mm per month). Annual averages are 16.5 °C for mean air temperature and 245 mm for total precipitation. According to recently interpolated graphics produced by the Consejo Insular de Aguas de Gran Canaria, maximum average annual rainfall in the island of Gran Canaria reaches 950 mm in areas where trade winds interact with topography to generate an almost continuous presence of clouds ("mar de nubes"; The climate of the Canary Islands is characterized by cool trade winds and cold surface waters [42], which create moderately temperate conditions during the warm season on land, despite the low latitude of the islands (~28 • N). In winter, the Azores anticyclone changes location, and the frequency of trade wind days decreases to about half of the total [43]. Large-scale climatic processes are translated into precipitation and temperature regimes on land that are heavily influenced by topography, especially elevation and exposure [43]. Overall, precipitation increases and air temperature decreases at higher elevations and on slopes directly exposed to north-easterly trade winds [44]. Monthly average total precipitation and mean air temperature for a representative station ( Figure S3) demonstrate limited seasonal differences in temperature, with August being the warmest month (mean of 21.1 • C) and January-February the coolest ones (mean of 12.4 • C).
Wet conditions are confined to the winter season (December-February), with monthly total precipitation being highest in December (38 mm) and lowest in June-August (6 mm per month). Annual averages are 16.5 • C for mean air temperature and 245 mm for total precipitation. According to recently interpolated graphics produced by the Consejo Insular de Aguas de Gran Canaria, maximum average annual rainfall in the island of Gran Canaria reaches 950 mm in areas where trade winds interact with topography to generate an almost continuous presence of clouds ("mar de nubes"; http://www.aguasgrancanaria.com/ cartografia/tematica/precipitaciones.php, accessed on 2 August 2023).

Field and Laboratory Measurements
Field measurements and samples were collected in 2001 ( Figure S4). Once a tree was chosen as the center of a plot, the 14 trees closest to it were identified. The height of these 15 trees was measured with a double-needle hypsometer (precision of 0.25 m), their stem was measured with a diameter tape (precision of 0.5 cm), and bark thickness was measured with an appropriate caliper (precision of 1 mm). A wood increment core was extracted at breast height (~1.3 m from the ground) on the north-facing side of each tree using a 40 cm Pressler borer. Existing information indicates that heartwood formation in Pinus canariensis is fairly uniform around the stem (Figures S1 and S2); hence, coring maximized the number of sampled trees without having to consider the core azimuth. Terrain slope and exposure were measured at the plot center, whose geographical coordinates were obtained using a GPS device. After measuring and coring the first 15 trees, the next 15 trees closest to the plot center were identified. The linear distance (m) from the plot center to tree 15 (D 15 ) and tree 30 (D 30 ), i.e., the furthest of the first 15 and the furthest of the last 15 trees, were measured and used to calculate a stand density index as follows [45]: 10, 000 2π where 10,000 is a conversion factor to express density as number of trees per hectare.
In the laboratory, the wood increment cores were analyzed to obtain the number of rings and the extent of the heartwood. When the core did not include the stem pith, its location was approximated using predetermined concentric circles [46], and the number of rings between the pith and the beginning of the core was then estimated using the average size of the innermost rings ( Figure S5). Visible (and estimated) rings were used to calculate cambial stem age. Cumulated ring widths were used to reconstruct stem radii at five-year intervals (R5, R10, R15, R20, ...) with a precision of 0.5 cm. The heartwood radius and the stem radius inside the bark were also measured on each core. We followed [40]'s statement that "The sapwood-heartwood boundary is sharply defined by heartwood colour in Pinus canariensis. Noncoloured heartwood is absent in this species, and thus, all heartwood traits refer to coloured heartwood."

Data Processing and Statistical Model
Because climate data were not available for the actual plot locations, we used topographical features to estimate indices of site moisture and insolation. A 30 m digital elevation model (DEM) was obtained using the OpenTopography DEM downloader plugin of QGIS [47] and used to produce a map of topographic wetness index (TWI) and a map of the Sky View factor. TWI is an index of where water will accumulate as a function of terrain slope and of the upstream contributing area [48]. Topographic variables (elevation, slope, and exposure, also in relation to the surrounding terrain) were further integrated into a Sky View factor, which varies from 1 for completely unobstructed land surfaces (such as horizontal surfaces or peaks and ridges) to 0 for completely obstructed land surfaces [49].
Frequency histograms, box plots, dot charts, and pairwise sample correlations (both Pearson's and Spearman's) were used for exploratory data analysis and to investigate relationships among plot and tree variables [50]. Whenever possible, relationships were examined by plot as well as for the entire dataset. As geographical coordinates were only available for the plot center, spatial autocorrelation of tree variables was geostatistically estimated [51,52] for the whole dataset by assigning the same location to all trees within a plot. The maximum pairwise Euclidean distance between the 30 plots was 17.1 km; 700 m distance intervals were used for variogram estimation and modeling because those intervals included more than 10 paired distances up to~13 km.
The statistical relationship between the "tea" radius and other measured variables was estimated using a generalized linear mixed model (GLMM; [53][54][55]). To account for the presence of nested data, i.e., the tree measurements that were taken within each plot, we used a random effect for the sampled plot and a fixed effect for the tree stem diameter, with an interaction with the plot that allowed both a random intercept and a random slope so that the "tea"-diameter relationship could change from place to place. The potential influence of additional tree variables was tested by adding them as fixed terms, one at a time, to this initial GLMM and then comparing the model formulations using the Akaike information criterion (AIC; [56]). Spatial autocorrelation was modeled by an exponential covariance structure, and each GLMM was fit using restricted maximum likelihood estimation (REML; [57]).
Pine heartwood was often absent; hence, the "tea" radius (Rtea) was characterized by a large frequency of zero values. Because no data transformation can ameliorate this issue, we also used a binary response (1 for Rtea >0, and 0 otherwise) in order to apply mixed-effects logistic regression [58]. Introducing a random effect for plot location into the standard logistic regression model [59] generates a mixed-effects logistic regression. The advantage of a random intercept model is that the probability of a tree having "tea" is correlated to other trees on the same plot. It is worth noting that there is no random slope in this model, and that the predictor (diameter) is a continuous, but truncated, variable because it cannot be <0. Because of this, prior to model fitting, the diameter was "centered", i.e., transformed into deviations from the mean [57]; otherwise, the intercept would represent the probability that a tree stem of zero diameter had "tea", even though there cannot be a stem of zero diameter. By centering the stem diameter, the intercept acquired the more meaningful interpretation of the probability that a tree of average diameter had "tea". All numerical analyses were performed in the R computing environment [60]. The mixed-effects models were fit using packages nlme [61], glmmML [62], and MASS [63].

Sampled Trees and Exploratory Data Analysis
Sampled plots (Table S1, Figure 2) were distributed over an 880 m elevation interval, from 995 to 1875 m above mean sea level. Both thinned (13) and unthinned (17) plots were characterized by terrain that varied from almost flat (4%-9% slope) to steep (60%-70% slope), with different exposures. The map of TWI ( Figure 2a) was plotted using a pseudocolor scheme, but it should be noted that because of the limited rainfall, there are no rivers in the island of Gran Canaria, only continuous streams in three of the island's ravines.
Stand density showed a very large outlier for one of the plots (Table S1), which was considered to be correct, but its influence on correlations was reduced by log-transforming (natural logs) these indices. Because the original plot numbers were found to be inversely correlated with elevation and with the "easting" coordinate (X_UTM , Table S1) but directly correlated with the "northing" coordinate (Y_UTM , Table S1), such geographical dependencies were eliminated by reassigning plot numbers at random in subsequent analyses.
Tree-level variables showed reasonably bell-shaped frequency distributions with the exception of heartwood radius (Rtea; Figures 3 and S6), which was highly skewed because of a large number of zeros, i.e., stems where heartwood could not be identified. Stem age for "tea" onset varied, with the earliest presence estimated at cambial ages of 19 years (Figure 4).
Among tree-level variables, stem diameter was the best predictor of heartwood radius (Figures 3 and S6). Additional significant correlations between Rtea and other tree-level variables were most likely influenced by their collinearity with stem diameter (Figures 3 and S6). Some spatial autocorrelation in heartwood radius was suggested by the variogram model ( Figure S7).

Mixed-Effect Models
The generalized linear mixed model (GLMM), which was estimated with the nlme R package using restricted maximum likelihood estimation (REML), included both a random intercept and a random slope to account for plot-related differences in the relationship between tree stem diameter and "tea". It also included an exponential covariance structure to account for spatial autocorrelation in the "tea". Based on the model output (Table  S2a, Figure 5), the fixed effects, i.e., the intercept and slope for the relationship between the stem diameter and the "tea" radius, were highly significant (p-value « 0). A high negative correlation (−0.969) existed between the random intercepts and slopes. The Akaike information criterion (AIC) of this model (884 ; Table S2a) was slightly higher than the AIC for the GLMM without spatially autocorrelated errors (880; detailed results for that model are not included). Stand density showed a very large outlier for one of the plots (Table S1), which was considered to be correct, but its influence on correlations was reduced by log-transforming (natural logs) these indices. Because the original plot numbers were found to be inversely correlated with elevation and with the "easting" coordinate (X_UTM , Table S1) but directly correlated with the "northing" coordinate (Y_UTM , Table  S1), such geographical dependencies were eliminated by reassigning plot numbers at random in subsequent analyses.
Tree-level variables showed reasonably bell-shaped frequency distributions with the Among tree-level variables, stem diameter was the best predictor of heartwood radius (Figures 3 and S6). Additional significant correlations between Rtea and other tree-level variables were most likely influenced by their collinearity with stem diameter (Figures 3 and S6). Some spatial autocorrelation in heartwood radius was suggested by the variogram model ( Figure S7).  Among tree-level variables, stem diameter was the best predictor of heartwood radius (Figures 3 and S6). Additional significant correlations between Rtea and other tree-level variables were most likely influenced by their collinearity with stem diameter (Figures 3 and S6). Some spatial autocorrelation in heartwood radius was suggested by the variogram model ( Figure S7).   After centering the stem diameter, a mixed-effects logistic regression model was fit to the data without any error covariance structure. The R package glmmML, which estimated the model parameters by maximum likelihood (Table S2b), showed that this model had a much lower AIC (571), thereby suggesting that the relationship between stem diameter and hardwood radius was more effectively modeled using the presence/absence of the "tea" rather than its measured size. Based on the model output (Table S2b), the probability (p) that the stem of pine j on plot i had "tea" is given by logit(p ij ) = −0.038 + 0.089 diam ij + a i with a i ∼ N(0, 0.976 2 ).

Mixed-Effect Models
The generalized linear mixed model (GLMM), which was estimated with the nlme R package using restricted maximum likelihood estimation (REML), included both a random intercept and a random slope to account for plot-related differences in the relationship between tree stem diameter and "tea". It also included an exponential covariance structure to account for spatial autocorrelation in the "tea". Based on the model output (Table S2a, Figure 5), the fixed effects, i.e., the intercept and slope for the relationship between the stem diameter and the "tea" radius, were highly significant (p-value « 0). A high negative correlation (−0.969) existed between the random intercepts and slopes. The Akaike information criterion (AIC) of this model (884 ; Table S2a) was slightly higher than the AIC for the GLMM without spatially autocorrelated errors (880; detailed results for that model are not included). The thick black line is the stem diameter-"tea" relationship on a typical plot as linear regression parameters in the GLMM were estimated with respect to an individual plot (thin black lines) due to the random intercept and slope (see text for details).
After centering the stem diameter, a mixed-effects logistic regression model was fit to the data without any error covariance structure. The R package glmmML, which estimated the model parameters by maximum likelihood (Table S2b), showed that this model had a much lower AIC (571), thereby suggesting that the relationship between stem diameter and hardwood radius was more effectively modeled using the presence/absence of the "tea" rather than its measured size. Based on the model output (Table S2b), the probability (p) that the stem of pine j on plot i had "tea" is given by

Random Intercept and Slope Model with Autocorrelated Errors
Stem diameter (cm) Heartwood radius (cm) Figure 5. Relationship between stem size and heartwood radius in Pinus canariensis plantations. The thick black line is the stem diameter-"tea" relationship on a typical plot as linear regression parameters in the GLMM were estimated with respect to an individual plot (thin black lines) due to the random intercept and slope (see text for details).
Because the random intercept a i was assumed to be normally distributed with mean 0 and variance 0.976 2 , 95% of a i values would be between −1.96 × 0.976 and 1.96 × 0.976, which resulted in the 95% confidence interval of (−1.913, 1.913). GLMM-predicted probabilities of "tea" presence with respect to (centered) stem diameter for pines at all plots are shown graphically in Figure 6 using the logistic curves for a "typical" plot and for the upper and lower limits of this 95% confidence interval. Based on the model results, going to a typical plot and sampling a pine of average stem diameter (the 0 value on the x-axis of Figure 6), the probability of finding "tea" would be~0.49, which is the y-axis value obtained from the population curve (thick red line in Figure 6). Therefore, there is essentially a 50-50 chance of finding "tea" on an average pine sampled on a typical plot. In addition, for 95% of plots, this probability can be anything between 0.12 and 0.87 (i.e., the y-axis values obtained from the thin red lines in Figure 6), pointing to a large degree of inter-site variation.

Further Testing of Fixed Effects
Further statistical testing focused on other potential predictors of stem heartwood radius. Briefly, no other fixed effect appeared to contribute significantly to the predictive power of the mixed-effects logistic regression model shown earlier, and we report here two examples of those additional analyses. First, we investigated the impact of thinning by adding a categorical dummy variable (either 0 or 1) as a second fixed effect, together with its interaction with the stem diameter. Based on the p-values, the treatment and its interaction with stem diameter were not significant (p-value > 0.4), either in the GLMM with spatial autocorrelation (Table S2c) or in the logistic GLMM (not shown). 0 and variance 0.976 2 , 95% of ai values would be between −1.96 × 0.976 and 1.96 × 0.976, which resulted in the 95% confidence interval of (−1.913, 1.913). GLMM-predicted probabilities of "tea" presence with respect to (centered) stem diameter for pines at all plots are shown graphically in Figure 6 using the logistic curves for a "typical" plot and for the upper and lower limits of this 95% confidence interval. Based on the model results, going to a typical plot and sampling a pine of average stem diameter (the 0 value on the x-axis of Figure 6), the probability of finding "tea" would be ~0.49, which is the y-axis value obtained from the population curve (thick red line in Figure 6). Therefore, there is essentially a 50-50 chance of finding "tea" on an average pine sampled on a typical plot. In addition, for 95% of plots, this probability can be anything between 0.12 and 0.87 (i.e., the y-axis values obtained from the thin red lines in Figure 6), pointing to a large degree of inter-site variation. Figure 6. Logistic curves showing the estimated probabilities of "tea" presence for a range of stem diameter anomalies (= centered values). The thick red line in the middle corresponds to a typical plot, i.e., the predicted probability for the "population of locations", or ai = 0. The other two red lines are obtained by adding and subtracting 1.913 from the random intercept to the predictor function. Hence, 95% of the plots have logistic curves between these two extremes, and the space between these two curves represents the variation of predicted values per plot.

Further Testing of Fixed Effects
Further statistical testing focused on other potential predictors of stem heartwood radius. Briefly, no other fixed effect appeared to contribute significantly to the predictive power of the mixed-effects logistic regression model shown earlier, and we report here two examples of those additional analyses. First, we investigated the impact of thinning by adding a categorical dummy variable (either 0 or 1) as a second fixed effect, together with its interaction with the stem diameter. Based on the p-values, the treatment and its interaction with stem diameter were not significant (p-value > 0.4), either in the GLMM with spatial autocorrelation (Table S2c) or in the logistic GLMM (not shown).
The second example reported here deals with the tree initial radial growth rate, represented by the cumulated size of the first five annual rings (R5 in Figure 3 and Figure  S6). As shown by smoothed scatterplots (red lines in Figure S6), this variable was positively correlated with stem diameter (linearly), height (nonlinearly), and the size of the heartwood (slightly). The (centered) R5 had a negative and barely significant (p-value = 0.0397) coefficient in the logistic GLMM (not shown), but its interaction with the stem diameter was not significant (p-value = 0.665). The mixed-effects logistic regression model was, however, slightly improved by adding this second fixed term, as shown by a minimally lower AIC (569.8). Because the interaction term was not significant, the model was refitted without it, and the AIC decreased a bit more (568; Table S2d). Because the Stem diameter anomaly (cm) Probability of 'tea' presence 10 Figure 6. Logistic curves showing the estimated probabilities of "tea" presence for a range of stem diameter anomalies (= centered values). The thick red line in the middle corresponds to a typical plot, i.e., the predicted probability for the "population of locations", or a i = 0. The other two red lines are obtained by adding and subtracting 1.913 from the random intercept to the predictor function. Hence, 95% of the plots have logistic curves between these two extremes, and the space between these two curves represents the variation of predicted values per plot.
The second example reported here deals with the tree initial radial growth rate, represented by the cumulated size of the first five annual rings (R5 in Figures 3 and S6). As shown by smoothed scatterplots (red lines in Figure S6), this variable was positively correlated with stem diameter (linearly), height (nonlinearly), and the size of the heartwood (slightly). The (centered) R5 had a negative and barely significant (p-value = 0.0397) coefficient in the logistic GLMM (not shown), but its interaction with the stem diameter was not significant (p-value = 0.665). The mixed-effects logistic regression model was, however, slightly improved by adding this second fixed term, as shown by a minimally lower AIC (569.8). Because the interaction term was not significant, the model was refitted without it, and the AIC decreased a bit more (568; Table S2d). Because the output from the glmmML package did not include the estimated correlations between fixed effects, we used the glmmPQL function from the MASS package to obtain that information ( Table S2e). The interaction term between the two fixed effects was confirmed to be nonsignificant (results not shown); hence, the model was refit without the interaction, and a high negative correlation (−0.607) was found between the two fixed effects, i.e. stem diameter and five-year cumulative growth. Based on the model output (Table S2e), the probability (p) that the stem of pine j on plot i had "tea" is given by with a i ∼ N(0, 1.034 2 ). A graphical representation of GLMM-predicted probabilities that pine j on plot i had "tea" would require a 3D plot to include the "stem diameter anomaly" axis ( Figure 6) as well as the new "R5 anomaly" axis. In addition, one can derive from the logit equation that there is still essentially a 50-50 chance of finding "tea" on an average pine sampled on a typical plot, with 95% of plots having probability between 0.11 and 0.88. Overall, introducing the R5 variable as a fixed effect improved the model relatively little, even though measuring R5 requires considerable extra effort both in the field (to collect increment cores) and in the laboratory (to analyze growth rings).

Discussion
In our study, the presence of heartwood or "tea" in Pinus canariensis stands that were planted in the island of Gran Canaria was most efficiently predicted by stem diameter measured at breast height. Based on the GLMM with exponential error covariance structure, some spatial autocorrelation occurred up to a linear distance of about 500 m, but there was no statistical advantage for including this model component. In other words, a random slope and intercept of the DBH-"tea" relationship already accounted for most of the spatially related variability in heartwood radius. Furthermore, the GLMM model became more effective when the relationship between stem diameter and heartwood was modeled using the presence/absence of the "tea" rather than its measured size.
Stem diameter was also the best estimator of heartwood radius for seven Chinese species in temperate climates [64]. The connection between heartwood formation and tree diameter was found to be stronger than the one with age in other species too, including Tectona grandis [65], Acacia melanoxylon R. Br. [66], and Eucalyptus globulus Labill. [67]. However, stem age was found by other authors to predict the amount of heartwood in Pinus canariensis [40,68] and in other species [16,27,28]. As the maximum estimated tree age at our study plots was 46 years (Figure 4), which is quite low due to the artificial regeneration of these stands, it is possible that the relationship between stem age and "tea" was masked by the connection between stem age and size, which is usually stronger at younger ages.
It should be noted that the onset of "tea" formation was found to be at about 30 years of age in naturally regenerated stands of Canary Island pine [40,68], whereas in the sampled plantations, it was about 2/3 of that at~19 years ( Figure 4). The earlier appearance of "tea" in these artificial stands could be linked to lower competition for resources because of understory removal at the time of planting and fairly regular spacing among pines. Pinus canariensis is a shade-intolerant species; hence, its natural regeneration occurs mostly after disturbances, such as wildfires or tree falls, which open forest gaps where regeneration is abundant, and growth is then slowed down by competitive interactions.
Initial stem growth, when cambial age is low, has been found to have a positive influence on the development of heartwood in natural stands of Pinus canariensis [35] and Pinus radiata [41,69]. In our study, cumulated wood increments during the first five years of cambial age (the "R5" variable), while statistically significant, added comparatively little to the performance of the GLMM. Furthermore, the sign of the fixed-effect coefficient for R5 was negative, which is the opposite of what was observed in previous studies. While we cannot exclude a numerical artifact caused by computational instabilities that are inevitable under collinearity of predictors [57,70], it is conceivable that larger tree rings would imply larger cell lumens, which in turn could favor sap flow, prevent lumen occlusion, and thereby delay the formation of heartwood. Future research should be aimed at investigating chemical and anatomical characteristics of Pinus canariensis heartwood to test such an intriguing hypothesis.
The lack of a significant difference in the relationship between stem diameter and heartwood formation between plots that were thinned and those that were not may be linked to several factors. It is interesting to note that a similar lack of reaction to pruning was reported for Pinus sylvestris [71], which experienced a negative defoliation impact in terms of stem growth rates but no significant connection to heartwood size. Because heartwood formation was not favored by thinning and firewood obtained from intermediate silvicultural treatments of Canary Island pine has no economic value, our study suggests that reforestation should be carried out at final densities or close to them so that thinning is not necessary.
While we estimated climatic conditions at the sampled plots using topographically derived variables, and found no significant correlations, other studies have indicated that climate can influence heartwood formation in both conifers [72,73] and angiosperms [74]. For teak (Tectona grandis) grown in India, wet sites yielded larger heartwood than dry sites [75], whereas the opposite case was found for teak planted in Ghana [76]. In the boreal forest, black spruce (Picea mariana (Mill.) BSP) heartwood was greater under wet conditions than under dry ones [77], whereas in natural stands of Pinus canariensis, humid climates corresponded to lower amounts of heartwood than dry climates [40]. Future research should further test the impact of climatic variability on heartwood formation because there is a possibility that our sampled plots were located in areas with higher moisture compared to the entire range of climatic conditions that are naturally experienced by Pinus canariensis in its native habitat.
The heartwood ("tea") of Pinus canariensis has been used for economic and cultural activities as a highly prized material because of its color, density, and durability. There is archeological evidence that even pre-Hispanic societies used "tea" for constructing barns and for funeral practices [78]. Intensive exploitation of the pine forest in the Canary Islands started in the 15 th century with Spanish colonization as quality wood was required to build furniture and ships and very large pines were present in the islands, most likely exceeding the height and diameter of any currently remaining pine [35].
The percentage of heartwood in mature trees is an essential characteristic for both wood quality and carbon sequestration [79]. Carbon absorbed from the atmosphere is deposited for a long time in the living tissue of woody plants, but there are differences in carbon amounts between heartwood and sapwood, which are usually inconsistent between species and within taxonomic groups [80]. As carbon stocks in forest ecosystems play a fundamental role in global management of natural resources in a changing world [81,82], future research should then be aimed at clarifying the factors that control heartwood formation in this species as well as others with the goal of designing the most carbonconscious forest policies.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/f14091719/s1: Figure S1: Cut stem showing the cross-sectional color difference between heartwood and sapwood in Pinus canariensis; Figure S2: Cut stem showing the cross-sectional decay resistance of heartwood compared to sapwood in Pinus canariensis; Figure S3: Walter-Lieth climate diagram for a representative station [83]; Figure S4: The first author coring a Pinus canariensis on one of the study plots; Figure S5: The Applequist estimator for average ring width of 2.5 mm; Figure S6: Paired-sample linear correlations between plot and tree variables; Figure S7: Omnidirectional sample variogram and fitted exponential model for heartwood radius; Table S1: Summary information for sampled plots; Table S2. Output of mixed-effects models estimated using R packages.