<figure><img src="image_1.png" /></figure>

# Atmospheric CO₂ forcing on Mediterranean biomes during the past 500 kyrs

Received: 15 August 2022

Accepted: 14 March 2023

Published online: 25 March 2023

<figure><img src="image_2.png" /></figure>

Check for updates

Andreas Koutsodendris Ⓔ<fnref n="1" />✉, Vasilis Dakos Ⓔ²,³, William J. Fletcher⁴, Maria Knipping⁵, Ulrich Kotthoff⁶,⁷, Alice M. Milner Ⓔ⁸, Ulrich C. Müller<fnref n="9" />,¹⁰, Stefanie Kaboth-Bahr<fnref n="1" />,¹¹, Oliver A. Kern<fnref n="1" />, Laurin Kolb<fnref n="1" />, Polina Vakhrameeva Ⓔ<fnref n="1" />, Sabine Wulf¹², Kimon Christanis¹³, Gerhard Schmiedl Ⓔ⁶ & Jörg Pross<fnref n="1" />,¹⁰

There is growing concern on the survival of Mediterranean forests under the projected near-future droughts as a result of anthropogenic climate change. Here we determine the resilience of Mediterranean forests across the entire range of climatic boundary conditions realized during the past 500 kyrs based on continuous pollen and geochemical records of (sub)centennial-scale resolution from drillcores from Tenaghi Philippon, Greece. Using convergent cross-mapping we provide empirical confirmation that global atmospheric carbon dioxide (CO₂) may affect Mediterranean vegetation through forcing on moisture availability. Our analysis documents two stable vegetation regimes across the wide range of CO₂ and moisture levels realized during the past four glacial-interglacial cycles, with abrupt shifts from forest to steppe biomes occurring when a threshold in precipitation is crossed. Our approach highlights that a CO₂-driven moisture decrease in the near future may bear an impending risk for abrupt vegetation regime shifts prompting forest loss in the Mediterranean region.

The Mediterranean region is characterized by a mosaic of subtropical and temperate forests thriving under water-stressed conditions<fnref n="1" />,². Besides representing vulnerable biodiversity hotspots³, these Mediterranean forests provide a wide range of critical ecosystem services such as support of soil, regulation of regional climatic and hydrological conditions, and supply of food and timber⁴. In light of yet unabated anthropogenic CO₂ release, there is increasing concern on the resilience of these forest biomes to the recurring droughts⁵,⁶ as they are projected for the Mediterranean region in the next 70 years⁴,⁷. Climate scenarios predict structural and compositional changes in

Mediterranean forest biomes unprecedented in the Holocene that will cause disruption of ecosystem services and trigger biodiversity loss⁸,<fnref n="9" />. To assess the resilience of Mediterranean forest biomes to future climate forcing, it is essential to determine the nature and rate of change, and forcing mechanisms that these biomes have experienced under the natural range of climate conditions during the Quaternary¹⁰.

Here, we analyze a unique, exceptionally long record of Mediterranean forest dynamics based on centennial-scale-resolution (-200 years) pollen data from drillcores from Tenaghi Philippon (Greece; Supplementary Fig. 1) covering the past 500 kyrs

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Institute of Earth Sciences, Heidelberg University, Heidelberg, Germany. ²Institute des Sciences de l'Évolution, Université de Montpellier, CNRS, IRD, EPHE, Montpellier, France. ³Institut d'Écologie et des Sciences de l'Environnement de Paris (iEES Paris), Sorbonne Université, Paris, France. ⁴School of Environment, Education and Development, The University of Manchester, Manchester, UK. ⁵Department of Molecular Botany, Institute of Biology, University of Hohenheim, Stuttgart, Germany. ⁶Center for Earth System Research and Sustainability, Institute of Geology, Hamburg University, Hamburg, Germany. ⁷Leibniz Institute for the Analysis of Biodiversity Change (LIB), Hamburg, Germany. ⁸Department of Geography, Royal Holloway University of London, London, UK.</footnote>

<footnote n="9">⁹Parlamentsstraße 32, 60385 Frankfurt am Main, Germany. ¹⁰Senckenberg Biodiversity and Climate Research Centre (SBiK-F), Senckenberg Gesellschaft für Naturforschung, Frankfurt am Main, Germany. ¹¹Institute of Geosciences, University of Potsdam, Potsdam-Golm, Germany. ¹²School of the Environment, Geography and Geosciences, University of Portsmouth, Portsmouth, UK. ¹³Department of Geology, University of Patras, Rio, Greece.</footnote>

<footnote>e-mail: andreas.koutsodendris@geow.uni-heidelberg.de</footnote>
continuously¹¹. Augmented by decadal-scale-resolution geochemical data (-30 years) that record precipitation variability, our datasets resolve vegetation and moisture dynamics at timescales that do justice to the tempo of anthropogenic climate change. The Tenaghi Philippon cores have excellent chronostratigraphical control based on the integration of radiocarbon dating, tephrochronology and orbital tuning (Methods), and cover the full range of climatic boundary conditions realized during the past four glacial/interglacial (G/I) cycles¹², including the highest (i.e., ~300 ppmv) atmospheric CO₂ concentrations recorded in ice cores¹³. Situated in an intermediate position between the temperate and tropical climate systems of the Northern Hemisphere (NH), and located at the moisture boundary separating the wetter western from the drier eastern Mediterranean region, the Tenaghi Philippon site is exceptionally sensitive to abrupt climate change¹¹,¹⁴. This is illustrated by the prompt response of tree populations in the Tenaghi Philippon area to even subtle climatic forcing as shown by pollen records for both glacial¹⁵⁻¹⁷ and interglacial¹⁸⁻²⁰ periods. Tenaghi Philippon was not a glacial refugium for trees in contrast to more humid intermontane locations²¹ such as Ioannina²² and Lake Ohrid²³ from which vegetation records spanning multiple glacial-interglacial cycles have also been recovered. Hence, our record from Tenaghi Philippon is representative for vegetation biomes that are close to their lower moisture-tolerance limits such as the 'subtropical dry forest' biome that extends across c. 76 million ha in the Mediterranean region (Supplementary Fig. 1)². Importantly, because our datasets are from a single site rather than being spliced from multiple sites with distinctly different environmental boundary conditions, they yield an internally consistent signal of climatic and ecological change.

# Results and discussion

## Vegetation dynamics at Tenaghi Philippon over the past 500 kyrs

Our palynological results document two dominant alternating vegetation regimes at Tenaghi Philippon under the full spectrum of climatic boundary conditions realized during the past 500 kyrs: an open vegetation dominated by steppe biomes (tree-pollen <35 %) during glacials and forest biomes (tree-pollen >65 %) during interglacials (Fig. 1). While the transitions from glacials to interglacials are persistently marked by rapid reorganizations from steppe to forest biomes, the transitions from interglacials to glacials are highly dynamic and exhibit increased frequencies of abrupt millennial-scale shifts between forest and steppe biomes. Our high-resolution dataset documents that abrupt shifts in vegetation regime at Tenaghi Philippon were not confined to the last G/I cycle¹⁵,¹⁹, but are a persistent feature of the past four G/I cycles. While the duration of the regime shifts is persistently <200 years across the record, our palynological data from core intervals with higher-than-average temporal resolution indicate that these shifts actually occur on decadal scales.

A pervasive characteristic of the terminal forest stages is the expansion of drought-tolerant, pine-dominated biomes coval with the decrease of cold-tolerant, but drought-sensitive montane and shrub biomes (Fig. 1). This structural change within forest biomes occurred across all past four interglacials, despite the fact that each of these interglacials exhibits a unique floristic composition¹¹,¹⁴. Based on our vegetation-independent age model, the demise of these drought-tolerant forest biomes persistently postdates the end of the warmest interglacial sub-stages (such as Marine Isotope Stages [MIS] 5e, 7e, 9e, and 11c) by about 5 to 20 kyrs (Supplementary Table 1) as determined in global temperature-proxy datasets such as the LR04 benthic δ¹⁸O stack and Antarctic ice-core records¹²,²⁴ (Fig. 2), as well as in regional Mediterranean climate records²⁵. The constant lag of the forest demise at Tenaghi Philippon with regard to global temperature-proxy datasets at times of glacial build-up implies a strong resilience of Mediterranean forest biomes to temperature decrease. This is supported by the continuous presence of cold-sensitive temperate and Mediterranean

