Deep Soil Conditions Make Mediterranean Cork Oak Stem Growth Vulnerable to Autumnal Rainfall Decline in Tunisia

: Tree rings provide fruitful information on climate features driving annual forest growth through statistical correlations between annual tree growth and climate features. Indices built upon tree growth limitation by carbon sequestration (source hypothesis) or drought-driven cambial phenology (sink hypothesis) can be used to better identify underlying processes. We used both analytical frameworks on Quercus suber , a sparsely studied species due to tree ring methodological issues, and growing on a favorable sub-humid Mediterranean climate and deep soil conditions in Tunisia (North Africa). Statistical analysis revealed the major role of autumnal rainfall before the growing season on annual tree growth over the 1918–2008 time series. Using a water budget model, we were able to explain the critical role of the deep soil water reﬁll during the wet season in affecting both the drought onset controlling growth phenology and the summer drought intensity affecting carbon assimilation. Analysis of recent climate changes in the region additionally illustrated an increase in temperatures enhancing the evaporative demand and advancing growth start, and a decline in rainfalls in autumn, two key variables driving stem growth. We concluded on the beneﬁts of using process-based indices in dendrochronological analysis and identiﬁed the main vulnerability of this Mediterranean forest to autumnal rainfall decline, a peculiar aspect of climate change under summer-dry climates.


Introduction
Tree rings are valuable climate archives that bear witness to the yearly climate conditions that effectively drive trunk expansion over the span of decades to centuries [1]. They have a seasonal resolution, are absolutely dated, and are available in a wide range of climatic settings [2]. They also provide fruitful information on climate features driving forest productivity, considered both as a resource and an ecosystem service through carbon sequestration, and allow for future assessment of climate change impacts and atmospheric feedbacks. However, statistical analyses relating annual tree growth to key climate variables provide poor information on ecophysiological processes driving tree growth. Recent advances balance the usually assumed hypothesis that tree growth is limited by carbon assimilation (source limitation) used in process-based models [3][4][5], by providing evidence of a stronger control by cambial activity (called the sink limitation) and its seasonal phenology [6][7][8][9][10]. The sink limitation hypothesis for Mediterranean climate states that annual tree growth is limited by phenological phases of cambial activity controlled by the dual role of cold temperatures in the winter and the onset/offset of drought [8], while the source limitation hypothesis states that stem growth is limited by carbon assimilation, highly affected by summer drought under Mediterranean climatic conditions. In turn, regional variations in soil field capacity within a homogeneous climate region that locally affect soil water balance and drought features [11] can lead to contrasted results in the statistical relationship between seasonal climate variable and tree rings.
Within this contextual framework, we explored the long-term tree growth/climate relationship for cork Oak (Quercus suber L.) by a dendrochronological study in Tunisia (on the North African southern bound of the Mediterranean basin) where this species is one of the most important forest species from an ecological, social, and economic value (covering 91,767 ha [12] among the 2.2 Mha [13] in the whole Mediterranean basin). Technical issues to perform such a study rely on the strict restriction on the felling of living trees and the difficulty of wood ring measurements on debarked trees for this species. The longest Quercus suber chronology in the literature is 25 years long [14] and contains a strong common signal using 14 trees never submitted to cork harvesting, and few dendrochronological studies managed to analyze tree ring growth on Quercus species in Tunisia [15][16][17][18]. Dendrochronological analyses are constrained by the difficulties in the identification of tree ring boundaries, which may lower the common growth pattern necessary for successful cross-dating and growth/climate comparisons [17].
We selected one of the most productive forests in Tunisia (Bellif) located on deep and almost stone-free soils [19,20] to get the thickest tree rings, and to find trees that are representative of the most productive forests used for their economic value in Tunisia, in an area long-ago identified as the fertile "granary of Rome" within the Mediterranean basin, but recently depicted as perishing [21]. The objectives of the study were: (i) to establish a statistically significant long-term tree-ring series on Quercus suber despite technical challenges on this species; (ii) identify climate control on tree growth by statistical analyses and test the process-based hypothesis of the source-and sink-limitation, in order to (iii) capture key aspects of climate changes to account for in the assessment of the sustainability of this ecosystem in the future.

Study Site
The study site is located at the Bellif national forest ( Figure 1) located in the north-western part of Tunisia (Nefza forest subdivision) on the North African southern bound of the Mediterranean basin where regional [22] and local scale [20,23] information on Quercus suber functioning have been performed. Two sites were chosen within a pure and natural stand of Quercus suber used for cork production and occasional sheep grazing. The two sites were approximately 3 km away from each other with the same climatic conditions and the same soil type but differing in disturbance history. The young site (YS, 37 • 2 41.7 N, 9 • 6 8.04 E) is a site with even-aged trees naturally regenerated after a tornado (1974). The mature site (MS, 37 • 2 12.44 N, 9 • 3 44.25 E) is a site without disturbance dominated by uneven-aged trees [20]. The single parent material of the forest soil is an oligocene sandstone interspersed with clay layers, yielding brown soils with strong biological activity [24].
The climate in the Bellif region is typically Mediterranean, characterized by hot and dry summers and mild winters ( Figure 2). Precipitations mainly occur from September to March and a drought period extended from May to August. The average annual precipitation is 1113 mm and the average annual temperature 19.3 °C, on the combined wettest and warmest side of the Mediterranean climate leading to high forest productivity [20]. Maximum and minimum air temperatures were 47.2 °C and 0.1 °C, respectively.

