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

Abstract

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 δ¹³C and δ¹⁸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 δ¹³C and δ¹⁸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 field 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 δ¹³C and δ¹⁸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.

1. 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 (δ¹³C) and oxygen (δ¹⁸O) isotope ratios provide valuable climate-related information [5,6]. Changes in δ¹³C and δ¹⁸O values are closely connected to plant physiological processes, particularly to processes associated with photosynthetic CO₂ uptake [7] and H₂O transpiration [8]. These processes are strongly coupled to natural variations in factors determining day-to-day fluctuations in weather and changes in climatic conditions over longer time periods. In comparison to growth rate, tree-ring δ¹³C and δ¹⁸O values reflect particularly annual differences in temperature, precipitation, humidity, and light intensity [9,10].

Moreover, the δ¹³C and δ¹⁸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 δ¹³C value is predominantly influenced by fluctuations in vapor pressure deficit [20]. In wet environments where the diffusion of CO₂ 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 (δ¹³C and δ¹⁸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.

2. Materials and Methods

2.1. 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].

2.2. 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₂ 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 δ¹³C measurements, the samples were combusted to CO₂ at 960 °C; δ¹⁸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₂ 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 δ¹³C values were referenced to caffeine (IAEA-600) and graphite (USGS24). The δ¹⁸O values were referenced to benzoic acid (IAEA-601 and IAEA-602). The δ¹³C and δ¹⁸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‰ (δ¹³C) and ≤0.09‰ (δ¹⁸O) for five consecutive measurements of the homogenized alpha cellulose sample. For details see Urban et al. [39]. The δ¹³C time series were corrected for the δ¹³C decrease in the atmosphere [5]. The correction for atmospheric δ¹³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 δ¹³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.

3. Results

3.1. 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 (1901–2018). 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.

3.2. Tree-Ring Stable Isotope Chronologies

Annually-resolved δ¹³C and δ¹⁸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 δ¹³C time series and the time series corrected for the increase in CO₂ concentration in the atmosphere are shown; hereafter, all presented δ¹³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 δ¹³C and δ¹⁸O, respectively) and EPS values (amounting to 0.80 and 0.83 for δ¹³C and δ¹⁸O, respectively).

The values of mean δ¹³C chronology range between −22.3‰ and −19.9‰ over the investigated period. The average δ¹³C value (±standard deviation) over the whole period is −21.0 ± 0.5‰; over the period 1901–2018, it is −21.2 ± 0.4‰. δ¹⁸O values show more variability, ranging between 24.8‰ and 30.2‰. Average values over the 1780–2018 and 1901–2018 periods are 27.7 ± 0.9‰ and 28.0 ± 0.9‰, respectively.

3.3. 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 non-significant. 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 δ¹³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 δ¹³C record, with a maximum correlation in August of the current year (r = −0.44).

The δ¹⁸O isotope record is generally more strongly correlated with climate parameters than the δ¹³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 δ¹⁸O and precipitation and between δ¹⁸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 δ¹⁸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 δ¹³C and δ¹⁸O chronologies.

3.4. Temporal Stability of the Climate Signal

The temporal stability of the climate signal was analyzed using bootstrap correlations between the δ¹³C and δ¹⁸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 δ¹³C and δ¹⁸O chronologies declines slightly through the analyzed period, especially after 1950 (Figure 5A). In contrast to δ¹³C, the climate signal of the δ¹⁸O chronology is not stable over the analyzed period for either temperature or precipitation (Figure 5A). As with the δ¹³C temperature signal, the δ¹⁸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 δ¹³C values in tree rings. We further tested the stability of the relationship between JJA precipitation and δ¹³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).

3.5. 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 δ¹³C appears over the whole Balkan Peninsula (Figure 6). The δ¹³C and δ¹⁸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 δ¹³C chronology is a stable, reliable proxy for JJA precipitation reconstruction over the eastern Mediterranean region.

4. 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 (δ¹³C and δ¹⁸O values) are better proxies for high-resolution 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 δ¹³C and δ¹⁸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 δ¹³C and summer precipitation. They found correlations between δ¹³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 δ¹⁸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 δ¹⁸O and the PDSI index and precipitation for the shorter JJA period. We observe a stronger relationship between PDSI and δ¹⁸O than between PDSI and δ¹³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 δ¹³C value is primarily modulated by the temperature-driven carboxylation rate and/or by the diffusive stomatal conductance of CO₂ into the leaves. Particularly at humid sites, where stomatal conductance and intercellular CO₂ concentration are not limiting factors of photosynthesis [52], higher temperatures result in higher δ¹³C values [53]. At dry sites, however, stomatal conductance may play a crucial role in ¹³C discrimination [54]. As shown by Farquhar et al. [55], low intercellular CO₂ concentrations lead to reduced discrimination of ¹³C, thereby resulting in an increase in δ¹³C values. Such a mechanism is likely responsible for the higher δ¹³C values at low sum of JJA precipitation observed in our study.

The ratio of photosynthetic CO₂ uptake to stomatal conductance, both of which are physiological processes contributing to carbon discrimination, is defined as plant water use efficiency (WUE). Accordingly, δ¹³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 δ¹⁸O values could potentially reflect changes in stomatal conductance due to a regulated rate of transpiration. It is assumed that δ¹⁸O values in plants are negatively correlated with stomatal conductance but independent of photosynthetic CO₂ uptake [58]. Therefore, the combined analysis of δ¹³C and δ¹⁸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 δ¹⁸O signal [54]. Among other factors, vapor pressure deficit (VPD) is thought to influence the δ¹⁸O signal most substantially. When VPD increases (i.e., environmental conditions become drier), more evaporation occurs and the remaining water becomes enriched in ¹⁸O. This is because water molecules containing the lighter ¹⁶O isotope evaporate more readily than water molecules containing the heavier ¹⁸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₂O from leaves and the diffusion of CO₂ into the leaves when VPD is high [52,59].

Indeed, several recent studies confirm that δ¹⁸O values can be a robust and reliable indicator of drought conditions and/or extreme hydroclimatic events [57,60,61]. We confirmed that δ¹⁸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 δ¹³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.

5. 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 δ¹³C and δ¹⁸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 δ¹⁸O and particularly δ¹³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.