Beyond width and density: stable carbon and oxygen isotopes in cork-rings provide insights of physiological responses to water stress in Quercus suber L

Study areas

In this study, we collected data from cork rings of selected trees in two distinct cork oak woodlands located in southwestern Portugal: Benavente (CL) (38°49′N, 8°49′W, 20 m a.s.l.) and ; Grândola (BS) (38°11′N, 8°37′W, 270 m a.s.l.) (Fig. 1). These cork oak woodlands have been described by the authors in their previous works (Costa et al., 2022; Costa et al., 2016).

Soil type, geologic formations and related biophysical site characteristics of these areas (CL and BS), fully described in previous work (Costa, Pereira & Oliveira, 2002; Costa et al., 2016; Costa et al., 2020; Mendes & Ribeiro, 2010; Mendes et al., 2016; Mendes et al., 2019), have distinctly affected soil water availability, specifically during the dry and hot summer season. At CL, trees have access to shallow groundwater, the main reliable water source for the maintenance of their physiologic activity during summer (Costa, Pereira & Oliveira, 2002). In contrast, trees at BS have a summer dormancy period (Oliveira et al., 1994) that is a drought stress-imposed rest period (Cherubini et al., 2003). At BS, tree growth must rely on groundwater that it is only used when the water table is reachable by tree roots, in early spring (Costa, Madeira & Oliveira, 2008; Costa et al., 2016).

Both study areas have a Mediterranean climate characterized by dry and hot summers and rainy winters fully described in Costa et al. (2022).

Precipitation and temperature raw data were obtained from the SNIRH (Sistema Nacional de Informação de Recursos Hídricos) database and were computed to obtain series of climate parameters for the period 1962-2012. Precipitation and temperature were calculated monthly, from February of the previous year to November of the current growth year. One (monthly), five- and twelve-months response functions were established for cork-rings δ¹³C, δ¹⁸O, width (C_RW) and maximum density (C_MD).

Cork ring δ¹³C and δ¹⁸O data collection

Cork samples were collected from four healthy trees at each site (CL and BS) during the cork harvesting seasons (June) of 2012 and 2013. Such cork samples are complex, multi-layered samples, encompassing reproduction cork-rings that were induced by multiple and consecutive cork harvests (Fig. 2). Thus, these samples allow addressing long-term chronologies of cork rings in individual trees (Costa et al., 2015; Costa et al., 2022). In each sample, cork rings were visually dated using image analysis techniques (Costa et al., 2015; Costa & Cherubini, 2021) for the period 1962–2013.

Generally, cork production cycles at the mountain site (BS) are longer (10 years) than at the peneplain site (CL) (9 years) (Costa et al., 2022). At BS, cork samples were harvested in 2012, all comprised five consecutive cork production cycles of 10 years, initiated in 1962 (first cork production cycle), 1972, 1982, 1992 and 2002 (fifth cork production cycle), with the final cork harvest in 2012. At CL, cork samples per tree were harvested in 2013, to which corresponded five consecutive cork production cycles of 9 years starting in 1968 (first), 1977, 1986, 1995 and 2004 (fifth), and harvested in 2012 to which corresponded five cork production cycles initiated in 1964 (first), 1973, 1982, 1991 and 2000 (exceptionally, the fifth cycle lasted 12 years) (Fig. 2).

Cork-ring width (C_RW) and maximum density (C_MD) chronologies were generated for each cork sample (Costa et al., 2022). Since C_RW is heavily influenced by cork harvesting, each cork-ring width series was detrended by fitting a polynomial function (Costa, Pereira & Oliveira, 2002). Cork-ring width indices (residuals) (IC_RW) were calculated by dividing the observed C_RW values by the predicted ones and IC_RW chronologies were then built. On the other hand, C_MD chronologies did not seem to contain any cork harvesting related trend (Costa et al., 2022) and therefore did not require any statistical detrending. Cork-ring maximum density indices (IC_MD) were calculated by dividing the observed C_MD values by their mean value in each cork production cycle and, this way the IC_MD chronologies were built.