tree taxa beyond the end of the warmest interglacial sub-stages (Fig. 1). On this basis, the resilience of Mediterranean forest biomes across glacial inceptions appears to be determined by moisture availability¹⁵,²⁶. However, validation of this hypothesis has yet been precluded by the absence of paired high-resolution (i.e., decadal- to centennial-scale) precipitation and vegetation proxy-record timeseries from the same locality.

To this end, we evaluate the log(Ca/Fe) fluctuations in the Tenaghi Philippon cores as an independent proxy for local precipitation. As the minerogenic content of the peat at Tenaghi Philippon is connected to an extensive karstic system in the mountain ranges surrounding the site, increased precipitation results in high Ca and low Fe contents in the Tenaghi Philippon sediments and vice versa (Methods). The log(Ca/Fe) values peak during interglacial maxima, indicating maximum precipitation at these times, and progressively decrease as global climate evolves gradually towards larger ice-sheets, colder conditions and lower greenhouse-gas concentrations of glacials (Fig. 2; Supplementary Fig. 2). The application of rolling correlations between the log(Ca/Fe) dataset from Tenaghi Philippon and the CO₂ and CH₄ records from Antarctic ice cores yields strong correlations between the records during peak interglacial conditions, whereas weak correlations prevail during the terminal interglacial phases (Supplementary Fig. 3). While this demonstrates that the hydrological record from Tenaghi Philippon responds sensitively to global climate dynamics, it also implies a systematic decoupling between the gradual moisture decrease and the abrupt forest contractions at times of glacial inception (Fig. 2).

## Assessing causal associations using convergent cross-mapping

To validate the hypothesis of a causal association between vegetation-regime shifts and hydroclimate change, we apply convergent cross-mapping (CCM) that can resolve causalities in complex systems with only partially coupled variables²⁷, hence overcoming the shortcomings of simple correlation methods. Specifically, the CCM analysis allows to distinguish between direct and indirect causation from simple correlation between timeseries by assessing the extent to which an affected variable reliably estimates states of a causal variable (termed 'CCM skill') (ref. ²⁷). In other words, if a variable A can be predicted using values of variable B, then we conclude that variable A is causally affecting variable B (Methods). As the cross-mapped estimates improve in 'CCM skill' with the length of the timeseries²⁷, our high-resolution palynological and geochemical datasets are exceptionally well suited for this approach.

Our results reveal increasing 'CCM skills' for the tree-pollen and log(Ca/Fe) datasets from Tenaghi Philippon (Fig. 3), suggesting a significant causal association between forest dynamics and precipitation. To confirm this unidirectional causality from bidirectional causality, which is often related to synchrony caused by forcing from a third independent factor, we evaluate the optimal time displacements of the 'CCM skills' (refs. ²⁸,²⁹) (Methods). Optimal time displacement means that we test the 'CCM skill' between tree-pollen and log(Ca/Fe) for different time lags. In that way, we can distinguish true causal associations from simple correlations and explore at which time lag the causal relationship between the two variables is strongest. Despite the bidirectional coupling between the tree-pollen and log(Ca/Fe) datasets as derived by the 'CCM skills', the peak for CCM estimates of the tree-pollen timeseries based on the log(Ca/Fe) timeseries, i.e., the effect of log(Ca/Fe) on tree-pollen, has an optimal lead of ~4–5 kyrs behind that for estimates of tree-pollen timeseries based on the log(Ca/Fe) timeseries (Fig. 3; Supplementary Table 2). Hence, the response of tree-pollen to log(Ca/Fe) change is significantly faster (<1 kyrs) than the response of log(Ca/Fe) to changes in tree-pollen (>4.5 kyrs). Collectively, the strong causality between the two timeseries and the rapid response of the tree-pollen signal to the log(Ca/Fe) variability testify that the abrupt vegetation-regime shifts at Tenaghi Philippon are determined by a crossing of precipitation thresholds.
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Overview of the pollen record from Tenaghi Philippon spanning the past 500 kyrs. a Tree-pollen percentages including (dark green) and excluding *Pinus* (olive green). b Temperate (light green) and Mediterranean (purple) tree taxa. c Montane tree taxa. d Pioneer tree taxa. e Shrubs. f Steppic taxa. g Other herbs.</figcaption></figure>

h Cyperaceae and Poaceae (light blue) and helophytes (dark blue). i Aquatic taxa. j Fern spores. k Algae. Marine Isotope Stages (MIS) adapted from ref. ⁶⁸, and warmest interglacial sub-stages (vertical dashed lines) after ref. ²⁴.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Vegetation and precipitation data from Tenaghi Philippon and selected climate records for the past 500 kyrs. **a** LR04 benthic $\delta^{18}$O stack$^{60}$. **b** and **c** atmospheric CH$_4$ (ref. $^{31}$) and CO$_2$ (ref. $^{13}$) concentrations from Antarctic ice cores, respectively. **d** log(Ca/Fe) and **e** tree-pollen percentages from Tenaghi Philippon including (dark green) and excluding *Pinus* (light green). **f** *Neogloboquadrina pachyderma* 'sinistral' (NPS) percentages from ODP Site 983 (ref. $^{32}$). **g** and **h** Ti/Al and Ba/Al ratios from ODP Site 967, respectively (ref. $^{33}$). **i** Orbital parameters$^{34}$.</figcaption></figure>

Horizontal dashed lines for the LR04 benthic $\delta^{18}$O stack (at 3.6 ‰) and the CO$_2$ concentration in Antarctic ice cores (at 260 ppmv) represent threshold values for interglacial conditions$^{12}$. Marine Isotope Stage (MIS) boundaries adapted from ref. $^{68}$, boundaries of warmest sub-stages (vertical dashed lines) of interglacials are from ref. $^{24}$. GHG greenhouse-gas concentrations, ITCZ Intertropical Convergence Zone, SST sea-surface temperature.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Convergent cross-mapping of the precipitation [log(Ca/Fe)] and vegetation (tree-pollen) datasets from Tenaghi Philippon. **a** Correlation of cross-mapped against observed values ('CCM skill') as a function of the timeseries (library) length 'L'. The blue line represents the effect of log(Ca/Fe) on tree-pollen, whereas the red line represents the effect of tree-pollen on log(Ca/Fe). Both lines are above the blue and red shaded areas, which represent the 5th to 95th percentiles of the 'CCM skill' for 100 surrogate timeseries from null models, suggesting that convergence is significant for both pairs. As higher 'CCM skill' indicates stronger causality (Methods), the CCM analysis demonstrates a strong bidirectional coupling between precipitation [log(Ca/Fe)] and forest cover [tree-pollen] at Tenaghi Philippon. The black dashed line represents the cross-correlation of the two timeseries. **b** Time displacements maximizing CCM skill corresponding to the causal relationship between log(Ca/Fe) and tree-pollen (marked with stars). As the maximum 'CCM skill' of the effect of log(Ca/Fe) on tree-pollen is negative and closer to 0 (blue line) than the effect of tree-pollen on log(Ca/Fe) (red line), the analysis suggests a real causal and rapid effect of precipitation change on forest cover. Shaded areas represent the 95% confidence interval for predictions from the smoothing model.</figcaption></figure>

