Stable Isotopes in Tree Rings of Pinus heldreichii Can Indicate Climate Variability over the Eastern Mediterranean Region

: A long-term context is important for understanding past climatic variability. Although tree-ring widths (TRWs) are widely used as a proxy for reconstructing past climate, the use of annually-resolved values of δ 13 C and δ 18 O tree-ring stable isotopes (TRSIs) is increasing and may provide further valuable information. Here, we present a 487-year-long TRW chronology and 240-year-long TRSI chronology for Bosnian pine ( Pinus heldreichii H. Christ) and compare them to each other. We demonstrate that both δ 13 C and δ 18 O values are better proxies for temperature, precipitation, and drought than TRW. The correlations between these climate parameters and TRSIs are strongest for the combined summer (JJA) period. The results of temporal and spatial ﬁeld correlation indicate that TRSI chronologies are stable, reliable proxies for JJA precipitation reconstruction over the whole Balkan Peninsula and surrounding eastern Mediterranean region. However, the stability of the temperature signal of the both δ 13 C and δ 18 O chronologies declines after the 1950s. Our work supports the emerging evidence that TRSI data track climate variability more accurately than a conventional TRW approach and can be subsequently used for the reconstruction of past climate.


Introduction
In addition to being geographically widespread, trees live for hundreds to thousands of years. Their longevity and sensitivity to their surrounding conditions makes trees valuable integrators of local environmental and climatic information [1][2][3]. These environmental changes are recorded on a yearly basis in tree rings. The growth and formation of individual tree rings reflect not only the internal physiology of wood formation, but also the influences of environmental conditions, particularly temperature and water availability. Tree rings therefore have enormous potential for the detection and reconstruction of past climatic conditions, including hydro-climatic extremes (e.g., [4]).
Besides conventional tree-ring width (TRW) analysis, tree-ring stable carbon (δ 13 C) and oxygen (δ 18 O) isotope ratios provide valuable climate-related information [5,6]. Changes in δ 13 C and δ 18 O values are closely connected to plant physiological processes, particularly to processes associated with photosynthetic CO 2 uptake [7] and H 2 O transpiration [8]. These processes are strongly coupled to natural variations in factors determining day-today fluctuations in weather and changes in climatic conditions over longer time periods.
In comparison to growth rate, tree-ring δ 13 C and δ 18 O values reflect particularly annual differences in temperature, precipitation, humidity, and light intensity [9,10].
Moreover, the δ 13 C and δ 18 O values reflect species-specific sensitivity to these environmental drivers, in addition to the overall habits of tree species, including the root system structure [11][12][13]. Previous studies have shown strong links between tree-ring stable isotopes (TRSIs) and many climatological phenomena including: drought/sunlight duration in the alpine region [14]; relative air humidity, cloud cover, and drought [15]; the Palmer drought severity index (PDSI; [16]); the standardized precipitation-evapotranspiration index (SPEI; [17]); and precipitation amounts [18]. However, the links between TRSIs and climate parameters are often stronger at wet sites than dry ones [19]. In dry environments, the δ 13 C value is predominantly influenced by fluctuations in vapor pressure deficit [20]. In wet environments where the diffusion of CO 2 across the stomata is not limited, solar radiation and the photosynthetic properties of the Rubisco enzyme are key factors [21].
Dendrochronological research in the Dinarides, a mountain range in southern and southeastern Europe that separates the Balkan Peninsula from the Adriatic Sea, has been limited to date. There is therefore huge potential in this region for the development of long, climate-sensitive tree-ring chronologies suitable for the reconstruction of sunshine duration [22], summer temperature [23,24], and drought [25].
In this study, we analyzed and compared the climate signals of TRWs and TRSIs in primaeval forests of Bosnian pine (Pinus heldreichii H. Christ), known for their exceptional longevity (e.g., [26]). Existing studies of the climate sensitivity of Bosnian pine TRWs indicate that although this species offers great potential for creating long-term chronologies, the response of Bosnian pine to climate is generally limited [27][28][29].
The main aim of this research was therefore to determine whether alternative proxies such as isotope data offer stronger climate reconstruction potential in comparison to conventional TRW chronologies [22,24,30]. Our hypothesis is that isotope data (δ 13 C and δ 18 O values) are better proxies for high-resolution reconstructions of past hydroclimatic variability and can significantly improve the climate signal in Bosnian pine tree-ring data.