Due to technical constrains to obtain enough cork mass for isotopic analysis, not all cork rings could be analyzed. In all samples, we selected the 1st and the 6th rings, respectively with the cork age (y_cork) equal to 1 and 6, of each of the considered five production cycles. The 1st ring is the largest one, representative of the trees response to cork harvesting and the 6th ring is a middle cycle ring, representative of an average cork growth, in the absence of stressful disturbances (Fig. 2). Data on carbon and oxygen isotopes were thus ascribed to both cork ages in all the consecutive cycles, regardless of the time series. By selecting fixed years (1st and 6th cork rings) of each cork-production cycle we expected to obtain information reflecting the intra- and inter-tree variability in tree physiological responses to drought conditions, under the influence of cork harvesting. A total of 78-records for (δ¹³C) and for (δ¹⁸O) were collected in the cork-rings of the study areas, 44 records for BS and 34 for CL.

Stable isotope ratio analyses were performed at the Stable Isotopes and Instrumental Analysis Facility (SIIAF) of the Centre for Ecology, Evolution and Environmental Changes (cE3c), University of Lisbon, Portugal. The stable-carbon isotope ratio (¹³C/¹²C) of cork was measured by continuous flow isotope mass spectrometry (CF-IRMS), on a Sercon Hydra 20-22 (Sercon, Hook, Hampshire, UK) stable isotope ratio mass spectrometer, coupled to a EuroEA (EuroVector, Italy) elemental analyzer for online sample preparation by Dumas-combustion. The stable-oxygen isotope ratio (¹⁸O/¹⁶O) of cork and suberin was also determined by continuous flow isotope mass spectrometry (CF-IRMS), coupled to a high temperature pyrolysis unit (HT-PyrOH, Eurovector, Italy) and a EuroEA3000 (EuroVector, Lombardy, Italy) elemental analyzer for online sample preparation by thermal degradation (thermolysis). The precision of the isotope ratio analysis, calculated using values from six to nine replicates of standard material interspersed among samples in every batch analysis, was ≤0.2‰. Results are given in the δ-notation, i.e., in the delta notation relative to the international standards: ${\displaystyle \delta \mathrm{sample}=\left(\mathrm{Rsample}/\mathrm{Rstandard}-1\right)\times \text{1,000},}$ where Rsample is the molar fraction of ¹³C/¹²C, or ¹⁸O/¹⁶O ratio of the sample and Rstandard, of the standards of Vienna Pee Dee Belemnite (VPDB) for carbon and Vienna Standard Mean Ocean Water (VSMOW) for oxygen.

Statistical analyses

The climate-isotope relationships were tested based on monthly Pearson correlations considering one month and five- and twelve-months windows, from February of the prior growth year till November of the current growth year. We included prior-year climate variables in the analysis because climate in the preceding growing season often influences stable isotope values in the following season due to carbon carry-over effects (Gessler et al., 2014).

Linear mixed-effects modelling was used to examine the influence of the covariate cork age (y_cork) at the two sites (non-water-limited area, CL and water-limited area, BS) on the isotope-derived physiological responses of cork oak and; to examine the variability of the cork δ¹³C and δ¹⁸O relationships with climate parameters, precipitation and temperature. Climate variables were thus considered as the fixed factors for the response functions δ¹³C and δ¹⁸O composition of cork rings. Also, linear mixed-effects models would explicitly integrate among-tree variability, considering trees as a random effect, related with tree responses to limited water availability. Model fitting was done by adding random effects of site Xy_cork (four classes) and tree (8 classes), through random intercept alone. We evaluated model performance using the Akaike information criterion (AIC). The model with the lowest AIC value was considered the best to predict response variables (δ¹³C and δ¹⁸O). We further evaluated the variation explained by the random and fixed effects together by calculating marginal R² to define the error or statistical uncertainty related with the inter-tree variability. All the analyses were conducted in R statistical environment (R Development Core Team, 2019) using the nlme package in R.

