Resin tapping influence on maritime pine growth depends on tree age and stand characteristics

Resin is a renewable forest resource that can increase the economic value of some forests in rural areas. Resin production is associated with climatic conditions. However, its impact on trees’ growth remains unclear. Here, we studied radial growth in six Portuguese Pinus pinaster forests that had been resin tapped for 5–7 years along a latitudinal and climatic gradient to understand whether resin tapping affects tree growth, and how it is affected by climate, stand and tree traits. Tree-ring width (TRW) on tapped and untapped trunk sides was compared before and during the tapping period. Tree-ring width decreased in the three youngest populations (< 30 years), with no changes in older populations (> 40 years), while TRW increased after resin harvesting began in the oldest stand (> 55 years). Annual resin-tapping impact (RTI), calculated as the ratio between TRW during tapping years and the 5-year average TRW before tapping, was below and above 1 in the younger and older stands, respectively. Among stand characteristics and across sites, RTI was negatively correlated with tree competition and positively correlated with TRW, cambial age, and tree height. Climatic conditions had a minor role on tree growth response to resin tapping. Our main conclusion is that the effect of resin extraction on growth is age-dependent. Our results encourage the co-production of resin and wood on maritime pine stands over 40 years old where resin extraction impact on growth is negligible.


Introduction
Resin tapping and harvesting is a traditional forestry activity that provides economical profit to rural populations and other services like fire surveillance and sustainable forest management (Pereira 2015;Soliño et al. 2018). Currently, there is a new demand for natural resin (Génova et al. 2014) leading to a price increase and to the resurgence of the activity in places where it was almost abandoned, such as in the Iberian Peninsula (Girón, 2021).
Pine resin extraction has been usually done alongside timber production in the Mediterranean basin (Soliño et al. 2018), especially in Portugal and Spain, where Pinus pinaster Aiton is the main species used for resin harvesting. Despite this, it has been hypothesized that secondary metabolism and growth may be competing for the available photosynthates (Züst & Agrawal 2017) and thus, repeated mechanical resin tapping could be stimulating resin biosynthesis at the expense of other carbon-dependent processes such as wood production (Du et al. 2021). Supporting this idea, lower growth rates and narrower tree-rings, as an estimate of biomass production, have been found in resin tapped pine stands (Papadopoulos 2013;Génova et al. 2014;Chen et al. 2015;Zeng et al. 2021). Papadopoulos (2013) showed a reduction of 14% in tree-ring width (6% of latewood) in tapped Pinus halepensis Mill trees, compared with non-tapped ones. Zeng et al. (2021) also found negative effects of resin harvesting on tree-ring width of Pinus tabuliformis Carr, although it only persisted for two years after the extraction ended. But, positive (Tomusiak & Magnuszewski., 2013;van der Maaten et al. 2017) or no resin tapping effects (Rodríguez-García et al. 2015;Du et al. 2021) on tree-ring formation have also been reported, suggesting low Communicated by Thomas Seifert.
1 3 carbon competition between growth and resin, which according to Hilty et al. (2021) can be considered as growth sensu latu-the combination of reversible and irreversible biomass or volume changes. Growth stimulation in tapped trees could also reflect reaction tissues if annual growth was measured too close to the tapping wound.
The contradictory results found in the literature concerning resin extraction effects on wood production might be also related to the studied species, the sampling design, and the resin tapping method and intensity. For example, Zeng et al. (2021) sampled both tapped and untapped sides to represent the overall tree-ring growth and observed that resin tapping in P. tabuliformis promotes asymmetrical radial growth, with narrower rings near the tapped face, while Garcia-Forner et al. (2021a, b) showed the opposite pattern in P. pinaster. It was also observed that wood reaction to tapping decreases when getting away from the wound upward as shown in Pinus sylvestris L (Tomusiak & Magnuszewski., 2013). Thus, resin extraction impact on growth will depend on the sampling design performed (Williams 2017) hindering comparisons between studies. But even when the same sampling procedure and species are used, other variables, e.g., environmental conditions, can affect growth sensitivity to resin harvesting (Garcia-Forner et al. (2021a, b).
Tree growth and resin production are both driven by environmental conditions and depend on carbon availability to produce new cells, which include the resin duct complex formed by the resin canal and the adjacent resin-producing epithelial cells that synthesize the defensive terpenoid mixture. In P. pinaster, some studies connect favorable temperatures, solar radiation and potential evapotranspiration with increasing tree growth and resin duct traits or resin yield (Rodríguez-García et al. 2014, variables that are phenotypically and genetically correlated (Vázquez-González et al. 2019). However, literature is still scarce when it comes to understand the resin tapping impact on wood production under a climatic gradient, the key factor to identify the environmental conditions that maximize both resin and wood production. Another variable that should be considered is the species provenance, as growth, as well as resin yield are subjected to strong genetic control varying among families and provenance regions (Vázquez-González et al. 2021). If resin harvesting depletes carbon stocks, wood production in fast-growing populations could be more susceptible to resin extraction than in slow-growing ones (Zeng et al. 2021).
Stand and tree characteristics and forest management are also important factors that should be considered since they may determine resources availability at individual level, and thus affect their growth rate and production of defenses. For instance, at tree level, resin yield increases at lower stand densities (McDowell et al. 2007;Rodríguez-García et al. 2014Zas et al. 2020). In high density stands, trees grow less, and smaller trees can be more affected (e.g., P. tabulaeformis in He & Wang 2021) or end-up being less resilient to drought (Fernández-de-Uña et al. 2015;Bottero et al. 2017). The age of the trees affects their vitality, growth, and their resin yield too. In this regard, Zas et al. (2020) found that even if the mechanisms remain elusive, resin production in P. pinaster increases as trees get older. Considering that in woody plants photosynthesis tends to decrease with age (Yoder et al. 1994;Bond 2000), it is necessary to understand whether resin yield increases as age increases due to carbon allocation shifts from growth to defenses, or because of higher carbon reserves in older trees.
Our aim is to understand whether resin tapping affects tree-growth rates, and how tree, stand characteristics and edaphoclimatic conditions modulate it. We have studied six populations of P. pinaster along a latitudinal gradient in Portugal and assessed tree-ring width before and during resin tapping. Our hypotheses are that: i) resin tapping reduces tree-ring width due to shifts in carbon allocation to sustain resin production, but the impact will depend on ii) local conditions, and iii) the characteristics of the trees, i.e., trunk girth and height, tree density and competition between neighbors, with younger trees being more affected due to a lower capacity to acquire carbon and/or lower carbon reserves.

