A common garden experiment was set up in 1997 using a strictly replicated experimental design at different sites in eastern and central Spain. The trial includes 52 P. halepensis populations from continental Spain, Balearic Islands (Spain), France, Italy, Greece and Tunisia, thus covering most of the species’ range (see Climent et al., 2008 for details). We chose a subset of 19 populations (Figure 1 and Supplementary Table S1) representing the main geographic and environmental gradients of origin, and with a good balance among them. Previous studies showed major significant differences among P. halepensis populations, with eastern populations displaying faster growth and lower investment in early reproduction (including low serotiny) and populations toward the southwest with the opposite trends (Climent et al., 2008; Santos-del-Blanco et al., 2013). One-year-seedlings were finally planted in a row-column design in four randomized complete blocks with four plants per population and block (16 trees per population). Height, diameter and cone number were periodically measured, enabling absolute and relative growth assessments and reproductive allocation estimations.

We selected two sites with contrasting environments allowing to separate population differentiation from phenotypic plasticity of bark thickness. Both sites differ in a range of environmental variables (Table 1), with drier site showing a 20% lower precipitation, a 100% higher drought intensity index (K) and more frequent and intense wind. Both sites contrast in soil structure, depth and water retention, which leads to a strong divergence in water availability. By contrast, soil nutrient availability between sites is very similar (Gil et al., 1996). Therefore, we refer to these differences as “drier” and “moister.” Both potential forest productivity (Sánchez-Palomares and Sánchez-Serrano, 2000) and survival were higher at the moister site (1.8-fold in productivity and 9% more survival compared to the drier site; Martín-Sanz et al., 2016; Figure 2 and Table 2). Moreover, site also affected significantly reproduction, degree of serotiny (long-lasting closed cones in the crown) and tree form.

Sampling was carried out in both experimental sites when trees were 18 years old (271 and 241 trees were sampled at the “moister” and “drier” sites, respectively). Stem diameter and bark thickness were measured basally (10 cm above the ground; hereafter referred as D10 and BT10, respectively) and at breast height (hereafter D130 and BT130, respectively) using a tree caliper and a standard bark gauge (Supplementary Figure S1). Since we were obliged to preserve the common garden experiment we abstained from more precise, but invasive or destructive bark measuring methods like extraction of bark portions (Jackson et al., 1999; Paine et al., 2010; Graves et al., 2014) or cutting stem discs to remove the bark (Rosell et al., 2014, 2015; Rosell, 2016). Therefore, by bark thickness (basal: BT10 or at breast height: BT130) we refer to total (inner and outer) bark. We performed a preliminary sampling of a few trees at each test site measuring bark thickness at three points surrounding the bole. We found that measuring bark thickness at a single point was accurate enough given the uniformity of P. halepensis bark at this age. Thus, for each height, we measured bark at a single point of the bole, always in the south orientation.

While simple and straightforward, we think that the linear relationship between bark thickness and diameter does not represent accurately enough the real resource allocation to bark. Our point is that allometry of bark allocation should be rather assessed through variables that are as-close-as-possible to the biomass of bark and sapwood (the functional part of the trunk). Therefore, focusing purely on the allometry of bark allocation, we looked at the allometric relationship between bark and sapwood volumes, basing the analysis on the relationship of percentage of bark volume and total volume from tree base to breast height (Martín-Sanz, 2018). Our sampled trees did not contain any heartwood (see Appendix A and Supplementary Figure A.1), so total volume is bark and sapwood (functional xylem) volumes.

Firstly, we calculated the total volume for the cone trunk from the tree base to breast height with the equation for a circular truncated cone (in dm³):

where h is the height difference between the base (10 cm above ground) and breast height (130 cm), Rb is the basal radius and Rbh is the radius at breast height.

Sapwood volume (dm³) was computed with the following equation:

where h is the height difference between the base (10 cm above ground) and breast height (130 cm), SRb is the sapwood basal radius and SRbh is the sapwood radius at breast height.

Finally, bark volume (dm³) and percentage of bark volume were estimated by Equations (3) and (4), respectively:

General linear mixed models (LMM) were fit with PROC MIXED procedure in SAS software (SAS Institute Inc., 2015). The incorporation of variables into the models was checked using the Akaike’s Information Criterion (AIC), selecting those models with lower AIC value. This criterion favors model fit and simplicity, based on the principle of parsimony (fewer parameters in the model). Residuals of each model were examined for normality using diagnostic plots and Kolmogorov-Smirnov normality tests.

We compared means of percentage of cone trunk from tree base to breast height taken up by bark volume (Equation 4) by site and populations using a LMM with environmental effect represented by trial sites and genetic effect by populations. In addition, we did the same model adding total volume of the cone trunk as a covariate in order to look for tendencies in percentage of bark volume. In both models, block within site was used as random factor.

Additionally, we wanted to compare our approach with the most commonly applied allometric method, thus we also used model II regression, i.e., standardized major axis (SMA) to examine the allometry of bark volume vs. the total volume of the cone trunk from the tree base to breast height (Henry and Aarssen, 1999; Niklas, 2006; Warton et al., 2006).

The classical allometric model is V_B = aV_T^b and is usually fit as logV_B = log₁₀a + b log₁₀V_T (Huxley, 1932; Niklas, 1994; Ter-Mikaelian and Korzukhin, 1997). Parameter a is the elevation or allometric coefficient and parameter b is the regression slope or allometric exponent. An exponent significantly different from 1 indicates an allometric relationship between the variables studied (increasing or decreasing with size, non-constant). The allometric exponents of each population at each site were compared with the isometric coefficient (b = 1) and to one another by multiple post hoc comparisons. These analyses were carried out with the smatr package (Warton et al., 2012) implemented in R software v3.3.2. (R Core Team, 2016). As the SMA analysis indicated that slopes differed among populations and sites, a general SMA test of elevation (allometric parameter) differences was not necessary.

We analyzed the effects of population (genetic effect), experimental site (environment) and their interaction on absolute bark thickness at breast height and at tree base with LMM, including block within site as random factor. We compared the population per site variation of absolute bark thickness with critical thickness values –the thickness above which the risk of cambium damage decreases considerably, therefore indicating a threshold for tree survival– both at breast height and at tree base. While our bark thickness data did not differentiate between inner and outer bark neither we compared further heat insulating properties as moisture content or density, we assumed that bark anatomy should not differ significantly among provenances of a single species, as found in other species (see Wesolowski et al., 2014). Due to the lack of experimental data of critical bark thickness for P. halepensis, we used data published for the genetically closest available tree species: a high thickness of 20 mm, postulated to allow survival to moderate fires (for Pinus pinea; Madrigal et al., 2019) and a low thickness of 10 mm that would allow survival only to low-intensity fires (van Mantgem and Schwartz, 2003; VanderWeide and Hartnett, 2011).

Looking for ecotypic trends in bark thickness, we tested for correlations among the mean values for bark thickness of our raw data (at breast and basal heights) and continuous environmental variables, as well as with fire frequency records from the 13 Iberian populations. The nineteen original bioclimatic variables for populations’ origin were obtained from Worldclim v.1.4 (Hijmans et al., 2005). We considered nine of these variables: annual precipitation, summer precipitation (June, July, and August), spring precipitation (March, April, and May), autumn precipitation (September, October, and November), precipitation of the driest month, annual mean temperature, mean temperature of the warmest month, mean temperature of the coldest month, and a continentality index (difference between mean temperature of the warmest month and mean temperature of the coldest month). Three spatial variables were also recorded for each population: longitude, latitude and altitude. We carried out a principal components analysis (PCA) with varimax rotation to reduce the number of environmental variables and allow a more synthetic interpretation. We decided the number of principal components to retain running a parallel analysis with 1000 iterations (Hayton et al., 2004) and selecting those principal components with eigenvalues for observed data larger than those from simulations. The analysis was performed with psych package (Revelle, 2017) on the R software v3.3.2. (R Core Team, 2016). Selected variables with loadings above 0.80 were used in Spearman correlation analysis with bark thickness. Natural fire frequency data (lightning fires, of relevance in the east part of Spain; Vázquez de la Cueva and Moreno, 1998) were defined as the number of fires in 90,000 ha of surface covered by P. halepensis forests each year during the period 1974–2010. A 10 × 10 km grid-unit was used to derive this fire frequency data for the 13 Iberian populations of P. halepensis. Unfortunately, we lacked data sources of, similarly, reliable and thorough fire information for the non-Iberian populations. Spearman correlations were done with Hmisc package (Harrell et al., 2018), implemented in R software v3.3.2. (R Core Team, 2016).