δ¹³C composition in cork rings and its correlations with climate

The mean isotopic compositions (δ¹³C) of corkrings ranged between −27.3‰ at CL and −26.8‰ at BS (Fig. 3). Average δ¹³C values at BS indicated that here cork was slightly more enriched in ¹³C (over 0.5‰) than at CL. However, there was no significant difference in cork-rings δ¹³C between the trees from BS and CL (nested ANOVA results showed no significant differences between the two groups (sites), F(1, 6) = 0.676, p-value = 0.442). Between trees within each study area, the variations were statistically significant (nested ANOVA results showed F(6, 70) = 6.951, p-value = 0.000). At CL, the cork of tree CL44 differed significantly from the other trees and; at BS, the cork of trees BS1 and BS8 differed significantly from the other two trees. CL44, BS1 and B8 were distinctly enriched in ¹³C (less negative δ¹³C values) when compared to the other trees (Fig. 3).

Furthermore, in all cork samples, except for the cork of tree CL44 (at CL), the mean isotopic compositions (δ¹³C) of 1st cork-rings (y_cork = 1) were less negative than 6th cork-rings (y_cork = 6) which indicates that, on average, cork formed immediately after cork harvesting is more enriched in ¹³C when compared to the cork formed later, hereafter the 6th year in the cork harvesting cycle.

The Pearson’s correlations (r) between δ¹³C of cork rings and mean temperature showed similarities among the trees of CL and BS, and were significantly negative (Figs. 4A–4B). This is a striking result which suggests that the δ¹³C (‰) of cork rings will be lower in warmer years, i.e., cork will be more ¹³C depleted.

Noticeably, besides a few significant monthly correlations (Fig. 4A), the relationships became more robust after grouping temperature variables in periods of 5- or 12-months, mainly at BS (Fig. 4B). The variation of cork δ¹³C at BS was much more sensitive to temperature than at CL, with a marked sensitivity to prior conditions since spring of the previous year throughout the growth year (Fig. 4B). The highest negative correlation coefficients were found for the mean temperatures previous to the growth year, in February–September (r =  − 0.423, p-value = 0.004) at BS and in April–May (r =  − 0.462, p-value = 0.006) at CL. Moreover, at BS, δ¹³C also correlated with temperature of the current growing season, in the spring period (March = −0.371, p-value = 0.013) (Figs. 4A–4B), in the onset of the cork growth flush.

Weaker correlations were found between δ¹³C of cork rings and precipitation in both study areas (Figs. 4C–4D). The precipitation of current-year March (at CL) and June (at BS) positively influenced cork δ¹³C (Fig. 4C). In contrast with temperature, no significant correlations were found with grouped precipitation values (periods of 5- to 12-months) (Fig. 4D). A relatively weak correlation was observed only at BS, in the summer of the current growth year in June–August (r = 0.350; p-value = 0.021). Also, in contrast with temperature, results revealed one irregularity of the relationship between cork δ¹³C and precipitation: although not statistically significant, correlations were either positive or negative correlations, regardless of the considered periods (Fig. 4D).

δ¹⁸O composition in cork rings and its correlations with climate

Cork-ring δ¹⁸O averaged 21.3 ± 0.9‰ (CL) and 20.9 ± 0.9‰ (BS). Average δ¹⁸O values indicate that at CL, cork was enriched in ¹⁸O, over 0.4‰ more than at BS (Fig. 5). However, there was no significant difference in cork-rings δ¹⁸O between the trees from BS and CL (nested ANOVA results showed no significant differences between the two groups (sites), F(1, 6) = 3.188, p-value = 0.121). Between trees within each study area, the variations were not statistically significant (nested ANOVA results showed F(6, 70) = 1.989, p-value = 0.079). Furthermore, mean values of suberin’s δ¹⁸O in cork rings ranged between 16.8‰ at CL and 16.3‰ at BS, in accordance with the values found for cork.