Study sites and resin harvesting
Maritime pine (Pinus pinaster Aiton) is a fast-growing heliophytic tree species with a great ecological plasticity and varied climatic spectrum (Farjon 2010). In Portugal, this evergreen conifer represents 22% of the total forest area (ICNF 2015) and is the main species used for resin tapping. According to Portuguese law (see Art. 4 th DL n.º 181/2015) resin tapping is performed in trees with more than 63 cm of trunk girth at 1.30 m from the ground, which is in line with previous studies that recommended resin tapping only in stands older than 30 years . Tapping starts on the trunk base at 20 cm above the ground by removing the bark and phloem every 15-20 days and continues upward for four years and up to 2 m high. Afterward, when tree size allows it, another vertical row is started 10 cm apart from the first wound and following the same procedure. Usually, the tapping season goes from March until November, coinciding with the radial growth of maritime pine.
For this study, we selected six stands of P. pinaster for resin production in the Northern and Central regions of Portugal, areas where this activity is traditionally done. Selection criteria of the stands was based on three main aspects: (1) to have a minimum of five years of tapping exploitation; (2) to have been exploited with only one resin wound at a time and thus harmonizing the resin tapping intensity between stands; (3) to have used the same resin extraction methodology. We selected paired stands relatively close climatically, but different in the average age to assess its role on the impact of resin tapping on growth. The selected populations were located at Paredes de Coura (PCO), Vila Pouca de Aguiar (VPA), Olho (OLH), Soure (SOU), Marinha Grande (MGR) and Pataias (PAT) (Fig. 1). Across study sites, resin tapping occurred between 2013 and 2019 ( Table 1). The number of tapped years was 5 in OLH and SOU; 6 in PCO, PAT and VPA and 7 in MGR (list of abbreviations in Table S1).