Tree-Ring Analysis
In summer 2009, 16 trees under cork production were cut and two cross sections were sampled at a height of 0.30 m and 1.30 m aboveground for dendrochronological analyses. Sampling trees were subjectively selected to represent the actual diameter class distribution on each site. Cross-sections were air dried in the laboratory and polished (40 to 500 grit). For the tree-ring width data, samples were processed using standard procedures of tree-ring analysis [25] and evaluated by "cross-dating" techniques [26]. Rings in each sample were dated to the calendar year of their formation. Tree-ring series of each dated sample were measured using a microcomputer with 0.01 mm accuracy (Lintab™-Rinntech™). After cross-dating and measuring all the samples, one tree (tree number 13) was The climate in the Bellif region is typically Mediterranean, characterized by hot and dry summers and mild winters ( Figure 2). Precipitations mainly occur from September to March and a drought period extended from May to August. The average annual precipitation is 1113 mm and the average annual temperature 19.3 • C, on the combined wettest and warmest side of the Mediterranean climate leading to high forest productivity [20]. Maximum and minimum air temperatures were 47.2 • C and 0.1 • C, respectively. The climate in the Bellif region is typically Mediterranean, characterized by hot and dry summers and mild winters ( Figure 2). Precipitations mainly occur from September to March and a drought period extended from May to August. The average annual precipitation is 1113 mm and the average annual temperature 19.3 °C, on the combined wettest and warmest side of the Mediterranean climate leading to high forest productivity [20]. Maximum and minimum air temperatures were 47.2 °C and 0.1 °C, respectively.

Tree-Ring Analysis
In summer 2009, 16 trees under cork production were cut and two cross sections were sampled at a height of 0.30 m and 1.30 m aboveground for dendrochronological analyses. Sampling trees were subjectively selected to represent the actual diameter class distribution on each site. Cross-sections were air dried in the laboratory and polished (40 to 500 grit). For the tree-ring width data, samples were processed using standard procedures of tree-ring analysis [25] and evaluated by "cross-dating" techniques [26]. Rings in each sample were dated to the calendar year of their formation. Tree-ring series of each dated sample were measured using a microcomputer with 0.01 mm accuracy (Lintab™-Rinntech™). After cross-dating and measuring all the samples, one tree (tree number 13) was

Tree-Ring Analysis
In summer 2009, 16 trees under cork production were cut and two cross sections were sampled at a height of 0.30 m and 1.30 m aboveground for dendrochronological analyses. Sampling trees were subjectively selected to represent the actual diameter class distribution on each site. Cross-sections were air dried in the laboratory and polished (40 to 500 grit). For the tree-ring width data, samples were processed using standard procedures of tree-ring analysis [25] and evaluated by "cross-dating" techniques [26]. Rings in each sample were dated to the calendar year of their formation. Tree-ring series of each dated sample were measured using a microcomputer with 0.01 mm accuracy (Lintab™-Rinntech™). After cross-dating and measuring all the samples, one tree (tree number 13) was eliminated from the study due to inconclusive dating. Thus, the number of sampled trees finally included in the analysis was 15 trees (10 from YS and 5 from MS). In total, 30 radii from 15 trees were successfully cross-dated (two radii per tree).
The individual ring width series were standardized by removing low and intermediate frequency variations using the appropriate ARMA (autoregressive moving-average) procedure with 3Pbase software [27]. Each ring-width value was replaced by a residual that represents only the unpredictable part of tree ring variation [16]. Standardization is necessary for the elimination of tree ring width decrease due to tree aging and in order to reinforce the climate signal in the tree rings [25,28].

Climate Data
Monthly climate data (maximum temperature (TMAX), minimum temperature (TMIN), and precipitation (P)) covering the whole dendrochronological time series  were derived from the updated CRU TS v 3.23 [29].
Daily climate data (Tmin, Tmax, and P) from the closest meteorological station (Nefza, Institut National de la Météorologie, Tunisia) were only available for the 1984-2008 period and were further used for the daily water budget model. Daily solar radiation was calculated from the following processing chain [30]: the theoretical clear sky solar radiation (R a ) was calculated from the daily timing of sunrise and sunset, the bi-hourly sun azimuthal angle, and the atmospheric transmittance according to longitude using R CRAN (v 3.3.1, R core Team, Vienna, Austria) packages "RAtmosphere" and "oce". Actual solar radiation (R s in MJ·m −2 ·day −1 ) was then calculated from maximum and minimum temperature and R a following the Hargreaves equation (Equation (1)): where the adjustment coefficient k was set to 0.18 [31]. Results of the Hargreaves equation were validated against solar radiation measured in the neighboring study site of Souk el Jema [32].