At CL, the cork of tree CL4 differed significantly from all other trees, and it was depleted in ¹⁸O when compared to the other trees (Fig. 5).

In all cork samples, except for the cork of tree CL49 (at CL), the mean isotopic compositions (δ¹⁸O) of 1st cork-rings (y_cork = 1) were slightly higher or equal than 6th cork-rings (y_cork = 6) which indicates that cork formed immediately after cork harvesting is somewhat more enriched in ¹⁸O when compared to the cork formed later, in the 6th year after the cork harvesting (Fig. 6).

The oxygen isotope signal in cork rings is more sensitive to precipitation than to temperature (Fig. 6). The correlations between δ¹⁸O and mean temperature were statistically significant only in trees at CL (Figs. 6A–6B). At this site, we found negative correlations with the temperature only in spring prior to the growth year. Besides the weak isolated monthly correlation with May (r =  − 0.345, p-value = 0.046) (Fig. 6A), this relationship became slightly more robust when after grouping temperature shows statistically significant correlations, in the summer prior to the growth year, between May–June (r =  − 0.373, p-value = 0.030) and May–August (r =  − 0.359, p-value = 0.037).

The correlations found between δ¹⁸O and precipitation were weak and negligible at BS but not at CL (Figs. 6C–6D). Here, δ¹⁸O relationship was positive and statistically significant with precipitation in the current spring, in March–May (r = 0.497, p-value = 0.003) and in spring-early summer, in March–June (r = 0.428, p-value = 0.012).

Correlations between width and density of cork-rings and climate

At the two study areas, none of the correlations between cork-ring isotopic compositions (δ¹³C and δ¹⁸O), cork-ring width (IC_RW) and cork-ring maximum density (IC_MD) were statistically significant.

The two study areas clearly differ in the relation between IC_RW and climate parameters. At BS, the results are inconsistent regarding precipitation: (i) within the year prior to the growth year, two distinct sub-periods: one in the spring (March–June), positively correlated with precipitation, and another in summer-autumn (June–October), negatively correlated and; (ii) in the current growth year, contrasting weaker positive correlations in summer-autumn (July–September) (Table 1). Here, despite the non-significant relation between IC_RW and mean temperature of the current growth year (in line with the non-significant relation with IC_MD, Table 1), the mean temperature within March and November prior to the growth year influenced cork growth much more than the cork density (Table 1).

At CL, an irregular pattern of correlations between ICRW and climate parameters was observed: negative correlations between precipitation and IC_RW, mainly in July–November of the current growth year, striking, strong and positive correlations between mean temperature in spring (April–May) of the year prior to the growth year, and negative correlation with mean temperature in the winter prior to the growth year (December–February) (Table 1).

Relationships between climate parameters (precipitation and temperature) and IC_MD were only relevant at BS; at CL the correlations were statistically not significant (Table 1). At former site, precipitation conditions prior to the growth year negatively affected IC_MD, from May to October. Conversely, the IC_MD correlations with precipitation of current growth year May and August were statistically significant and positive. Moreover, between temperature and IC_MD, only a relative weak significant correlation was found, with mean temperature of May prior to the growth year (r = 0.278, p-value = 0.046) (Table 1).

The correlations found between mean temperature and IC_RW and IC_MD at BS, and between precipitation and IC_RW and IC_MD at CL, were considerably weaker than the correlations found with δ¹³C (BS) and with δ¹⁸O (CL).

Influence of cork harvesting and site on the relation between climate and cork- rings isotope composition (δ¹³C and δ¹⁸O)