Research Area
The study area (42 • 34 N, 18 • 32 E, elevation 1894 m a.s.l.) is situated in the Orjen Mountains, a transboundary Dinaric Mediterranean limestone mountain range that stretches for about 25 km through Montenegro and Bosnia and Herzegovina. Orjen Peak is the highest peak in the sub-Adriatic Dinarides. The sub-Adriatic range contains evergreen deciduous forests and the vegetation transitions from temperate forests in its lower elevation belts to conifers and tundra in its northern and upper regions. The belt of primaeval Bosnian pine forest is located at the highest elevations above 1800 m a.s.l. The underlying forest soils are characterized by shallow layers of rendzina soil atop limestone and dolomite bedrock [31,32].
The Orjen Mountains are situated within the Mediterranean subtropical belt, at the intersection between Mediterranean and continental climates (Figure 1). At the study site, summers are hot and sunny with mean temperatures up to 19 • C in August, and autumn, winter, and spring constitute the rainy season and temperatures can drop to a minimum of −10 • C. Average annual precipitation (rain, snow) is about 5000 mm. Located as it is in the littoral Dinarides, the study area experiences strong precipitation contrasts throughout the year. November thunderstorms can deposit 2000 mm of rainwater over the course of just a few days, whereas August is frequently completely dry [33]. Forest fires are a frequent disturbance during the dry season [34].