Trees at the drier site showed a lower percentage of bark volume (43%) compared to the moister site (47%; Table 3), an apparently modest difference but statistically highly significant [F_{(1, 6)} = 23.87, P = 0.003]. Differences among populations were also significant [F_{(18, 468)} = 1.87, P = 0.016], with average percentages of bark volume ranging from 42 to 47% (Supplementary Table S2). Site by population interaction was not significant [F_{(18, 468)} = 0.88, P = 0.608; Supplementary Table S3], indicating a similar population ranking for percentage of bark between sites. Total volume as a covariate for percentage of bark was highly significant [total volume: F_{(1, 430)} = 112.98, P < 0.001], as well as volume × site interaction [F_{(1, 430)} = 7.61, P = 0.006]. The rest of factors and interactions were not significant (Supplementary Table S4). The allometric effect of total volume on percentage of bark volume was negative (Figure 3). While at lower total volumes, this allometric effect was similar between sites, at higher total volumes differences between sites were magnified.

Regarding the allometric (SMA) analysis, half of the populations showed constant investment in bark with tree size, i.e., isometric allometric exponents (b = 1), at the drier site and the other half, decreasing investment with tree size (b < 1). At the moister site, most populations had isometric allometric exponent (b = 1), except three populations showing decreasing investment in bark (b < 1; Supplementary Table S5). In general, populations at the moister site showed slopes closer to 1 (mean b = 0.92) than the same populations at the drier site (mean b = 0.86; differences between slopes were significant: LRT = 5.56, df = 1, P = 0.018).

Absolute bark thickness was significantly lower (P < 0.001) at the site with lower water availability both at breast and basal heights (Table 4 and Supplementary Table S6). We performed this analysis for all 19 populations taken together and specifically for the 13 Iberian populations (for which we have fire record data). Results for the Iberian populations followed the same pattern as those for all populations, so we only show the results for all populations.

At the tree base, all populations showed bark thicker than the low critical threshold (10 mm) at both experimental sites (Figure 4A–D). By contrast, while all populations reached the high threshold thickness of 20 mm at the moister site (21.3–26.3 mm, Supplementary Table S7), not all of them did at the drier site (17.3–24.3 mm, Supplementary Table S7 and Figure 4). As for the bark thickness at breast height, only three populations exceed the 10 mm thickness threshold at the drier site. At the site with higher water availability, most of the populations showed a bark thickness higher than 10 mm, but below 20 mm- only two populations did not reach the low critical value of 10 mm-. At breast height, the range of differences among populations within sites was more than twofold. Moreover, the ranking of populations for bark thickness at breast height remained rather stable: populations with higher growth rates (from Greece and Italy) had higher BT130, while the populations with lower growth rates (Tunisia and southern Spain) showed lower BT130 at both test sites. However, at the tree base the populations ranking varied considerably between the drier and moister sites, without clear patterns. Only one of the faster-growing populations (P211) had high BT10 at both test sites and one population from southern Spain remained below the average at both sites (population P172, Supplementary Table S7).

Bark thickness at breast height was closely linked to sapwood diameter among populations, whereas at the tree base bark thickness was little correlated to sapwood diameter, especially at the drier site (Figure 4A–D). To help understanding the differences in plasticity of bark thickness between sites, we depicted the sapwood-bark trends between sites for three extreme provenances: one that did not reach the minimum critical bark thickness even at the moister site (P152), another one in which the plasticity conducted to a different critical bark thickness between sites (P185), and finally a third population that achieved the critical bark thickness even at the drier site (P211, Figure 4E,F).