Selected mixed-effects linear models (M0–M5) were plotted against the site Xy_cork and trees random-effects (Table 2). The fixed-effect variables were the climate variables which displayed the highest bootstrap correlation coefficients. For δ¹³C were considered the mean temperatures of: February₍₋₁₎–June₍₋₁₎ (M1) (r =  − 0.355, p-value = 0.039 for CL and r =  − 0.397, p-value = 0.008 for BS); April₍₋₁₎–May₍₋₁₎ (M2) (r =  − 0.462, p-value = 0.006 for CL and r =  − 0.303, p-value = 0.013 for BS) and; March (M3) (r =  − 0.302, p-value = 0.082 for CL and r =  − 0.371, p-value = 0.013 for BS). For δ¹⁸O were considered the precipitation of: August₍₋₁₎ (M4) (r = 0.087, p-value = 0.624 for CL and r = 0.362, p-value = 0.016 for BS) and; March–June (M5) (r = 0.428, p-value = 0.012 for CL and r =  − 0.220, p-value = 0.157 for BS).

In all the models, δ¹³C in cork rings decreased with mean temperature and confirmed our previous results that δ¹³C of cork rings was lower in warmer years, i.e., ¹³C depleted (Figs. 4A–4B). In the best model for δ¹³C (M2—with the lowest AIC of 205.08 and explaining about 20% of the total variability), δ¹³C decreased on average about −0.267‰ per °C of monthly spring temperature (April₍₋₁₎–May₍₋₁₎) (Table 2). Despite the general negative correlation, in this model the predicted values of δ¹³C in the cork-rings formed immediately after cork harvesting (y_cork = 1) at a drier site (BS) were clearly higher (less negative δ¹³C), with an α₀ = −22.482, than in the later cork-rings (y_cork = 6) or at the CL wetter site (with either y_cork of 1 or 6), with an α₀ = −23.004, for the same negative slope (α₁ = −0.267, p < 0.001).

Model M20 for δ¹³C, with trees as random effect, had the lowest AIC, 181.12 and explained about 50% of the total variation (Table 2). This suggests a most strong influence of between-tree variability on the variations of cork-rings δ¹³C with mean temperature, mainly when compared to model M2 that explained a non-negligible 20% of total variation and that the site Xy_cork random effect explained up to 15% of total variability.

In the models with δ¹⁸O as the response function, δ¹⁸O increased with the monthly precipitation, in accordance with previous results on significant correlations (Figs. 6C–6D). The models M4 (with an AIC of 218.99 and explaining 6% of the total variation) and M5 (with an AIC of 220.29 and explaining 1% of the total variation), poorly explained the total variation. Furthermore, no variation was explained by the site Xy_cork random effect (Table 2). However, in the random intercept model (M5), the predicted values of δ¹⁸O in cork-rings formed at a drier site (BS) were clearly different and slightly depleted in ¹⁸O (lower δ¹⁸O, with an α₀ = 20.772), than in the wetter site (CL), with higher δ¹⁸O (α₀ = 21.039). These results would suggest the relative higher influence of site on the depth of the preferred source water for cork rings formation, independently on the effect of the cork harvesting.

Influence of water availability conditions on the δ¹³C of cork rings

Cork oak has the adaptive capacity to cope with variations in soil water availability during the growing season and can adjust its annual growth rhythm to prevailing environmental constrains, under the highly seasonal Mediterranean climate conditions (Cherubini et al., 2003; De Micco et al., 2016). Trees at Grândola (BS) typically present a bimodal growth pattern, as growth is subjected to a double constraint, by low temperatures in winter, and by drought in summer as in temperate climates (Cherubini et al., 2003). Generally, the onset of cork (and tree) growth occurs in spring, after a relative cold and wet winter, when soil water reserves have been replenished. In the summer, with limited water availability, trees experience a drought stress period and stop their stem (and cork) radial growth in July–August (summer rest period), resuming growth in autumn. However, trees at Benavente (CL) are able to reach the groundwater (Mendes et al., 2016), and to maintain their growth uninterruptedly from March to November; in fact, it is during the summer (June–July) that cork oaks’ radial growth (that it is mostly cork growth) reaches a maximum peak at this site (Costa, Pereira & Oliveira, 2003). So, different physiological growth patterns would be expected between the two sites.