Water Budget Model and Drought Features
Soil water storage integrated over the 4.5 m rooting depth, was simulated at a daily time-step with a soil water balance model developed for Mediterranean ecosystems [33,34] and further used in the studied region [22,32]. The drainage curve relating deep drainage to soil water storage depends on the stone content over the whole-soil profile. The model was driven by daily values of incoming solar radiation, mean temperature, and rain amount. Potential evapotranspiration was computed using the Priestley-Taylor equation [35]. Soil water storage and soil water potentials were related by a Campbell-type retention curve [36], whose parameters are strongly dependent on soil texture [37], to deduce the daily plant predawn water potential (see Annex A for model description and parameters). Values of leaf predawn water potential used for model validation were measured during the 2008 dry season on a sample of six trees (pressure chamber PMS instrument Co. model 1000). The daily simulations of predawn water potential were performed with the climatic data from the Nefza meteorological station for the periods 1984-2008. Annual or seasonal water deficit was quantified using the water stress integral (WSI), defined by Myers [38] as the yearly sum of simulated daily predawn water potentials, and proposed as a generic drought index explaining annual carbon assimilation in water-limited Mediterranean ecosystems [39,40]. For each year we estimated the dates of stem growth phenology that bounds the growing period: the DOY (day of the year) when stem growth starts (hereafter t0), the DOY when stem growth stops in early summer (hereafter t1), and the DOY when growth restarts after the dry period (hereafter t2). In the absence of an automated dendrometer on the site to evaluate the start of the growing season, we used the relationship between t0 and the January to March sum of temperatures T JFM found for Quercus ilex (Equation (2)) (t0 = 849.2 × exp (−0.6436 × TJFM) + 121) [8] ( as the species exhibit close seasonal wood phenology with Quercus suber [41]. t1 was predicted by the DOY when the simulated plant water potential reached a threshold of −1 MPa, a generic critical threshold for cambial elongation for this species [8]. t2 corresponds to the DOY when predawn water potential reaches −1 MPa after the dry period, and allows for secondary growth as previously observed on twig elongation in the study site [20].

Statistical Analysis
Response functions were conducted using a bootstrap method [42]. Bootstrapped response functions were calculated, for each site, using the residual values as a dependent variable and the 24 monthly climatic parameters as explanatory variables; precipitation and temperature of the 12 month period from October of the year prior to the tree ring formation (t − 1) to September of the year of the tree ring formation (t). This period is the most widely used for dendroclimatological and dendroecological purposes in the Mediterranean region [8,43,44].
The bootstrap process consists of calculating the regression on years drawn by lots, some 1000 simulations being involved in total. For each iteration, the regression coefficients and the multiple correlation of climatic variables are computed on randomly selected years (calibration years). An independent verification is done on the observations omitted from the subsample (verifications years). The significance of each bootstrapped regression coefficient is provided by the ratio between the mean value (MR) calculated from the results of these 1000 simulations and its standard deviation (S). When the ratio ranges from 1.65-1.95, 1.96-2.57, 2.58-3.29, and >3.29, the significance of the corresponding regression coefficients respectively reach 90%, 95%, 99%, and 99.9% [16]. The global significance of the response function is defined by the ratio between (MR) and (S) according to the same probability levels [45]. MRV and SV then correspond to MR and S of the verification years, and MRC and SC to the calibration years. Owing to the low number of samples, climatic variables were gathered to have significant response functions [46][47][48].
We assessed the links between tree growth and phenological variables by testing correlations with the "treeclim" R package built upon the DENDROCLIM2002 statistical tool dedicated to tree ring analysis [49,50]. Pearson's correlation coefficients and their significance were computed based on the standard bootstrap method with 1000 samples taken from the original distribution of phenological variables and tree ring data.
Temporal trends in time series were estimated using the Theil-Sen test [51] within the "openair" R package. The Theil-Sen test calculates the temporal slope for time series and is commonly associated with the Mann-Kendall test, this last test detecting only the significance of the trend without slope estimation. To ensure the consistency between the calculation of the slope and its significance (p-value), and to account for temporal autocorrelation, both forms of information were derived from a 1000 block bootstrap simulation, randomly evaluating the trend slopes on sampling one third of the length of the whole time series, and offering a quantification of uncertainty and significance level [52].

Cross-Dating and Chronologies
Our samples did not show any missing ring, but we faced extremely narrow rings at the external area of the oldest trees. To solve this problem, cross sections were obtained at stump height (basal sections) and at 1.30 m to increase the sample size with each tree. Our attempts to measure and cross-date the complete ring width series of 15 out of the 16 trees under cork production were successful. Our measurements were performed by avoiding the mechanical injury caused by cork removal. These chronologies were validated by the historical management data (cork removal).
The time span of our chronologies was 1974-2008 for YS and 1918-2008 for MS, with an average ring width of 3.71 mm and 2 mm, respectively. Wood ring widths showed a decreasing age trend for both sites ( Figure 3A,B). They ranged between a minimum of 1.13 mm and a maximum of 8.01mm for the 34 year-old trees in YS and between 0.78 mm and 3.30 mm, respectively, for the 71-102 year-old trees in MS. The quality of our chronology in the MS was not always stable and shows a lower statistical quality during the earlier period (juvenile phase of trees), therefore, we kept only the 1937-2008 period for the response function.

Cross-Dating and Chronologies
Our samples did not show any missing ring, but we faced extremely narrow rings at the external area of the oldest trees. To solve this problem, cross sections were obtained at stump height (basal sections) and at 1.30 m to increase the sample size with each tree. Our attempts to measure and crossdate the complete ring width series of 15 out of the 16 trees under cork production were successful. Our measurements were performed by avoiding the mechanical injury caused by cork removal. These chronologies were validated by the historical management data (cork removal).
The time span of our chronologies was 1974-2008 for YS and 1918-2008 for MS, with an average ring width of 3.71 mm and 2 mm, respectively. Wood ring widths showed a decreasing age trend for both sites ( Figure 3A,B). They ranged between a minimum of 1.13 mm and a maximum of 8.01mm for the 34 year-old trees in YS and between 0.78 mm and 3.30 mm, respectively, for the 71-102 yearold trees in MS. The quality of our chronology in the MS was not always stable and shows a lower statistical quality during the earlier period (juvenile phase of trees), therefore, we kept only the 1937-2008 period for the response function.

Climate-Growth Relationship
The correlation between annual wood ring growth and monthly precipitations was neither high nor statistically significant for individual months. However, when considering precipitation accumulated during some selected periods (Table 1), the correlation with tree ring width was high and significant for some cases, namely the periods before the start of the growing season at both sites. In YS, tree-ring width and precipitation from October of the previous year (October t − 1) to February of the current year (February t) revealed a strong positive correlation. However, in MS, the highest and most statistically significant correlation coefficients were obtained for the cumulated precipitation for the periods including autumn and early winter rainfall of the previous year (from October t − 1 to December t − 1). Positive correlations were obtained with March precipitations (95% significance) in YS and with May precipitations (99% significance) in MS. We also obtained a negative correlation between wood growth and precipitation levels of the current summer in MS.  Regarding the effect of temperature, the global response function was not statistically significant for both sites (Table 2). However, significant positive relationships between tree ring width and the minimum and maximum temperature of September (t) were observed in YS. In general, between October and February, there was a negative correlation between wood ring growth and temperature. This inverse relationship was more pronounced in MS than in YS. For MS, minimum and maximum temperatures from the previous December (t − 1) to February (t) were negatively correlated with growth. The latter response functions were significant at the 90% threshold. In both sites, minimum and maximum temperatures of the growing season, with the exception of the hot season, exhibited a positive influence on wood growth but without statistical significance. The joint effects of precipitation and temperature were analyzed (further referenced as P-TMAX for the joint effect of Precipitation P and maximum temperature TMAX, and P-TMIN for the joint effect of precipitation P and minimum temperature TMIN) in order to emphasize the combined response functions. For YS, the high autocorrelation regression coefficients of P-TMAX (RMV/SV P-TMAX (ON) = 3.93 for the period October-November (t-1) (ON) and RMV/SV P-TMAX (DJF) = 2.30 for the period December (t-1) to February (t)) (DJF) and P-TMIN (RMV/SV P-TMIN (ON) = 3.70 fort he period October-November (t-1) (ON) and RMV/SV P-TMIN (DJF) = 2.29 for the period December (t-1) to February (t)) (DJF) suggested that cumulated precipitations before the growing season from October (t − 1) to February (t) have a prevailing role on tree growth under minimum and maximum temperatures ( Figure 4). For MS, Figure 5 showed a positive correlation coefficient for P-TMAX for the period extending from October (t − 1) to December (t − 1) (OND) (RMV/SV P-TMAX (OND) = 4.45) and P-TMIN (RMV/SV P-TMIN (OND) = 4.44). During the summer period (June-August), tree ring width was negatively affected by the combined effect of precipitation and temperature in MS ( Figure 5).

Soil Water Balance and Wood Phenology-Growth Relationship
We performed a daily soil water balance simulation and the corresponding simulated daily leaf predawn water potentials over the 1985-2008 period for which daily climate variables were available. Figure 6 illustrates the simulated daily predawn water potential (MPa) for the year 2008, during which predawn water potentials were actually measured (black dots). Predawn water potentials varied between −0.5 MPa in winter when soil water content was at field capacity and −2.85 MPa at the end of the dry season in late summer/early autumn.
The 1985-2008 simulated soil water content is presented in Figure 7A for a field capacity of 890 mm, corresponding to the actual soil conditions on the study site (see Annex 1 for site characteristics). We observe the typical recurrent seasonal variation of soil water content between field capacity in winter and close to the wilting point in summer. However, we can notice peculiar years (1990,1995,1997,2000,2002) where autumnal rainfall does not refill the whole soil profile, and winter soil water content does not reach field capacity.

Soil Water Balance and Wood Phenology-Growth Relationship
We performed a daily soil water balance simulation and the corresponding simulated daily leaf predawn water potentials over the 1985-2008 period for which daily climate variables were available. Figure 6 illustrates the simulated daily predawn water potential (MPa) for the year 2008, during which predawn water potentials were actually measured (black dots). Predawn water potentials varied between −0.5 MPa in winter when soil water content was at field capacity and −2.85 MPa at the end of the dry season in late summer/early autumn.
The 1985-2008 simulated soil water content is presented in Figure 7A for a field capacity of 890 mm, corresponding to the actual soil conditions on the study site (see Annex 1 for site characteristics). We observe the typical recurrent seasonal variation of soil water content between field capacity in winter and close to the wilting point in summer. However, we can notice peculiar years (1990,1995,1997,2000,2002) where autumnal rainfall does not refill the whole soil profile, and winter soil water content does not reach field capacity.     Correlations between annual tree growth and phenological features are presented in Table 3. Annual tree growth is significantly negatively correlated with WSI with a Pearson's correlation coefficient = −0.467 for MS and −0.428 for YS (p-value < 0.05). Correlations with t0 driven by winter temperature are either uncorrelated or weakly correlated (not significant for MS and p < 0.1 for YS) with negative correlation coefficients. t1 appears as a significant explanatory variable for MS (p < 0.05) but not for YS. The length of the spring growing season t1 − t0 appeared as a better indicator of tree growth and is significantly positively correlated to tree growth for both sites (p < 0.1). Better correlations were even obtained when accounting also for autumnal growth restart t2, either through the drought duration t2 − t1 (negative correlations with p < 0.05) for both sites, or for the whole Correlations between annual tree growth and phenological features are presented in Table 3. Annual tree growth is significantly negatively correlated with WSI with a Pearson's correlation Correlations with t0 driven by winter temperature are either uncorrelated or weakly correlated (not significant for MS and p < 0.1 for YS) with negative correlation coefficients. t1 appears as a significant explanatory variable for MS (p < 0.05) but not for YS. The length of the spring growing season t1 − t0 appeared as a better indicator of tree growth and is significantly positively correlated to tree growth for both sites (p < 0.1). Better correlations were even obtained when accounting also for autumnal growth restart t2, either through the drought duration t2 − t1 (negative correlations with p < 0.05) for both sites, or for the whole growing season additionally accounting for t0t1 and autumnal growth restart (p < 0.1 for MS and p < 0.01 for YS), indicating the significant role of the secondary growth in the annual wood ring increment. Table 3. Pearson's correlations (coefficient correlation (Corr) and p-value) between annual tree ring width and annual water stress integral WSI, start of the growing season (t0), end of the growing season (t1), length of the spring growing season (t1 − t0), length of the whole year growing including autumnal secondary growth (t1 − t0 + t2), and length of the summer growth break (t2 − t1) for the mature stand (MS) and the young stand (YS). NS: non-significant, *: significant at p < 0.1, **: significant at p < 0.05, ***: significant at p < 0.01. To better understand the relationship between wood phenological phases and climate drivers, we correlated t1 with the significant seasonal precipitation amounts tested in the climate/growth relationship (Table 4). We observed that t1 is highly correlated to both autumnal rainfalls for October-November P_ON and spring rainfalls for March-April (P_MA) and Avril-May (P_AM) for a field capacity fixed at 890 mm (p-value < 0.01, Table 4), illustrating a combined control of autumnal and spring rainfall on drought onset under deep soil conditions, in accordance with the direct climate/tree growth relationship. We compared this result with t1 obtained for simulations performed with field capacity = 220 mm (shallow soil hypothesis) to test whether the deep soil value used in the model (field capacity at 890 mm) is actually a key contributing parameter in these relationships. Table 4 (FC = 220 mm) then presents the correlations between seasonal precipitations and the simulated drought onset t1 under shallow soil conditions. We observe in these conditions that t1 is only correlated to spring precipitations P_MA and P_AM (p-value < 0.01), so that the contribution of P_ON on t1 and annual tree growth is only observed for deep soil conditions. Similar results regarding the effects of soil depth on the seasonal contribution of precipitation to WSI were observed, with a significant effect of P_ON on WSI for deep soil conditions becoming non-significant when considering a shallow soil. In addition, for WSI, only summer precipitations from June to August (P_JJA) remained significant for shallow soil conditions. Figure 7B presents the daily soil water content for the 1985-2008 period under shallow soil conditions, and illustrates the systematic yearly refill of soil water content to field capacity during autumn and winter so that drought onset is no more dependent on this refill compared to deep soil conditions. Table 4. Correlation coefficients and p-value between drought onset (t1) and WSI (simulated from a water budget model with field capacity FC = 890 mm and FC = 220 m) with the sum of precipitations for selected months. NS: non-significant, *: significant at p < 0.1, **: significant at p < 0.05, ***: significant at p < 0.01.

Impact of Recent Climate Changes on Key Drivers of Tree Growth
We finally tested how recent climate changes may have affected the most significant key drivers of tree growth obtained in the climate/tree growth relationship. The slopes and their p-values of the temporal trends (Theil Sen test) over the period 1985-2008 are presented in Table 5, and significant trends are presented in Figure 9. We observed for the region a high temperature increase for summer months (slope = +0.104 • C·year −1 , p-value < 0.01) and weaker temperature increase for winter months (slope = +0.027 • C·year −1 , p-value = 0.06). No significant trend for WSI was observed. We also observed a significant trend toward earlier drought onset t1 (slope = −4.82 days·year −1 , p-value = 0.1). As we showed that t1 is due to the previous fall and current rainfall component of the soil water budget, we tested for temporal trends on P_ON and spring precipitations, which were only significant for P_ON (−4.82 mm·year −1 , p-value = 0.023). At the same time, the temporal trend for growth start t0 was also observed to be significantly decreasing (slope = −0.41, p-value = 0.06) as a result of winter temperature increase, overall leading to a non-significant trend in t0t1 as both t1 and t0 decrease with the same amplitude. In turn, we conclude that the role of precipitation from the fall period before the growing season in driving drought onset and tree growth observed in the climate/tree growth relationship for Quercus suber growing on deep soil in Tunisia is currently the most affected by recent climate changes. However, due to compensation processes, the growth period t0t1 and the yearly drought intensity WSI do not exert any temporal change. Table 5. Temporal trend in climate variables (temperatures for June, July, and August (T_JJA), October, November, December, and January (T_ONDJ)), precipitations for October and November (P_ON), March and April (P_MA), April and May (P_AM), and June, July, and August (P_JJA), and the simulated (field capacity 890 mm) wood phenology start date (t0), drought onset (t1), drought offset (t2), length of the wood growth period (t1 − t0), and the annual water stress integral (WSI) for the 1985-2008 period obtained from a bootstrapped Theil Sen test. NS: non-significant, *: significant at p < 0.1, **: significant at p < 0.05, ***: significant at p < 0.01.   Table 5. Temporal trend in climate variables (temperatures for June, July, and August (T_JJA), October, November, December, and January (T_ONDJ)), precipitations for October and November (P_ON), March and April (P_MA), April and May (P_AM), and June, July, and August (P_JJA), and the simulated (field capacity 890 mm) wood phenology start date (t0), drought onset (t1), drought offset (t2), length of the wood growth period (t1 -t0), and the annual water stress integral (WSI) for the 1985-2008 period obtained from a bootstrapped Theil Sen test. NS: non-significant, *: significant at p < 0.1, **: significant at p < 0.05, ***: significant at p < 0.01.

Challenges in Cork Oak Tree Ring Analysis
Although the reading of wood rings on evergreen oaks (Quercus ilex and Quercus suber, the predominant tree species in Mediterranean forests [53]) presents methodological issues under cork extraction management [54], we were able to build a long (90 year) chronology for Quercus suber. Vessel arrangement is usually altered after cork removal, leading to irregular semi-ring porosity, and faint growth ring boundaries [55]. Tree ring distinction and dating is therefore easier before cork extraction [54][55][56][57]. In the present study, the annual ring width ranged between 1.13 mm and 8.01 mm

Challenges in Cork Oak Tree Ring Analysis
Although the reading of wood rings on evergreen oaks (Quercus ilex and Quercus suber, the predominant tree species in Mediterranean forests [53]) presents methodological issues under cork extraction management [54], we were able to build a long (90 year) chronology for Quercus suber. Vessel arrangement is usually altered after cork removal, leading to irregular semi-ring porosity, and faint growth ring boundaries [55]. Tree ring distinction and dating is therefore easier before cork extraction [54][55][56][57]. In the present study, the annual ring width ranged between 1.13 mm and 8.01 mm for young trees (YS) and between 0.78 mm and 3.30 mm for mature trees (MS). These values were similar to the few studies previously performed for similar ecosystems, and reporting wood ring increments of 2 mm·year −1 [58] and 2.61 mm·year −1 [14] for young trees as well as values ranging from 1 to 4 mm·year −1 [54] and from 0.68 to 3.20 mm·year −1 [14] for trees under cork production.

A Weak Role of Temperature on Cork Oak Growth
From the tree ring time series and the monthly climate variables, the global response function was always statistically non-significant for temperatures, although some conclusions may arise from the monthly analysis (the ratio between monthly partial autoregression coefficient and its standard deviation). Our observed negative correlation between October (t − 1) temperature and tree growth was already reported [59,60], but is in disagreement with the results on other deciduous Quercus species which have shown a positive effect of October (t − 1) temperature on oak growth for sessile oaks [61], red oaks [62], and pedunculate oaks [63].
We can hypothesize that higher September temperatures lead to an extended period of physiological activity at the end of the growing season, which allows for more time to complete relocation of transportable assimilates from leaves to perennial parts of the tree [63]. Moreover, a negative correlation was observed in both sites between the wood ring width and the temperature from December (t − 1) to February (t). This correlation was more pronounced (statistically significant) in trees in the MS than in the YS. The negative correlations seem to confirm that high temperatures during winter, when precipitations are higher, are not favorable to diameter growth, as already reported for cork growth [59,64], and may result from the respiratory loss of carbohydrate reserves at high temperatures [65]. In addition, increasing potential evapotranspiration (PET) during winter might enhance soil water loss and lead to earlier drought start.
Finally, our results showed that tree ring growth was enhanced by an increase in temperature from the onset of the physiological activity until the beginning of autumn (September) but it was repressed by the high temperature in June-August of the current year, which is typical for trees growing in a Mediterranean climate. Similar results were previously reported for cork oak growth [59,60,66]. We might hypothesize here an enhanced potential evapotranspiration rate during summers with high temperatures leading to a belated growth restart in the fall, as the offset of the dry season (t2) was identified as a significant driver of annual tree growth in our system.

A Dual Role of Previous Autumnal and Current Spring Precipitation on Tree Growth
Our statistical analysis on monthly climate variables showed that inter-annual variations on tree ring width were mostly affected by precipitations compared to temperatures. In YS, this period extended from October (t − 1) to February (t) and from October (t − 1) to December (t − 1) for MS. There was no significant overall correlation between precipitation in January and February precipitations and wood ring width. This result actually agrees with studies on Quercus species in neighboring Quercus suber ecosystems of Algeria with similar soil and climate conditions [47].
The direct relation between annual tree growth and the autumnal precipitations of the previous year could be explained by the benefits from these favorable conditions to synthesize and store carbohydrates reserves, since initial cambial divisions in the growing season are fueled from energy reserves held through the previous dormant season [61,67]. For the evergreen Quercus species in general, carbohydrate storage mainly occurs in autumn (October to December) and spring mobilization is especially high [61,68]. Spring mobilization of carbohydrate reserves is also related to its hydraulic properties, so the large vessels observed for this species are very sensitive to winter embolism [69]. A large part of the previous year's earlywood vessels are embolized by frost events in winter, and require the use of stored carbohydrates for production of new earlywood large vessels before leaf expansion for the spring recovery of hydraulic conductivity [70]. Additionally, roots usually continue to grow at the end of the growing season, if they have a sufficient supply of nutrients and water, until the soil temperatures become too low [71]. The growth of the trunk and the roots is, indeed, very related to the quantity of reserves accumulated before the growing season [28,44]. In turn, for Quercus suber, a higher growth rate is observed when the carbon reserves accumulated during the previous season are likely to support a strong vegetative development [72].
Regarding precipitation effects during the dry season, our results are in accordance with several studies showing a negative effect of summer water deficit on oak stem growth [16,63,[72][73][74]. A positive correlation between summer rainfall and radial growth can be the consequence of soil water being used immediately during the dry period for carbon assimilation [64], taking advantage of the warm temperature and of the still large photoperiod, and leading to a correlation between Net Ecosystem Exchange (NEE) and tree growth [75]. This relationship suggests a regulation of the evapotranspiration processes involved in the metabolism, in direct connection with respiration and photosynthesis [16]. However, the weak correlation between growth and precipitation in our study might be explained by the characteristic deep roots of oak species [76] which increases the maximum extractable soil water reserve and leads to a minimization of the effects of summer drought [57,64,70,77]. At the MS, however, a negative correlation was observed between wood ring widths and summer precipitation levels. This response could be related to the stormy and short events of summer precipitation that can break the dormancy which constrains the thickness of the cellular walls and the maturation of the latewood [17]. Another explanation could be that the stored carbohydrates can be diverted to support a high fructification instead of the radial growth [17,44,78].

Soil Water Budget Defining Drought Intensity and Wood Growth Phenology as A Key Tool for Understanding Tree Growth/Climate Relationship under Constrasted Soil Depth Conditions
To strengthen the statistical tree ring/climate analysis and explore more deeply the underlying processes driving tree growth, we analyzed the relationship between annual tree ring increment and the drought-driven phenological phases as well as the annual water stress integral limiting carbon assimilation under the theoretical framework of the "source" or "sink" limitation hypotheses. We found the best significant relationship between drought onset/offset and tree growth, suggesting that annual tree growth is related to the length of the growing season. This relationship was enhanced when additionally considering a secondary phenological phase in the fall, suggesting a significant contribution of the secondary growth in the annual tree growth in our site. This result then extends early results found for Quercus ilex [8,79] to Quercus suber, and highlights the potential role of the secondary growth under a mild Mediterranean climate where secondary growth is common, as previously observed on this site for leaf production and twig elongation [20].
In addition, we also identified the critical role of soil water holding capacity in determining the phenological phases and WSI and explaining the major contribution of the previous fall precipitation on tree growth. In our case, where the study site is characterized by a deep and stone-free soil, drought onset was highly correlated to winter precipitation levels only when considering a field capacity of 890 mm, allowing for different interannual soil water refill during the wet season, which was not observed when using soil parameters for shallow soils in the model. Our results then support and explain previous observations under deep soil conditions in Spain (Extremadura), where precipitation during November (t − 1) and December (t − 1) exerted a large positive influence on cork growth, suggesting the importance of the water holding capacity of the soil [59,60]. This importance of soil characteristics was also observed in southwestern Portugal [64], where the rainfall in November (t − 1) and December (t − 1) was not significant for the tree-diameter growth, when the sandy soil of the stand induces a low capacity to retain water. Similar results were observed for different Oak species growing on contrasted soil types across the Mediterranean basin [16]. The deep root system of Quercus suber, observed in caves at depths of up to 10 m [76], then contributes to secure the water uptake from the soil to maintain transpiration and tree growth [80,81], and in summer, Quercus suber continues to extract sufficient water from the soil [82] to maintain high leaf and tissue hydration for cell elongation [83]. Under deep soil conditions and mild temperatures, autumn and winter precipitation refills water reserves, a key aspect of the climate regime for tree growth during the following growing season [84][85][86]. We could not conclude from our analysis about the dominant role of either the growth period phenology t0t1t2 or the summer drought WSI affecting carbon assimilation. However, both variables were highly correlated to the previous fall precipitation P_ON under deep soil conditions, which was the key climate variable in the tree ring/climate statistical analysis. In turn, whatever the hypothesis (the sink limitation driven by drought onset or the source limitation driven by summer drought intensity), stem growth then appeared to be driven by the soil water refill during the previous fall under our deep soil conditions, an unexpected variable to consider in summer dry ecosystems.

Quercus Suber Vulnerability to Recent Climate Changes under Deep Soil Conditions: An Unstable Balance between Increasing Winter Temperatures and Earlier Prolonged Summer Drought
Our main result showed the statistical relation between annual tree growth and the autumnal rainfall before the growing season, which we could explain by the important contribution of the deep soil water refill during that period. Without clearly disentangling the predominant role of the length of the growing period (t1 − t0 + t2) or the intensity of the dry season WSI, as both variables were correlated to the same autumnal rainfall amount, we could show that both variables did not exert any significant trend over the recent period in allowing for the ecosystem sustainability. Actually, we observed that the start of the drought period (t1) was significantly earlier in the last decades in the region, as a result of this decreasing winter precipitation (P_ON) as previously observed in the region [87], so that both the length of growing period and the drought period should be affected. However, we also showed that the start of the growing season (t0) also moved toward an earlier start date as a consequence of warmer winter temperatures, compensating for the earlier (t1) and leading to constant (t1 − t0) growth period, a result recently pointed out for Quercus ilex in southern France [88]. The same conclusion is also reached for the summer drought intensity WSI, for which the earlier drought start is compensated by a slightly earlier (although non-significant) drought offset (t2). We then warn here, that despite no significant major trends in the two key drivers of tree growth (growth period or carbon assimilation), these two drivers are controlled by climate variables that are highly affected by recent climate changes, with compensating effects until now, but which could potentially result in the switching to a shrinkage of the growth phenological period, as recently observed with the winter temperature warming pause in southern France [88], or a prolonged drought at the end of the dry season as a consequence of enhanced potential evapotranspiration due to temperature increase as forecasted in the region [87,89,90].

Conclusions
Our study could provide a rare long-term wood-ring-width chronology for the Quercus suber species under cork oak production with strong common climatic signals. In our site, which is characterized by a favorable water balance on deep soil conditions, we identified that wood growth is mainly controlled by the water stored before the growing season and depends less on the spring and summer rainfall. Regarding the recent decline of cork oak forest in the region during the last few decades [21], management strategies should be based on fall/winter rainfall projections for the forthcoming decades, a neglected aspect of climate change in summer dry ecosystems such as the Mediterranean climate, but yet, already identified as a key control in some forested ecosystems [84]. The use of a functional process-based index for tree growth integrated both the climate variability and soil conditions, hardly captured by the statistical analysis on monthly climate variables alone. Our correlation coefficients were significant but weaker than previous studies [8] using a similar approach, illustrating the high uncertainty in the calculation of phenological phases, but nevertheless a conservation of the statistical significance despite this weakness. This study then illustrates the potential benefits of this functional index for further use in dendrochronological analysis, even under data-limited conditions for water budget calibrations. Process-based indices as phenological phases or drought indices driven by soil water budget, appear to represent efficient tools for integrating both the climatic and edaphic conditions and potentially leading to contrasted trends and patterns within uniform climate regimes at the regional scale [11].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Description of the Soil Water Balance Model
Daily variations in Soil water content (SWC) are calculated used the following equation (Equation (A1)): where P is daily precipitation based upon its interception by the canopy, E is soil evaporation, T is transpiration, and D is deep drainage, all being expressed in mm. Runoff and subsequent lateral water fluxes are not considered. Each day a fraction of incident precipitation is directly intercepted and evaporated at the canopy level, as a function of leaf area index (LAI) (Equation (A2)).
In (mm) = 0.5 LAI (A2) Throughfall water vertically infiltrates into the successive layers from the surface layer to the deeper ones following a three layer bucket model of different depths (0 cm-20 cm, 20 cm-1 m, 1 m-4.5 m). Soil depth is set constant, but a rock fragment content parameter is used to modify stand water availability among soil types. On a daily basis, the difference between actual soil water content and soil water content at field capacity (FC) is calculated inside each layer. This excess water is daily exported to the deeper layer. Excess water in the deepest layer is finally exported as deep drainage. Vegetation ability to extract water from soil is related to the predawn leaf water potential (Ψ), successfully used in vegetation functional models on Mediterranean species. Thus, each soil layer i is characterized by its water content SWCi (mm) and its water potential Ψi, these two variables being related by the power function model for the retention curve [36] (see Equation (A3)): where SC is the saturated water content (mm), Ψ e is the air entry potential (MPa), k is the rock fragment content (%), and b the power parameter. These parameters are derived from soil texture [37]. Vegetation transpiration and soil evaporation are then calculated as a function of PET. Daily PET is estimated using the Priestely-Taylor equation [35], which is well adapted for large scale simulations where local fine climate data are lacking, and suitable for Mediterranean biomes [91]. The PT coefficient a is set to 1.14. Thus, the maximum transpiration ET max of the day is related to PET and LAI according to the Beer Lambert (BL) law (see Equation (A4)). BL = 1 − e (−k·LAI) (A4) with k being the Beer Lambert extinction coefficient and LAI being the leaf area index. Actual transpiration (AET) is then calculated as a function of the maximum transpiration (ETmax) when soil water content is at field capacity, and is a variable simulating the stomatal closure of vegetation due to water stress. This limitation is calculated according to soil water potential (Ψ) for values ranging between field capacity and the specific minimum water potential for water extraction (Ψ lim ), beyond which uptake is zero (see Equation (A5)).
AET/ETmax= 1 − Ψ soil /Ψ lim (A5) According to the diffusion theory, cumulative soil evaporation (E) is considered proportional to the square root of the elapsed time after the start of evaporation, and is distributed across the profile according to an exponential decay. A list of inputs, species parameters of the model, and stand characteristics (from in situ measurements [20]) are provided below in Tables A1 and A2.