Parallel analysis associated with PCA for environmental variables at the populations’ origins revealed two principal components (variance explained: PC1 = 45%, PC2 = 36%; Martín-Sanz, 2018). Considering loadings above 0.80, the first principal component was positively related to spring and summer rainfalls, as well as to the precipitation of the driest month. PC2 was negatively related to the continentally index and positively to autumn rainfall (see Supplementary Table S8). Pooling together data of both sites, BT10 was correlated (at 90% confidence, P = 0.10) with summer rainfall and precipitation of the driest month (Spearman rho = 0.46, P = 0.049, rho = 0.42, P = 0.072, respectively). BT10 just at the drier site was correlated with summer and spring rainfall, and precipitation of the driest month (rho = 0.50, P = 0.031, rho = 0.42, P = 0.074, and rho = 0.53, P = 0.020, respectively, Figure 5). At the moister site we did not find any association between bark thickness and environmental variables (Supplementary Table S9). Natural fire frequency at the original populations’ habitat was positively correlated with bark thickness at the tree base, but only at the drier site (Spearman rho = 0.62, P = 0.024).

This study highlights the existence of phenotypic plasticity for bark traits in a conifer, responding to different resource availability among environments. Using a replicated common garden experiment at two highly contrasting sites we could differentiate environmental, allometric and genetic effects (Merilä and Hendry, 2014).

Following our expectations, at the drier site our sampled P. halepensis trees allocated fewer resources to the bark, even considering the large difference in size among sites (see Figure 3 where the percentage of bark volume is compared between tree equal sizes -equal total volumes-). Noteworthy, the decreasing relationship between the percentage of bark volume and the total volume was more pronounced at the drier site, confirmed also by the allometric analyses.

Instead of using the usual ratio with stem diameter (see for example, Paine et al., 2010; Pausas, 2015; Midgley and Lawes, 2016) we followed a new methodological approach intending to be closer to the real allometry between the bark and the rest of the tree. We considered this better defined by comparing the percentage of bark volume vs. the total volume of the trunk, considering also that all stem was sapwood in our sampled P. halepensis trees (Supplementary Figure A.1). As expected, due to the well-documented relationship between bark thickness and stem diameter (Adams and Jackson, 1995; Pinard and Huffman, 1997; Lawes et al., 2013; Poorter et al., 2014), bark volume was closely related to total volume. This allometric effect, which has an adaptive significance, implies that the evolutionary forces that act on plant size will produce changes in bark thickness (Rosell, 2016; Rosell et al., 2017). Due to this relationship, when the goal is to compare differential allocation to the bark among species or populations, plant size must be taken into account (Hempson et al., 2014; Poorter et al., 2014; Rosell and Olson, 2014). However, while differential relative allocation to bark is essential under a functional-evolutionary perspective, survival to fire clearly depends on absolute bark thickness (Midgley et al., 2010; Lawes et al., 2011b). Confirming our main hypothesis, the lower allocation to bark in the growth-limiting environment led to lower absolute bark thickness bark, such that at 18 years of age (close to an average-low fire return interval for P. halepensis) trees did not achieve the critical thickness of 20 mm. Moreover, our results confirmed a strong population genetic effect on absolute bark thickness, with eastern Europe populations showing thicker bark, and populations from North Africa and Southern Spain displaying thinner bark. Difference of bark thickness at breast height among extreme populations was more than twofold (7.2 mm -population 152– and 15.8 mm –population 211–). This difference among populations is comparable to that found for trees of similar height in P. pinaster (Tapias et al., 2004), one of the pine species with best known strong ecotypic variation (see for example Alía et al., 1995). Moreover, plasticity of bark absolute thickness was rather homogeneous among provenances, hence the site-related scarcity of resources led to lower bark thickness affecting all populations alike.

While it has been shown that inner living bark is involved in the storage and translocation of water and photosynthates (Srivastava, 1964; Scholz et al., 2007; Romero, 2014), thermal protection depends on the outer dead bark thickness (Romero and Bolker, 2008; Pásztory and Ronyecz, 2013), rather than on the total bark thickness (Graves et al., 2014). Moreover, it seems that the decreasing rate in outer bark thickness with tree height is greater in fire-resistant species than in fire-sensitive ones, in contrast to inner bark thickness which diminishes with tree height or total bark thickness which is maintained along the bole (Graves et al., 2014). Therefore, differentiating between inner and outer bark thicknesses in future studies will help us to better understand the ecology and functional role of bark.