Under low water availability (at BS), the higher cork rings δ¹³C, i.e., cork-rings enriched in ¹³C (on average, −26.8‰ at BS against −27.3‰ at CL) (Fig. 3) would be mostly explained by a strong regulation (reduction) of stomatal conductance, possibly affecting photosynthesis, confirming our first hypothesis on that the cork rings δ¹³C variation may reflect the physiological and structural adjustments (e.g., stomatal conductance) to cope with drought. The former values are slightly lower but in line with the range of values between −26.3‰ and −26.1‰ found in tree rings of the evergreen oak Q. ilex (Zalloni et al., 2018), where wood formation is determined by a double climate stress (summer drought and winter cold) and usually results from two phases of cambial activity separated by the summer rest period. Since cork mass is 40–50% suberin, a lipid-based polymer (Graça & Pereira, 2000; Graça, 2015), our results on ¹³C depleted cork-rings in comparison with tree-rings are in accordance with previous studies reporting that the metabolic pathway involved in lipid synthesis is much more ¹³C depleted than the one involved in (wood) cellulose synthesis. This occurs due to the ¹³C depletion of primary assimilates in the leaves to produce lipids, in contrast with ¹³C enrichment to produce cellulose (De Niro & Epstein, 1977; Melzer & Schmidt, 1987; Gessler & Ferrio, 2022). Moreover, our results on the correlations of cork rings δ¹³C with climate variables confirmed the critical role of temperature. This agrees with the fact that trees are forced to concentrate their photosynthetic activity in the warmer months of spring and autumn at both sites but mainly at the drier one (BS) (Fig. 4 and Table 2).

Noticeably, there is an apparent paradox in our results as trees at both sites: CL and BS, showed a general decreasing trend of δ¹³C in cork rings (i.e.,¹³C depleted cork rings) with increasing temperatures, and significant negative correlations between δ¹³C and mean temperature of the previous-year spring (Figs. 4A–4B). However, this depletion of ¹³C in cork rings can be explained by the possibility of the occurrence of post-photosynthetic processes associated with the remobilization and use of previous year(s) synthetized and stored carbon for cork formation, due to the shortage, or even the ¹³C discrimination of recent photo-assimilates in current (early) spring which are effectively ¹³C enriched (with higher δ¹³C) due to reduction of stomatal conductance. In cork oaks, a transfer of carbon reserves from the non-conducting to the conducting phloem tissue, close to the zone where phellogen is active, was reported by Knapic et al. (2007). These authors speculated that a displacement of these reserves in the form of extractives from heartwood to the outer sapwood for cork formation occurs at expenses of wood growth. It is possible that isotopic fractioning occurs during transport of metabolites As according to Cerasoli et al. (2004), the mobilization and utilization of stored carbon in cork oaks implies the hydrolysis of starch and the synthesis of sucrose, both enriched in ¹³C (Eglin et al., 2010) with a reduction of its concentration in stem wood. One can suggest that the preferential carbon reserves used for cork growth would be depleted in ¹³C, perhaps formed in less water stressed periods, and/or that the fractioning at the biochemical pathways of suberin synthesis would lead to compounds with very low δ¹³C (Damesin & Lelarge, 2003).