Climatic data collection
Climatic data for each location was downloaded from the Climate Explorer of the Royal Netherlands Meteorological Institute (https:// clime xp. knmi. nl) using grid points 0.5 or 0.25. Due to their close proximity, MGR and PAT were located at the same grid point, and thus shared climatic data. Monthly time series of mean air temperatures (average, maximum and minimum; °C) and total precipitation (mm month −1 ) were obtained from E-OBS 21.0e 0.25. Self-calibrating Palmer Drought Severity Index (PDSI) calculated from Climatic Research Unit Time Series (CRU TS) data was obtained for 1901-2017 0.5° Global 3.26 early data set, and vapor pressure (hPa) from CRU TS 4.04 0.5° data set. Vapor pressure deficit (VPD, hPa) was calculated from average temperature and vapor pressure using the following formula: where T and VP are the means of the monthly average air temperature and vapor pressure, respectively. Maximum and minimum VPD were also calculated for each location using the respective temperatures and vapor pressure.
(1) VPD = 6.11 * exp ((17.27 * T)∕(T + 273.3)) −VP Fig. 1 Location and climatic diagrams of the study sites. The color in the map represents the mean annual precipitation minus potential evapotranspiration (P-PET, in mm) for the period 1990-2019. Negative P-PET values indicate drier sites. For each site, an ombrothermic diagram is shown representing the average monthly air temperature (T, in red) and precipitation (P, in blue) for the period 1990-2019, and the gray area indicates the drought period (2 T > P). The annual average of air temperature and precipitation for a period of 30 years is displayed at the top of each climatic diagram

Stand and tree characteristics
Twenty P. pinaster trees per location were selected to investigate growth responses and the structure of the stand. We first evaluated the trees' canopy and then selected healthy trees with average trunk girth for each stand. Trunk girth at breast height (TG), total height and tree height until the first branches with green needles were measured. The difference between the two last variables was used as an estimate of the tree canopy size (trunk foliated length).
To quantify the degree of competition between trees and its potential effect on radial growth of selected trees, we calculated two competition indices based on trees' height and trunk girth (CI H and CI TG , respectively). This information is not available for PCO because trees were felled before collecting all data. We first identified all living trees in a circumference of 5-m radius from the trunk of each target tree, and then measured their trunk girth at breast height and the two described heights. The number of present neighbors was used to calculate tree density of the location (trees ha −1 ). CI H and CI TG were calculated as the sum of the ratio between each competitor tree in the circumference and the target tree as follows: (2) where tg i is the trunk girth at breast height of i th competitor and tg is the trunk girth at breast height of the target tree; and where h i is the trees' height of i th competitor and h is the height of the target tree. Two soil samples per target tree were also collected at 50 cm (0-10 cm depth) from the north and south side of the trunk. Pooled samples were dried at 60 ºC for four days and sifted with a 2-mm sieve. Soil pH was measured with a pH meter (JENWAY Model 3505) in a 1:2.5 (mass/volume, soil/ water) suspension. The percentage of organic matter (OM) was assessed by loss on ignition at 450 °C for 4 h in a muffle furnace (Nelson & Sommers, 2018;Roberston 2011).

Dendrochronological sampling and processing
Wood increment cores were collected on the 20 target trees per site in 2019 at PCO and VPA and in 2020 on the remaining sites. Two cores were taken from each tree at breast height (ca. 1.3 m), with an increment borer of 12-mm diameter. Core extraction was done between resin-tapping wounds, as all studied trees had been tapped for more than four consecutive years. The core between tapping wounds represents the resin tapping face and the one from the opposite side represents the untapped face (Fig. 2). The orientation of each core was recorded to verify if the exposition of the tapping wound influenced growth response to tapping within and between sites. After harvesting, each core was immersed in 99% acetone for 24 h at 4 ºC to clean the samples from resin, dirt, and debris. Then, they were air dried, glued on wooden mounts and progressively polished (150, 400, 800 and 1200G) with a belt sander (Bosch GBS 100 AE). Polished samples were scanned at 2400 dpi and 48 bits in RGB using a photo-scanner (EPSON Perfection 4990 Photo) and the EPSON Scan software (2020 Seiko Epson Corporation). Increment core images were analyzed using the xRing package for R to automatically detect tree-ring borders and to measure tree-ring width (TRW ) (Campelo et al. 2019; R Core Team 2020). After visually cross-dating by comparing TRW within tree (tapped and untapped sides) and site, the number of annual rings from the pith was used to determine the cambial age. The pith was not present in 1 tree at VPA, PCO, MGR and SOU and in 3 trees at OLH. Trees with correlations below 0.4 were discarded for further analysis. Final chronologies were computed with 15 trees for OLH, 18 for PCO and SOU, 20 for MGR, PAT and VPA. To remove size and age-related trends, a cubic smoothing spline was applied to each individual ring-width time-series. Due to the large variation in the average length of ring-width time series between sites, the recommendations of Klesse (2021) were followed to define the length of the spline. Thus, the grand mean of the mean of the time-series length for each site was used. In doing so, each ring-width time-series was detrended using a spline with a frequency response cut-off of 0.5 and a length of 26 years to enhance the climatic signal and to allow intersite comparisons. Autoregressive modeling was performed on each standardized series to remove temporal autocorrelation and to maximize the climatic signal. For each site, the obtained residual ring-width indices (RWI) were averaged using a bi-weighted mean to produce a site chronology. For the common interval to all sites, the statistical quality of each chronology was evaluated using the following variables: the expressed population signal (EPS) that indicates the degree to which a given chronology approaches an infinitely replicated chronology (Wigley et al. 1984a, b), the r.tot which is the mean of all correlation between different cores, r.wt which is the mean of correlations between series from the same tree over all trees and r.bt which is the mean inter-series correlation between all series from different trees. The chronologies and detrending were developed using dplR (Bunn 2010; R Core Team 2020) and detrendeR (Campelo et al. 2012) packages for R version 3.6.2 (R Core Team 2020).

Statistics
We characterized the climatic differences between sites by conducting a principal component analysis (PCA) of the main climatic variables between 1999 and 2019, the common period for all the chronologies. Variables with low contribution percentage among principal components were removed from the PCA. The variables included in the final PCA were: mat (mean temperature), pp (mean annual precipitation), ts (temperature seasonality), maxvpd (maximum vapor pressure deficit) and mpdsi (mean Palmer Drought Severity Index). All variables were previously scaled, and the analysis was computed using the function PCA from the FactoMineR package (Lê et al. 2008) for R version 3.6.2 (R Core Team 2020).
The parameters of the studied trees (TG, tree height, trunk foliate length) and the structure of the stand (tree density, core exposition, soil pH and OM, CI TG and CI H ) between sites were compared using linear models. All response variables were normally distributed, except soil pH and OM that were square root transformed to meet the assumptions of parametric analyses. The respective tree/stand parameters were the dependent variable, site the fixed factor and tree the random factor. Bonferroni method was used for the pairwise post hoc tests and alpha < 0.05 as a significant level. These analyses were performed using the function from the nlme package v.3.1-14 (Pinheiro et al. 2020) and the cld function Fig. 2 Representation of an increment core sampled on the trunk of a tree with two adjacent tapping wounds. Concentric circles represent annual tree-rings. The green area represents the wound area where resin tapping-bark and phloem removal-was performed, and the red arrow shows the direction of the increment core, which starts between resin tapping wounds (tapped side) and ends on the opposite side (untapped side) of the trunk. The distance between tapped faces is around 10 cm according to the Portuguese law regulating the harvesting and marketing of pine resins (see Art. 4th DL n. º 181/2015) from the multcomp during pairwise post hoc tests (Hothorn et al. 2008) function for R (R Core Team 2020).
To comprehend if resin tapping affects tree-ring width at site level, we compared the TRW during the tapping years and the same number of years before the start of resin harvesting, e.g., 2006-2012 before tapping vs. 2013-2019 during tapping in MGR. Then, linear mixed models were performed using TRW as dependent variables, square root transformed to meet a normal distribution. Site, tapping (before vs. during), and trunk side (tapped vs. untapped) were considered as fixed factors, while tree within site was considered a random factor. Pairwise post hoc tests with Bonferroni corrections were carried out to compare differences between trunk side, and period of tapping for each site. The same procedure was applied to RWI to assess the results after removing age influence on growth trends.
To compare resin harvesting impact between years and stands with different growth rates, we standardized data calculating the resin tapping impact (RTI), i.e., the ratio between TRW of each tapping year and the mean TRW for the 5 years before tapping started. The RTI corresponds to the resistance index proposed by Lloret et al. (2011), which estimates the performance during and before a disturbance, i.e., resin tapping in this study. We have used 5 years averages for the pre-disturbance period to reduce the weight of any natural trend, and to account for inter-annual climatic variability. Values > 1 and < 1 indicate that wider and narrower rings were produced during tapping relative to the previous period, respectively. To compare the RTI between tapping years and sites, we used linear mixed models where RTI was the dependent variable, site, and the tapping year the fixed factors and tree the random factor. More complex models with autocorrelation were considered but did not improve model performance according to Akaike and Bayesian information criteria. When significant effects were detected, a pairwise post hoc test with Bonferroni corrections was applied to account for multiple comparisons between years and sites. For these analyses, it was used the lme function from the nlme package v.3.1-14 (Pinheiro et al. 2020) and the cld function from the multcomp during pairwise post hoc tests (Hothorn et al. 2008).
To assess the effect of tree traits, stand characteristics, and climate on RTI, the Spearman's rank correlation coefficient was calculated within and across sites. The parameters used for the tree and stand characteristics were TRW , TG, tree height, tree foliated height, tree density, the competition indices (CI TG and CI H ), soil pH, soil OM and cambial age. For the PCO site, where trees were felled before we could collect data for radial growth. Only TRW , TG, and the soil pH and OM variables were used. The climatic parameters used were based on the scores of the first, second, third, and fourth principal components of the PCA (mat, ts, pp, maxvpd and mpdsi). These analyses were conducted using rcorr.adjust function of the Hmisc package for R version 3.6.2 (R Core Team 2020).

Climate and stand characteristics
The sampling sites showed a latitudinal climatic gradient with the coldest location, VPA, with an average temperature of 12.1 ºC, also located at the highest altitude, 791 m ( Table 1, Abbreviations in Table S1). The warmest location was SOU, whereas MGR and PAT were the driest ones, with negative P-PET values (Fig. 1). PCO showed the highest precipitation levels (1471 mm on average) and highest average VPD (6.88 hPa) (Fig. 1, Table S2). The average Palmer Drought Severity Index (PDSI) was close to 0 and similar between sites, but slightly more negative in the lower latitudes (Table S2).
The PCA retained 5 climatic variables: mat, pp, ts, maxvpd and mpdsi. The two first principal components (henceforth described as PC1 and PC2) explained 38.9% and 23.2% of the overall variance, respectively (Fig. 3). The variables ts and maxvpd were positively related to PC1, while mat, pp and mpdsi were negatively correlated. Precipitation (pp) was negatively correlated with PC2, while the other variables were positively related. The PC1 and PC2 aggregated the sites in three main groups: i) PCO, ii) VPA and iii) OLH, SOU, MGR and PAT (Fig. 3). PCO was positively related to pp and mpdsi, and negatively to mat, while VPA had a positive correlation with ts. The rest of the locations correlated positively with PC2, with mat, maxvpd and mpdsi (Fig. 3).
Soil OM and pH were different among sites (Table 1) with PCO and VPA (Table 1) showing higher OM and lower pH values. Across sites, average cambial age varied from 58 to 29 years-old TG and height between 88 and 108 cm and 18 and 22 m, respectively (Table 1). The biggest and tallest trees were found in MGR, the oldest stand; however, TG in PAT, the stand located at the same climatic grid point, was similar, even if the latter was the youngest sampled population (Table 1). Marinha Grande (MGR) and SOU showed, on average, less than two neighbors, in contrast with the 3.4 per target tree identified in OLH; however, competitor trees were smaller in the latter (Table 1). Olho (OLH) was the densest stand with the highest inter-tree competition (554 trees/ha and CI TG 2.30), and MGR the one with the lowest density and competition indices (357 trees/ha and CI TG 1.60). Across locations, the tapped side of the trunk (tapped increment core) tended to be exposed to south, varying between southeast and southwest (Table 1).