Tree Core Sampling, Tree-Ring Width, and Stable Isotope Measurement
In total, 23 Bosnian pine trees were sampled in the primaeval forest. One core per tree was extracted at breast height (1.3 m) using a Pressler borer (Haglof Company Group, Sweden) with a 5-mm inner diameter. To avoid the compression of the wood, the cores were sampled in the direction parallel to the slope [35]. After careful surface preparation to maximize the visibility of the tree rings, all core samples were measured using a VIAS TimeTable measuring system (SCIEM, Vienna, Austria) with a measuring length of 78 cm and a resolution of 1/100 mm). The tree-ring width (TRW) was measured to an accuracy of 0.01 mm, and each series was cross-dated using PAST4 [36]. The series were statistically controlled using the COFECHA program version Cofecha (Laboratory of Tree-Ring Research, University of Arizona, Tucson, AZ, USA) [37].
Four core samples were selected for the isotopic analysis. Each of the precisely dated annual increments was separated with a scalpel under a stereomicroscope and packed into Teflon filter bags (F57; Ankom Technology, Macedon, NY, USA). The alpha-cellulose was then extracted via the modified Jayme-Wise isolation method [38]. Teflon filter bags were washed twice for two hours each using 5% NaOH solution at 60 • C, followed by washing with 7% NaClO 2 solution (pH 4-5) for 30 hours at 60 • C. The samples were subsequently dried at 50 • C for 24 h, sealed in Eppendorf microtubes, and stored in the dark at room temperature (21 • C) before analysis. The samples of alpha-cellulose (0.2-1.0 mg) were weighed into tin boats and silver capsules (Elementar Analysensysteme, Langenselbold, Germany) for the determination of carbon and oxygen isotopes, respectively. For the δ 13 C measurements, the samples were combusted to CO 2 at 960 • C; δ 18 O samples were pyrolyzed to CO at 1450 • C using a high-temperature combustion cell of an elemental analyzer (EA) varioPYRO cube (Elementar Analysensysteme, Germany). The stable isotopes in the CO 2 and CO gases were then determined by a continuous flow isotope ratio mass spectrometer (irMS), ISOPRIME100 (Isoprime, Manchester, UK).
The system was calibrated using certified reference materials with known isotopic ratios from the International Atomic Energy Agency (IAEA, Vienna, Austria). The δ 13 C values were referenced to caffeine (IAEA-600) and graphite (USGS24). The δ 18 O values were referenced to benzoic acid (IAEA-601 and IAEA-602). The δ 13 C and δ 18 O values (‰) were given with respect to the deviation from the Vienna Pee Dee Belemnite (VPDB) and Vienna Standard Mean Ocean Water (VSMOW) standards, respectively. The long-term reproducibility of these standards, evaluated as standard deviation, was ≤0.05‰ (IAEA-600 and USGS24), ≤0.08‰ (IAEA-601), and ≤0.11‰ (IAEA-602). Standard deviations were ≤ 0.04‰ (δ 13 C) and ≤0.09‰ (δ 18 O) for five consecutive measurements of the homogenized alpha cellulose sample. For details see Urban et al. [39]. The δ 13 C time series were corrected for the δ 13 C decrease in the atmosphere [5]. The correction for atmospheric δ 13 C depletion was based on the compilation of ice core and direct measurements from Mauna Loa (https://www.esrl.noaa.gov/gmd/dv/data/). No additional physiological corrections of δ 13 C [18] were applied.
Because recent results show that stable isotopes manifest as a constant spread versus a level relationship over the lifespan of the tree [40], only the TRW datasets were standardized to suppress non-climatic factors. Negative exponential curves (NegExp), cubic smoothing splines (with a 50% frequency response cut-off at 150 years; spline150yr) and regional curve standardization (RCS) method were applied to remove age-related growth trends using the ARSTAN software version ARS41d_xp (Tree-Ring Laboratory, Lamont Doherty Earth Observatory of Columbia University, Palisades, NY, USA) [41]. TRW indices were calculated as residuals from estimated growth curves after applying an adaptive power transformation to the raw measurement series [42]. The final chronologies from each of the three detrending techniques were calculated using robust bi-weighted means. The expressed population signal (EPS; [43]) and inter-series correlation (Rbar) were calculated to assess the quality of each chronology. Because no significant differences were observed between the different chronology versions (Figure 2A), the residual chronology, after applying the 150-year smoothing spline, was selected to calculate the correlations with the monthly mean temperatures, precipitation sums, and Palmer drought severity index (PDSI; [44]). The correlations were calculated using the R package "treeclim" version 2.0.5.1 [45].
Temporal and spatial correlations between TRW and TRSI chronologies and seasonal climatic parameters were analyzed using monthly gridded data (0.5 × 0.5 • grid) from the CRU TS3.24.01 database, available via the Royal Netherlands Meteorological Institute (KNMI) Climate Explorer platform (http://climexp.knmi.nl) [46]. Local climate observation data for the study area are very rare and time-limited, as is characteristic of the Dinaric region as a whole. For that reason, gridded CRU TS3.24.01 climatic data for the period of 1901−2018 were used. In our case, data from a nearby station (Crkvice) were compared with the CRU TS3. 24.01 data to determine how well the downloaded data approximates local climate conditions ( Figure 1D). Pearson's correlation coefficients were calculated from April of the previous year to October of the current year for each climatic factor. The temporal stability of the climate signal was analyzed using moving window correlations with a 45-year interval plus 1 year. Analyses were performed with the most significant seasonal variables for the 1901-2018 period using CRU TS3.24.01 climate data.

Figure 2. (A)
Raw tree-ring width chronology with a trend approximated by an exponential function, sample replication, and three slightly different pine residual chronologies after applying three different standardization techniques: cubic smoothing splines with 50% frequency cut-off at 150 years, negative exponential functions, and regional curve standardization (RCS) method. The upper left inset provides insight into data characteristics. Basic statistical parameters (MSL-mean segment length; AGR-average growth rate; SD-standard deviation; AC1-first-order autocorrelation) are shown. (B) Expressed population signal (EPS) and inter-series correlation (Rbar) of the indexed TRW series. EPS and Rbar statistics were calculated over 30-year windows lagged by 25 years.

Tree-Ring Width Chronology
The 487-year TRW chronology covers the period of 1531-2018 ( Figure 2). The minimum length of the TRW series (160 years) and the mean segment length (297 years) indicate that the chronology represents old trees. Although replication of the chronology decreases backwards in time, it does not drop below 23 TRW series during the studied period . The expressed population signal (EPS) fluctuates between 0.87 and 0.94. The high first-order autocorrelation (0.80) of the raw TRW chronology indicates large temporal memory.

Tree-Ring Stable Isotope Chronologies
Annually-resolved δ 13 C and δ 18 O chronologies obtained from four separate cores ( Figure 3A) are highly synchronous over the 1780-2018 period as shown by EPS and Rbar ( Figure 3B). Both the uncorrected δ 13 C time series and the time series corrected for the increase in CO 2 concentration in the atmosphere are shown; hereafter, all presented δ 13 C time series are atmospheric-corrected data. The strength of the common signal over the whole period is confirmed by the consistently high mean Rbar (amounting to 0.51 and 0.55 for δ 13 C and δ 18 O, respectively) and EPS values (amounting to 0.80 and 0.83 for δ 13 C and δ 18 O, respectively).

TRW and TRSI Responses to Climate
Correlations between TRW and TRSI chronologies and seasonal climatic parameters were obtained using the monthly gridded data (0.5 × 0.5 • grid) from the CRU TS3.24.01 database. The sum of the current year summer precipitation is most strongly correlated with TRW ( Figure 4). Summer (JJA) precipitation and PDSI are positively correlated with TRW, whereas monthly mean temperature over the same period is negatively correlated with TRW. Generally, correlations between radial growth and precipitation sum and PDSI of the previous growing season and the winter months are low and statistically nonsignificant. The correlation coefficients substantially increase during summer months of the current year. However, the mean monthly temperature of the previous growing season has no significant effect on TRW. The positive effect of winter temperatures on growth (up to April) gradually becomes negative in summer. This pattern suggests that the availability of water during the current summer is positively related to radial growth. Precipitation during the previous growing season is less important for pine growth. The δ 13 C isotopic signal shows a significant positive correlation with temperature in June (r = 0.22) and August (r = 0.23) of the current year, and an even stronger negative correlation with precipitation during the summer months, particularly June (r = −0.45). For precipitation, the highest correlation is found when the whole summer period (JJA) is included (r = 0.60). Regarding the PDSI index, almost all months of both previous and current years are significantly negatively correlated with the δ 13 C record, with a maximum correlation in August of the current year (r = −0.44).
The δ 18 O isotope record is generally more strongly correlated with climate parameters than the δ 13 C record. Significant positive correlations with temperature are found in July, August, and September of the previous year and during the spring-summer period (May-August) of the current year. For individual months, the strongest correlation was found for August (r = 0.43); the correlation is even stronger for the combined JJA period (r = 0.48). Significant negative correlations with precipitation were found from May to August. As with temperature, the highest correlations between δ 18 O and precipitation and between δ 18 O and PDSI were observed for the JJA period (r of −0.48 and −0.43, respectively). Generally, the extension of JJA period and/or selection of other combined period does not improved the correlations of TRW and TRSI records with climate parameters (data not shown). Robust correlations between δ 18 O and PDSI index were also observed for separate months from May through October (Figure 4). Correlations of individual time series with climate parameters have shown similar results as the correlations of averaged δ 13 C and δ 18 O chronologies.

Temporal Stability of the Climate Signal
The temporal stability of the climate signal was analyzed using bootstrap correlations between the δ 13 C and δ 18 O chronologies and the most highly correlated climatic factors, which in our case are the JJA precipitation sum and the mean JJA temperature, respectively. The stability signal of JJA precipitation is prominent throughout the analyzed period (1901-2018; Figure 5A). In contrast to precipitation, the stability of the temperature signal of the both δ 13 C and δ 18 O chronologies declines slightly through the analyzed period, especially after 1950 ( Figure 5A). In contrast to δ 13 C, the climate signal of the δ 18 O chronology is not stable over the analyzed period for either temperature or precipitation ( Figure 5A). As with the δ 13 C temperature signal, the δ 18 O signal decreases over time, especially after the 1950s. Our results confirm that the JJA precipitation sum is the most dominant and temporally stable climatic factor influencing the δ 13 C values in tree rings. We further tested the stability of the relationship between JJA precipitation and δ 13 C values using separate 60-year calibration and verification periods ( Figure 5B). The calibration-verification model shows positive and significant (p < 0.01) predictive capability, regardless of which period is used for the calibration (Figure 5B,C).

Spatial Field Correlations
The JJA precipitation sum and mean JJA temperatures downloaded from the CRU TS3.24.01 gridded dataset were further used to explore spatial field correlations. The strongest correlation between JJA precipitation sum and δ 13 C appears over the whole Balkan Peninsula (Figure 6). The δ 13 C and δ 18 O chronologies and mean JJA temperature only correlate strongly for the period 1901-1958, but diverge after the 1950s (Figure 6). The results of spatial field correlation indicate that the δ 13 C chronology is a stable, reliable proxy for JJA precipitation reconstruction over the eastern Mediterranean region. Figure 6. Spatial field correlations between δ 13 C and δ 18 O data and summer (JJA) precipitation sums and mean JJA temperatures for two subperiods. Subperiods are determined based on moving correlation analysis and split calibrationverification results. Spatial correlations between tree-ring stable isotope values and climatic variables were analyzed using monthly gridded data (0.5 × 0.5 • grid) from the CRU TS3.24.01 database.

Discussion
In the past decade, a number of studies have been published regarding the development of multi-century temperature-and precipitation-sensitive tree-ring chronologies from primaeval pine forests in the Balkans [28,47]. These studies show that the climate signal in the TRWs of pine species can be fairly weak [28,48], and more strongly dependent on local site conditions than on regional climate patterns [48]. However, the climatic signal in the width of late wood seems to be stronger (spring precipitation; [49]); maximum late wood density exhibits an even stronger climate signal (summer temperature; [28,47,49]).
In this study, we present a 487-year-long TRW chronology, thus extending the chronology for this region by 50 years [24,30]. Because the climate signal in TRW chronologies of P. heldreichii is relatively limited for robust climate reconstruction (also shown by [24]), we tested the hypothesis that TRSIs (δ 13 C and δ 18 O values) are better proxies for highresolution reconstructions of hydroclimatic variability. Indeed, our results show higher correlations between TRSIs and climate parameters compared to TRW. The highest correlations were found when temperature, precipitation, and PDSI of the whole summer period (June-August; JJA) are included. These results are in general agreement with an earlier study by Hafner et al. [50] that confirm a robust positive relationship between δ 13 C and δ 18 O values and summer temperature, and a negative relationship with precipitation in the southeastern European Alps. Compared to the results of Levanič et al. [24], we observe a somewhat stronger relationship between δ 13 C and summer precipitation. They found correlations between δ 13 C and JJA temperature up to 0.6, although it was only −0.51 for JJA precipitation over the Balkan Peninsula. Even results from geographically distant locations of south-eastern China [16] confirmed that δ 18 O is a strong indicator of both precipitation and regional PDSI index reaching negative correlation coefficient below −0.6 for the June to October period of the current year. Similarly, we observe significant negative correlations between δ 18 O and the PDSI index and precipitation for the shorter JJA period. We observe a stronger relationship between PDSI and δ 18 O than between PDSI and δ 13 C, which is in accordance with the results of Esper et al. [19] and Rybníček et al. [51].
From an eco-physiological perspective, the tree-ring δ 13 C value is primarily modulated by the temperature-driven carboxylation rate and/or by the diffusive stomatal conductance of CO 2 into the leaves. Particularly at humid sites, where stomatal conductance and intercellular CO 2 concentration are not limiting factors of photosynthesis [52], higher temperatures result in higher δ 13 C values [53]. At dry sites, however, stomatal conductance may play a crucial role in 13 C discrimination [54]. As shown by Farquhar et al. [55], low intercellular CO 2 concentrations lead to reduced discrimination of 13 C, thereby resulting in an increase in δ 13 C values. Such a mechanism is likely responsible for the higher δ 13 C values at low sum of JJA precipitation observed in our study.
The ratio of photosynthetic CO 2 uptake to stomatal conductance, both of which are physiological processes contributing to carbon discrimination, is defined as plant water use efficiency (WUE). Accordingly, δ 13 C values of tree-ring cellulose are thought to provide particularly important insights into the tree's WUE and other related physiological processes [7,56,57]. However, it is not possible to determine whether a change in WUE is caused by carboxylation, stomatal conductance, and/or a combination of the two.
Based on the theory, variability in δ 18 O values could potentially reflect changes in stomatal conductance due to a regulated rate of transpiration. It is assumed that δ 18 O values in plants are negatively correlated with stomatal conductance but independent of photosynthetic CO 2 uptake [58]. Therefore, the combined analysis of δ 13 C and δ 18 O values represents a substantial improvement compared to previous work. It should be noted that both the isotopic signal of source water and post-photosynthetic and post-evaporative oxygen atom exchange processes could affect the final tree-ring δ 18 O signal [54]. Among other factors, vapor pressure deficit (VPD) is thought to influence the δ 18 O signal most substantially. When VPD increases (i.e., environmental conditions become drier), more evaporation occurs and the remaining water becomes enriched in 18 O. This is because water molecules containing the lighter 16 O isotope evaporate more readily than water molecules containing the heavier 18 O isotope. VPD has been shown to be one of the most effective drivers of changes in stomatal aperture, leading to reductions in transpiration of H 2 O from leaves and the diffusion of CO 2 into the leaves when VPD is high [52,59].
Indeed, several recent studies confirm that δ 18 O values can be a robust and reliable indicator of drought conditions and/or extreme hydroclimatic events [57,60,61]. We confirmed that δ 18 O values are tightly correlated with JJA precipitation sum over the substantial part of the Balkan Peninsula. The results of spatial field correlation, however, indicate that the δ 13 C chronology is an even stronger and more stable proxy for summer precipitation reconstruction. However, no spatio-temporal correlation between TRSIs and temperature during 1959-2018 was detected. This might be related to the widely documented phenomena known as "divergence", where higher instrumental temperatures are not reflected in the TRWs proxies (e.g., [62]). We hypothesize that an increasing summer temperature over recent decades is no longer a limiting factor for growth of pines at high altitudes with very short vegetation periods. Moreover, lower amounts of summer precipitation, together with an increased evapotranspiration due to higher temperature, provide a stronger and stable precipitation signal over the past century. More research including an extended dataset and detailed analyses is necessary to explain the causality of the divergence phenomenon.
There are several hypotheses trying to explain this discrepancy, but probably combinations of different environmental and anthropogenic factors are involved (e.g., [63,64]). Stine and Huybers [63] showed that no divergence is observed in arctic ecosystems with sufficient (non-limiting) light availability (i.e., no volcanic/dimming/pollution effects); these authors recommend an isotopic signal analysis to account for the anomaly. Brownlee et al. [64] have shown that one of the suggested causes of divergence, drought, is not likely an explanatory factor in the Alaskan forests. A recent overview of divergence problem in the climate reconstruction is provided by Wilmking et al. [3], in which the need for correct statistical treatment of the data is strongly emphasized to avoid improper interpretation of the past climate changes and extremes. As many as two-thirds of published studies do not test for the stationarity of tree growth vs. climate reconstructions, whereas more than half of those which did found non-stationarity in the data. This suggests a dynamic nature of the tree-environment relationship; accordingly, all results should be interpreted with a caution, especially those not performing the stationarity test.

Conclusions
In this study, we present a 487-year-long TRW chronology, thus extending the chronology for this region by 50 years [24,30]. This is also the first study to use stable carbon and oxygen isotopes for climate reconstruction in this geographical region. By comparing TRW records with δ 13 C and δ 18 O records, we show that tree-ring stable isotopes from P. heldreichii growing at high elevations correlate strongly with factors reflecting the climate variability during the summer months (JJA). These correlations are stronger than those found for tree-ring widths. Values of δ 18 O and particularly δ 13 C are stable and reliable proxies of JJA precipitation and can be further used for paleoclimatic reconstructions of hydroclimatic conditions over the whole Balkan Peninsula and surrounding eastern Mediterranean region.

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