Trees at Benavente (CL) are with unlimited water availability and thus less water stressed. Cork growth should rely mainly on current photosynthates, relatively more depleted in ¹³C, in the onset of cork formation, in spring. Moreover, in autumn and warmer winters, trees use and probably store current photosynthates for the next season, which are less enriched in ¹³C (lower δ¹³C) given the lower tree water stress. In BS, relatively higher mean values of ¹³C in cork rings reflect limited water availability and trees water stress. Spring cork growth flush should rely not only on ¹³C enriched current photosynthates but on enriched ¹³C carbon reserves, which were stored during the previous dormant (winter) season (Ghalem et al., 2018), similarly to spring xylem (early wood) formed in tree-rings of other Mediterranean evergreen oaks, e.g., Q. ilex (Ferrio et al., 2003; Gea-izquierdo, Cherubini & Cañellas, 2011; Shestakova et al., 2014; Zalloni et al., 2018), and even in deciduous Q. petraea (Barbaroux & Bréda, 2002).

Furthermore, suberin is the main structural compound of the cork cell walls and it has a very special biochemical composition, a polyester lipid build up from long chain (C16–C24) fatty acids, which have carboxylic acid and/or hydroxyl functionalities at both chain ends (Graça, 2015). When phellogen is active, the phellem cells undergo a relative quick normal sequential differential stage: expansion, massive deposition of suberin layers in the internal side of cell walls, and programmed cell death (Soler et al., 2007; Teixeira & Pereira, 2010; Inácio et al., 2018). Such construction and integration of suberin in cork could lead to a ¹³C depletion and a concomitant enrichment of ¹³C respired CO₂, strongly and directly related with temperature, during cork formation as it is well known that respiratory processes discriminate against ¹³C.

Another possible explanation for the δ¹³C negative correlation with temperature is that, while the tight stomatal regulation is very responsive to water stress, as in a true drought-avoiding species, meaning that trees effectively display high stomatal resistance during drought, even reducing photosynthetic metabolism which may result in the growth impairment (Kurz-Besson et al., 2006), the rate of carboxylation and net photosynthesis would be much less sensitive or more complacent to water stress, which would decisively influence δ¹³C of carbohydrates, and lipids.

Overall, the particularly unexpected result on ¹³C isotopic signature in cork might reflect the dynamics of both processes, primary CO₂ fixation (stomatal regulation) and the downstream metabolic processes, carbon photosynthetic assimilation or post-assimilation of carbon reserves, including the metabolites used for export and correspondent fractionation processes, at different times and under different environmental conditions. These are some alternative hypotheses which can potentially explain the varying rates of incorporation of relatively heavy carbon (¹³C) in suberin during cork formation, and which merits future research.

The examination through mixed modelling of cork rings δ¹³C variations between trees (random intercept model, M2) might help to understand the range of varying rates of ¹³C in cork. According to our results, cork rings at the drier site (BS), in the first year after cork harvesting (y_cork = 1) were clearly different and had less negative δ¹³C (¹³C enriched) than other cork-rings (y_cork = 6, at BS; or y_cork = 1 or 6, at CL) much more ¹³C depleted, despite sharing common variation in response to temperature (Table 2, Fig. 3). Recently harvest trees in BS are under severe water stress due to water loss from the cork harvested stem (and branches) which much exceeds the water lost by leaf transpiration during the summer season (Correia et al., 1992; Oliveira & Costa, 2012). This suggests that under these conditions, trees would close stomata, with the consequent reduction of CO₂ concentration in the intercellular space, limiting CO₂ availability for the carboxylation, which is then forced to less discrimination and fix a relative greater proportion of the heavier isotope (¹³C) (Farquhar et al., 1989). Also, water stress trigger suberin biosynthesis and suberization (Aguado et al., 2012; Boher et al., 2018), similarly to the promotion of xylogenesis in tree rings (Castagneri et al., 2018), and trees need to urgently regenerate their cork layers to survive. Under these conditions, only in the first year after cork harvest (y_cork = 1) and only at the drier site (BS), trees likely used mixing carbon pools in their metabolic pathways of suberin synthesis: from current assimilates (already relatively enriched in ¹³C) and from carbon reserves, which will be relatively more enriched in ¹³C, producing ¹³C enriched cork rings.