While variations in atmospheric CH₄ content, high-northern-
latitude temperature, and the position of the Intertropical Con-
vergence Zone (ITCZ) are considered the main driving forces of pre-
cipitation variability in the Mediterranean region over the past
500 kyrs²⁵,²⁶,³⁰, it remains unclear whether these processes are perti-
nent to the projected precipitation shifts resulting from anthro-
pogenic CO₂ release³. In addition, any influence of these processes on
vegetation dynamics other than through the modulation of moisture
availability remains largely unexplored. To identify which of these
climatic processes may play a critical role in near-future climate sce-
narios, we applied the CCM analysis to the log(Ca/Fe) and tree-pollen
datasets from Tenaghi Philippon and representative climate-proxy
timeseries, including global atmospheric greenhouse-gas (CO₂ and

CH₄) concentrations from Antarctic ice cores¹³,³¹, Neogloboquadrina pachyderma 'sinistral' (NPS) abundances from Ocean Drilling Program (ODP) Site 983 in the North Atlantic as an indicator for high-latitude surface-ocean temperature³², and the Ti/Al and Ba/Al records from ODP Sites 967/968 as indicators for dust input and sapropel formation in the Eastern Mediterranean Sea associated with low-latitude migration of the ITCZ over Africa³³. Moreover, we included orbital precession and obliquity³⁴ in our analysis, although orbital configuration affects climate dynamics on millennial timescales and thus has an insignificant role in near-future climate change¹². Our working hypothesis is that if orbital configuration is the predominant force of moisture variability and vegetation dynamics at Tenaghi Philippon over the past 500 kyrs, it can be assumed that Mediterranean forest biomes will be marginally affected by CO₂-induced climate change in the future. To achieve a mechanistic understanding of the interacting processes, we further assessed the optimal time displacements of the 'CCM skills' to distinguish true causal associations from simple correlations and differentiate direct from indirect causality²⁸ (Supplementary Table 2; Supplementary Fig. 4). The chronological limitations inherent in all the investigated climate-proxy timeseries may bear an effect on the statistical output. Hence, to explore potential effects of age-model uncertainties on the statistical results we also ran the 'CCM' analysis for all the timeseries pairs considering the maximum error range by shifting the Tenaghi Philippon chronology by +/- 5 kyrs (Supplementary Fig. 5).

The CCM analysis yields a mixture of causal associations among the investigated variables, allowing to decipher the extent of high- and low-latitude climate forcing on precipitation and vegetation dynamics at Tenaghi Philippon. Specifically, the combined evaluation of the 'CCM skills' and optimal time displacements suggests that only the global atmospheric CO₂ record has a causal effect on both the tree-pollen and precipitation records from Tenaghi Philippon (Fig. 4; Supplementary Fig. 4). Importantly, because the statistical output of the CCM runs remains stable even when accepting chronological errors of up to 5 kyrs, this effect appears to be independent of age-model uncertainties (Supplementary Fig. 5). Our results further document a strong causality between CH₄ concentration and precipitation at Tenaghi Philippon, and no causal associations of the Tenaghi Philippon datasets with precession, obliquity, ITCZ migration, and dust. Again, these associations appear to be independent of age-model uncertainties. Our data further point to a causal association between precipitation at Tenaghi Philippon and North Atlantic temperature (as reflected by the log(Ca/Fe) and NPS timeseries, respectively; Supplementary Fig. 4); however, our sensitivity tests indicate a potential age-model effect on the 'CCM skills' (Supplementary Fig. 5), leading us to not further considering this correlation. Overall, it appears that only atmospheric CO₂ variability has a rapid causal effect on vegetation dynamics at Tenaghi Philippon. Notably, this effect is realized by means of two discrete pathways. First, there is a direct causal effect of CO₂ variability on vegetation (Fig. 4) through plant physiology processes³⁵. Second, the causal associations document an indirect effect related to climatic processes: the CO₂ variability dictates moisture change at Tenaghi Philippon (Fig. 4), which in turn exerts control on vegetation dynamics (Figs. 3 and 4). Collectively, the CCM analysis suggests an important effect of global atmospheric CO₂ variability on Mediterranean forest biomes during multiple G/I cycles of the Quaternary.

Effect of atmospheric CO₂ on biome shifts

The causality between the gradual variability in atmospheric CO₂
concentration and the abrupt vegetation-regime shifts at Tenaghi
Philippon over the past 500 kyrs attests to the crossing of critical
thresholds that determine the state of Mediterranean forest biomes
that are close to their lower moisture-tolerance limits. Cross-plot
analysis of the tree-pollen and log(Ca/Fe) datasets reveal a non-linear
<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Summary of causal associations of the climate-proxy records as revealed by convergent cross-mapping (CCM) and time displacements of the 'CCM skills'.** Precipitation and forest cover at Tenaghi Philippon based on the log(Ca/Fe) and tree-pollen percentage datasets, respectively; orbital parameters from ref. ³⁴ atmospheric CH₄ (ref. ³¹) and CO₂ (ref. ¹³) concentrations from Antarctic ice cores; position of the Intertropical Convergence Zone (ITCZ) and dust influx based on the Ba/Al and Ti/Al ratios from ODP Site 967, respectively (ref. ³³); North Atlantic surface-water temperature based on Neogloboquadrina pachyderma 'sinistral' percentages from ODP Site 983 (ref. ³²). Significance of causal associations was determined when 'CCM skills' of the proxy timeseries were outside the 5th and 95th percentiles computed from the randomly generated timeseries and the area under curve (AUC) value was higher than 0.8 (Supplementary Fig. 4; Supplementary Table 2). Bold arrows indicate significant causality (AUC > 0.9), whereas dashed arrows indicate weak causality (0.9>AUC > 0.8); crossed-out variables have no causal associations with the datasets from Tenaghi Philippon.</figcaption></figure>

relationship of the two timeseries, unequivocally demonstrating a bimodal distribution of tree-pollen percentages for a wide range of precipitation levels. Forest biomes prevail above a value of 0.8 in the log(Ca/Fe) record, whereas steppe biomes dominate below a log(Ca/Fe) value of 0.25 (Fig. 5). Our analysis shows that shifts between the two distinct vegetation regimes occur at a threshold value of -0.5 in log(Ca/Fe). We explain this to represent the -400 mm/year boundary in annual precipitation that is required for subtropical dry forest biomes to thrive in the Mediterranean region¹,². The results further indicate that both forest and steppe biomes can occur within a range of ±0.3 from the log(Ca/Fe) threshold (Fig. 5). This effect portrays vegetation multiple equilibria with prevailing climate conditions related to migration processes, longevity, and community dynamics³⁶. On balance, our results demonstrate that regime shifts from forest to steppe biomes over the past 500 kyrs have occurred abruptly with the crossing of a precipitation threshold rather than as a gradual decrease in forest cover with precipitation. Most notably, peak precipitation levels and maximum forest cover at Tenaghi Philippon have prevailed persistently under CO₂ concentrations above 260 ppmv as they are realized during full interglacial conditions¹² (Fig. 5). Consequently, our analysis highlights that a precise determination of when forests at Tenaghi Philippon lost resilience within the natural CO₂ range in the past can be achieved by examining the precipitation deficit as portrayed in the decadal-scale-resolution log(Ca/Fe) geochemical data.

Owing to the anthropogenic CO₂ increase to non-analog levels of the past 500 kyrs, it remains unclear whether the relationship of atmospheric CO₂, precipitation and vegetation dynamics as realized during the past four G/I cycles at Tenaghi Philippon can be employed to forecast near-future vegetation-regime shifts in the Mediterranean region. Plants can acclimate rapidly to increasing CO₂ through the reduction of the diffusive stomatal conductance of their leaves until they reach the limits of their phenotypic plasticity, which is twice as high as the CO₂ concentrations predicted for the next century³⁷. Moreover, higher CO₂ concentrations may allow plants to maintain carbon uptake with lower stomatal conductance and, to some extent,

help alleviate drought stress³⁸. As plant physiology hence can compensate to a certain degree the moisture decrease under increasing CO₂, determining the lower limit of moisture availability for maintaining forest resilience is the most critical factor for assessing future ecological stress in the Mediterranean region. Despite the non-analog states between the Quaternary and future CO₂ concentrations, our paired palynological and geochemical record from Tenaghi Philippon can provide an indication for the minimum amount of precipitation required to preserve subtropical dry forest biomes in the Mediterranean region. At Tenaghi Philippon, the warmest interglacial conditions of the past 500 kyrs are persistently marked by log(Ca/Fe) values above 1.5 (Supplementary Fig. 2), and forest biomes are only maintained above a threshold value of 0.8 in the log(Ca/Fe) ratio (Fig. 5). Hence, a minimum decrease in annual precipitation of 40-45 % from full interglacial levels could set off a regime shift from forest to steppe biomes under natural CO₂ concentrations. Climate scenarios for a CO₂-induced temperature rise of 2 °C relative to present in the coming decades predict a decrease of up to 30 % in annual precipitation in the Mediterranean region⁵. This adds up to the natural reduction of water availability in the Mediterranean region in comparison to the wettest interval of the Holocene³⁹. Climate scenarios also predict decreased surface runoff and summer soil moisture in the Mediterranean region during the 21st century⁴⁰. Specifically, the water deficiency will stem from an increase in canopy water consumption with additional vegetation growth along with longer and warmer growing seasons, despite the increased surface resistance to evapotranspiration and vegetation total water use efficiency⁴⁰. Against the background of declining resilience of arid forest biomes to climate variability over the past three decades⁴¹,⁴², the combined natural and human-induced water deficit in the near-future along with the increasing human pressure on Mediterranean ecosystems⁴³ raises concerns for the crossing of the moisture threshold required for the persistence of subtropical dry forest biomes. Our 500-kyr long record from Tenaghi Philippon indicates that an abrupt collapse of forest communities and a rapid switch into a steppe environment would almost be inevitable, urging for sustainable conservation of Mediterranean forests to bolster their resilience to near-future climate change.

