1 2 3 4 5 6 7 8 9 10 11 12 13 Mediterranean winter rainfall in phase with African monsoon during past 1.36 million years Bernd Wagner1*†, Hendrik Vogel2†, Alexander Francke1,3, Tobias Friedrich4, Timme Donders5, Jack H. Lacey6, Melanie J. Leng6,7, Eleonora Regattieri8,9, Laura Sadori10, Thomas Wilke11, Giovanni Zanchetta8, Christian Albrecht11, Adele Bertini12, Nathalie CombourieuNebout13, Aleksandra Cvetkoska5, 11, Biagio Giaccio14, Andon Grazhdani15, Torsten Hauffe11, Jens Holtvoeth16, Sebastien Joannin17, Elena Jovanovska11, Janna Just1,18, Katerina Kouli19, Ilias Kousis20, Andreas Koutsodendris20, Sebastian Krastel21, Niklas Leicher1, Zlatko Levkov22, Katja Lindhorst21, Alessia Masi10, Martin Melles1, Anna M. Mercuri23, Sebastien Nomade24, Norbert Nowaczyk25, Konstantinos Panagiotopoulos1, Odile Peyron17, Jane M. Reed26, Leonardo Sagnotti27, Gaia Sinopoli10, Björn Stelbrink11, Roberto Sulpizio28, 29, Axel Timmermann30,31, Slavica Tofilovska22, Paola Torri32, Friederike Wagner-Cremer5, Thomas Wonik33, Xiaosen Zhang34 14 1 17 18 3 21 22 5 Institute of Geology and Mineralogy, University of Cologne, Cologne, Germany. 15 16 2 19 20 4 International Pacific Research Center, University of Hawaii at Manoa, Honolulu, Hawaii, USA. 23 6 Institute of Geological Sciences & Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland. School of Earth, Atmospheric, and Life Science, University of Wollongong, Wollongong, Australia. Palaeoecology, Department of Physical Geography, Utrecht University, Utrecht, The Netherlands. National Environmental Isotope Facility, British Geological Survey, Nottingham, UK. 24 25 7 26 8 29 10 27 28 Centre for Environmental Geochemistry, School of Biosciences, University of Nottingham, UK. Dipartimento di Scienze della Terra, University of Pisa, Pisa, Italy. 9 Institute of Earth Sciences and Earth Resources-Italian National Research Council (IGGCNR), Pisa, Italy. Dipartimento di Biologia Ambientale, Università di Roma "La Sapienza", Rome, Italy. 30 31 11 32 12 35 14 Istituto di Geologia Ambientale e Geoingegneria – CNR, Rome, Italy. Department of Animal Ecology & Systematics, Justus Liebig University Giessen, Giessen, Germany. Dipartimento di Scienze della Terra, Università di Firenze, Firenze, Italy. 13 33 34 CNRS UMR 7194, Muséum National d’Histoire Naturelle, Institut de Paléontologie Humaine, Paris, France. 15 Faculty of Geology and Mineralogy, University of Tirana, Albania. 37 16 School of Chemistry, University of Bristol, Bristol, UK. 40 18 36 38 39 17 CNRS UMR 5554, Institut des Sciences de l’Evolution de Montpellier, Université de Montpellier, Montpellier, France. Fachbereich Geowissenschaften, Universität Bremen, Bremen, Germany. 1 19 41 42 Faculty of Geology and Geoenvironment, National and Kapodistrian University of Athens, Athens, Greece. 45 21 Institute of Geosciences, Christian-Albrechts-Universität zu Kiel, Kiel, Germany. 46 22 University Ss Cyril and Methodius, Institute of Biology, Skopje, FYROM. 49 50 24 53 26 43 44 20 47 48 23 51 52 25 Paleoenvironmental Dynamics Group, Institute of Earth Sciences, Heidelberg University, Heidelberg, Germany. Dipartimento di Scienze della Vita, Laboratorio di Palinologia e Paleobotanica, Università di Modena e Reggio Emilia, Modena, Italy. Laboratoire des Sciences du Climat et de l’Environnement, UMR 8212, CEA/CNRS/UVSQ et Université Paris-Saclay, Gif-Sur-Yvette, France. Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany. The Department of Geography, Geology and Environment, University of Hull, Hull, UK. 54 27 Istituto Nazionale di Geofisica e Vulcanologia, Rome, Italy. 55 28 Dipartimento di Scienze della Terra e Geoambientali, University of Bari, Bari, Italy. 29 IDPA-CNR, Milan, Italy. 57 30 Center for Climate Physics, Institute for Basic Science, Busan, South Korea 58 31 Pusan National University, Busan, South Korea 61 33 Leibniz Institute for Applied Geophysics (LIAG), Hannover, Germany. 62 34 Institute of Loess Plateau, Shanxi University, Taiyuan, China. 64 *Correspondence to: wagnerb@uni-koeln.de. 65 † these authors contributed equally to this work. 56 59 60 63 32 Dipartimento di Scienze della Vita, Laboratorio di Palinologia e Paleobotanica, Università di Modena e Reggio Emilia, Modena, Italy. 66 67 2 68 Precipitation is a key factor for socioeconomic development in densely populated and summer 70 difficult to project accurately. While current climate model simulations indicate a progressive 72 well constrained1. Only a few continental proxy records capable of capturing hydroclimate 74 underlying orbital geometries, necessary to validate climate model data on Quaternary time 76 transient climate model hind cast, to show that high winter precipitation anomalies occur 78 activity. While this is counter-intuitive at first sight, our data suggest that increased sea- 80 pressure systems entering the Mediterranean during phases characterized by low continental 69 dry regions such as the Mediterranean realm. Seasonal and regional changes are critical, but 71 summer drying over the next century, precipitation changes during winter months are less 73 change cover multiple Northern Hemisphere summer insolation maxima2,3 with different 75 scales. Here we use a 1.36 million year proxy time series from Lake Ohrid, coupled to a long 77 during phases with strong seasonal contrast in insolation and high African summer monsoon 79 surface temperatures amplify local cyclogenesis while also refuelling North Atlantic low 81 ice volume and high atmospheric CO2 concentrations. Comparison with modern reanalysis 83 to those driving the reconstructed precipitation increases. Our extended record covers multiple 85 performance. 82 data shows that current drivers of rainfall amount in the Mediterranean share some similarities 84 insolation maxima and therefore is an important benchmark for testing climate-model 86 87 88 89 90 91 92 Mediterranean climates are characterized by strong seasonal contrasts between dry and warm summers, and wet and mild winters. The amount and temporal extent of precipitation during the winter half-year (October through March) determines the prevailing type of vegetation and water availability for agrarian land-use in the Mediterranean borderlands. In recent decades, reduction of winter precipitation has become a regular phenomenon in this region, 3 93 with anthropogenic greenhouse gas (GHG) and aerosol forcing identified as potential 95 Pathway (RCP) 4.5 and 8.5 scenarios, predict a progressive summer drying over the next 97 well constrained, with different simulation runs showing trends both towards wetter and drier 99 current modelling approaches are useful for decision makers5,6. 94 contributors4. Current climate model simulations, using the Representative Concentration 96 century1. Precipitation changes during the Northern Hemisphere (NH) winter months are less 98 conditions. The uncertainty in winter precipitation projections limits the extent to which 100 Long-term, empirical baseline data are helpful to constrain uncertainties in climate 101 modelling proxy records. Proxy records and modelling experiments suggest that enhanced 103 intertropical convergence zone (ITCZ) and increase in African monsoon strength during 105 winter insolation (NHWI) minima2,7,8,9. However, most continental records that are capable of 107 orbital geometries. In fact, the majority of records are limited to the Holocene10,11, yet the 109 interglacials, due to lower eccentricity. Terrestrial proxy time series covering multiple NHSI 111 Mediterranean Sea provide continuity throughout the Plio-Pleistocene and capture cessations 113 layers12,13. While multiple factors contribute to sapropel formation, increased freshwater 115 considered the most important14,15. Hence, the Mediterranean sapropel record is thought to be 117 a direct indicator of precipitation in, and runoff from, the entirety of the Mediterranean realm. 102 precipitation in the Mediterranean region is in phase with the northward shift of the 104 precession minima causingNorthern Hemisphere summer insolation (NHSI) maxima and 106 capturing hydroclimate change do not cover multiple NHSI maxima with different underlying 108 Early Holocene NHSI maximum was relatively weak compared to most other Quaternary 110 maxima from the Mediterranean region are scarce2,3. Sediment records from the 112 of deep-water ventilation associated with the formation of prominent, organic-rich sapropel 114 input, particularly from the African continent during NHSI-forced monsoon maxima, is 116 an excellent indicator of the relative timing of increased African monsoon strength rather than 118 Reconstructed precipitation increases in the northern Mediterranean borderlands during 4 119 sapropel formation have been interpreted to be a product both of intensified summer and 121 stronger wintertime storm tracks2 or air-sea temperature difference, and locally induced 123 Alternatively, conceptual models based on proxy time series have suggested increases in the 120 122 124 winter precipitation15,16. Modelling experiments explain increased winter precipitation by convective precipitation that dominate freshwater budget changes on obliquity time scales17. frequency and intensity of low-pressure systems evolving in the Mediterranean region, mostly 125 during fall and early winter7,8,16. Hence, a well-dated proxy record covering multiple glacial- 127 addressing long-standing questions regarding the underlying mechanisms, timing, and 129 concentration, orbital geometries, continental ice sheet volume and extent). 131 Myr sedimentary record from Lake Ohrid (Fig. 1, Extended Data Fig. 1). Climate variations at 133 We compare our sedimentary proxy time series with transient climate simulation data and 135 variability and seasonality, as well as phase relationships to orbital forcing. 137 primarily fed by an extensive karst aquifer system, which supplies ions (mainly Ca2+ and 139 m-long composite sediment succession from the lake centre, comprised of fine-grained hemi- 141 with no evidence of unconformities or erosion surfaces. Independent age control from 16 143 Fig. 2, Extended Data Table 1, Extended Data Table 2) provides a robust chronological 126 interglacial cycles and being sensitive to changes in Mediterranean hydroclimate is key to 128 amplitude of precipitation variability under different climate boundary conditions (GHG 130 Here, we assess precipitation variability in a continuous, independently dated 1.36 132 this site represent broader climate variability across the northern Mediterranean borderlands18. 134 prominent monsoon records, to provide a mechanistic understanding of precipitation 136 Lake Ohrid is of tectonic origin and 293 m deep. The lake is hydrologically open and 138 HCO3-) to the lake and filters particulate matter19. Scientific drilling in 2013 resulted in a 584- 140 pelagic muds in the upper 447 m18,20. Sedimentation is thought to have been uninterrupted, 142 interspersed tephra layers in combination with magnetostratigraphy (Fig. 1, Extended Data 5 144 framework. This framework allows us to match changes in orbital parameters with our proxy 146 spans the last 1.36 Myr (Fig. 1). 148 excluding pine (AP-P), deciduous oaks), and hydrological variability (total inorganic carbon 150 precessional (~21 ka) component (Fig. 2; Extended Data Figs 3, 4, and 5). During periods of 152 hydrological and vegetation proxy data (Fig. 2). We interpret these peaks in TIC (mainly from 154 activity of, and ion supply from, the karst aquifers combined with higher aquatic productivity 145 147 149 data to refine the age-depth relationships. The data demonstrate that the Lake Ohrid record Indicators for detrital input (quartz, potassium), catchment vegetation (arboreal pollen (TIC), Ca/K, δ18Ocalcite, δ13Ccalcite) show clear orbital-scale cyclicity, also characterized by a 151 global ice volume minima and NHSI maxima, we observe prominent peaks in the 153 endogenic calcite) and Ca/K (a proxy for the concentration of calcite) to result from enhanced 155 due to warmer conditions19. Pollen show a simultaneous increase in vegetation cover, 156 particularly deciduous oaks, during early phases of interglacials. Deciduous oaks benefit from 158 suggest greater soil development, while lower δ18Ocalcite (Extended Data Fig. 3) indicate more 160 suggest higher temperatures along with maxima in annual precipitation amount and potential 157 a limited length of the summer dry season21. Lower δ13Ccalcite values during these periods 159 positive precipitation/evaporation (P/E) balance18. Thus, aquatic and terrestrial datasets 161 shorter summer aridity during interglacials (Extended Data Fig. 4). 163 Lake Ohrid record in a regional context, we analysed climate data time series derived from a 162 To provide a better understanding of the observed precipitation variability from the 164 transient 784 kyr simulation using the earth system model LOVECLIM22,23 (Extended Data 166 period 1979–2017. Temperature time series of the 5°x5° Lake Ohrid grid cell simulated by 168 (Extended Data Fig. 3), such as the LR04 benthic oxygen isotope stack24 (r=–0.8737 or 165 Fig. 6) as well as NOAA reanalysis precipitation data of the Lake Ohrid region for the time 167 the LOVECLIM earth system model closely resemble records of first-order global ice volume 6 169 r2=0.76 based on 1000-year averages of both data sets). The close match to changes in the 170 amount of detrital siliciclastics and tree pollen (AP-P) confirms the sensitivity of the Lake 172 highest amplitudes in precipitation time series occur during phases of reduced ice volume, 174 simulated precipitation and our precipitation proxy time series (r2=0.38) and the persistence of 176 response recorded at Lake Ohrid also captures changes in regional hydroclimate back to 1.36 171 Ohrid record to global-scale climate fluctuations (Fig. 2; Extended Data Figs 3 and 4). The 173 with prominent peaks during NHSI maxima. The significant positive relationship between 175 the relationship with the orbital parameters (Extended Data Fig. 4) suggest that the local 177 Ma (Fig. 2). 179 monsoon systems during precession minima and NHSI maxima is a prominent example of 178 Seen both in paleo records and in climate model simulations, the intensification of NH 180 orbitally-forced changes in precipitation variability14,15,25. Iconic records of monsoon strength, 182 foraminifera oxygen isotope records14,15,27, show a positive phase relationship with Lake 184 northward displacement of atmospheric circulation systems, including the position of the 186 subsidence over, and persistence of, high-pressure systems in the Mediterranean region, 188 Reduced NHWI has highest impact on tropical and subtropical latitudes2 and leads to low 190 Furthermore, this cooling results in a reduced meridional temperature gradient leading to a 192 relationship between the Lake Ohrid precipitation record (Fig. 2, Extended Data Figs 3 and 4) 194 region when NHWI is low. 181 such as the Chinese speleothem26, eastern Mediterranean sapropel12,13,26 and planktonic 183 Ohrid hydrological proxy time series (Fig. 2). Strengthening of NH monsoons results from a 185 Hadley cells and the ITCZ during NH summer. The shift of the Hadley cell amplifies 187 leading to warmer and drier summers17, and higher sea-surface temperatures (SST)16,28. 189 latitude cooling and a southward shift of the ITCZ and the NH Hadley and Ferrel cells. 191 weakening of the westerlies based on the thermal wind relationship. The observed 193 and the monsoon archives suggests increased precipitation during the winter half-year for this 7 195 The Lake Ohrid record, in combination with the transient simulation time series and 196 the NOAA reanalysis data, may provide fundamental insights into the mechanisms invoked 198 last 39 years show high precipitation anomalies (defined as above two standard deviations) to 200 atmospheric pattern associated with these precipitation events exhibits a trough in the Gulf of 197 by orbital forcing on Mediterranean precipitation. The monthly NOAA reanalysis data of the 199 occur between the months of September and December (Extended Data Fig. 7a,b). The 201 Genoa region (Extended Data Fig. 7c), pointing to either increased cyclogenesis over or 202 advection of North Atlantic low pressure into the western Mediterranean region. 204 agreement with the reanalysis data; the model, however, underestimates the annual mean 203 205 206 The annual cycle of simulated Lake Ohrid precipitation in LOVECLIM is in good precipitation (Extended Data Fig. 8b). Maxima in our simulated precipitation time series (defined as above two standard deviations) indicate a positive anomaly from September to 207 November (SON) in agreement with the reanalysis data (Fig. 3, Extended Data Fig. 8b). 209 both the NOAA and LOVECLIM data show pronounced troughs in the central Mediterranean 211 observations support previous modelling experiments suggesting that weakened atmospheric 213 increased contrast between warm SST and lower continental air temperatures17, fuel 215 at the beginning of the fall, when the stronger thermal inertia of the sea relative to the land 217 in the NH atmospheric circulation cells during the winter half-year, which also favours a more 219 lead to increased winter rainfall in the Mediterranean mid-latitudes. 208 Despite important differences in the geographical expansion of geopotential height anomalies, 210 area and an increase of rainfall during winter half-year in our focus region (Fig. 3). Our 212 stratification and reduced hemispheric temperature contrasts2, in combination with an 214 precipitation increase in the Mediterranean. Such a preconditioning is particularly pronounced 216 promotes local cyclogenesis17,29. Local cyclogenesis in combination with the southward shift 218 southerly trajectory for storm tracks across the North Atlantic and into the Mediterranean2, 8 220 Owing to the significant positive correlation between the simulation and our proxy 221 time series (Extended Data Fig. 4), in terms of timing and amplitude, we infer that this 223 similar to the NH summer monsoon records, we observe a strong influence of NHSI and a 225 series, suggesting persistence of the mechanism during different climate boundary conditions. 227 sapropel records (Fig. 2) indicates a strong coherence of African summer monsoon strength 222 mechanism primarily controlled precipitation at Lake Ohrid for the last 1.36 Myr. Indeed, 224 reduced winter temperature contrast in the NH throughout the entirety of our multiproxy time 226 The positive phase relationship between the Lake Ohrid precipitation proxy time series and 228 229 230 231 232 and widespread Mediterranean winter half-year precipitation. Some peaks in our precipitation proxy time series, which are not represented by sapropel layers (Fig. 2), may indicate lower monsoon strength and reduced runoff from the African continent or that the general setting required for sapropel deposition and preservation was not established in the Mediterranean Sea during these periods15. During colder and drier glacial periods3 with increased global ice 233 volume, lower atmospheric CO2 concentrations, and stronger mid-latitude westerlies, 235 the sensitivity simulations conducted to disentangle the individual effects of orbital forcing, 234 insolation forcing on precipitation appears suppressed in our record. This is in agreement with 236 NH ice sheets, and CO2 on Lake Ohrid precipitation (Extended Data Fig. 6). 238 also exerts a strong control on precipitation variability in the Mediterranean mid-latitudes 240 well-constrained information on precipitation maxima during phases of lower 242 equivalence of the past regional key drivers of precipitation extremes to those produced by 244 simulation uncertainties and makes these results also relevant to predictions for the future 237 Precessional forcing on insolation is not only the key driver of the NH monsoons, it 239 during the Quaternary. Lake Ohrid sediment cores record highly resolved and chronologically 241 intrahemispheric temperature contrast and peak SST’s over the last 1.36 Myr. The apparent 243 continued anthropogenic increase of atmospheric GHG concentrations may help to reduce 245 evolution of Mediterranean climate. 9 246 247 248 References: 1. IPCC, 2013: Annex I: Atlas of Global and Regional Climate Projections [van 249 Oldenborgh, G. J., M. Collins, J. Arblaster, J. H. Christensen, J. Marotzke, S. B. Power, 251 Basis. Contribution of Working Group I to the Fifth Assessment Report of the 253 Tignor, S. K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex and P. M. Midgley (eds.)]. 255 (2013). 250 M. Rummukainen and T. Zhou (eds.)]. In: Climate Change 2013: The Physical Sci¬ence 252 Intergovernmental Panel on Climate Change [Stocker, T. F., D. Qin, G.-K. Plattner, M. 254 Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA 256 2. Kutzbach, J. E., Chen, G., Cheng, H., Edwards, R. & Liu, Z. Potential role of winter 257 rainfall in explaining increased moisture in the Mediterranean and Middle East during 259 (2014). doi:10.1007/s00382-013-1692-1 periods of maximum orbitally-forced insolation seasonality. Clim. Dyn. 42, 1079−1095 258 260 3. Philippon, revised chronostratigraphy and long-term vegetation trends. Quat. Sci. Rev. 261 262 25, 3416−3430 (2006). doi:10.1016/j.quascirev.2006.09.002 263 4. 265 5. 267 6. 264 266 268 269 Tzedakis, P. C., Hooghiemstra, H. & Pälike H. The last 1.35 million years at Tenaghi Hoerling, M. et al. On the increased frequency of Mediterranean drought. J. Climate 25, 2146−2161 (2012). doi:10.1175/JCLI-D-11-00296.1 Weisheimer, A. & Palmer, T. N. On the reliability of seasonal climate forecasts. J. Royal Soc., Interface 11, 20131162 (2014). doi.org/10.1098/rsif.2013.1162 Totz, S., Tziperman, E., Coumou, D., Pfeiffer, K. & Cohen, J. Winter Precipitation Forecast in the European and Mediterranean Regions Using Cluster Analysis. Geophys. Res. Lett. 44, 12,418–12,426 (2017). doi.org/10.1002/2017GL075674 10 270 7. 272 8. 271 early part of the Last Interglacial. Geology 40, 919−922 (2012). doi:10.1130/G33204.1 275 276 277 278 279 280 281 282 283 284 285 286 Toucanne, S. et al. Tracking rainfall in the northern Mediterranean borderlands during sapropel deposition. Quat. Sci. Rev. 129, 178−195 (2015). 273 274 Milner, A. M. et al. Enhanced seasonality of precipitation in the Mediterranean during the doi:10.1016/j.quascirev.2015.10.016 9. Stockhecke, M. et al. Millennial to orbital-scale variations of drought intensity in the Eastern Mediterranean. Quat. Sci. Rev. 133, 77−95 (2016). doi: 10.1016/ j.quascirev.2015.12.016 10. Roberts, N. et al. Stable isotope records of Late Quaternary climate and hydrology from Mediterranean lakes: the ISOMED synthesis. Quat. Sci. Rev. 27, 2426−2441 (2008). doi:10.1016/j.quascirev.2008.09.005 11. Magny, M. et al. North–south palaeohydrological contrasts in the central Mediterranean during the Holocene: tentative synthesis and working hypotheses. Clim. Past 9, 2043−2071 (2013). doi:10.5194/cp-9-2043-2013 12. Emeis K.-C., Camerlenghi A., McKenzie J. A., Rio D.& Sprovieri R., The occurrence and significance of Pleistocene and Upper Pliocene sapropels in the Tyrrhenian Sea. Mar. Geol. 100, 155−182 (1991). doi:10.1016/0025-3227(91)90231-R 287 13. Kroon, D. al. Oxygen isotope and sapropel stratigraphy in the Eastern Mediterranean 289 Scientific results, A. H. F. Robertson, K.-C. Emeis, C. Richter, A. Camerlenghi. Eds. 288 290 291 292 293 during the last 3.2 million years, in Proceedings of the Ocean Drilling Program. (College Station, Texas, 1998), vol. 160, pp 181−190 (1998). 14. Rossignol-Strick, M. Mediterranean Quaternary sapropels, an immediate response of the African monsoon to variation of insolation. Palaeogeogr. Palaeoclimatol. Palaeoecol. 49, 237−263 (1985). doi:10.1016/0031-0182(85)90056-2 11 294 15. Rohling, E. J., Marino, G. & Grant, K. M. Mediterranean climate and oceanography, and 295 the periodic development of anoxic events (sapropels). Earth Sci. Rev. 143, 62−97 297 16. Tzedakis, P. C. Seven ambiguities in the Mediterranean palaeoenvironmental narrative. 299 17. Bosmans, J. H. C. et al. Precession and obliquity forcing of the freshwater budget over 296 298 300 301 302 303 304 305 306 307 308 309 310 311 312 (2015). doi:10.1016/j.earscirev.2015.01.008 Quat. Sci. Rev. 26, 2042−2066 (2007). doi:10.1016/j.quascirev.2007.03.014 the Mediterranean. Quat. Sci. Rev., 123, 16−30 (2015). doi:10.1016/j.quascirev.2015.06.008 18. Wagner, B. et al. The environmental and evolutionary history of Lake Ohrid (FYROM/Albania): Interim results from the SCOPSCO deep drilling project. Biogeosciences 14, 2033−2054 (2017). doi:10.5194/bg-14-2033-2017 19. Vogel, H., Wagner, B., Zanchetta, G., Sulpizio, R. & Rosén, P. A paleoclimate record with tephrochronological age control for the last glacial-interglacial cycle from Lake Ohrid, Albania and Macedonia. J. Paleolimnol. 44, 295−310 (2010). doi:10.1007/s10933-009-9404-x 20. Francke, A. et al. Sedimentological processes and environmental variability at Lake Ohrid (Macedonia, Albania) between 637 ka and the present. Biogeosciences 13, 1179−1196 (2016). doi:10.5194/bg-13-1179-2016 21. Forner, A. et al. Extreme droughts affecting Mediterranean tree species’ growth and 313 water-use efficiency: the importance of timing. Tree Physiol. 38, 1127−1137 (2018). 315 22. Friedrich, T., Timmermann, A., Tigchelaar, M., Timm, O. E. & Ganopolski, A. Nonlinear 314 316 317 doi:10.1093/treephys/tpy022 climate sensitivity and its implications for future greenhouse warming. Sci. Adv. 2 (2016), p. e1501923. doi:10.1126/sciadv.1501923 12 318 23. Timmermann, A. & Friedrich, T. Late Pleistocene climate drivers of early human 320 24. Lisiecki, L. E. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed 319 321 322 migration. Nature 538, 92−95 (2016). doi:10.1038/nature19365 benthic δ18O records. Paleoceanography 20, PA1003 (2005). doi:10.1029/2004PA001071 323 25. Cheng, H. et al. The Asian monsoon over the past 640,000 years and ice age 325 26. Konijnendijk, T. Y. M., Ziegler, M. & Lourens, L. J. Chronological constraints on 324 326 327 328 329 330 331 332 333 334 335 336 337 338 terminations. Nature 534, 640−646 (2016). doi:10.1038/nature18591 Pleistocene sapropel depositions from high-resolution geochemical records of ODP Sites 967 and 968. Newslett. Stratigr. 47, 263−282 (2014). doi:10.1127/0078-0421/2014/0047 27. Colleoni, F., Masina, S., Negri, A. & Marzocchi, A. Plio–Pleistocene high–low latitude climate interplay: a Mediterranean point of view. Earth Planet. Sci. Lett. 319–320, 35−44 (2012). doi:10.1016/j.epsl.2011.12.020 28. Martrat, B., Jimenez-Amat, P., Zahn, R. & Grimalt J. O., Similarities and dissimilarities between the last two deglaciations and interglaciations in the North Atlantic region. Quat. Sci. Rev. 99, 122−134 (2014). doi:10.1016/j.quascirev.2014.06.016 29. Trigo, R. M., Osborne, T. J. & Corte-Real, J. M. The North Atlantic Oscillation influence on Europe: climate impacts and associated physical mechanisms. Clim. Res. 20, 9−17 (2002). doi:10.3354/cr020009 30. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the earth. Astron. Astrophys. 428, 261−285 (2004). doi:10.1051/0004-6361:20041335 339 340 13 341 Acknowledgments: The Hydrobiological Institute in Ohrid (S. Trajanovski and G. Kostoski) 343 logistic support for site surveys and the scientific drilling campaign. Drilling was carried out 345 Skinner provided logistic and technical advice prior and during drilling operation. The 347 project was funded by the International Continental Scientific Drilling Program (ICDP), the 349 University of Cologne, the British Geological Survey, the INGV and CNR (both Italy), and 351 tephra, which was 40Ar/39Ar dated with funding from LEFE "INTERMED" grant (CNRS- 342 and the Hydrometeorological Institute in Tirana (M. Sanxhaku and B. Lushaj) provided 344 by Drilling, Observation and Sampling of the Earth’s Continental Crust (DOSECC). A. 346 Scientific Collaboration on Past Speciation Conditions in Lake Ohrid (SCOPSCO) drilling 348 German Ministry of Higher Education and Research, the German Research Foundation, the 350 the governments of the republics of North Macedonia and Albania. V. Scao collected the V5 352 INSU) to S. Nomade. 354 Author Contributions: B Wagner and H Vogel designed the study and contributed equally. 356 conceived major scientific ideas of this study. A Francke (sedimentology, chronology), T 358 and L Sadori (palynology) contributed and oversaw key datasets used in the study. They 353 355 BW initiated and coordinated the SCOPSCO drilling project and drilling campaign. HV 357 Friedrich (LOVECLIM modelling), T Donders (palynology), J Lacey (isotope geochemistry), 359 360 361 362 coordinated together with F Cremer-Wagner, M Leng, E Regattieri, T Wilke and G Zanchetta discussion and interpretations of proxy data groups and model results. Specific data were provided by A Bertini (pollen, MIS 19–21, MIS 25–28, MIS 42–43), N Combourieu-Nebout (pollen, MIS 1–4, MIS 8, MIS 14–15), B Giaccio (tephrostratigraphy), S Joannin (pollen, 363 MIS 1–4, MIS 13–16, MIS 30), J Just (paleomagnetic data), K Kouli (pollen, MIS 6–8, MIS 365 Koutsodendris (pollen, MIS 11–12, MIS 15), N Leicher (tephrostratigraphy), A Masi (pollen, 364 10, MIS 16–19, MIS 28–30, MIS 33), I Kousis (pollen, MIS 11–12, MIS 15), A 14 366 MIS 5–6, MIS 20–25, MIS 31–32), A M Mercuri (pollen, MIS 6, MIS 34), S Nomade 368 MIS 35–43), O Peyron (pollen, MIS 1–4, MIS 13–16, MIS 30), L Sagnotti (paleomagnetic 370 6, MIS 34). S Krastel, K Lindhorst, and T Wonik coordinated the seismic survey of Lake 372 correlation. A Grazhdani, M Melles, J Reed, and Z Levkov contributed to the conception of 374 micropaleontological and organic geochemistry data, which confirmed that the sediment 376 provided model infrastructure and resources. All authors contributed to the discussion and 377 interpretation of the data and provided comments and suggestions to the manuscript. 378 379 Author Information: Reprints and permissions information is available at 367 (tephrochronology), N Nowaczyk (paleomagnetic data), K Panagiotopoulos (pollen, MIS 7–8, 369 data), G Sinopoli (pollen, MIS 5–6), R Sulpizio (tephrostratigraphy) and P Torri (pollen MIS 371 Ohrid, the selection of the coring location and the geophysical measurements needed for core 373 the work. A Cvetkoska, J Holtvoeth, E Jovanvoska, S Tofilovska, and X Zhang provided 375 succession from the DEEP site covers the entire history of Lake Ohrid. A Timmermann 380 www.nature.com/reprints. Authors declare no competing interests. Correspondence and 381 requests for materials should be addressed to wagnerb@uni-koeln.de. Data are available in 383 https://doi.pangaea.de/10.1594/PANGAEA.896848. Data used for LOVECLIM are available 382 the main text, in the supplementary materials and in the Pangaea database at 384 at https://climatedata.ibs.re.kr/grav/data/loveclim-784k. 386 Figure legends: 388 is based on tephrostratigraphic correlation of 16 tephra layers to their radiometrically dated 385 387 Fig. 1. Chronology and location of the Lake Ohrid DEEP site record. (a) The age model 389 proximal deposits (red, first-order tie points) (b) tuning of total organic carbon (TOC) minima 390 in the DEEP site record vs. inflection points in insolation and winter season length (blue, 15 391 second-order tie points), and cross evaluation of two paleomagnetic age reversals (a; dashed 393 247 meters composite depth (mcd) of the record20 (see Methods). For the ages and errors of 395 ±2,000 years. (c) The insert shows the location of Lake Ohrid and the approximate position of 397 Fig. 2. Lake Ohrid precipitation indicators and global monsoon records for the last 1.4 399 interval/oxidized sapropel, violet = ghost sapropel)12,13,27; (b) Chinese Speleostack δ18O 26 in 392 lines). The age model was calculated following the methodological approach for the upper 394 the tephra layers, see Extended Data Table 1. The tuning points (green) include an error of 396 the intertropical convergence zone (ITCZ) in summer and winter. 398 million years. (a) Eastern Mediterranean (EM) Sapropel ages (green = sapropel, red = red 400 ‰ relative to VPDB; (c) Medstack δ18O planktonic28 in ‰ relative to VPDB; SST=sea- 402 relative to VPDB; (e) Lake Ohrid deciduous oaks pollen percentage; (f) Lake Ohrid total 404 between the tropic of cancer and the arctic circle30; (h) annual mean precipitation amount for 406 excluding Pinus pollen (AP-P) percentages. Tenaghi Philippon arboreal pollen (AP) 401 surface temperature, SSS=sea-surface salinity; (d) Lake Ohrid δ13C endogenic calcite in ‰ 403 inorganic carbon (TIC) concentrations; (g) Northern Hemisphere winter insolation difference 405 the Lake Ohrid grid cell from the LOVECLIM simulation; (i) Lake Ohrid arboreal pollen 407 percentages3 (k) and LR04 benthic δ18O stack25 in ‰ relative to VPDB with odd numbers for 409 radiometrically dated tephra layers, blue and white diamonds the position of reversals of 411 Fig. 3. Simulated Lake Ohrid precipitation and atmospheric anomaly pattern associated 413 cell. Data based on 1,000-year averages. Dashed line indicates two standard deviations above 415 Methods for details on the model simulations. (b) Composite anomalies of September- 408 interglacials (l) are shown for comparison. Red and white diamonds indicate the position of 410 Earth’s magnetic field in the Lake Ohrid sediment record. 412 with precipitation maxima. (a) Simulated precipitation (cm yr-1) for the Lake Ohrid grid 414 the mean. Red shading highlights precipitation values exceeding two standard deviations. See 16 416 417 November (SON), 800 hPa geopotential height (m, shading) and wind (m s-1, vectors) associated with precipitation maxima shown in (a). 418 419 420 421 422 17 423 Methods: 424 Lake and lake hydrology 426 climate zone with average monthly air temperature ranging from +26°C during summer to - 428 yr-1 with increasing altitude and occurs primarily during winter months31. The lake is ~30 km 425 Lake Ohrid (41°02’N, 20°43’E, 693 m a.s.l.; Fig. 1c) is located in the sub-Mediterranean 427 1°C during winter. Precipitation in the Lake Ohrid watershed increases from 698 to 1,194 mm 429 long, ~15 km wide, and has a maximum water depth of 293 m (Extended Data Fig. 1). 431 water input. Due to an oligotrophic state, bottom waters remain partly oxygenated for several 430 Sublacustrine karst springs (55%), direct precipitation, and river inflow (45%) constitute the 432 years, although the lake is oligomictic and a complete overturn occurs only every few years at 433 present32. 434 435 Sediment cores 437 Deep Lake Drilling System (DLDS) of Drilling, Observation and Sampling of the Earth’s 439 interdisciplinary Scientific Collaboration on Past Speciation Conditions in Lake Ohrid 441 Program (ICDP). The composite sediment record is based on 6 parallel boreholes that reached 443 99.8%. Small gaps occur between 204.719 and 204.804 mcd (8.5 cm) and between 447.89 mcd, and between 55 and 50 mcd. Subsampling in the upper 447.12 mcd excluded mass 436 Sediment cores from the Lake Ohrid DEEP site were recovered in spring 2013, using the 438 Continental Crust (DOSECC) and within the framework of the multinational and 440 (SCOPSCO) project that was co-sponsored by the International Continental Scientific Drilling 442 a terminal depth of 568 m33. Sediment recovery from 0 to 456.1 m composite depth (mcd) is 444 and 448.19 mcd (30 cm)33. Mass movement deposits (<3 cm) occur between 117 and 107 445 446 movement and tephra deposits. 447 448 Scanning-X-ray fluorescence (XRF) analysis 18 449 Scanning-XRF analysis was performed at the University of Cologne, Germany, on split core 451 Analytics) equipped with an energy dispersive silicon drift detector and a Cr-tube set to 30 453 integrated in Q-spec (Cox Analytics). 455 Elemental analysis 457 following freeze-drying and homogenization at the University of Cologne. For total carbon 459 sample material was dispersed in 10 ml deionized water. TC was determined at combustion of 461 100 and a DIMATOC 200 (DIMATEC Corp., Germany). The total organic carbon (TOC) 450 surfaces at 2.5 mm increments and 10 s dwell time using an ITRAX XRF core scanner (Cox 452 kV/30 mA. Raw data were processed and element-specific photon energy peaks were 454 456 Elemental analysis was performed on 16-cm-spaced samples (2794 samples, ~480 yr) 458 (TC) and total inorganic carbon (TIC) measurements, an aliquot of 40 mg of the homogenized 460 900°C and TIC was measured after treatment with 40% H3PO4 at 160°C using a DIMATOC 462 content was calculated as the difference between TC and TIC. 464 Fourier Transform Infrared Spectroscopy (FTIRS) 466 32 cm (1462 samples, ~1,000 yr). Measurements were performed using a Bruker Vertex 70 468 and a HTS-XT accessory unit (multisampler) in an air-conditioned laboratory at the 463 465 Relative concentration changes for quartz were assessed using FTIRS, on samples spaced at 467 equipped with a lN2-cooled MCT (mercury-cadmium-telluride) detector, a KBr beam splitter, 469 University of Bern, Switzerland. For this purpose, 11 mg of each sample and 500 mg of oven- 470 471 dried spectroscopic grade KBr (Uvasol®, Merck Corp.) were homogenized and scanned 64 times at a resolution of 4 cm-1 (reciprocal centimetres) for the wavenumber range from 3,750 472 to 520 cm-1 in diffuse reflectance mode. Data processing encompassed a linear baseline 474 zero (3,750 and 2,210–2,200 cm-1). Peak areas diagnostic for symmetric stretching of SiO4 in 473 correction to remove baseline shifts and tilts by setting two points of the recorded spectrum to 19 quartz (778 and 798 cm−1), and representative for relative abundance34,35 were integrated 475 476 using the OPUS (Bruker Corp.) software package. 477 478 Palynology processing and analysis 480 following processing, identification, and counting approaches as described in36. Dry sediment 479 Pollen analysis was carried out on sediment samples spaced at 64 cm (697 samples, ~2000 yr) 481 (1.0–1.5 g) samples were treated with cold HCl (37%vol), cold HF (40%vol), and hot NaOH 482 483 (10%vol) to dissolve carbonates, silicates, and humic acids, respectively. Glycerin-mounted residues were analysed by transmitted light microscopy to a mean of ~533 (incl. Pinus) and 484 ~250 (excl. Pinus) grains/sample. Relative abundances are based on the total terrestrial pollen 486 taxon36. Deciduous oak abundances represent the combined percentages of Quercus robur and 485 sum excl. Pinus due to overrepresentation and potential long-distance transport of this 487 Q. cerris types37, which is commonly used as an indicator for mid-elevation, relatively humid 488 forest across the Mediterranean38,39,40,41. 489 490 Isotope analysis 492 16 cm through zones of higher TIC (>0.5%), comprising a total of 1309 sediment samples. 494 sediment and oxidize reactive organic material. Potential biogenic carbonate was removed by 496 to a fine powder in an agate mortar. CO2 was evolved from 10 mg CaCO3 powders by 491 Oxygen and carbon isotopes were analysed on bulk carbonate (calcite)42 in samples spaced at 493 The samples were immersed in 5% NaClO solution for 24 h to gently disaggregate the 495 sieving and the <64 μm fraction washed with deionized water, dried at 40°C, and then ground 497 reaction with anhydrous H3PO4 overnight inside a vacuum at a constant temperature of 25°C. 498 The liberated CO2 was cryogenically purified under vacuum and collected for analysis on a 500 standard delta notation (δ18Ocalcite and δ13Ccalcite, respectively) in per mille (‰) calculated to 499 VG Optima dual inlet mass spectrometer. Oxygen and carbon isotope values are reported in 20 501 the Vienna Pee Dee Belemnite (VPDB) scale using a within-run laboratory standard (MCS) 502 503 standard was <0.1‰ (±1σ) for δ18O and δ13C. 504 505 Magnetostratigraphic analyses 507 demagnetization (10 steps up to 100 mT) was measured on ~900 discrete cube (6.3 cm3) 509 GeoForschungsZentrum, Potsdam, Germany, using a 2G Enterprises cryogenic 511 principle component analysis (PCA) after removal of low-coercivity magnetic overprints. calibrated against international NBS standards. Analytical reproducibility for the within-run 506 Remanent magnetization in its natural state (NRM) and after step-wise alternating field 508 samples with an average 48-cm-spacing at the Paleomagnetic Laboratory at the 510 magnetometer. Paleomagnetic directions (declination and inclination) were calculated using 512 After identification of geomagnetic polarity transitions, ~500 additional samples were taken at 513 2 to 3-cm-spacing across these transitions for high-resolution analysis at the Istituto Nazionale 514 di Geofisica e Vulcanologia, Rome, Italy, using the same analytical set up and routine as in 516 overprints the primary paleomagnetic signal43, paleomagnetic transitions are faithfully 515 517 Potsdam. As glacial intervals of the core contain diagenetically formed greigite, which preserved only in interglacial intervals, at the base of the Jaramillo sub-Chron (373.8 mcd) 518 and at the Matuyama/Brunhes (M/B) boundary (287.6 mcd). 520 Tephrostratigraphic analysis 522 mcd of the record44,45,46. Two additional tephra layers from the lower (>247 mcd) part of the 524 based on geochemical fingerprinting of single glass shards using Wavelength Dispersive 519 521 Eleven tephra and three cryptotephra layers have successfully been identified in the upper 247 523 DEEP site record are introduced here. The tephrostratigraphic correlation of these tephras is 525 Electron Microprobe Analysis (WDS-EPMA) as described in 44. 21 526 40 Ar/39Ar dating was performed at the LSCE facility (CEA, UVSQ and University 527 Paris-Saclay). V5 tephra (=OH-DP-2669 layer) was collected in Montalbano-Jonico 529 fraction 0.6-1.0 mm, were extracted from V5 and irradiated for 2 h in the Cd-lined, in-core 528 530 531 532 533 (Southern Italy, N40°17’32.8’’; E16°33’27.4’’). Twenty pristine sanidine crystals, of the CLICIT facility of the Oregon State University TRIGA reactor (Irradiation CO 001). Subsequently, 14 crystals were individually loaded in a copper sample holder and put into a double vacuum Cleartran window. Each crystal was individually fused using a Synrad CO2 laser at 10-15% of nominal power (~50 W). The extracted gas was purified for 10 min by two 534 hot GP 110 and two GP 10 getters (ZrAl). Ar isotopes (36Ar, 37Ar, 38Ar, 39Ar and 40Ar) were 536 217 SEV SEN coupled to an ion counter. Neutron fluence J for each sample is calculated 538 the total decay constant of 48. J-values computed from standard grains is 0.00053220 ± 540 analytical period, and was relative to a 40Ar/36Ar ratio of 298.56 49. 542 boundaries comprising up to 500 µm large platy glass shards and minor elongated 544 position between the M/B boundary (287.6 mcd) and OH-DP-2060 (Tufo di Bagni Albula, 535 analysed by mass spectrometry using a VG5400 equipped with an electron multiplier Balzers 537 using co-irradiated Alder Creek Sanidine (ACs-2) standard with an age of 1.1891Ma 47 and 539 541 543 0.00000160. Mass discrimination was estimated by analysis of air pipette throughout the Tephra OH-DP-2669 is a 2.5 cm thick, yellowish layer with sharp upper and lower micropumices. Its distinct trachytic composition (Extended Data Fig. 2) and the stratigraphic 545 524.84 ka 44; Extended Data Table1) narrow potential tephrostratigraphic equivalents. Tephra 547 with a similar trachytic composition50,51 for this interval (Extended Data Fig. 2, Extended 549 majority of the SC1-35.30/SUL2-1 and OH-DP-2669 analyses correlate well with the more 551 field of the less evolved group V5a (Extended Data Fig. 2, Extended Data Table 2). Tephra 546 layer SC1-35.30/SUL2-1 from the Sulmona basin in the Italian Apennines is the only tephra 548 Data Table 2). SC1-35.30/SUL2-1 was correlated with tephra V5 from the MJS 52,53. The 550 evolved group of V5 (V5b: SiO2 >63% wt.; CaO <1.5% wt.). Only few analyses plot in the 22 552 layer SUL2-1 and V5 were 40Ar/39Ar dated at 722.8±2.4 ka50 and 719.5±12.6 ka53, 554 The previous proposed correlation of SUL2-1/V5 with the Parmenide ash found in the 556 Parmenide ash (710±5 ka) 54,55,56 and the differences in the geochemical data to OH-DP-2669 558 Tephra OH-DP-2898 is a ~0.8 cm thick, whitish-yellowish band of lenses comprising 553 respectively. 555 Crotone basin50,52 is not considered here due to a slightly younger 40Ar/39Ar age of the 557 (Extended Data Fig. 2, Extended Data Table 2). 559 fine-grained glass shards with a high degree of vesicularity and a phonolitic composition 561 indicative for interglacial conditions20. The comparison of OH-DP-2898 glass composition 560 (Extended Data Fig. 2). It is located ~2 m below the M/B boundary, in calcareous sediments 562 with those of Sulmona tephra SUL2-19, -20, -25, -29 and -31 in a similar 564 tephra close to the M7B transition, SUL2-22, -23, and -27, have a composition similar to OH- 566 (Extended Data Fig. 2, Extended Data Table 2). SUL2-27 is geochemically indistinguishable 568 geochemically indistinguishable from OH-DP-2898 and shares a similar stratigraphic position 570 OH-DP-2898 with tephra V4 from the MJS is not possible due to differences in the 572 age of 773.9±1.3 ka of V4 52, quasi-synchronous position during the 10Be peak or M/B 563 magnetostratigraphic position exclude a correlation (Extended Data Fig. 2). Other Sulmona 565 DP-2898, but SUL2-23 has slightly lower alkali and higher CaO, FeO, TiO2 concentrations 567 from OH-DP-2898, but deposited in glacial sediments of the MIS 20 57. SUL2-22 is also 569 below the M/B boundary58,59 and at the transition from MIS 20 to MIS 19 57. A correlation of 571 compositional range (Extended Data Fig. 2, Extended Data Table 2) and a younger 40Ar/39Ar 573 transition60. Also a correlation of OH-DP-2898/SUL2-22 with tephra V3 of the MJS 575 Data Fig. 2, Extended Data Table 2) and deposition of V3 during glacial conditions of MIS 20 577 climatostratigraphic position 55,61,62, but differs geochemically from OH-DP-2898/SUL2-22. 574 (801.2±19.5 ka) is excluded due to differences in the geochemical composition (Extended 576 60 . The Pitagora ash from the Crotone basin is found in a similar magneto- and 23 578 Therefore, we regard a correlation of OH-DP-2898 with SUL2-22 as most robustly and use its 579 40 581 (Extended Data Table 1). This update includes the Campanian Ignimbrite (Y-5/OH-DP- 583 from the Sulmona section64. The tephrostratigraphy of the Fucino record65 improved and 585 of TF-17, correlated to OH-DP-0624, yielded a much more precise age of 158.8 ± 3.0 ka, 580 582 Ar/39Ar age of 791.9±1.9 ka 58 for our chronology. In addition to the new tephra correlations, we updated ages for the upper tephra layers 0169)63 and tephra layers OH-DP-0404/POP2 and OH-DP-0435/X-6, based on new results 584 reassessed the correlations established for OH-DP-0617 and OH-DP-0624 44. 40Ar/39Ar dating 586 which replaced the age of Vico B/OH-DP-0617 (162±6 ka)66. 588 provided a new chronological tie-point at 410±2 ka 67. The previously established correlation 587 589 Furthermore, the correlation of cryptotephra OH-DP-1700.6 with the Vico β eruption45 of tephra layer OH-DP-1955 with tephra layer SC-5 from the Mercure basin44 was rejected in 590 the light of its large uncertainty (±10.9 ka) and the new tephrostratigraphic data. 592 DP-2669, by updating the value of the atmospheric Ar-composition (40Ar/36Ar: 298.5 instead 594 Data Table 1) using the decay constant of 48 and an age of 1.1891 Ma for the ACs-2 flux 596 undistinguishable within uncertainty and thus used for our chronology. All other 40Ar/39Ar 597 598 ACs-2 (1.1891 Ma) and Fish Canyon sanidines (FCs) ages of 28.294 Ma. 599 600 Chronology 602 of the DEEP site sediment succession down to 447.12 mcd uses tephrochronological data44,45, 591 593 595 Reassessment of the raw Ar-isotope data of SC1-35.30/SUL2-1, the equivalent to OH- of 295.5 originally) and removing xenocrysts58 yielded a new age of 715.02±5.4 ka (Extended standard47. Our new 40Ar/39Ar age of V5 (716.2±5.4 ka; MSWD = 0.8, P = 0.7) is used were recalculated using the software ArAR68 with a given decay constant and age for 601 Following the methodological approach for the upper 247 mcd of the record20, the chronology 603 46 as 1st-order tie points and tuning of climate-sensitive proxy data (TOC; ~480 yr resolution) 24 604 against orbital parameters as 2nd-order tie points considering that maxima in TIC represent 606 DEEP site since the penultimate glacial period (Y-5, X-6, P-11, and A11/12) occur at depths 608 passages in March correspond to the inflection points of increasing local summer insolation 610 equinoxes) at the latitude of Lake Ohrid (41°N; Fig. 1). Increasing summer insolation 612 organic matter (OM) supply to the sediments. An extended winter season improves lake-water 614 Thus, minima in TOC result from moderate OM supply to the sediments and improved 616 resolution in the DEEP site record used for tuning purposes. 618 cryptotephra layers and 66 2nd-order tie points obtained from orbital tuning were cross 605 607 609 interglacial periods19,20. Some chronologically well-constrained tephra layers deposited at the where TOC shows minima at times of the perihelion passage in March20. These perihelion (21st June) and winter-season length (number of days between the September and March 611 promotes high summer temperatures, primary productivity in the water column and increases 613 mixing which enhances oxidation of OM in the water column and the surface sediments20. 615 oxidation of OM at the sediment surface and are due to their available high temporal 617 The independent chronological information obtained from the 16 tephra and 619 evaluated by the two paleomagnetic age constraints (base of the Jaramillo sub-Chron and 621 considering overall uniform (mem.strength=60, mem.mean=0.9, thick=80 cm) sedimentation 620 Matuyama/Brunhes M/B; Fig. 1). The age model was calculated using Bacon 2.2 69, 622 rates (acc.shape=1.5, acc. mean=20) at the DEEP site33. An error of ±2,000 years was applied 623 to the 2nd-order tie points to account for tuning inaccuracy. The 95% confidence intervals of 624 ages for specific depths produced by the Bacon Bayesian age modelling average at ±5,500 626 447.12 m of the DEEP site record covers the last 1.364 Myr, continuously. 628 Cave speleothem record70 and found agreement within errors of the chronologies. Arboreal 625 627 years with a maximum of ±10,680 years. The resulting chronology implies that the upper We20 evaluated the DEEP site’s chronology against the 0-160 ka U/Th dated Soreq 25 629 pollen (AP) percentages in the DEEP site record are also in agreement with those from the 630 orbitally-tuned Tenaghi Philippon record3 back to 1.364 Ma (Fig. 2). 632 Model simulations and forcing 634 impacts of orbital forcing, Northern Hemisphere (NH) ice sheets, and variations in 631 633 Transient simulations with the Earth system model LOVECLIM were conducted to study the 635 atmospheric greenhouse gases (GHGs) on glacial-interglacial climate change. 637 atmospheric component of LOVECLIM is the spectral T21, three-level model ECBilt72 based 639 ice component of LOVECLIM consists of a free-surface Ocean General Circulation Model 641 Atmosphere and ocean components are coupled through the exchange of freshwater and heat 636 LOVECLIM is a coupled ocean-atmosphere-sea ice-vegetation model71. The 638 on quasi-geostrophic equations extended by estimates of ageostrophic terms. The ocean-sea 640 with a 3°x3° horizontal resolution coupled to a dynamic-thermodynamic sea-ice model73. 642 fluxes. The vegetation model VECODE74 computes the evolution of terrestrial vegetation 644 The transient simulations of the last 784,000 years were forced by time-dependent 643 cover based on annual mean surface temperature and precipitation. 645 boundary conditions for orbital parameters, atmospheric GHG concentrations, NH ice sheet 647 calculated according to76. Atmospheric GHG concentrations were prescribed according to 646 orography, and albedo following the methodology described in75. The orbital forcing was 648 reconstructions from EPICA Dome C for CO277 as well as CH4 and N2O78. Orbital forcing and 649 atmospheric GHG concentrations were updated every model year. The effects of NH ice 651 with an acceleration factor of 5, which compresses 784,000 forcing years into 156,000 model 653 model simulation is an updated version of the one presented in75 and uses a higher climate 650 652 sheets on albedo and land topography were prescribed according to79. The forcing was applied years. This acceleration factor is appropriate for quickly equilibrating surface variables. The 26 654 sensitivity resulting in a better representation of the glacial-interglacial surface temperature 655 amplitude23. 657 described above (Extended Data Fig. 6). The sensitivity simulations cover the last four glacial 659 ice sheets and orbital parameters to glacial-interglacial climate change. The first sensitivity 661 atmospheric GHG concentrations. The “GHG effect” can then be calculated as the difference 663 simulation uses transient forcing as described above but constant PI NH ice sheets (extent and 665 simulation and this simulation. Two simulations were designed to study the role of orbital 667 used. However, one simulation was run under constant PI atmospheric CO2 concentration of 656 Four sensitivity simulations were conducted in addition to the full-forcing simulation 658 cycles (408,000 years) and aim at exploring the individual effects of atmospheric GHGs, NH 660 simulation uses transient forcing as described above but constant preindustrial (PI) 662 between the simulation using the full forcing and this simulation. The second sensitivity 664 albedo). The “NH ice sheet effect” is calculated as the difference between the full-forcing 666 forcing under warm and cold climate. For both simulations, transient orbital parameters are 668 280 ppm, whereas the second simulation uses a constant atmospheric CO2 concentration of 669 200 ppm resulting in a colder background climate. 671 Data analysis 673 deciduous oak pollen percentage data, a wavelet power spectrum was computed for the 675 interpolation) at 0.3 kyr (TIC) and 1.0 kyr (pollen), and subsequently submitted to continuous 677 approach by 81. Results of the CWT show persistent presence of 100 kyr and ~21 kyr orbital 670 672 To assess the temporal evolution of dominant periodicities in the DEEP site TIC and 674 respective time series. The time series were resampled at regular intervals (linear 676 wavelet transform (CWT, Morlet window) using PAST v.3.21 software80 following the 678 frequencies, and a clear presence of 41 kyr in the early half of the pollen record. Relative to 27 679 the pollen, the CWT results of the TIC show a more pronounced 100 kyr cyclicity over the 680 entire record, and less pronounced 21 kyr signals. 682 maxima against precession forcing, the bandpass-filtered 18–25 kyr component of the proxy 684 Partial least squares regression (PLSR) was used to test the correlation of TIC and 681 To quantitatively test the observed correlation between deciduous oak and TIC 683 data was regressed against precession based on the La2004 orbital solution30. 685 deciduous oaks as predictive variables with LOVECLIM temperature and precipitation output 687 filtered using a frequency centred at 0.05 and a bandwidth of 0.02 prior to multivariate 686 data. PLSR was performed using SIMCA 14 (Sartorius Stedim Biotech). All datasets were 688 statistical analysis to accommodate for slight age offsets between proxy and simulation data. 690 Methods and Extended Data files references: 689 691 31. Popovska, C. & Bonacci, O. Basic data on the hydrology of Lakes Ohrid and Prespa. 693 32. Matzinger, A., Spirkovski, Z., Patceva, S. & Wüest A. Sensitivity of ancient Lake Ohrid 692 694 695 Hydrol. Process. 21, 658–664 (2007). doi:10.1002/hyp.6252 to local anthropogenic impacts and global warming, J. Great Lakes Res. 32, 158–179 (2006). doi:10.3394/0380-1330(2006)32[158:SOALOT]2.0.CO;2 696 33. Wagner, B. et al. The SCOPSCO drilling project recovers more than 1.2 million years of 698 34. Farmer, V. C. The infrared spectra of minerals, edited by: V. C. Farmer, Mineralogical 700 35. Chukanov, N. V. Infrared Spectra of Mineral Species. Springer, Dordrecht, Heidelberg, 697 699 701 history from Lake Ohrid. Sci. Drill. 17, 19−29 (2014). doi:10.5194/sd-17-19-2014 Society Monograph 4, 227 pp, Adlard & Son, Dorking, Surrey, (1974). New York, London (2014). doi:10.1007/978-94-007-7128-4 28 702 703 704 36. Sadori, L. et al. Pollen-based paleoenvironmental and paleoclimatic change at Lake Ohrid (south-eastern Europe) during the past 500 ka. Biogeosciences 13, 1423–1437 (2016). doi:10.5194/bg-13-1423-2016 705 37. Beug, H.-J. Leitfaden der Pollenbestimmung für Mitteleuropa und angrenzende Gebiete. 707 38. Cheddadi, R. et al. Imprints of glacial refugia in the modern genetic diversity of Pinus 706 708 709 710 711 712 Verlag Dr. Friedrich Pfeil, München, Germany (2004). sylvestris. Global Ecol. Biogeogr. 15, 271–282 (2006). doi:10.1111/j.14668238.2006.00226.x 39. Rossignol-Strick, M. The Holocene climatic optimum and pollen records of sapropel 1 in the Eastern Mediterranean, 9000–6000 BP. Quat. Sci Rev. 18, 515–530 (1999). doi:10.1016/S0277-3791(98)00093-6 713 40. Langgut, D., Almogi-Labin, A., Bar-Matthews, M. & Weinstein-Evron, M. Vegetation 715 Interglacial cycle (86 ka): new marine pollen record. Quat. Sci Rev. 30, 3960–3972 717 41. Combourieu-Nebout, N. et al. Climate changes in the central Mediterranean and Italian 714 and climate changes in the South Eastern Mediterranean during the Last Glacial– 716 (2011). doi:10.1016/j.quascirev.2011.10.016 718 719 720 721 722 723 724 725 vegetation dynamics since the Pliocene. Rev. Palaeobot. Palynol. 218, 127–147 (2015). doi:10.1016/j.revpalbo.2015.03.001 42. Lacey, J. H. et al. Northern Mediterranean climate since the Middle Pleistocene: a 637 ka stable isotope record from Lake Ohrid (Albania/Macedonia). Biogeosciences 13, 1801−1820 (2016). doi:10.5194/bg-13-1801-2016 43. Just, J. et al. Environmental control on the occurrence of high-coercivity magnetic minerals and formation of iron sulfides in a 640 ka sediment sequence from Lake Ohrid (Balkans). Biogeosciences 13, 2093–2109 (2016). doi:10.5194/bg-13-2093-2016 29 726 727 728 729 730 731 732 733 734 735 736 737 44. Leicher, N. et al. First tephrostratigraphic results of the DEEP site record from Lake Ohrid (Macedonia and Albania). Biogeosciences 13, 2151–2178 (2016). doi:10.5194/bg13-2151-2016 45. Kousis, I. et al. Centennial-scale vegetation dynamics and climate variability in SE Europe during Marine Isotope Stage 11 based on a pollen record from Lake Ohrid. Quat. Sci. Rev. 190, 20−38 (2018). doi:10.1016/j.quascirev.2018.04.014 46. Francke, A. et al. Sediment residence time reveals Holocene shift from climatic to vegetation control on catchment erosion in the Balkans. Global Planet. Change 177, 186–200. 2019. doi:10.1016/j.gloplacha.2019.04.005 47. Niespolo, E. M., Rutte, D., Deino, A. L. & Renne, P. R. Intercalibration and age of the Alder Creek sanidine 40Ar/39Ar standard. Quat. Geochronol. 39, 205−213 (2017). doi:10.1016/j.quageo.2016.09.004 738 48. Renne, P. R., Balco, G., Ludwig, K. R., Mundil, R. & Min, K. Response to the comment 740 the Fish Canyon sanidine standard, and improved accuracy for 40Ar/39Ar geochronology” 742 doi:10.1016/j.gca.2010.06.017 739 by W. H. Schwarz et al. on “Joint determination of 40K decay constants and 40Ar∗/40K for 741 by P. R. Renne et al. (2010). Geochim. Cosmochim. Acta 75, 5097−5100 (2011). 743 49. Lee, J. Y. et al. A redetermination of the isotopic abundances of atmospheric Ar. 745 50. Giaccio, B. et al. Revised Chronology of the Sulmona Lacustrine Succession, Central 747 51. Giaccio, B. et al. Tephra layers from Holocene lake sediments of the Sulmona Basin, 744 746 748 Geochim. Cosmochim. Acta 70, 4507−4512 (2006). doi:10.1016/j.gca.2006.06.1563 Italy. J. Quat. Sci. 28, 545−551 (2013). doi:10.1002/jqs.2647 Central Italy: implications for volcanic activity in Peninsular Italy and tephrostratigraphy 30 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 in the Central Mediterranean area. Quat. Sci. Rev. 28, 2710−2733 (2009). doi:10.1016/j.quascirev.2009.06.009 52. Petrosino, P. et al. The Montalbano Jonico marine succession: An archive for distal tephra layers at the Early–Middle Pleistocene boundary in southern Italy. Quat. Internat. 383, 89−103 (2015). doi:10.1016/j.quaint.2014.10.049 53. Ciaranfi, N. et al. Integrated stratigraphy and astronomical tuning of Lower–Middle Pleistocene Montalbano Jonico section (Southern Italy). Quat. Internat. 219, 109−120 (2010). doi:10.1016/j.quaint.2009.10.027 54. Massari, F. et al. Interplay between tectonics and glacio-eustasy: Pleistocene succession of the Crotone basin, Calabria (southern Italy). Geol. Soc. Am. Bull. 114, 1183−1209 (2002). doi:10.1130/0016-7606(2002)114<1183:IBTAGE>2.0.CO;2 55. Capraro, L. et al. Climatic patterns revealed by pollen and oxygen isotope records across the Matuyama-Brunhes Boundary in the central Mediterranean (southern Italy). Geol. Soc., London, Spec. Publ. 247, 159−182 (2005). doi:10.1144/GSL.SP.2005.247.01.09 56. Capraro, L. et al. Chronology of the Lower–Middle Pleistocene succession of the southwestern part of the Crotone Basin (Calabria, Southern Italy). Quat. Sci. Rev. 30, 1185−1200 (2011). doi:10.1016/j.quascirev.2011.02.008 766 57. Giaccio, B. et al. Duration and dynamics of the best orbital analogue to the present 768 58. Sagnotti, L. et al. Extremely rapid directional change during Matuyama-Brunhes 767 769 770 interglacial. Geology 43, 603−606 (2015). doi:10.1130/G36677.1 geomagnetic polarity reversal. Geophys. J. Internat. 199, 1110−1124 (2014). doi:10.1093/gji/ggu287 31 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 59. Sagnotti, L. et al. How fast was the Matuyama–Brunhes geomagnetic reversal? A new subcentennial record from the Sulmona Basin, central Italy. Geophys. J. Internat. 204, 798−812 (2016). doi:10.1093/gji/ggv486 60. Simon, Q. et al. Authigenic 10Be/9Be ratio signature of the Matuyama–Brunhes boundary in the Montalbano Jonico marine succession. Earth Planet. Sci. Lett. 460, 255−267 (2017). doi:10.1016/j.epsl.2016.11.052 61. Rio, D. et al. Reading Pleistocene eustasy in a tectonically active siliciclastic shelf setting (Crotone peninsula, southern Italy). Geology 24, 743−746 (1996). doi:10.1130/00917613(1996)024<0743:RPEIAT>2.3.CO;2 62. Macri, P., Capraro, L., Ferretti, P. & Scarponi, D. A high-resolution record of the Matuyama-Brunhes transition from the Mediterranean region: The Valle di Manche section (Calabria, Southern Italy). Phys. Earth Planet. Inter. 278, 1−15 (2018). doi:10.1016/j.pepi.2018.02.005 63. Giaccio, B., Hajdas, I., Isaia, R., Deino, A. & Nomade, S. High-precision 14C and 40 Ar/39Ar dating of the Campanian Ignimbrite (Y-5) reconciles the time-scales of climatic-cultural processes at 40 ka. Sci. Rep. 7, 45940 (2017). doi:10.1038/srep45940 64. Regattieri, E., et al. A last Interglacial record of environmental changes from the Sulmona Basin (central Italy). Palaeogeogr. Palaeoclimatol. Palaeoecol. 472, 51−66 (2017). doi:10.1016/j.palaeo.2017.02.013 65. Giaccio, B. et al. First integrated tephrochronological record for the last ~190 kyr from the Fucino Quaternary lacustrine succession, central Italy. Quat. Sci. Rev. 158, 211−234 (2017). doi:10.1016/j.quascirev.2017.01.004 66. Laurenzi, M. A. & Villa, I. 40Ar/39Ar chronostratigraphy of Vico ignimbrites. Period. Mineral. 56, 285−293 (1987) 32 795 796 797 798 799 800 801 802 803 804 67. Karner, D. B., Marra, F. & Renne, P. R. The history of the Monti Sabatini and Alban Hills volcanoes: groundwork for assessing volcanic-tectonic hazards for Rome. J. Volcanol. Geotherm. Res. 107, 185−219 (2001). doi: 10.1016/S0377-0273(00)00258-4 68. Mercer, C. M. & Hodges, K.V. ArAR — A software tool to promote the robust comparison of K–Ar and 40Ar/39Ar dates published using different decay, isotopic, and monitor-age parameters. Chem. Geol. 440, 148−163 (2016). doi:10.1016/j.chemgeo.2016.06.020 69. Blaauw, M. & Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. Bayes. Analys. 6, 457−474 (2011). doi:10.1214/ba/1339616472 805 70. Grant, K. M. et al. Rapid coupling between ice volume and polar temperature over the 807 71. Goosse, H. et al. Description of the Earth system model of intermediate complexity 806 808 809 810 811 812 813 814 815 816 817 818 past 150,000 years. Nature 491, 744–747 (2012). doi:10.1038/nature11593 LOVECLIM version 1.2. Geosci. Model Dev. 3, 603–633 (2010). doi:10.5194/gmd-3603-2010 72. Opsteegh, J. D., Haarsma, R. J., Selten, F. M. & Kattenberg A. ECBILT: a dynamic alternative to mixed boundary conditions in ocean models. Tellus, Ser. A, Dyn. Meterol. Oceanogr. 50, 348–367 (1998). doi:10.3402/tellusa.v50i3.14524 73. Goosse, H. & Fichefet, T. Importance of ice-ocean interactions for the global ocean circulation: A model study. J. Geophys. Res. 104, 23337–23355 (1999). doi:10.1029/1999JC900215 74. Brovkin, V., Ganopolski, A. & Svirezhev, Y. A continuous climate-vegetation classification for use in climate-biosphere studies. Ecol. Modell. 101, 251–261 (1997). doi:10.1016/S0304-3800(97)00049-5 33 819 75. Timmermann, A. et al. obliquity and CO2 effects on Southern Hemisphere climate during 821 76. Berger, A. Long-term variations of daily insolation and Quaternary climate change. J. 820 822 823 the past 408 ka. J. Clim. 27, 1863–1875 (2014). doi:10.1175/JCLI-D-13-00311.1 Atmos. Sci. 35, 2362–2367 (1978). doi:10.1175/15200469(1978)035<2362:LTVODI>2.0.CO;2 824 77. Lüthi, D. et al. High-resolution carbon dioxide concentration record 650,000-800,000 826 78. EPICA community members. Eight glacial cycles from an Antarctic ice core. Nature 429, 828 79. Ganopolski, A. & Calov, R. The role of orbital forcing, carbon dioxide and regolith in 830 80. Hammer, O. PAleontological Statistics (PAST) Version 3.21 reference manual, Natural 832 81. Torrence, C. & Compo, G. P. A practical guide to wavelet analysis. Bull. Am. Meteorol. 834 82. Lindhorst, K. et al. Sedimentary and tectonic evolution of Lake Ohrid 836 83. Melard, G. Algorithm AS 197: A fast algorithm for the exact likelihood of 825 827 829 831 833 835 837 838 years before present. Nature 453, 379–382 (2008). doi:10.1038/nature06949 623–628 (2004). doi:10.1038/nature02599 100 kyr glacial cycles. Clim. Past, 7, 1415–1425 (2011). doi: 10.5194/cp-7-1415-2011 History Museum, University of Oslo (2018). https://folk.uio.no/ohammer/past/ Soc. 79, 61−78. (1998). doi:10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2 (Macedonia/Albania). Basin Res. 27, 84−101 (2015). doi:10.1111/bre.12063 autoregressive-moving average models. Appl. Stat. 33,104–114 (1984). doi:10.2307/2347672 839 84. Bar-Matthews, M., Ayalon, A., Gilmour, M., Matthews, A. & Hawkesworth, C. J. Sea- 841 Eastern Mediterranean region and their implication for paleorainfall during interglacial 840 land oxygen isotopic relationships from planktonic foraminifera and speleothems in the 34 842 843 844 intervals. Geochim. Cosmochim. Acta 67, 3181−3199 (2003). doi: 10.1016/S00167037(02)01031-1 85. Zanchetta, G. et al. Aligning and synchronization of MIS5 proxy records from Lake 845 Ohrid (FYROM) with independently dated Mediterranean archives: implications for 847 2016 846 848 849 850 851 852 853 854 855 856 857 DEEP core chronology. Biogeosciences, 13, 2757–2768 (2016). doi:10.5194/bg-13-2757- 86. Le Bas, M. J. Le Maitre, R. W. Streckeisen, A. & Zanettin, B. A Chemical Classification of Volcanic Rocks Based on the Total Alkali-Silica Diagram. J. Petrol. 27, 745–750, (1986). doi:10.1093/petrology/27.3.745 87. Wagner, B. et al. The last 40 ka tephrostratigraphic record of Lake Ohrid, Albania and Macedonia: a very distal archive for ash dispersal from Italian volcanoes. J. Volcanol. Geotherm. Res. 177, 71–80 (2008). doi:10.1016/j.jvolgeores.2007.08.018. 88. Zanchetta, G. et al. Tephrostratigraphy, chronology and climatic events of the Mediterranean basin during the Holocene: An overview. Holocene 21, 33–52 (2011). doi:10.1177/0959683610377531 89. Siani, G., Sulpizio, R., Paterne, M. & Sbrana, A. Tephrostratigraphy study for the last 858 18,000 C-14 years in a deep-sea sediment sequence for the South Adriatic. Quat. Sci. 860 90. Albert, P. G. et al. Revisiting the Y-3 tephrostratigraphic marker: a new diagnostic glass 859 861 862 863 864 865 Rev. 23, 2485–2500 (2004). doi:10.1016/j.quascirev.2004.06.004. geochemistry, age estimate, and details on its climatostratigraphical context, Quat. Sci. Rev. 118, 105–121 (2015). doi:10.1016/j.quascirev.2014.04.002 91. Satow, C. et al. A new contribution to the Late Quaternary tephrostratigraphy of the Mediterranean: Aegean Sea core LC21, Quat. Sci. Rev. 117, 96–112 (2015). doi:10.1016/j.quascirev.2015.04.005 35 866 867 868 869 92. Giaccio, B. et al. Isotopic (Sr-Nd) and major element fingerprinting of distal tephras: an application to the Middle-Late Pleistocene markers from the Colli Albani volcano, central Italy. Quat. Sci. Rev. 67, 190–206 (2013). doi: 10.1016/j.quascirev.2013.01.028 93. Petrosino, P., Jicha, B. R., Mazzeo, F. C. & Russo Ermolli, E. A high resolution 870 tephrochronological record of MIS 14–12 in the Southern Apennines (Acerno Basin, 872 doi:10.1016/j.jvolgeores.2014.01.014 871 873 Italy). J. Volcanol. Geotherm. Res. 274, 34–50 (2014). 94. Marra, F., Karner, D. B., Freda, C., Gaeta, M. & Renne, P. Large mafic eruptions at 874 Alban Hills Volcanic District (Central Italy): Chronostratigraphy, petrography and 876 doi:10.1016/j.jvolgeores.2008.11.009 875 eruptive behavior. J. Volcanol. Geotherm. Res. 179, 217–232 (2009). 877 878 879 Data Availability 881 database at https://doi.pangaea.de/10.1594/PANGAEA.896848. Data used for LOVECLIM 880 Data are available in the main text, in the supplementary materials and in the Pangaea 882 are available at https://climatedata.ibs.re.kr/grav/data/loveclim-784k. 884 Code Availability 886 the IBS Center for Climate Physics: https://climatedata.ibs.re.kr/grav/data/loveclim-784k. 883 885 Model data produced by the LOVECLIM simulations are available through the data centre of 887 Additional data are available upon request made to Tobias Friedrich (tobiasf@hawaii.edu). 889 Extended Data Legends 888 36 890 Extended Data Figure 1 | Map of Lake Ohrid and its surrounding area. Geology, 891 topography, and bathymetry compiled from19,82 and geological maps of Albania and North 893 of 293 m. The water depth at the DEEP drill site is 240 m. 895 Extended Data Figure 2 | Correlation of tephra layers at the DEEP site with tephra 897 Al2O3, (c) CaO vs. TiO2, (d) Na2O vs. K2O, and (e) total alkali vs. silica (TAS) diagram86 892 894 Macedonia. The lake is located at an altitude of 693 m a.s.l. and has a maximum water depth 896 layers found in mid-distal records. Bi-oxide plots of (a) CaO vs. FeOtotal, (b) CaO vs. 898 show the correlation of OH-DP-2669 with the tephra layers SC1-35.30/SUL2-1/V5 and the 900 CaO vs. TiO2, (i) Na2O vs. K2O, and (k) TAS diagram show the correlation of OH-DP-2898 Error bars of the Parmenide Ash refer to 54.Tephra ages, geochemical data, tephrostratigraphic 899 differences to the Parmenide ash. Bi-oxide plots of (f) CaO vs. FeOtotal, (g) CaO vs. Al2O3, (h) 901 with tephra SUL2-22 and the differences to SUL2-23, -27, -31, V4, V3, and the Pitagora ash. 902 903 discussion and references are provided in Extended Data Tables 1 and 2 and in Methods. 904 905 Extended Data Figure 3 | Lake Ohrid LOVECLIM simulation data and sedimentary 907 for the Lake Ohrid grid cell from the LOVECLIM simulation; (b) simulated precipitation 909 organic carbon (TOC) concentrations; (d) Lake Ohrid δ13C endogenic calcite in ‰ relative to 911 sedimentary quartz content; (g) Lake Ohrid K intensities in kilo counts and displayed using a 913 mean; (i) Lake Ohrid Ca intensities in kilo counts and displayed using a 11 pt running mean; 906 paleoclimate and paleoenvironment proxies. (a) Simulated surface-air temperature (SAT) 908 amount for the Lake Ohrid grid cell from the LOVECLIM simulation; (c) Lake Ohrid total 910 VPDB; (e) Lake Ohrid δ18O endogenic calcite in ‰ relative to VPDB; (f) Lake Ohrid relative 912 914 915 11pt running mean; (h) Lake Ohrid ratio of Ca/K intensities displayed using a 11 pt running (k) Lake Ohrid total inorganic carbon (TIC) concentrations; (l) Lake Ohrid deciduous oaks pollen percentages; (m) Lake Ohrid arboreal pollen excluding Pinus pollen (AP-P) 37 916 percentages; red and white diamonds indicate the position of radiometrically dated tephra 918 Ohrid sediment record. (b), (d), (e), (K), (l) and (m) are from Fig. 2. 920 Extended Data Figure 4 | Data analysis. Continuous wavelet transform results for 922 Ohrid DEEP where colour represents the signal amplitude at a given time and spectral period 924 test according to 81) against a red-noise background spectrum with autocorrelation coefficient 926 in PAST (80 based on 83). Thick grey line indicates the “cone of influence” outside of which 928 pass-filtered 18-25 ky component of (c) % TIC and (d) the % deciduous oak against 917 919 layers, blue and white diamonds the position of reversals of Earth’s magnetic field in the Lake 921 percentages of total inorganic carbon (TIC; a) and deciduous oak pollen (b) time series from 923 (yellow highest, red lowest power). Black contour is the 5% significance level (chi-squared 925 of 0.95, estimated through an autoregressive–moving-average (ARMA) model implemented 927 boundary effects can influence the results. Least squares regression (red line) between band 929 precession at 1 ky resolution. Blue lines indicate 95% bootstrapped (n=1999) confidence 931 response (steeper slope) of the deciduous oaks. Partial least squares regression (PLSR) using 933 precipitation output data as observations. PLSR was performed using SIMCA 14 (Sartorius 935 of 0.02 prior to multivariate statistical analysis to accommodate for slight offsets in age 937 correlations of simulated temperatures (e) and of simulated precipitation (f) data to proxy 939 compared to temperature. 930 intervals. Results show significant negative relationships for both proxies, with a stronger 932 TIC and deciduous oaks as predictive variables and LOVECLIM (e) temperature and (f) 934 Stedim Biotech). All datasets were filtered using a frequency centred at 0.05 and a bandwidth 936 differences between proxy and simulation data. Results show highly significant positive 938 data, with a higher sensitivity of TIC and deciduous oaks towards changes in precipitation 940 38 941 Extended Data Figure 5 | Lake Ohrid precipitation indicators and global monsoon 943 based on Soreq Cave speleothem δ18O data and U/Th chronology84; (b) simulated 945 Ohrid deciduous oaks pollen percentage; (d) Lake Ohrid total inorganic carbon (TIC) 942 944 records during MIS 5. (a) Ages of sapropels and humid phases in the Eastern Mediterranean precipitation amount for the Lake Ohrid grid cell from the LOVECLIM simulation; (c) Lake 946 concentrations; (e) Chinese Speleostack δ18O 25 in ‰ relative to VPDB; red and white 948 record. The chronology of the MIS 5 interval in the Lake Ohrid DEEP site record is based on 947 diamonds indicate the position of radiometrically dated tephra layers in the Lake Ohrid 949 85 951 Extended Data Figure 6 | Simulated Lake Ohrid precipitation for full-forcing run and 953 (black) and a simulation using only orbital forcing under a warm background climate (red). 955 climate (blue). (c) Black line as in (a) and a simulation using full-forcing except for a constant 950 . 952 sensitivity simulations. (a) Lake Ohrid precipitation (cm yr-1) for full-forcing simulation 954 (b) Black line as in (a) and a simulation using only orbital forcing under a cold background 956 preindustrial NH ice sheet. (d) Black line as in (a) and a simulation using full-forcing except 958 only cover the last 408 kyr. Please see Methods for details on the sensitivity simulations. 960 Extended Data Figure 7 | NOAA reanalysis data for the Mediterranean region. (a) 962 monthly means. Dashed line indicates two standard deviations above the mean. (b) Composite 964 maxima shown in (a) and referring to the months shown in (c). (c) Monthly distribution of 957 959 961 for constant preindustrial GHG concentrations. Please note that the sensitivity simulations Reconstructed precipitation (cm yr-1) for the Lake Ohrid reanalysis grid cell. Data based on 963 anomalies of 850 hPa geopotential height (m) associated with Lake Ohrid precipitation 965 precipitation maxima shown in (a). 966 39 967 Extended Data Figure 8 |. Mean seasonal cycle of Lake Ohrid precipitation - model 969 precipitation (cm yr-1) for all model years (green) and model years with annual-mean 968 simulation and NOAA reanalysis data. (a) Mean seasonal cycle of simulated Lake Ohrid 970 precipitation exceeding two standard deviations (magenta). Please see also Fig. 3a. (b) Mean (blue) and simulated for the 1–0 kyr period (red). The annual means were removed for better 971 seasonal cycle of Lake Ohrid precipitation (cm yr-1) derived from NOAA reanalysis data 972 973 comparison and are provided in the panel. 974 975 Extended Data Table 1 | Selected tephra layers from Lake Ohrid and their correlation 976 with tephra layers of other records. 40Ar/39Ar ages from literature were recalculated using a 978 (FCs) at 28.294 Ma73. Tephra ages in bold are used for age-depth modelling in Fig. 1. Age 979 uncertainties are provided according to the original reference (Reference age). 980 981 Extended Data Table 2 | Average compositions of OH-DP-2669 and OH-DP-2898 and 983 SC1-35.50 from 50; V5, V4, V3, Pitagora ash from 52 and the Parmenide ash from 54. x̅ = 977 982 984 decay constant73 and Alder Creek sanidine (ACs-2) at 1.1891 Ma74 or Fish Canyon sanidine potential equivalent correlations. Data of SUL2-1, SUL2-22, SUL2-23, SUL2-27 from 51; mean; S = standard deviation; n= number of analysis. 985 40 986 987 Fig. 1 988 989 41 Fig. 2 EM Sapropels a Speleostack b -10 + SE Asian monsoon strength - -8 -6 δ13Ccalcite (VPDB ‰) 0 2 2 L. Ohrid -2 c d 0 20 0 + temp./ karst runoff - TIC (%wt) L. Ohrid 8 4 f insolation difference g 260 LOVECLIM 270 30 + simulated precipitation - h 20 100 60 40 0 100 k 5 1 l 4 50 7 9 11 13 15 17 19 23 3 27 0 31 25 21 29 33 35 37 39 41 43 45 5 0 200 400 600 Age (ka) 42 800 1,000 1,200 1,400 AP (%) LR04 20 T. Philippon δ18Obenthic (VPDB ‰) + tree cover - L. Ohrid 80 i AP-Pinus (%) precipitation (cm yr-1) higher N. Hemisphere winter insol. 250 Tropics - Arctic (W m-2) 0 lower L. Ohrid 40 e deciduous oaks (%) 60 + soil moisture - 992 993 994 Medstack +(-) SST (SSS) -(+) δ18Oplanktonic (VPDB ‰) δ18O(VPDB ‰) 990 991 ice volume + Fig. 3 Lake Ohrid annual-mean precipitation (cm yr -1) 995 996 40.0 36.0 32.0 28.0 24.0 20.0 16.0 (a) 100 0 200 300 400 Age (ka) 500 700 600 60°N (b) 50°N 11 40°N 2 30°N 20°N -7 -16 1 m s-1 GPH-anomaly and wind-anomaly (800 hPa, SON) for precipitation above 2xSTD 10°N 997 20 30°W 15°W 0° 15°E 43 30°E 45°E 60°E -25 m