Tree-ring width before and during tapping
In general, the tree-ring width index series considering the common interval within (Table S3) and between sites (Table S4;   showed EPS values above 0.85, except in OLH and SOU that showed slightly lower values from 1999 to 2019 and including the tapped and untapped side of the trunk. PAT, the youngest stand, was the site with the highest average TRW values (5.53 mm), followed by OLH (4.22 mm), PCO (3.73 mm), VPA (2.54 mm), MGR (2.27 mm) and SOU (1.91 mm). In general, all chronologies showed the age growth trend producing narrower tree-rings with increasing tree age (Fig. 4). TRW of the untapped side was smaller in PCO and OLH; however, this trend was present before the resin tapped period (Fig. 4). To understand if TRW was affected by resin harvesting, we compared the mean TRW between the resin tapped years and the same number of years before tapping with the two principal components explaining 62.1% variance overall. Different sites are depicted by different symbols and colors. Due to its proximity, MGR and PAT share the same climatic data calculated separately for the tapped and untapped trunk sides (Fig. 5). Pine trees in two of the youngest populations, PCO and OLH, showed a significant decrease of the TRW during the tapping period on both sides of the trees, while PAT-the youngest and fastest growing population-showed a TRW decrease only in the untapped side (Fig. 5). All the other sites (VPA, SOU and MGR) showed no significant growth decrease during tapping, regardless of the side of the tree (Fig. 5). In MGR, the stand with the oldest tree population, TRW increased during the tapping period on the tapped side of the trunk (Fig. 5). The same analysis was performed using the treering-width index (RWI), where the TRW age trend was removed (Fig. S1). In this case, RWI increased during the tapping period in VPA and MGR, regardless of the side of the tree (Fig. S1). OLH showed the opposite trend, that is, RWI was higher before than after the tapping started, while PCO, PAT and SOU showed no significant differences (Fig. S1).

Fig. 5
Boxplots of the tree-ring widths (TRW , mm) for each side of the tree (tapped in blue and untapped in red) before and during resin tapping. Different letters indicate statistically significant differences between groups at each location (p < 0.05)
The Spearman's rank correlation coefficients between RTI, tree, and stand characteristics across sites were significant and positive for TRW (+ 0.38), cambial age (+ 0.53), total height (+ 0.24), and soil OM (+ 0.11), and negative with CI TG (− 0.22) and soil pH (− 0.14) ( Table S5, Fig. 7).  (Table S5). In general, TG was strongly and positively associated with RTI in VPA and PAT, to total height in PAT and OLH, but to soil OM only in PCO. Negative RTI were shown for cambial age in MGR and PAT (Table S5).
Spearman correlations were carried out between the RTI and the scores of the first, second, third and fourth components of the PCA (Fig. 8). The scores of the principal components were not significantly associated with RTI, except for the PC4 that preserved only the 10.6% of the climatic variance and was mainly explained by mat (32.8%). In the latter, with the increasing mat, the RTI decreases to values below 1 (Fig. 8d). Fig. 7 Spearman's rank correlations coefficients between RTI with cambial age a, total height b, trunk girth c and CI TG d across sites. Each gray dot represents individual trees (average RTI during the tapping period and the characterization variable measured during the sampling campaign) and the blue dots represent the average value per site with its respective standard error bar (SE) and regression line (in red). Relationships between variables are represented by dashed lines when non-significant and by solid lines when significant. The Spearman's correlation coefficient (R) and the p-values are shown (for p < 0.05)

Discussion
The aim of this study was to understand if resin tapping affects wood growth of P. pinaster, and whether tree, stand characteristics and edaphoclimatic conditions can modulate it or not. Our results show that radial growth is not systematically reduced on maritime pine trees exploited for resin production and hint that negative impacts are more frequent in younger (< 30 years-old) than older populations as hypothesized. This result indicates that, as well as being a profitable forest activity, resin and timber coproduction in maritime pine may be sustainable, especially when performed in populations older than 40-year-old. Fig. 8 Spearman's rank correlation coefficients between the resin tapping impact (RTI) and the scores of the first a, second b, third c and fourth d climatic principal components (PC1, PC2, PC3 and PC4; see Fig. 3) for all studied locations. Each dot represents each tapping year for each site. Relationships between variables are represented by dashed lines when non-significant and by solid lines when significant (p < 0.05). The PC1 is more weighted by the variable ts and the PC2 by pp. PC3 has mpdsi has the highest contributed variable and PC4 is with mat. The Spearman correlation coefficient (R) and the p-values are shown In general, when the TRW before and during resin tapping is compared, the youngest populations (PCO, OLH and PAT) showed a decreasing growth trend during resin extraction, while the oldest populations showed no significant differences (VPA, SOU) or even wider rings during this period (MGR). Tomusiak & Magnuszewski (2013) and van der Maaten et al. (2017) also showed positive influence of tapping on radial growth on Scots pine and no systemic reductions in growth as we found in the older stands. On the contrary, Chen et al. (2015) showed that younger tapped Masson pine trees were able to enlarge and repair their trunks faster than older trees, as their tree-ring index only declined by 4.8%. However, they only investigated two sites, and no details about stand structure were given. Discrepancies can also be partially explained by the use of non-detrended series. Studying resin tapping impact through tree-ring width may be confusing because the biological trend associated with tree age is still present. Yet, our results consistently showed that growth resistance increases in older stands when RWI instead of the TRW was considered. In this case, two of the oldest stands (VPA and MGR) showed higher growth indices during tapping than the predisturbance period. Although resin production was not measured here, Zas et al. (2020) showed that resin production increases considerably with tree age in maritime pine stands growing under similar climatic conditions in northern Portugal. Determining a minimum age to start resin harvesting in this species becomes, therefore, a key factor to design sustainable timber and resin forest exploitations, and future studies should shed light on this.
Besides the stand age, other parameters must be taken into account when trying to understand growth sensitivity to resin harvesting, for instance the design of the increment cores' sampling. Wood formation in the trunk of tapped pines is not uniform as growth is suppressed in the region due to the mechanical injury of the bark and phloem tissues, and asymmetrical rings can be produced in the remaining surface depending on their distance to the wound (Grissino-Meyer et al. 2001). Although our study did not measure the effect of tapping on wood properties, other studies have shown that resin tapping results in higher wood deformation and elasticity (Silva et al. 2018), but also increases wood density (especially on the tapped side), which modifies its mechanical properties (García-Iruela et al. 2016;Silva et al. 2018;Wu et al., 2022). Furthermore, if the wood core is extracted too close to the wound, the increase in TRW can just be related to the wood reaction to the injury. To take this into account, Zeng et al. (2021) extracted six increment cores per tree, and Papadopoulos (2013) measured TRW and latewood in three radial directions in Aleppo pines. Nonetheless, impact assessment, as well as comparisons between studies, is usually difficult as most of the sampling designs usually include the extraction of a single increment core per tree and height taken from the tapped or untapped face (Tomusiak & Magnuszewski., 2013;Hood & Sala 2015;Zas et al. 2020;Du et al. 2021). In this study, we harvested stands subjected to the same tapping intensity (number of wounds at a time and tapping frequency), and collected wood cores using the same methodology, so that the comparison between sites can disclose the effect of stand age and edaphoclimatic characteristics among the studied sites.
Across sites, the locations characterized by taller trees producing wider tree-rings, low inter-tree competition at TG level, and greater content in soil organic matter tended to show higher growth rates during tapping than the previous period (RTIs above 1). On the one hand, larger and taller trees may have better access to light and can potentially photosynthesize more, acquiring more carbon that can be used for either resin or wood production. Furthermore, taller trees usually have deeper roots and are less sensitive to variations in precipitation because they can tap deeper water sources and thus, keep high transpiration and photosynthesis rates for longer under drought (Brando 2018; Giardina et al. 2018). On the other hand, the negative CI TG correlation with RTI suggests that competition for resources could affect the capacity to respond to harvesting, as well as shown in response to other disturbances (Rozendaal et al. 2020;Marqués et al. 2021). In addition, it has been shown that at least in young Pinus oocarpa Schiede ex Schltdl, high competition levels lead to decreasing resin yield, although the opposite can occur in adult trees (Egloff 2020). Therefore, and specially in young stands, it is necessary to carry out forest management practices such as thinning to reduce inter-tree competition, as well as the associated growth and resin production constraints.
Within sites, only PCO-the tree population growing in more acidic soils-showed significant correlations between soil organic matter, pH and RTI, positive and negative, respectively. Organic matter is an important source of nitrogen, a macronutrient for plants and central for photosynthesis, with leaves with higher amounts of nitrogen showing a higher photosynthetic capacity and increasing the growth potential of the trees. The negative association between soil pH and RTI could simply reflect the effect of this variable on the mineralization of the organic matter, and/or to a higher nutrient availability under more acidic soils (Curtin et al. 1998).
Selected pine stands in this study could represent different P. pinaster provenances (ICNF 2012) and thus, their annual growth could be strongly influenced by climatic conditions and, to a lesser extent, by genetic plasticity (Rozas et al. 2020). Despite this, the climatic sensitivity of tapped and untapped trees appears to be similar, as shown by van der Maaten et al. (2017). In their study, tapped and non-tapped P. sylvestris trees presented comparable sensitivities even in the occurrence of an extreme climatic event. In line with this, our results showed that the climatic conditions had a marginal effect on RTI across stands. Only 10% of RTI variance was explained by PC4 indicating that growth sensitivity across sites increases with increasing average temperatures. Without taking into account the potential effect of temperature on resin production and considering the studied region-where resin is traditionally harvested-warmer places seem to be less profitable for timber production when the activity is combined with resin harvesting, a trend that could be exacerbated by rising temperatures.
In conclusion, this study showed that the age of maritime pine stands had a strong effect on the impact of resin extraction on wood growth, with younger stands (< 30 years-old) showing a tendency to form narrower rings-reduced wood production, and a lower capacity to buffer the resin extraction disturbance compared with older stands. Nonetheless, we did not have access to the resin yield of the different stands, an important information to deepen our knowledge on the impact of resin extraction on the carbon economy of maritime pine. Based on our results, to maximize resin and wood co-production, resin extraction should be performed in populations with more than 40 years of age, and on managed stands with controlled tree density to reduce competition and thus provide greater growth resistance to resin extraction.