## Methods

### Study site and regional setting

The sedimentary archive preserved in the subsurface of the Tenaghi Philippon peatland in the Drama Basin of NW Greece (see ref. ¹¹ for an in-depth review) yields a key terrestrial archive for deciphering global climate dynamics¹²,⁴⁴. The Tenaghi Philippon peatland is situated at an elevation of c. 40 m above sea level surrounded by mountain ranges with altitudes of up to 2232 m (Supplementary Fig. 1). From the early Pleistocene onwards, limnotelmatic conditions prevailed across large parts of the Drama Basin, resulting in a c. 172-m-thick, stratigraphically continuous succession in the southeastern part of the basin that is overwhelmingly composed of fen peat; this succession represents the thickest known peat deposit worldwide¹¹,⁴⁵.

Climatically, Tenaghi Philippon is characterized by a typical Mediterranean seasonal regime modified by continental influence associated with lower winter temperatures and wetter summers as they are characteristic for the NE Mediterranean region. Meteorological data from the nearby town of Drama (Supplementary Fig. 1) record a mean annual temperature of 16.1 °C, and mean July and January temperatures of 27.2 and 5.2 °C, respectively, for the period 1975-2008 (ref. ⁴⁶). Owing to the region's pronounced topography, katabatic flow of cold air masses from the surrounding mountains into the Drama Basin plain causes much lower winter temperatures on the basin floor than at comparable elevations along the nearby coast⁴⁷. Winter temperature is also influenced by anomalously cold and windy episodes connected to southward outbreaks of polar air masses from the Russian High⁴⁸ often causing the development of frost. Modern
<figure><img src="image_9.png" /></figure>

**Fig. 5 |** Forest cover as a function of precipitation at Tenaghi Philippon for the past 500 kyrs. **a** Cross-plot of tree-pollen percentages and log(Ca/Fe) values illustrating the persistence of two distinct and stable vegetation regimes for a wide range of precipitation levels. Rather than a gradual increase in tree cover with increasing precipitation at Tenaghi Philippon, non-linear regime shifts between steppe (low tree-pollen %) and forest (high tree-pollen %) biomes occur above a precipitation threshold (as suggested by the step-like blue line from a generalized additive model [GAM] fit). Purple circles show values at periods of atmospheric CO₂

concentrations above 260 ppmv (ref. ¹³) representing interglacial conditions¹². The gray shaded areas show the density of the data points (darker color for higher density). **b** Histograms of tree-pollen percentages at different log(Ca/Fe) values demonstrating a bimodal distribution with steppe biomes (low tree-pollen %) prevailing below a value of -0.2 in log(Ca/Fe) and forest biomes (high tree-pollen %) dominating above a value of -0.8. Regime shifts between the two biomes occur at a threshold value of -0.5 in log(Ca/Fe).

precipitation in the Drama Basin is connected to the penetration of westerly storm tracks across southern Europe and Mediterranean cyclogenesis (ref. ⁴⁹ and references therein). Precipitation in the basin shows a strong seasonal distribution with a maximum in autumn and winter, and moderately dry summers. Based on the Drama meteorological station, mean annual precipitation is 520 mm, with August and September being the driest (mean: 28 and 29 mm, respectively) and December (mean: 67 mm) the wettest months; 40 % of the annual precipitation occur from November to February, and only 25 % occur from June to August (period 1975–2008; ref. ⁴⁶). The annual precipitation in the surrounding mountain ranges is up to 1100 mm (period 1966–1991; ref. ⁵⁰) and thus twice as high as in the basin itself. A 100 m increase in elevation is associated with a precipitation increase between -30 and -55 mm/year in the western and eastern mountain ranges, respectively⁵⁰.

Owing to its complex geomorphology and broad range of climatic conditions, the Drama Basin and its surroundings are home to a highly diverse flora⁵¹⁻⁵⁴ comprising four main ecozones, i.e., temperate mountain systems, temperate continental forests, subtropical mountain systems, and subtropical dry forests² (Supplementary Fig. 1). Dominant conifer trees are *Abies alba* Mill. and *Pinus nigra* J.F. Arnold, whereas the most common deciduous trees include *Acer spp.*, *Carpinus orientalis* Mill., *Castanea sativa* Mill., *Cornus mas* L., *Corylus avellana* L., *Crataegus monogyna* Jacq., *Fagus sylvatica* L., *Ostrya carpinifolia* Scop., *Quercus frainetto* Ten., *Q. pubescens* Willd., *Q. petraea* (Matt.) Liebl., and *Tilia tomentosa* Moench. Sclerophyllous evergreen shrubs and trees comprise *Arbutus unedo* L., *Erica arborea* L., *Juniperus* sp., *Pistacia terebinthus* L., *Quercus coccifera* L., and *Q. ilex* L. Noteworthy, the mountain ranges surrounding the Drama Basin harbor the Balkan endemic trees *Abies borissii-regis* Mattf., *Acer heldreichii* Boiss., *Acer hyrcanum* Fisch. & C.A. Mey., *Pinus peuce* Griseb., *Salix xanthicola* K. I. Chr. as well as the southernmost occurrences of *Betula pendula* Roth and *Picea abies* (L.) H. Karst. on the Balkan Peninsula⁵⁴.

## Tenaghi Philippon sediment cores