The environmental correlations found in this work suggest the existence of ecotypic patterns in this species that deserve further investigation to understand their evolutionary meaning. We found a highly expected significant correlation between higher frequency of natural forest fires and thicker basal bark. Among the environmental characteristics at the provenances’ seed source, the most significant correlation occurred between higher rainfall and thicker bark. These results combined may appear counter-intuitive, but it is necessary to bear in mind that in Mediterranean ecosystems higher rainfall leads to high fuel accumulation that may increase fire hazard during yearly summer droughts. Actually, permanently dry forests are usually less likely to burn due to lower fuel loads, even when trees may be physically more flammable (Fuhlendorf and Smeins, 1997; Kitzberber et al., 1997; Fernandes, 2013). A more accurate estimation of the fire regime including not only fire frequency data but also intensity, spread patterns or seasonality, and amount and type of fuel, would indeed allow us to better understand the role of fire on the variation of bark thickness in P. halepensis.

As an obligate seeder, the vital strategy of P. halepensis is the maintenance of an aerial seed bank in serotinous cones, enough for the persistence of the population in case of fire (Tapias et al., 2001; Ne’eman et al., 2004). Therefore, it is essential that trees can survive low or moderately intense fires until accumulating this “sufficient” aerial seed bank. This is usually achieved when P. halepensis trees are between 15 and 20 years old (Moreira et al., 2011), but both intraspecific variation and plasticity can largely modify that reference age (Martín-Sanz et al., 2016). Lacking specific data on the critical minimum bark thickness for this species, we chose two extreme values based on the literature (see van Mantgem and Schwartz, 2003; Wesolowski et al., 2014; Pausas, 2015; Madrigal et al., 2019). We found that our trees generally reached the critical bark thickness at the base, but not at breast height. This is consistent with a greater thermal insulation at the stem base and a steep bark tapering along the bole found in several Mediterranean pine species (De Ronde, 1982; Pageaud, 1991; Jackson et al., 1999). This can be regarded as an adaptive solution to reduce immaturity risk (Lamont et al., 1991; Keeley et al., 1999), i.e., ensuring individual survival until reaching a sufficient aerial seed bank enabling recruitment after lethal fires. Importantly, our studied trees did not reach the critical breast height bark thickness of 20 mm and, at the base of the tree, this value was only achieved at the moister site. Therefore, as predicted, the plasticity associated to different resource availability (mostly water availability during the vegetative period, since light, nutrients and CO₂ are not limiting in our study sites) clearly affects the probability of P. halepensis to survive fires; the lesser resources, the lower expected survival facing fires. Interestingly, depending on the species and habitats, variable and even opposite relationships have been found between bark thickness and productivity or soil fertility (Climent et al., 2004; Schubert, 2014; Jager et al., 2015; Richardson et al., 2015). However, as far as we know, our study is the first to assess the strict environmental effect separated from the genotype-environment interaction; thus, the comparison with earlier works is not straightforward.

The maximum life expectancy of P. halepensis due to fire return interval is usually estimated between 20 and 50 years, while the average is about 25 years (Agee, 1998; Vázquez de la Cueva and Moreno, 1998; Arianoutsou, 2001; Tessler, 2012). However, fire interval could be as short as 6 years in some areas of its natural distribution (Tessler, 2012). When fire interval is shorter than 15 years, P. halepensis recruitment would be limited, and its populations may totally disappear (Roitemberg and Ne’eman, 2000; Eugenio et al., 2006; Herman, 2009; Tessler et al., 2014). Therefore, the age of our sampled trees -18 years- is highly meaningful considering the species’ life-history.

Preceding studies in P. halepensis have shown relevant ecotypic trends for growth, reproduction and water use efficiency such that provenances from dry continental environments showed early and intense reproduction and serotiny, as well as conservative water use in detriment of vegetative growth (Ferrio et al., 2003; Voltas et al., 2008, 2015; Santos-del-Blanco et al., 2013). Therefore, adaptive geographic variation of this species seems to be derived from differential resource allocation among key life-history processes (growth, reproduction, fire tolerance, constitutive and induced defenses, and water stress tolerance; Sampedro et al., 2014; Climent et al., 2015).