Despite the absence of significant correlations between cork ring parameters width (IC_RW) and maximum density (IC_MD) and δ¹³C isotopic composition, the way that trees generate cork cells or increase their walls density in response to climate conditions could be different, depending on the way they respond to water stress conditions and control stomata and/or photosynthetic activity. The IC_RW and IC_MD correlations with mean temperature (Table 1) were considerably weaker than the correlations found for δ¹³C (Fig. 4). This suggests that cork is formed, or its density is increased as a result of an overall tree performance on their biomass production, while ¹³C signatures reflect tree adaptive mechanisms of response to water stress and drought.

Influence of water availability conditions on the δ¹⁸O of cork rings

Under a Mediterranean type climate, with water as the main limiting factor for plant growth, distinct cork δ¹⁸O should primarily reflect the variability in the water source used during growing season between trees. Furthermore, cork δ¹⁸O should also be responsive to the stomatal conductance, but not sensitive to photosynthetic activity, as the δ¹³C (Farquhar & Lloyd, 1993; Roden & Farquhar, 2012).

Cork rings δ¹⁸O revealed some differences between the two study areas. At drier site (BS) cork was more depleted in ¹⁸O than at the wetter site (CL), with higher δ¹⁸O (Fig. 5). Here, these δ¹⁸O higher values indicate mostly the use of a shallower, more isotopically enriched water source, but also probably some (leaf) water ¹⁸O enrichment, given the unlimited water availability suggesting a higher stomatal conductance (Fig. 5). Moreover, the spring March-May precipitation only affected the δ¹⁸O of cork rings at wetter site (CL) (Fig. 6C) suggesting some dependence on growing season precipitation for cork formation in the long and uninterrupted cork growth period (Costa, Pereira & Oliveira, 2003). Thus, relatively more enriched ¹⁸O cork-rings are expected.

In contrast, at the drier site (BS), cork-rings are more ¹⁸O depleted (Fig. 5). Moreover, no significant correlations were found with precipitation. These results might suggest that cork oak obtain their water from deeper soil horizons (Sarris, Christodoulakis & Körner, 2007), trees are less dependent on growing season precipitation (or even precipitation in winter previous to growth year) for cork formation and have a strong reliance on storage carbon reserves (Ferrio et al., 2015; Szymczak et al., 2019). Other possible explanation for the non-significant correlations with precipitation is that water reserves in deeper groundwater can represent a long-term mean isotopic value of winter precipitation and hence do not only reflect the amount of winter precipitation prior to the growth season (Szymczak et al., 2019). Furthermore, in contrast to CL, trees at BS have a stronger stomatal conductance regulation that probably would constrain the (leaf) water ¹⁸O enrichment.

Cork oak is a deep-rooted species, with a very effective water uptake from the seasonally recharged groundwater, and with a dual root system which allow to use of soil water, even in the drier summer season (Vaz et al., 2010; David et al., 2007). Moreover, the preferential water source of cork oak is deeper groundwater sources (Altieri et al., 2015). However, at BS, a direct comparison between δ¹⁸O of cork rings and seasonal fluctuations of precipitation can be hampered in several ways. First, seasonal cork (suberin) synthesis is decoupled from seasonal precipitation and distinct information can be encoded in the ¹⁸O of the cork rings; second, tree uses mixing pools of carbon which are formed with distinct water pools, for suberin synthesis in distinct growth seasons; third, depending on rooting depth, trees have to access to different water pools with varying isotopic signatures and it can be assumed that groundwater replenishment in winter previous to the growth season cannot be sufficient and trees depend also on precipitation during the growth season as other main water source.

Altogether, mixed modelling results points to the existence of a relationship between cork δ¹⁸O signatures and trees preferential water sources, confirming our second hypothesis, without the complexity of the effect of climate drivers and of the cork harvesting on the variability of δ¹⁸O in cork rings (Fig. 5; Table 2). However, further investigations are needed in order to clarify the interpretation of our results on the δ¹⁸O signature in cork rings.