The study is based on a composite sequence developed by the alignment of two parallel drillcores retrieved in a distance of 50 meters from each other (TP-2005 Phi1 and Phi2 cores; 40° 58' 24.0" N, 24° 13' 25.2" E) and one drillcore (TP-2009; 40° 57' 39.5" N, 24° 16' 03.1" E) retrieved 4.4 km south from the TP-2005 core location (Supplementary Fig. 1). Correlations between the cores were established based on a set of independent proxies, i.e., cryptotephra identification⁵⁵, X-ray fluorescence (XRF) core-scanning data (Supplementary Fig. 6), and palynostratigraphy.

We here present information from the interval 4–91 m composite depth of the Tenaghi Philippon archive spanning the past c. 500 kyrs continuously (see ‘Chronology’). The uppermost 4 m of the composite record have been excluded from this study because the peat sequence is disturbed by anthropogenic activities related to the drainage of the peatland and agriculture. Lithologically, the succession consists predominantly of fen peat. Specifically, the intervals 4.0–58.5 m and 67.0–82.0 m are made up of fen peat, the intervals 58.5–67.0 m and 84.0–91.0 are dominated by organic-rich mud, and the interval 82.5–84.0 m consists of lake marl.

## Peat geochemistry

The geochemical composition of the inorganic components in the Tenaghi Philippon peat cores is predominantly influenced by mineral-rich groundwater and, to a lesser extent, by spring and surface waters flowing into the Drama Basin from the surrounding mountain ranges⁴⁵,⁵⁶. The subsurface of the Drama Basin contains a shallow unconfined aquifer with a depth of ~1.5 m at Tenaghi Philippon⁴⁵,⁵⁶. Based on hydrological investigations, the groundwater recharge predominantly occurs from carbonate terrains (i.e., limestones and marbles) in the surrounding mountain ranges (Supplementary Fig. 1) accounting for >50 % of the total precipitation⁵⁰. Rainwater from the surrounding mountains reaches the aquifer in the Drama Basin within few weeks⁵⁰. To a lesser extent, groundwater recharge also occurs
within the Drama Basin due to the positive precipitation/evapo-
transpiration balance³⁰ even during drier-than-average conditions as
depicted in meteorological data since 1975 (ref. ⁴⁶). As such, the
groundwater in the Drama Basin is strongly enriched in Ca as a result of
carbonate dissolution. In addition, relatively high Fe and Si con-
centrations are contained in the groundwater at the basin margins
originating from alluvial fan sediments deposited close to the moun-
tain ranges⁵⁰,⁵⁶. Considering that a dominance of mechanical over
chemical erosion of carbonates in the Eastern Mediterranean region
only emerges under hyper-arid conditions at a mean annual pre-
cipitation <200 mm (ref. ⁵⁷), the carbonate content in the drillcores
from the Drama Basin can be firmly attributed to result from carbonate
dissolution related to the activity of the karst system in the sur-
rounding mountain ranges.

Detailed coal petrographic and geochemical characterization of the uppermost -10 m of the Tenaghi Philippon peatland suggests distinct differences in inorganic-geochemical composition in the peat deposited during the Last Glacial *versus* that deposited during the Holocene⁵⁸⁻⁶⁰: the Holocene peat is rich in Ca, whereas peat from the Last Glacial is depleted in Ca, but enriched in Fe, K, Ti, and Si; moreover, the degree of peat humification and gelification of the organic matter, which is enhanced under low-oxygen conditions and during periods of high water content in the sediments⁶¹, is high in the Holocene and low in the Last Glacial sediments. In addition, mineralogical analysis across the Last Interglacial interval of the core material shows highest contents of Ca-bearing minerals in the peat during the time of sapropel formation in the Eastern Mediterranean Sea, which is related to increased precipitation in the northern Mediterranean borderlands²⁵,⁶²,⁶³. Highest contents of siliciclastics are registered during the transition phase from full interglacial conditions to the subsequent glacial, resulting from a higher erosion of surface sediments due to increased forest openness⁶⁴. Together, these results suggest that (i) increased precipitation during interglacials leads to a highly active karstic system, which in turn results in a high groundwater table, enhanced telmatic conditions at Tenaghi Philippon, and high Ca content of the peat; (ii) decreased precipitation during glacials leads to a less active karstic system causing a lowering of the groundwater table, and enhanced siliciclastic sedimentation from the basin and the peatland-exposed margins related to surface soil erosion due to low vegetation cover.

To confirm the observation from the uppermost -10 m of the core material that the interglacial peat at Tenaghi Philippon is enriched in Ca, whereas the glacial peat is enriched in Fe, K, Ti, and Si, we examined the variability in elemental composition over the past four G/I cycles by applying a PCA analysis to the elemental dataset derived by XRF core scanning. The PCA yields two major components that account for 73.1% of the total variance of the elemental dataset (Supplementary Fig. 7). The first component (PC1) explains 55.7% of the total variance and is marked by positive loadings of Al, Fe, K, Mn, Si, and Ti, and negative loadings of Ca and S. The second component (PC2) explains 17.4% of the total variance. Its positive pole is primarily driven by S and Fe, and to a lesser extent by Ti, whereas its negative pole is primarily driven by Ca and Mn. The two principal components can be explained to differentiate between intervals of high groundwater table rich in carbonates (negative scores in both PCs) and intervals of low groundwater table rich in siliciclastics due to enhanced clastic sedimentation from the basin margins (positive scores in PC1) and in sulfur due to organic matter decomposition under anaerobic conditions (positive scores in PC2). In light of the significance of Ca and Fe in reconstructing karstic system dynamics at Tenaghi Philippon, and considering that these elements yield the highest XRF counts throughout the studied sequence, we use the log(Ca/Fe) ratio as an indicator of moisture variability with high(low) values corresponding to high(low) groundwater table and, by extension, high(low) precipitation.

Chronology

The age model is based on combined acceleration mass spectrometry (AMS) ¹⁴C dating, tephrochronological and cyclostratigraphic techniques. Firstly, 67 available ¹⁴C dates⁶⁵,⁶⁶ and ten independently dated (crypto)tephra layers out of a total of 57 detected in the uppermost 91 m composite depth of the Tenaghi Philippon cores⁵⁵,⁶⁵,⁶⁷ were employed to establish chronostratigraphic anchor points and calculate average sedimentation rates. Subsequently, cyclostratigraphic analysis was performed on the log(Ca/Fe) record, which yielded pronounced cyclic variability in the 12–16 m range and an additional, albeit minor, cycle within the 5–6 m period band (Supplementary Figs. 8 and 9). Considering an average sedimentation rate of ~15 cm/kyr for the Philippi peatland during the past 500 kyrs, the 12–16 m cycle corresponds to a period of 80–106 kyrs, and the 5–6 m cycle to a period of 33–40 kyrs. These periodicities are strongly reminiscent of orbital (short) eccentricity and obliquity, respectively (ref. ³⁴), suggesting that the cycles in the log(Ca/Fe) record represent G/I variability paced by orbital eccentricity during the Middle and Late Pleistocene⁴⁴,⁶⁸, with interglacials corresponding to high log(Ca/Fe) values and vice versa (see also ‘Peat geochemistry’).

Subsequently, a continuous chronology was developed by tuning the XRF-based log(Ca/Fe) record to the ice-sheet model of ref. ⁶⁹, assigning maxima (minima) in the log(Ca/Fe) record to minima (maxima) in global ice volume (Supplementary Fig. 9). The tuning approach presupposes an in-phase behavior between changes in the log(Ca/Fe) record and the moisture budget at Tenaghi Philippon during G/I cycles. This assumption is confirmed by the in-phase relationship between the log(Ca/Fe) record and the ice-sheet model during the last ~40 kyrs (corresponding to the uppermost 13 m in the Tenaghi Philippon sequence), which is robustly dated by ¹⁴C and tephrochronology⁶⁵,⁶⁶. Considering the error estimates of the ice-sheet model towards radiometric ages for the past ~650 kyrs⁶⁹, the age-model uncertainties for the Tenaghi Philippon chronology are <5 kyrs. The robustness of the age model is confirmed by the independently dated (crypto)tephra layers detected in the Tenaghi Philippon cores⁵⁵,⁶⁵,⁶⁷ all of which fall within the 5-kyrs error uncertainty (Supplementary Fig. 9).

Palynological analysis

Our study is based on 2347 palynological samples from the 4–91 m composite core depth interval of the Tenaghi Philippon archive. Sampling resolution is between 2 and 6 cm, which based on the age model for the composite core yields a mean temporal resolution of 209 years. For the core intervals corresponding to MIS 6–7 and MIS 9–13, a total of 1183 new samples were counted. For MIS 1–5 and MIS 8, our analysis is based on the raw data of ref. ¹⁵,¹⁷,¹⁹ that were augmented by additional pollen counts in order to increase their temporal resolution and statistical fidelity.

Palynological preparation comprised sediment freeze-drying, weighing, spiking with *Lycopodium* spores, treatment with HCl (10 %), NaOH (10 %), HF (when necessary; 40 %), heavy-liquid separation with Na₂WO₄ × 2H₂O (when necessary), acetolysis, sieving through a 7 µm mesh, and slide preparation using glycerine jelly. Between 53 and 3661 (mean: 500) pollen grains were counted per sample excluding pollen from aquatic plants and fern spores. Pollen percentages were calculated on the number of pollen grains from terrestrial plants excluding Cyperaceae and Poaceae because of their natural over-representation at Tenaghi Philippon¹⁴. The percentages of Cyperaceae and Poaceae, aquatic taxa, fern spores and algae were calculated on the basis of the main sum plus the counts for the individual taxa.

X-ray florescence core scanning

The TP-2005 and TP-2009 cores were scanned with an AVAATECH (GEN-4) XRF core scanner equipped with an OXFORD 'Neptune 5200' series 100 W X-ray source with a Rhodium anode at the Institute of Earth
Sciences, Heidelberg University. Prior to scanning, all core sections were warmed to room temperature, and the core surfaces were smoothed and covered with 4-µm-thick Ultralene® foil to avoid contamination of the detector window and desiccation of the cores. Scanning was carried out at a 10 kV energy level without filter using a current of 200 mA, a counting time of 10 s, and slit sizes of 5 mm downcore and 5 or 10 mm crosscore. Data acquisition was via a RAYSPEC SiriusSD 65 mm² Silicon Drift Detector (Model 878-0616B), and a BRIGHTSPEC Topaz-X Multi-channel Analyzer. Processing of the X-ray spectra was performed using the BRIGHTSPEC bAxilBatch software (Version 1.4). The individual element counts were normalized to the total counts of all processed elements for each measurement excluding Ag and Rh because these elements are biased by the signal generation as they are included in the beam collimator of the detector and the X-ray source, respectively. In addition, log elemental ratios were used to eliminate non-linear matrix effects and constant-sum constraints<sup>70,71</sup>. Based on the age model, the average temporal resolution of the XRF data is 32 years.

## Convergent cross-mapping

To examine potential causalities between the Tenaghi Philippon pollen and log(Ca/Fe) datasets and other paleoclimate proxy records, we applied convergent cross-mapping (CCM) that allows us to distinguish causality from simple correlation in multivariate timeseries of deterministic dynamical systems<sup>27</sup>. Representing a powerful method to infer causality among complex Earth system processes<sup>72</sup>, the CCM analysis has been successfully applied to recent ecological<sup>73</sup> and climate timeseries<sup>74</sup>, as well as Antarctic paleoclimate datasets<sup>28</sup>. The CCM analysis is based on Takens' Theorem<sup>75,76</sup>, which states that if time-series variable A influences timeseries variable B in a dynamical system, then the state A(t) can be reliably predicted from the state of B(t). On this basis, CCM allows us to assess the extent to which the record of the affected variable B reliably estimates states of the causal variable A. This is achieved by comparing Pearson's correlation coefficient (ρ) between the predicted and observed values of A(t) (termed 'CCM skill' after ref. 27). As convergence is used to distinguish between direct and indirect causation from simple correlation, the cross-mapped estimates improve in 'CCM skill' with the length of the timeseries that is used to predict A from B.

CCM was performed using the ccm function of the rEDM package (Version 1.2.3.; ref. 77) in R (Version 3.6.3; ref. 78). All datasets were interpolated linearly to equidistant 1 kyr time steps prior to CCM analysis. Optimal CCM prediction parameters (i.e., the combination of embedding dimension [*E*] and time delay [*τ*] from a range of *E* [1–10] and *τ* [1–4] that yielded the highest forecast skill *ρ*) are derived from simplex projection using the *simplex* function included in the rEDM package, utilizing the entire timeseries as both library and predictor variable. Significance was determined by cross-mapping timeseries A against 100 bootstrapped and randomly phase-shifted surrogate iterations of timeseries B, and vice versa<sup>79</sup>. The sensitivity of the method was assessed by calculating the 'CCM skill' for every timeseries against 100 surrogate timeseries, utilizing five combinations of *E* and *τ* that yielded the highest forecast skill *ρ* as derived from the simplex projection. The results were evaluated using receiver operating characteristics (ROC) curves and, more specifically, the area under curve (AUC). Both ROC and AUC calculations were performed using the ROCR package (Version 1.0-11; ref. 80) in R (Version 3.6.3; ref. 78).

To determine the true direction of causality and to identify a potential occurrence of dynamic synchrony (i.e., bidirectional CCM convergence despite unidirectional causation)<sup>75</sup> that frequently occurs in strongly coupled dynamic systems such as the Earth's climate system, we also applied time-delayed CCM, employing a prediction horizon of 10 time steps (1 kyr each) in both directions<sup>28,29</sup>. This approach is based on the principle that in dynamic systems with delayed responses to extrinsic forcing, the optimal 'CCM skill' for a variable should occur with an asymmetrical time delay for both variables. In contrast, dynamic

synchrony (i.e., variables A and B are forced by a third variable Z) should occur nearly instantaneous and with symmetrical time displacement for both variables<sup>28,29</sup>. The time displacements can be utilized further to reveal delayed onsets of cause and effect in dynamic systems. A maximum 'CCM skill' close to a time displacement of 0 (or to the embedding lag *τ*) implies a (nearly) instantaneous effect. The further the maximum 'CCM skill' diverges from 0, the slower is the response time of said effect. Negative time displacements that lead to higher CCM skill would imply a real causal effect as long as there is a difference in the optimal lag between the CCM skills of one variable to the other. In contrast, a positive time displacement would imply that the affected variable is better at predicting future values of the driving variable and can be used to rule out causality<sup>29</sup>. Time-delayed CCM was performed in R (Version 3.6.3) using the *ccm* function by setting the *tp* parameter to non-zero (e.g., ±10) (ref. 29). The time displacement yielding the maximum 'CCM skill' was determined by smoothing (*x* - *y*<sup>2</sup>, *n* = 100, *span* = 1) the 'CCM skill' of 100 bootstrapped library sets of size 100 as a function of time displacement (±10) (ref. 28). Values towards the limits of the time displacement were not considered due to high uncertainties of the smoothing function.

## Generalized additive models

Generalized additive models (GAM) is a statistical modeling technique utilizing non-linear smoothing functions to model non-linearities in datasets<sup>81</sup>. GAM values were calculated using the *gam* function from the *mgcv* package (Version 1.8-31; ref. 82) in R (Version 3.6.3; ref. 78) using cubic regression splines (*cr*) and a smoothing dimension of *k* = 25.

## Principal component analysis

PCA was applied to the elemental data derived through XRF core scanning in order to examine the relationships among selected elements with robust signals across the entire sequence (i.e., Al, Ca, Fe, K, Mn, P, S, Si, and Ti). The analysis was performed using the function *prcomp* implemented in the R base package (Version 3.6.3; ref. 78), which utilizes a singular value decomposition function to introduce unit variances.

## Rolling correlation

Rolling correlations were used to identify changes in the correlation between two times series through time. The method utilizes a rolling window in which the correlation is calculated independently of the surrounding data, with the adjusting window size yielding information on shorter (small window size) or longer (big window size) timeframes. The analysis was performed in R (Version 3.6.3; ref. 78) using the *surrogateCor* function from the *astrochron* package (Version 1.0; ref. 83), with window sizes ranging between 20 and 50 time steps, increasing in increments of 10, and with every time step representing 1 kyr.

## Spectral analysis

To identify orbital frequencies in the log(Ca/Fe) record we applied REDFIT analyses<sup>84</sup> as implemented in PAST software. The significance levels at 80%, 90%, and 95% were tested using a Monte-Carlo simulation. To trace the spatial development of the identified dominant frequencies in the depth domain we further applied evolutive harmonic spectral analysis (EHA) using the multitaper method<sup>85</sup> as implemented in the function *eha* in the *astrochron* R software package<sup>83</sup>. The log(Ca/Fe) dataset was resampled at 0.01 m levels prior to the application of the EHA.

## Data availability

The pollen and XRF data generated in this study have been deposited in the PANGAEA database under accession codes (https://doi.org/10.1594/PANGAEA.943592, https://doi.org/10.1594/PANGAEA.943593). Source data are provided with this paper.
## Code availability

The fully reproducible code for the convergent cross-mapping (CCM) analysis is available on 'GitHub'(https://zenodo.org/badge/latesttdoi/452267408).

## References

1. San-Miguel-Ayanz, J. et al. (eds.). *European Atlas of Forest Tree Species* (Publication Office of the European Union, Luxembourg, 2016).

2. Food and Agriculture Organization of the United Nations. *Global Forest Resources Assessment 2000: Main Report* (FAO, 2000).

3. Nolan, C. et al. Past and future global transformation of terrestrial ecosystems under climate change. *Science* **361**, 920–923 (2018).

4. Schröter, D. et al. Ecosystem service supply and vulnerability to global change in Europe. *Science* **310**, 1333–1337 (2005).

5. Cramer, W. et al. Climate change and interconnected risks to sustainable development in the Mediterranean. *Nat. Clim. Change* **8**, 972–980 (2018).

6. Gazol, A. et al. Forest resilience to drought varies across biomes. *Glob. Change Biol.* **24**, 2143–2158 (2018).

7. Tramblay, Y. et al. Challenges for drought assessment in the Mediterranean region under future climate scenarios. *Earth Sci. Rev.* **210**, 103348 (2020).

8. Guiot, J. & Cramer, W. Climate change: The 2015 Paris Agreement thresholds and Mediterranean basin ecosystems. *Science* **354**, 465–468 (2016).

9. Newbold, T., Oppenheimer, P., Etard, A. & Williams, J. J. Tropical and Mediterranean biodiversity is disproportionately sensitive to land-use and climate change. *Nat. Ecol. Evol.* **4**, 1630–1638 (2020).

10. Willis, K. J., Bailey, R. M., Bhagwat, S. A. & Birks, H. J. B. Biodiversity baselines, thresholds and resilience: testing predictions and assumptions using palaeoecological data. *Trends Ecol. Evol.* **25**, 583–591 (2010).

11. Pross, J. et al. The 1.35-Ma-long terrestrial climate archive of Tenaghi Philippon, northeastern Greece: Evolution, exploration, and perspectives for future research. *Newsl. Stratigr.* **48**, 253–276 (2015).

12. Past Interglacials Working Group of PAGES. Interglacials of the last 800,000 years. *Rev. Geophys.* **54**, 162–219 (2016).

13. Bereiter, B. et al. Revision of the EPICA Dome C CO₂ record from 800 to 600 kyr before present. *Geophys. Res. Lett.* **42**, 542–549 (2015).

14. Tzedakis, P. C., Hooghiemstra, H. & Pälike, H. The last 1.35 million years at Tenaghi Philippon: revised chronostratigraphy and long-term vegetation trends. *Quat. Sci. Rev.* **25**, 3416–3430 (2006).

15. Müller, U. C. et al. The role of climate in the spread of modern humans into Europe. *Quat. Sci. Rev.* **30**, 273–279 (2011).

16. Tzedakis, P. C. et al. Ecological thresholds and patterns of millennial-scale climate variability: the response of vegetation in Greece during the last glacial period. *Geology* **32**, 109–112 (2004).

17. Fletcher, W. J., Müller, U. C., Koutsodendris, A., Christanis, K. & Pross, J. A centennial-scale record of vegetation and climate variability from 312 to 240 ka (Marine Isotope Stages 9c-a, 8 and 7e) from Tenaghi Philippon, NE Greece. *Quat. Sci. Rev.* **78**, 108–125 (2013).

18. Pross, J. et al. Massive perturbation in terrestrial ecosystems of the Eastern Mediterranean region associated with the 8.2 kyr B.P. climatic event. *Geology* **37**, 887–890 (2009).

19. Milner, A. M. et al. Vegetation responses to abrupt climatic changes during the Last Interglacial Complex (Marine Isotope Stage 5) at Tenaghi Philippon, NE Greece. *Quat. Sci. Rev.* **154**, 169–181 (2016).

20. Ardenghi, N. et al. Temperature and moisture variability in the eastern Mediterranean region during Marine Isotope Stages 11–10 based on biomarker analysis of the Tenaghi Philippon peat deposit. *Quat. Sci. Rev.* **225**, 105977 (2019).

21. Médail, F. & Diadema, K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. *J. Biogeogr.* **36**, 1333–1345 (2009).

22. Tzedakis, P. C. Vegetation change through glacial-interglacial cycles: a long pollen sequence perspective. *Philos. Trans. R. Soc. Lond. Ser. B: Biol. Sci.* **345**, 403–432 (1994).

23. Donders, T. et al. 1.36 million years of Mediterranean forest refugium dynamics in response to glacial-interglacial cycle strength. *Proc. Natl Acad. Sci. USA* **118**, e2026111118 (2021).

24. Tzedakis, P. C. et al. Can we predict the duration of an interglacial? *Clim. Past* **8**, 1473–1485 (2012).

25. Wagner, B. et al. Mediterranean winter rainfall in phase with African monsoons during the past 1.36 million years. *Nature* **573**, 256–260 (2019).

26. Tzedakis, P. C. Seven ambiguities in the Mediterranean palaeoenvironmental narrative. *Quat. Sci. Rev.* **26**, 2042–2066 (2007).

27. Sugihara, G. et al. Detecting causality in complex ecosystems. *Science* **338**, 496–500 (2012).

28. van Nes, E. H. et al. Causal feedbacks in climate change. *Nat. Clim. Change* **5**, 445–448 (2015).

29. Ye, H., Deyle, E. R., Gilarranz, L. J. & Sugihara, G. Distinguishing time-delayed causal interactions using convergent cross mapping. *Sci. Rep.* **5**, 14750 (2015).

30. Tzedakis, P. C., Pälike, H., Roucoux, K. H. & de Abreu, L. Atmospheric methane, southern European vegetation and low-mid latitude links on orbital and millennial timescales. *Earth Planet. Sci. Lett.* **277**, 307–317 (2009).

31. Louergue, L. et al. Orbital and millennial-scale features of atmospheric CH₄ over the past 800,000 years. *Nature* **453**, 383–386 (2008).

32. Barker, S. et al. Early interglacial legacy of deglacial climate instability. *Paleoceanogr. Paleoclimatol.* **34**, 1455–1475 (2019).

33. Grant, K. M. et al. A 3 million year index for North African humidity/aridity and the implication of potential pan-African Humid periods. *Quat. Sci. Rev.* **171**, 100–118 (2017).

34. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. *Astron. Astrophys.* **428**, 261–285 (2004).

35. Hetherington, A. M. & Woodward, F. I. The role of stomata in sensing and driving environmental change. *Nature* **424**, 901–908 (2003).

36. Svenning, J. C. & Sandel, B. Disequilibrium vegetation dynamics under future climate change. *Am. J. Bot.* **100**, 1266–1286 (2013).

37. de Boer, H. J. et al. Climate forcing due to optimization of maximal leaf conductance in subtropical vegetation under rising CO₂. *Proc. Natl Acad. Sci. USA* **108**, 4041–4046 (2011).

38. Zhu, Z. et al. Greening of the Earth and its drivers. *Nat. Clim. Change* **6**, 791–795 (2016).

39. Mauri, A., Davis, B. A. S., Collins, P. M. & Kaplan, J. O. The climate of Europe during the Holocene: a gridded pollen-based reconstruction and its multi-proxy evaluation. *Quat. Sci. Rev.* **112**, 109–127 (2015).

40. Mankin, J. S., Seager, R., Smerdon, J. E., Cook, B. I. & Williams, A. P. Mid-latitude freshwater availability reduced by projected vegetation responses to climate change. *Nat. Geosci.* **12**, 983–988 (2019).

41. Forzieri, G., Dakos, V., McDowell, N. G., Ramdane, A. & Cescatti, A. Emerging signals of declining forest resilience under climate change. *Nature* **608**, 534–539 (2022).

42. Smith, T., Traxl, D. & Boers, N. Empirical evidence for recent global shifts in vegetation resilience. *Nat. Clim. Change* **12**, 477–484 (2022).

43. Kéfi, S. et al. Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems. *Nature* **449**, 213–217 (2007).

44. Lang, N. & Wolff, E. W. Interglacial and glacial variability from the last 800 ka in marine, ice and terrestrial archives. *Climate* **7**, 361–380 (2011).
45. Christianis, K. Ein Torf erzählt die Geschichte seines Moores. *Telma* **13**, 19–32 (1983).

46. Greek National Meteorological Service. *Meteorological data from station WMO No 16607 (Drama)*, accessed on April upon request (Greek National Meteorological Service, 2020).

47. Flocas, A. A. & Angouridakis, V. E. Extreme values analysis of air temperature over Greece. *Arch. für. Meteorologie, Geophysik und Bioklimatologie, Ser. B* **27**, 47–57 (1979).

48. Saaroni, H., Bitan, A., Alpert, P. & Ziv, B. Continental polar outbreaks into the Levant and eastern Mediterranean. *Int. J. Climatol.* **16**, 1175–1191 (1996).

49. Lionello, P. (ed.) *The Climate of the Mediterranean Region: from the past to the Future* (Elsevier, 2012).

50. Panilas, G. *Ground Water Problems during Lignite Exploitations in the Open Pit Mines: the Case of Drama's Lignite Deposit*. PhD thesis (University of Patras, Greece, 1998).

51. Karagiannakidou, V. & Kokkini, S. Vegetation and altitudinal zones of Mount Menikion, NE Greece. *Bot. Hel.* **98**, 149–160 (1988).

52. Karagiannakidou, V., Konstantinou, M. & Papademetriou, K. Floristic and phytogeographical research on the upper montane and the subalpine grassland flora of East Macedonia, Greece. *Feddes Repert.* **106**, 193–213 (1995).

53. Tsiripidis, I. & Athanasiadis, N. Contribution to the knowledge of the vascular flora of NE Greece: floristic composition of the beech (*Fagus sylvatica* L.) forests in the Greek Rodopi. *Willdenowia* **33**, 273–297 (2003).

54. Phitos, D., Constantinidis, T. & Kamari, G. (Eds.) *The Red Data Book Of Rare And Threatened Plants Of Greece* (Hellenic Botanical Society, Patras, Greece, 2009).

55. Vakhrameeva, P. et al. Eastern Mediterranean volcanism during marine isotope stages 9 to 7e (335–235 ka): insights based on cryptotephra layers at Tenaghi Philippon, Greece. *J. Volcanol. Geotherm. Res.* **380**, 31–47 (2019).

56. Georgakopoulos, A. et al. The Drama basin water: quality and peat/lignite interaction. *Environ. Geol.* **41**, 121–127 (2001).

57. Ryb, U. et al. Styles and rates of long-term denudation in carbonate terrains under a Mediterranean to hyper-arid climatic gradient. *Earth Planet. Sci. Lett.* **406**, 142–152 (2014).

58. Christianis, K., Georgakopoulos, A., Fernández-Turiel, J. L. & Papazisimou, S. Geological factors influencing the concentration of trace elements in the Philippi peatland, eastern Macedonia, Greece. *Int. J. Coal Geol.* **36**, 295–313 (1998).

59. Kalaitzidis, S. & Christianis, K. Mineral matter in the Philippi peat in relation to peat/lignite-forming conditions in Greece. *Energy Sources* **24**, 69–81 (2002).

60. Kalaitzidis, S., Christianis, K., Georgakopoulos, A., Fernández-Turiel, J. L. & Papazisimou, S. Influence of geological conditions during peat accumulation on trace element affinities and their behavior during peat combustion. *Energy Fuels* **16**, 1476–1482 (2002).

61. Teichmüller, M. The genesis of coal from the viewpoint of coal petrology. *Int. J. Coal Geol.* **12**, 1–87 (1989).

62. Kotthoff, U. et al. Climate dynamics in the borderlands of the Aegean Sea during formation of sapropel S1 deduced from a marine pollen record. *Quat. Sci. Rev.* **27**, 832–845 (2008).

63. Rohling, E. J., Marino, G. & Grant, K. M. Mediterranean climate and oceanography, and the periodic development of anoxic events (sapropels). *Earth Sci. Rev.* **143**, 62–97 (2015).

64. Milner, A. M. et al. Enhanced seasonality of precipitation in the Mediterranean during the early part of the Last Interglacial. *Geology* **40**, 919–922 (2012).

65. Wulf, S. et al. The marine isotope stage 1–5 cryptotephra record of Tenaghi Philippon, Greece: Towards a detailed tephrostratigraphic framework for the Eastern Mediterranean region. *Quat. Sci. Rev.* **186**, 236–262 (2018).

66. Staff, R. A. et al. Reconciling the Greenland ice-core and radio-carbon timescales through the Laschamp geomagnetic excursion. *Earth Planet. Sci. Lett.* **520**, 1–9 (2019).

67. Vakhrameeva, P. et al. The cryptotephra record of the Marine Isotope Stage 12 to 10 interval (460–335 ka) at Tenaghi Philippon, Greece: Exploring chronological markers for the Middle Pleistocene of the Mediterranean region. *Quat. Sci. Rev.* **200**, 313–333 (2018).

68. Lisiecki, L. E. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed benthic $\delta^{18}$O records. *Paleoceanography* **20**, PA1003 (2005).

69. Imbrie, J. & Imbrie, J. Z. Modeling the climatic response to orbital variations. *Science* **207**, 943–953 (1980).

70. Weltje, G. J. et al. Prediction of geochemical composition from XRF core scanner data: A new multivariate approach including automatic selection of calibration samples and quantification of uncertainties. In: *Micro-XRF Studies Of Sediment Cores. Developments In Paleoenvironmental Research*. (eds. Croudace, I. W. & Rothwell, R. G.) Vol. 17, 507–534 Springer, Dordrecht (Springer, Dordrecht, 2015).

71. Kern, O. A. et al. X-ray fluorescence core scanning yields reliable semiquantitative data on the elemental composition of peat and organic-rich lake sediments. *Sci. Total Environ.* **697**, 134110 (2019).

72. Runge, J. et al. Inferring causation from time series in Earth system sciences. *Nat. Commun.* **10**, 2553 (2019).

73. Ushio, M. et al. Fluctuating interaction network and time-varying stability of a natural fish community. *Nature* **554**, 360–363 (2018).

74. Dima, M., Nichita, D. R., Lohmann, G., Ionita, M. & Voiculescu, M. Early-onset of Atlantic Meridional Overturning Circulation weakening in response to atmospheric $CO_2$ concentration. *npj Clim. Atmos. Sci.* **4**, 27 (2021).

75. Takens, F. Detecting strange attractors in turbulence. In: *Dynamical Systems and Turbulence Warwick 1980. Lecture Notes in Mathematics* (eds. Rand, D. & Young, L.S.). Vol. **898**, 366–381 (Springer, Berlin/Heidelberg, 1981).

76. Deyle, E. R. & Sugihara, G. Generalized theorems for nonlinear state space reconstruction. *PLoS ONE* **6**, e18295 (2011).

77. Ye, H. et al. rEDM: Applications of Empirical Dynamic Modeling from time series. *Zenodo* https://doi.org/10.5281/zenodo.1935847 (2018).

78. R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, https://www.r-project.org/ (R Core Team, 2020).

79. Ebisuzaki, W. A method to estimate the statistical significance of a correlation when the data are serially correlated. *J. Clim.* **10**, 2147–2153 (1997).

80. Sing, T., Sander, O., Beerenwinkel, N. & Lengauer, T. ROCR: visualizing classifier performance in R. *Bioinformatics* **21**, 3940–3941 (2005).

81. Hastie, T. J. & Tibshirani, R. J. *Generalized Additive Models* (Chapman & Hall, CRC Press, 1990).

82. Wood, S. N. *Generalized Additive Models: An Introduction With R, Second Edition* (Chapman & Hall, CRC Press/Taylor & Francis Group, 2017).

83. Meyers, S. R. *astrochron: An R Package for Astrochronology*. https://cran.r-project.org/package=astrochron (2014).

84. Schulz, M. & Mudelsee, M. REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. *Comput. Geosci.* **28**, 421–426 (2002).

85. Thomson, D. J. Spectrum estimation and harmonic analysis. *Proc. IEEE* **70**, 1055–1096 (1982).

## Acknowledgements

We thank I. Glass, S. Kalaitzidis, S. Liner, A. Neugebauer, and G. Siavalas for technical assistance. This study was supported by the German Science Foundation (DFG; grants PR651/3 and PR651/19 to JP, and KO4960/3 and KO4960/5 to AK), the Biodiversity and Climate Research Center
(BIK-F) of the Hessian Initiative for Scientific and Economic Excellence (LOEWE), and the Wilhelm Schuler Foundation.

**Author contributions**

A.K. and J.P. designed the study, coordinated the discussion and interpretation of the proxy datasets, and wrote the paper; W.J.F., M.K., U.K., A.K., A.M.M., and U.C.M. contributed palynological data; A.K. developed the XRF record; V.D., O.A.K., and L.K. carried out the statistical analyses; P.V. and S.W. provided tephrochronological data; S.K.-B. and A.K. developed the age model; L.K. generated the maps; K.C., A.K., U.C.M., J.P., and G.S. contributed to the planning and implementation of the coring campaigns. All co-authors provided comments and suggestions to the manuscript.

**Funding**

Open Access funding enabled and organized by Projekt DEAL.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-37388-x.

**Correspondence** and requests for materials should be addressed to Andreas Koutsodendris.

**Peer review information** *Nature Communications* thanks Eelco J Rohling and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023