# Moist and warm conditions in Eurasia during the last glacial of the Middle Pleistocene Transition

María Fernanda Sánchez Goñi, Thomas Extier, Josué Polanco-Martínez, Coralie Zorzi, Teresa Rodrigues, André Bahr

► To cite this version:

María Fernanda Sánchez Goñi, Thomas Extier, Josué Polanco-Martínez, Coralie Zorzi, Teresa Rodrigues, et al.. Moist and warm conditions in Eurasia during the last glacial of the Middle Pleistocene Transition. Nature Communications, 2023, 14 (1), pp.2700. 10.1038/s41467-023-38337-4. hal-04374190

HAL Id: hal-04374190

https://hal.science/hal-04374190v1

Submitted on 5 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Moist and warm conditions in Eurasia during the last glacial of the Middle Pleistocene Transition

Received: 21 November 2022

Accepted: 26 April 2023

Published online: 10 May 2023

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

Check for updates

María Fernanda Sánchez Goñi<fnref n="1" />,²✉, Thomas Extier ②,
Josué M. Polanco-Martínez ③,⁴, Coralie Zorzi ②, Teresa Rodrigues⁵ &
André Bahr ⑥

The end of the Middle Pleistocene Transition (MPT, ~ 800-670 thousand years before present, ka) was characterised by the emergence of large glacial ice-sheets associated with anomalously warm North Atlantic sea surface temperatures enhancing moisture production. Still, the direction and intensity of moisture transport across Eurasia towards potential ice-sheets is poorly constrained. To reconstruct late MPT moisture production and dispersal, we combine records of upper ocean temperature and pollen-based Mediterranean forest cover, a tracer of westerlies and precipitation, from a subtropical drill-core collected off South-West Iberia, with records of East Asia summer monsoon (EASM) strength and West Pacific surface temperatures, and model simulations. Here we show that south-western European winter precipitation and EASM strength reached high levels during the Marine Isotope Stage 18 glacial. This anomalous situation was caused by nearly-continuous moisture supply from both oceans and its transport to higher latitudes through the westerlies, likely fuelling the accelerated expansion of northern hemisphere ice-sheets during the late MPT.

The strong non-linear behaviour of the climate system in response to the Milankovitch astronomical forcing is particularly well exemplified during the Middle Pleistocene Transition (MPT, ~1.2–0.67 Ma) when 100,000-year ice age cycles emerged without any significant change in orbital cyclicities dominated by the combination of 23,000 and 41,000-year cycles<fnref n="1" />. Imbrie et al.'s conceptual model² shows that the growth and decay of ice-sheets amplify, through albedo, winds and ocean mass properties, the system's initial modest responses to orbital forcing by changing the transport of heat to boreal latitudes and atmospheric CO₂ concentration. At the end of the MPT (Marine Isotope Stage-MIS 19 to 17, ~800–670 ka), which is marked by the manifestation of large continental ice-sheets, paradoxical responses are observed with warm sea surface temperatures (SST) in the North

Atlantic mid-latitudes generating excess moisture that fuelled ice-sheet growth in Europe and further north³⁻⁵. This climatic evolution culminated in the maximum southern extent of the Eurasian inland ice-sheets during MIS 16 at ~650 ka⁶ contemporaneously with a shift to greater ice accumulation in North America⁷. Within the end of the MPT, MIS 18 (centred at ~730 ka) is a warm glacial in the North Atlantic Ocean³,⁸ and represents an anomaly during this period. However, the cause of the mild expression of MIS 18 remains unclear yet. Its glacial limits are unknown and geomorphological data suggests that the glaciation was limited to the high latitudes, especially in the Southern Hemisphere⁹. Furthermore, little is known about the Eurasian environments and atmospheric climate during the final stages of the MPT. We hypothesise that the combination of unusual warming in the mid/

<footnote n="1">¹Ecole Pratique des Hautes Etudes (EPHE, PSL University), Paris, France. ²Univ. Bordeaux, CNRS, Bordeaux INP, EPOC, UMR 5805, 33600 Pessac, France. ³Unit of Excellence GECOS, IME, University of Salamanca, 37007 Salamanca, Spain. ⁴Basque Centre for Climate Change (BC3), 48940 Leioa, Spain. ⁵Divisão de Geologia e Georecursos Marinhos, Instituto Português do Mar e da Atmosfera, Rua Alfredo Magalhães Ramalho, 6, 1495-006 Lisboa, Portugal. ⁶Institute of Earth Sciences, Heidelberg University, Im Neuenheimer Feld, 234, 69120 Heidelberg, Germany. ✉ e-mail: maria.sanchez-goni@u-bordeaux.fr</footnote>
high-latitude North Atlantic, strong SST gradient in the tropical West
North Pacific Ocean and intermediate glacial background conditions
should have enhanced long-term moisture production leading to
interglacial climates in the subtropical latitudes of Eurasia during MIS
18, and thus ice accumulation in the southern sectors of Eurasia.

Here, we document the Eurasian environments and climates
during the final part of the MPT using high-resolution pollen-based
vegetation data, a proxy for precipitation, and alkenone-based SST
from IODP Site U1385 located in the South-West Iberian margin
(37°34.285′N, 10°7.562′W, 2578 m depth) combined with published sea
subsurface temperature records from the same site³ (Fig. 1). These
records are compared with geochemical¹⁰ and grain size¹¹,¹² loess
records accounting for the evolution of the East Asian monsoon in
North-East China, and alkenone-based SST for the tropical and sub-
polar West North Pacific Ocean¹³,¹⁴ (Fig. 1). Robust statistical analysis
techniques applied to the paleoclimatic records combined with iLO-
VECLIM model simulations were used to identify the main atmospheric
and oceanographic processes underlying the warm ocean-maximum
ice growth paradox at the end of the MPT (Methods).

**Results and discussion**

Interglacial climate in South-West Europe during MIS 18 glacial

We qualitatively estimate forest cover changes using arboreal versus
non-arboreal pollen percentages. This simple approach that acknowl-
edges the non-linear relationship between plant abundances and pol-
len percentages¹⁵, is supported by a wide array of studies on the
modern relationship between pollen assemblages and vegetation (e.g.,
refs. 16,17) and, specifically, by the recent work of the pollen repre-
sentation of the vegetation in the Tagus basin (South-West Iberia)¹⁸,
the most important source of pollen preserved at IODP Site U1385. We
infer the strong development of the Mediterranean forest cover,
mainly composed of deciduous *Quercus* (-40–60 pollen %) and scler-
ophyllous trees and shrubs (evergreen *Quercus*, *Olea*, *Phillyrea*, *Pistacia*
and *Cistus*, -10 pollen %), during the MIS 18d-b interval, 740 to 725 ka,
with a brief forest setback at around 733 ka (MIS 18c) (Fig. 2, Supple-
mentary Fig. S1). This interval is bracketed by the dominance of semi-
desert landscapes and heathlands indicating cold-dry and cold-humid
climates during MIS 18e and MIS 18a, respectively. The forest pollen
percentage reached 60 % at c. 730 ka, a similar value characterising

other interglacials in South-West Iberia such as MIS 13 (-500 ka¹⁹), the last interglacial (MIS 5, -130 ka²⁰) and the Holocene²¹. Moreover, our data record all the phases characterising the succession of vegetation during interglacial periods in the Mediterranean region (*Juniperus* and *Betula* pioneer woodlands/deciduous and evergreen *Quercus-Olea-Phillyrea-Pistacia* forest/deciduous and evergreen *Quercus* woodlands) (Fig. S1)²². Site U1385 indicates that the expansion of the deciduous oaks parallels the development of the Mediterranean sclerophyllous taxa, with a synchronous maximum expansion. The pollen assemblages recorded at this site during the warm phases are very similar to the modern ones inferred from samples collected in the deciduous oak woods of the Tagus basin characterised by -10–20% of sclerophyllous pollen taxa from the thermomediterranean belt and -40–60% of deciduous oak pollen taxa from the mesomediterranean belt¹⁸. No arboreal trees extinct at present in Europe are recorded during this interval at Site U1385. The vegetation inferred from the pollen assemblages at Site U1385 is, therefore, very close to the present-day Mediterranean forest of South-West Iberia, dominated by broadleaf trees, such as the oak and mixed sclerophyll forests²³.

Based on transient model simulations with time-varying insolation
and atmospheric CO₂ concentrations²¹, the strong Mediterranean for-
est cover development reflects high amount of winter precipitation,
and thus relatively zonal and weak westerlies directed towards
southern Europe. This interpretation is supported by the high statis-
tical correlation between the NAO index, controlling the intensity and
direction of the westerlies, i.e., winter precipitation, and the present-
day Mediterranean forest cover changes²⁴. The moderate values (-10%)
of sclerophyllous plants (evergreen *Quercus*, *Olea*, *Pistacia*, *Phillyrea*
and *Cistus*) reflect furthermore a weak seasonal climate characterised
by relatively cold and wet summers, compared to other interglacials
and, particularly, to MIS 5, marked by sclerophyllous plant values
peaking at 20%. Furthermore, the record of a few pollen grains of the
summer-drought intolerant tree *Castanea*²⁵ at c. 725 ka highlights the
temperate and humid climate during MIS 18, with the westerlies
affecting southern Iberia all over the year. SST at the South-West
Iberian margin also reached the highest values, -18 °C, (Fig. 2) and no
substantial freshwater fluxes affected this margin during this interval
as indicated by the %C₃₇.₄ record (Fig. 2). Upper ocean warming off
South-West Iberia during MIS 18 was caused by enhanced subsidence

<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Location of the sites discussed in the text. IODP Site U1385 (this study), ODP Sites 980, 981, 983, 984, U1313, U1314 and DSDP 607 (ref. 4), Lake Ohrid⁴⁰, Tenaghi Philippon³⁹, Chinese Loess Plateau selected sequences: Gulang¹⁰ and Pingliang, Jingshuan, Baoji, Lingtai, Puxian¹¹,¹², Jingyuan⁴¹, ODP Site 882 (ref. 14), ODP Site 1143 (ref. 13). NAC North Atlantic Current, STG Subtropical Gyre, SPG: Subpolar Gyre,</figcaption></figure>

ENACW: Eastern North Atlantic Central Waters, KC Kuroshio Current, EASM East Asia summer monsoon, EAWM East Asia winter monsoon, SH Siberian High, after 78. Green arrow: present-day westerlies during a low index of the North Atlantic Oscillation (centred at -40°-45°N⁴⁹), reaching more southern latitudes during the glacial periods as shown by data and model simulations for the Last Glacial Maximum⁷⁹.
<figure><img src="image_5.png" /></figure>

**Fig. 2 |** *Paleoclimatic records of Marine Isotope Stage (MIS) 18 from the Iberian margin. a* insolation at 65°N in July (black), obliquity (blue) and precession (red)²⁶, summer energy (dashed black line), T350 defines the number of summer days in which daily insolation is above 350 W/m² (ref. 27), **b** δ¹⁸O of benthic foraminifera from IODP Site U1385 (blue)³⁰, sub-stages of MIS 18 follow⁸⁰, **c** alkenone-based sea surface temperature (SST, red) and percentages of C₃₇:₄, tracer of freshwater pulses

in the Iberian margin, from IODP Site U1385 (black, this study), **d** temperature gradient between T thermocline from IODP Site U1385 (37°N) and SST from IODP Site U1313 (41°N)³, **e** total pollen concentrations (this study), **f** semi-desert (orange) and heathlands (Ericaceae, blue) pollen percentages, Mediterranean forest (green) pollen percentages mainly composed of deciduous oaks and sclerophyllous (evergreen *Quercus*, *Pistacia*, *Olea*, *Phillyrea and Cistus*, yellow) (this study).

of warm East North Atlantic Central Waters (ENACW) resulting in the generation of relatively warm glacial thermocline waters off Iberia³ (Fig. 2). The resultant extensive upper ocean heat reservoir provided a substantial source of heat and moisture at that time off Iberia that likely favoured the strong development of the regional forest cover during the MIS 18 glacial.

## A warm and wet increasing trend in South-West Europe from MIS 19 to MIS 17

The entire period from MIS 19 to MIS 17 was characterised by an enhanced subduction of warm mid-latitude surface waters to the thermocline off the Iberian margin giving rise to sustained high SST at Site U1385 (ref. 3, Fig. 3). The mid-latitude sourcing of thermocline waters at Site U1385 is illustrated by the gradient between thermocline temperature (T<sub>therm</sub>) at Site U1385 and SST at mid-latitude Site U1313 (41°N) (thereafter "T<sub>therm</sub> SST gradient"), with a low gradient indicating enhanced subduction of warm mid-latitude waters in contrast to a high gradient indicative of a subpolar source³. Gradients of ~6 °C during

warm phases of MIS 19-17 are among the lowest for the entire interval MIS 44-14 (ref. 3) and provide evidence for the high production of warm thermocline waters off Iberia. Atmospheric moisture provided by the accumulation of warm subtropical waters in the mid-latitude eastern North Atlantic could have contributed through the westerlies to the recorded extended mild-humid and forested conditions in central and South-East Europe from MIS 21 to 17³,⁵. Surprisingly, the subduction of warm waters to thermocline level during MIS 18 was similar and even stronger compared to both interglacials MIS 17 and 19, respectively, despite that insolation²⁶, summer energy²⁷, atmospheric CO₂ concentrations²⁸ and relative sea level²⁹,³⁰ were lower during MIS 18 (Fig. 3).

The comparison of the glacial MIS 18 forest expansion with those during the previous and succeeding interglacials is striking (Fig. 3). We observe that the strongest Mediterranean forest development occurred during MIS 17, centred at 700 ka, reflecting a high amount of regional winter precipitation³¹. In contrast, MIS 19 (~785 ka), under the influence of both similar ice volume and higher atmospheric CO₂
<figure><img src="image_6.png" /></figure>

**Fig. 3 | Global and Eastern North Atlantic paleoclimatic records for the interval Marine Isotope Stage (MIS) 19-17.** a insolation at 65°N in July (black), obliquity (blue) and precession (red)²⁶, b winter Latitudinal Insolation Gradient (LIG) between 15°N and 65°N (this study), c atmospheric CO₂ concentrations²⁸, d δ¹⁸O of benthic foraminifera from IODP Site U1385 (blue)³⁰, a-e: sub-stages of MIS 18⁸⁰, and relative sea level changes from ODP Site 1123 (ref. 29) compared to present-day (black), e temperature gradient between T thermocline from IODP Site U1385

(37°N) and sea surface temperature (SST) from IODP Site U1313 (41°N)³, f alkenone-based SST from IODP Site U1385 (this study), g simulated South-West Iberia winter rainfall (blue), Mediterranean forest pollen percentages (green) from IODP Site U1385, and simulated South-West Iberia tree fraction percentages (light green) (this study). Panels e, g straight lines indicate ordinary least squares fits. Panel e thick line indicates weighted average fit with a 5-sample window. Grey bands indicate the forested phases that define the terrestrial interglacials.

concentration, is marked by limited forest expansion³² indicating lower winter precipitation compared to MIS 17. More interestingly, the MIS 18 glacial was more forested, reflecting stronger winter rainfall, compared to the preceding MIS 19 despite that the latter interglacial was characterised by higher insolation, sea level, atmospheric CO₂ concentrations and similar warm SST compared to MIS 18 (Fig. 3). Previous paleoceanographic studies have shown that prior to MIS 16 the North Atlantic deep-water formation, releasing heat and humidity to the atmosphere and controlling the global climate, was located further south and west (Boreal heat pump) compared to the period after MIS 16 (Nordic heat pump)³³. This position allowed the arrival of a substantial amount of moisture to southern Europe during the interval MIS 19-17 but does not explain the wetter conditions during the MIS 18 glacial compared to the MIS 19 interglacial.

To ground-truth the proxy-based evidence, we utilised simulated annual and seasonal air temperatures and precipitation as well

as tree fraction from MIS 19 to 17 generated with the iLOVECLIM model solely forced with changes in insolation, ice-sheet reconstruction and greenhouse gas concentrations (Figs. 3, 4 and S6). The numerical model data in fact show high levels of winter rainfall for the three MISs and similar tree fraction percentages during MIS 18 and MIS 19 as inferred from the pollen data. In contrast to the proxy data, the simulated tree fraction is the weakest during MIS 17 (Fig. 3). The simulated winter rainfall is the highest during MIS 17 but shows an overall good agreement with the simulated tree fraction during MIS 18 (Fig. 3). This mismatch between model and proxy reconstructions could be explained by the difficulty in quantitatively estimating the forest cover from pollen data¹⁵. However, our approach identifies substantial qualitative forest cover and winter rainfall changes that tightly parallel changes in the Tₜherm SST gradient (Fig. 3). Thus, the mismatch could be as well the result of a feedback process that is not well reproduced in iLOVECLIM such as the poor
<figure><img src="image_7.png" /></figure>

**Fig. 4 |** Paleoclimatic records from Eastern North Atlantic margin and the Chinese Loess Plateau compared to subpolar and tropical West North Pacific sea surface temperature (SST) records for the interval Marine Isotope Stage (MIS) 19-17. a insolation at 65°N in July (black), and precession (red)²⁶, b SST record from ODP Site 882 located in the subpolar West North Pacific Ocean¹⁴, c SST record from ODP Site 1143 located in the tropical West North Pacific Ocean¹³, d Grain Size Loess (GSL) record from the Chinese sequences of Pingliang, Jingshuan, Baoji, Lingtai, Puxian indicating changes in East Asia winter monsoon (EAWM)¹¹, e Fe/K record from the Chinese Loess Plateau Gulang sequence indicating changes in East

Asia summer monsoon (EASM) (red)¹⁰, and iLOVECLIM simulated summer rainfall (grey, this study—the data have been smoothed using a 10-year running mean), f temperature gradient between T thermocline from IODP Site U1385 (37°N) and SST from IODP Site U1313 (41°N)³, g Mediterranean forest pollen percentages (green) from IODP Site U1385, and simulated South-West Iberia winter rainfall (blue) (this study). Panels **b**, **c**, **f**, **g**: straight lines indicate ordinary least squares fits. Panel **e**, **f**: thick lines indicate weighted average fit with a 23- and 5-sample windows, respectively. Grey bands indicate the MIS interglacials 19 and 17.

prediction of the ocean thermal gradient despite a robust SST
simulation.

In order to identify the main driver that may explain the increasing
trend of forest cover and winter rainfall from MIS 19 to MIS 17, and
particularly more rainfall during MIS 18 glacial compared to MIS 19
interglacial, we investigate four large-scale reconstructed climate
variables characterising the eastern North Atlantic (Fig. 3): (a) the
amount of winter precipitation using the Mediterranean forest pollen
percentage record (this study), (b) the position and intensity of the
westerlies using the winter latitudinal insolation gradient as proposed
by Davis and Brewer³⁴, (c) the position of the North Atlantic moisture
source using the T<sub>therm</sub> SST gradient³, and (d) the growth and decay of
the northern hemisphere ice-sheets using the benthic foraminifera
δ¹⁸O record from Site U1385 (ref. 30) that strongly parallels the LR04
δ¹⁸O stack curve³⁵. Since large-scale reconstructed climate time series

are unevenly spaced, we used the binned correlation technique³⁶,³⁷ to
convert the unevenly spaced paleoclimate time series to evenly spaced
time series (Fig. S4). To estimate the time-dependence of the strength
of association between the large-scale reconstructed time series and
the Mediterranean forest pollen (MF), we used a robust Pearson cor-
relation approach that takes into account the serially dependence of
the time series under study³⁸ (Fig. S4). The Pearson's correlation ana-
lysis shows a high indirect correlation, with a correlation coefficient
("r") statistically significant of -0.538 [-0.735; -0.258], between the
Mediterranean forest and the T<sub>therm</sub>-SST gradient record for the whole
studied interval, MIS 17, 18, and 19. The correlation is also statistically
significant and slightly higher, -0.586 [-0.781; -0.288] when MIS 18 is
taken individually. This indirect correlation indicates that the north-
ward displacement trend of the ENACW, weakening the T<sub>therm</sub>
SST gradient trend from MIS 19 to MIS 17, correlates well with the
long-term Mediterranean forest increase, i.e., winter rainfall, over this interval. An indirect correlation also exists between Mediterranean forest and ice-volume changes with a *r* of -0.623 [-0.829; -0.268]. No statistical significant correlation was found (*r* = -0.178 [-0.525; 0.220]) between Mediterranean forest cover and winter latitudinal insolation gradient. The progressive increase in warmer and wetter conditions in South-West Iberia, induced by the northward latitudinal shift of the oceanic moisture source to the Iberian margin, appears to be a distinctive feature of South-West Europe at the end of the MPT. Simulated relative humidity, i.e., the moisture content in the lowest atmospheric layer (at 800 hPa), between Sites U1313 and U1385 shows the northward migration of the moisture source from MIS 19 to MIS 18 although it does not reproduce the northern position of the moisture source during MIS 17 (Fig. S5) suggested from marine and pollen data. Pollen records from Albania/North Macedonia and North-East Greece show that, while high temperate forest cover characterised South-East Europe during MIS 18³⁹,⁴⁰, no linear increasing forest cover trend is identified in this region from MIS 19 to MIS 17 (Fig. S6).

## Contemporaneous strong moisture production in the North Atlantic and West North Pacific during the MIS 18 glacial

Interestingly, the long-term trend towards wetter winter conditions recorded in South-West Iberia from MIS 19 to 17 (Fig. 4) parallels the progressive increase in the intensity of the East Asian summer monsoon (EASM) (Fig. 4), as inferred from the Fe/K ratio, a proxy for summer monsoon-induced changes, in the Gulang loess section (37° 30' N, 102° 53' E, North-West margin of the Chinese Loess Plateau, Fig. 1)¹⁰, both regions being located at similar subtropical latitudes (-35-37°N). The 8¹³C of the loess carbonate from the Jingyuan sequence (36.35°N, 104.6°E) also shows summer precipitation during MIS 18 as high as during MIS 19⁴¹, also indicating the interglacial character of MIS 18 glacial. The EASM was as strong or stronger during the weaker summer insolation and more glaciated MIS 18 period than during the MIS 19 interglacial, and weaker compared to MIS 17 (Fig. 4). The intensity of the East Asian winter monsoon (EAWM), marked by cold and dry outbursts in the Chinese Loess Plateau related to the development of Siberian Highs⁴², was similar or more intense during MIS 18 compared to MIS 19 as shown by the mean grain size of quartz particles (MGSQ)¹² and grain size loess (GSL) stack records from the Chinese Loess Plateau¹¹ (Fig. 4 and S6). Simulated summer rainfall in East Asia shows an increase in EASM during MIS 18 although weaker compared to the reconstructed loess-based EASM (Fig. 4). The amount of oceanic moisture production in the tropical North Pacific across the MIS 19-17 interval could explain the observed long-term wetting trend in East Asia. The intense EAWM recorded during MIS 18 was associated with SST cooling in the high latitudes¹⁴ and strong warming in the tropical sector¹³ (Fig. 4 and S6) causing an increase of the meridional SST gradient in the West North Pacific Ocean. This SST gradient increase may lead to the strong production of water vapour in the Indian and Pacific tropical sectors enhancing the EASM over North-East China⁴³ as reflected by the Gulang and Jingyuan sections.

At present, a strong SST gradient exists off the coast of East Asia during the early summer monsoon as the result of the previous winter's cold air outbreaks increasing rainfall over southern Taiwan (24°N)⁴⁴. During the midsummer, the peak westerlies migrate north of the Plateau, the extratropical northerlies weaken, leaving only the monsoon low-level circulation that penetrates North-East China, which, coupled with stronger monsoonal southerlies, leads to the northward migration of the rain band⁴⁴⁻⁴⁶. The stronger EASM during MIS 18 compared to MIS 19 in North China would be the result of the penetration of higher amounts of oceanic moisture into East Asia channelled by the northward migration of the westerlies. The simulated relative humidity between ODP Sites 882 and 1143 shows indeed a northward position of the moisture source during both MIS 19 and MIS 18 (Fig. S5). During the MIS 18 glacial, the strong summer precipitation

in East Asia was contemporaneous with the high moisture production and winter rainfall recorded in the subtropical latitudes of the North Atlantic, an atmospheric configuration that is compatible with the strong development of the Siberian Highs⁴².

During MIS 17, interglacial conditions associated with higher insolation would produce the intensification of the SST gradient in the West North Pacific and Indian tropical sectors and, therefore the increase of the EASM amplified by the contemporaneous increase of cold outbursts (Fig. 4), compared to MIS 18. The stronger winter rainfall in South-West Europe during MIS 17 compared to MIS 18 (Fig. 3 and S6) from both pollen and climate simulation would be the result of the nearest position of the ENACW subduction centre to South-West Iberia caused by the progressive increase of the northward transport of heat and moisture by the North Atlantic Current into the Nordic Seas associated with the establishment of the Nordic heat pump after ~700 ka⁴,³³.

These results reveal for the first time the paradoxical interglacial character of MIS 18 glacial in the Eurasian subtropical latitudes, and the role of the oceanic moisture production in explaining the strong precipitation during a glacial period. Not only the direction and the intensity of the westerlies but the amount of oceanic moisture production should be taken into account to explain the recorded evolution of South-West European precipitation and EASM. Furthermore, sensitive model experiments show that during the intervals marked by northern ice-sheets of intermediate size, such as our studied interval, the ice growth rate will potentially increase with the strengthening of the North Atlantic Current as summer temperature remains at 0 °C at the south of the major ice-sheets while precipitation increase over this area at least in the Eurasian sector⁴. Therefore, the nearly continuous moisture production in the North Atlantic and West North Pacific subtropical sectors and its progressive northward transport by the westerlies during the end of the MPT may have substantially contributed to the strongest glaciation of the last millions of years in Eurasia and North America during MIS 16, leading to the strong dominant 100,000-year ice age cycles. This process could also be at work in previous moderate glacial/interglacial periods of the MPT. The long-term decreasing T<sub>therm</sub> SST gradient trend in the North Atlantic from MIS 25 to 22 is similar to that observed from MIS 19 to 17³ and could have led to the remarkable expansion of the Northern Hemisphere ice during MIS 22.

## Methods

### Environmental setting of IODP Site U1385

Site U1385 was recovered during IODP Expedition 339 "Mediterranean Outflow". The site is located on a spur, the Promontorio dos Principes de Avis, along the continental slope of the South-West Iberian margin, which is elevated above the abyssal plain avoiding the influence of turbidites³⁰. The water depth of Site U1385 places it under the influence of Northeast Atlantic Deep Water today, although it was influenced by southern sourced waters during glacial periods⁴⁷.

The South-West Iberian margin is located in the north-eastern edge of the subtropical gyre, under the influence of Eastern North Atlantic central water (ENACW). The surface water column is affected by the Portugal current (PC), which brings cold nutrient-rich water from the northern latitudes and forms the ENACW of subpolar origin (ENACWsp), and by the Azores current (AC) which brings warm water from the Azores front generating the ENACW of subtropical origin (ENACWst)⁴⁸. The general distribution of water masses is influenced by the seasonal migration of the Azores anticyclonic cell and its associated large-scale wind pattern.

Climate in South-West Iberia is directly affected by the intensity and direction of the westerlies that are, in turn, controlled by the North Atlantic Oscillation (NAO)⁴⁹,⁵⁰. During winter the North Atlantic westerlies bring moisture to South-West Iberia, while a high-pressure cell develops in the North Atlantic during summer, which generates strong northerly trade winds inducing coastal upwelling⁵¹. This climate
seasonality is characterised by wet and mild winters (Tmin: 5–1 °C; Tmax: 13–8 °C) and hot and dry summers, annual precipitation <600 mm⁵², and leads to the development of a Mediterranean vegetation in the adjacent landmasses dominated by deciduous oak at middle elevation, and evergreen oak, olive tree, *Pistacia*, *Phillyrea* and rockroses (*Cistus*) at lower elevations¹⁸. Ericaceae (heathlands) are restricted to rainy mountains (annual precipitation >600 mm), strong oceanic conditions and locations with wet soils⁵³.

## Stratigraphical framework

The sedimentary section recovered at Site U1385 (150 m-long compo-
site record) shows hemipelagic continental margin sediments depos-
ited under normal marine conditions with a fully oxygenated water
column and average sedimentation rates of 10 cm/ky⁵⁴. The stratigraphy
of Site U1385 was built upon a combination of chemo-stratigraphic
proxies³⁰. Ca/Ti ratio measured every cm in all holes by core scanning
X-ray fluorescence (XRF) was used to construct a composite section,
and the low resolution (20 cm) benthic foraminifera oxygen isotopic
record (δ¹⁸Obenth) was correlated to the marine δ¹⁸Obenth stack of LR04³⁵
to provide the age model that we present here³⁰.

## Pollen-based vegetation and climate reconstruction

Sediment subsamples of 1-cm thickness and 1.5–5 cm³ volume were prepared for pollen analysis using the standard protocol for marine samples (https://www.epoc-u-bordeaux.fr/index.php?lang=fr&page=eq_paleo_pollens), employing coarse-sieving at 150 µm, successive treatments with cold HCl, cold HF at increasing strength and micro-sieving (10 µm mesh). Known quantities of *Lycopodium* spores in tablet form were added to permit the calculation of pollen concentrations. Slides were prepared using a mobile mounting medium to allow rotation of the grains and counted using a Primo Star light microscope at 400 and 1000-fold magnifications for routine identification of pollen and spores, respectively. One hundred twenty-four samples were analysed every 4 cm (450-yr average temporal resolution). Pollen counts oscillate between 100 and 152 terrestrial pollen grains excluding *Pinus* (main pollen sum), aquatics and spores. The main pollen sum plus *Pinus* oscillates between 205 and 1438. The number of pollen morphotypes in most of the samples, 98 samples out from 124, ranges from 20 to 27, and from 16 to 19 morphotypes in the remaining samples. Pollen percentages for terrestrial taxa were calculated against the main sum of terrestrial grains, while percentages for *Pinus* were calculated against the main sum plus *Pinus*. Aquatic pollen and spore percentages are based on the total sum (Pollen + spores + indeterminables + unknowns) (Fig. S1). The uncertainties of the pollen percentages of the main ecological groups at 95% have been calculated using the “exactci” function in the PropCIs v. 0.3.0 package (Package “PropCIs”). The calculation of Cls is based on the ‘exact’ Clopper-Pearson method assuming a binomial proportion⁵⁵,⁵⁶. The average uncertainty of the calculated pollen percentage of the Mediterranean forest values in our analysis is less than 8% (Fig. S2). Total spores and pollen concentrations oscillate between 2000 and 28,000 grains.cm⁻³ (Fig. 2). Changes in grain concentrations do not parallel changes in pollen percentages and, therefore, these latter changes indicate actual variations in forest cover and composition. The interpretation of the pollen diagram was assisted by a constrained hierarchical cluster analysis based on Euclidean distance between samples (Fig. S1). Analysis was performed in the R environment v. 3.6.3 using the “chclust” function from the *Rioja* package⁵⁷. This analysis identifies three main pollen zones that we interpret as a strong expansion of the forest cover bracketed by two open vegetation phases.

## Sea surface temperature reconstruction and identification of freshwater pulses

Sea Surface Temperature (SST) was reconstructed using di- and tri-
unsaturated alkenones of 37 carbons atoms, which are organic

compounds synthesised by marine coccolithophore algae in a
temperature-related proportion. Alkenones, in particular the tetra-
unsaturated compounds (C₃₇.₄), can also be used to track episodes of
massive cold freshwater input (from iceberg melting or river dis-
charges), which are responsible for decreasing salinities in the surface
water masses⁸,⁵⁸. Alkenones are part of the total lipid extracted
(TLE) fraction which can be extracted from 2 g of sediment by soni-
cation with dichloromethane and hydrolysed with 6% potassium
hydroxide in methanol. After derivatization with bis(trimethylsilyl)tri-
fluoroacetamide, the TLE was identified using a Bruker mass spectro-
meter detector and quantified with a Varian gas chromatograph Model
3800 equipped with septum programmable injector and a flame
ionisation detector with a CPSIL-5 CB column. As a gas carrier was used
hydrogen at 2.5 ml.min⁻¹. Alkenone concentrations were determined
using n-hexatriacontane as an internal standard. Reproducibility tests
showed that uncertainty in the alkenone index Uк₃₇ determinations
were lower than 0.015 (ref. 59). The Uк₃₇ was calculated based on the
di- and tri- unsaturated concentrations⁶⁰, and converted into tem-
perature values using the global core top calibration of annual SST⁶¹.
The average temporal resolution between samples is 528 years and the
uncertainty of SST reconstruction is ±0.5 °C (ref. 62).

## Data analyses

To correlate the paleoclimate records, we used first the binned cor-
relation technique³⁶,³⁷ to convert the unevenly spaced paleoclimate
time series under analysis to evenly spaced time series. This technique
is based on a novel estimation approach proposed by ref. 36 for esti-
mating the correlation between two paleoclimate time series with
different timescales. The idea is that autocorrelation means that
memory enables values obtained on different time points to be cor-
related. Binned correlation is performed by resampling the unevenly
spaced paleoclimate time series under analysis into time bins on a
regular grid, assigning the mean values of the variable under scrutiny
within those bins. This method was recently implemented in the R
package BINCOR³⁶, which is freely available on CRAN. To estimate the
time-dependence of the strength of association between the large-
scale reconstructed time series and the Mediterranean forest pollen
(MF), we used a robust correlation approach that takes into account
the serially dependence of the time series under study. We have fol-
lowed the method developed by³⁸ to estimate the Pearson's correla-
tion coefficients that takes into account the serially dependence. This
method also includes an estimation of a confidence interval (95%)
obtained through a nonparametric stationary bootstrap technique
with an average block length proportional to the maximum estimated
persistence time of the data under analysis. This statistical technique is
implemented in the software PearsonT3, which is freely available from
http://www.climate-risk-analysis.com.

## iLOVECLIM model

We used the intermediate complexity climate model iLOVECLIM in version 1.1.5 (ref. 64), which is a development branch of the LOVECLIM model in version 1.2 described in Goosse et al.⁶³. iLOVECLIM does not present major differences in the physics of the atmosphere and ocean compared to LOVECLIM. The based components of this version are the atmosphere ECBilt⁶⁶, the ocean CLIO⁶⁷ and the vegetation and land surface VECODE⁶⁸. The modelled tree fraction in VECODE includes shrub plants such as Ericaceae⁶⁸. Oxygen isotopes (¹⁸O and ¹⁶O) have been implemented in the coupled model and evaluated against observed data in water samples and carbonates⁶⁹. An equilibrium run was first performed for 5000 years under climatic conditions at 800 ka. We then ran a transient run between 800 and 670 ka. The insolation⁷⁰, greenhouse gas concentrations²⁸,⁷¹ and prescribed ice-sheets reconstruction⁷² have been updated with an acceleration factor of 10 for the transient simulation. The model has been run at T21 spatial resolution (5.6° in longitude and latitude) and the output
are computed with an annual timestep and with a monthly timestep for the precipitation only. iLOVECLIM was previously successfully applied in the Asian monsoon region to investigate the oxygen isotopic composition of precipitation and calcite in association with monsoonal precipitation changes over the last 150 ka⁷³. The δ¹⁸O in precipitation computed in the model is extracted over the East Asia region (-25°N-41.5°N and -90°E-118°E). The mean δ¹⁸Ocalcite is then computed from δ¹⁸Oprecipitation, corrected from ice-sheet contribution to the global seawater δ¹⁸O (ref. 74) and from atmospheric temperature at 2 m height using the following Eq. (1) for synthetic calcite under equilibrium conditions⁷⁵.

$$1000 \ln \alpha (\text{calcite} - \text{H}_2\text{O}) = 18.03 \left( \frac{10^3}{T} \right) - 32.42, \quad (1)$$

where α is the fractionation factor and T the temperature in Kelvin. The values are then converted from SMOW scale to VPDB scale using the equation of Coplen et al. ⁷⁶ (2), and smoothed using a 10-year running mean:

$$\delta^{18}\text{O}_{\text{calciteVPDB}} = 0.97002 \delta^{18}\text{O}_{\text{SMOW}} - 29.98. \quad (2)$$

with 0.97002 the isotopic fractionation factor between VPDB and SMOW, and 29.98 expressed in per mil. The air temperature at 2 m height, annual precipitation and tree fraction are also computed over South-West Iberia region (-36°N-41.5°N and -5.6°W-11.2°W), which includes most of the entire Tagus basin from where the majority of pollen grains arrive to the South-West Iberia margin sediments¹⁸.

## Data availability

The elemental data that support the findings of this research are provided in PANGAEA database (https://doi.org/10.1594/PANGAEA.957610).

## Code availability

Codes for the statistical analyses (Bincor and Pearson's correlations) are available in the Supplementary Information. The iLOVECLIM source code and developments are hosted at http://forge.ipsl.jussieu.fr/ludus (ref. 77) but are not publicly available due to copyright restrictions. Access can be granted on demand by request to D. M. Roche (didier.roche@lsce.ipsl.fr) to those who conduct research in collaboration with the iLOVECLIM user group.

## References

1. Mudelsee, M. & Stattegger, K. Exploring the structure of the mid-Pleistocene revolution with advanced methods of time-series analysis. Geol. Rundsch. 86, 499-511 (1997).

2. Imbrie, J. et al. On the structure and origin of major glaciation cycles 2. The 100,000-year cycle. Paleooceanography 8, 699-735 (1993).

3. Bahr, A. et al. Oceanic heat pulses fueling moisture transport towards continental Europe across the mid-Pleistocene transition. Quat. Sci. Rev. 179, 48-58 (2018).

4. Barker, S. et al. Strengthening Atlantic Inflow Across the Mid-Pleistocene Transition. Paleooceanogr. Paleoclimatol. 36, e2020PA004200 (2021).

5. Catunda, M. C. A. et al. Subsurface heat channel drove sea surface warming in the high-latitude north atlantic during the mid-pleistocene transition. Geophys. Res. Lett. 48, e2020GL091899 (2021).

6. Batchelor, C. L. et al. The configuration of Northern Hemisphere ice sheets through the Quaternary. Nat. Commun. 10, 3713 (2019).

7. Hodell, D. A., Channell, J. E. T., Curtis, J. H., Romero, O. E. & Röhl, U. Onset of "Hudson Strait" Heinrich events in the eastern North Atlantic at the end of the middle Pleistocene transition (~640 ka)? Paleooceanography 23, PA4218 (2008).

8. Rodrigues, T. et al. A 1-Ma record of sea surface temperature and extreme cooling events in the North Atlantic: a perspective from the Iberian Margin. Quat. Sci. Rev. 172, 118-130 (2017).

9. Ehlers, J. & Gibbard, P. L. in Quaternary Encyclopaedia. 1023-1031 (Elsevier, 2007).

10. Sun, Y. et al. Persistent orbital influence on millennial climate variability through the Pleistocene. Nat. Geosci. 14, 812-818 (2021).

11. Ding, Z. L. et al. Stacked 2.6-Ma grain size record from the Chinese loess based on five sections and correlation with the deep-sea δ18O record. Paleooceanography 17, 5-1-5-21 (2002).

12. Sun, Y., Clemens, S. C., An, Z. & Yu, Z. Astronomical timescale and palaeoclimatic implication of stacked 3.6-Myr monsoon records from the Chinese Loess Plateau. Quat. Sci. Rev. 25, 33-48 (2006).

13. Li, L. et al. A 4-Ma record of thermal evolution in the tropical western Pacific and its implications on climate change. Earth Planet. Sci. Lett. 309, 10-20 (2011).

14. Martínez-Garcia, A., Rosell-Melé, A., McClymont, E. L., Gersonde, R. & Haug, G. H. Subpolar link to the emergence of the modern equatorial pacific cold tongue. Science 328, 1550-1553 (2010).

15. Zanon, M., Davis, B. A. S., Marquer, L., Brewer, S. & Kaplan, J. O. European Forest Cover During the Past 12,000 Years: A Palynological Reconstruction Based on Modern Analogs and Remote Sensing. Front. Plant Sci. 9, https://doi.org/10.3389/fpls.2018.00253 (2018).

16. Huntley, B. & Birks, H. J. B. An Atlas of Past and Present Pollenmaps for Europe: 0-13,000 B.P. years ago. (Cambridge University Press, 1983).

17. Prentice, I. C. Modern pollen spectra from lake sediments in Finland and Finnmark, north Norway. Boreas 7, 131-153 (1978).

18. Morales-Molino, C., Devaux, L., Georget, M., Hanquiez, V. & Sánchez Goñi, M. F. Modern pollen representation of the vegetation of the Tagus Basin (central Iberian Peninsula). Rev. Palaeobot. Palynol. 276, 104193 (2020).

19. Oliveira, D. et al. Combination of insolation and ice-sheet forcing drive enhanced humidity in northern subtropical regions during MIS 13. Quat. Sci. Rev. 247, 106573 (2020).

20. Sánchez Goñi, M. F., Eynaud, F., Turon, J.-L. & Shackleton, N. J. High resolution palynological record off the Iberian margin: direct land-sea correlation for the Last Interglacial complex. Earth Planet. Sci. Lett. 171, 123-137 (1999).

21. Oliveira, D. et al. Unraveling the forcings controlling the vegetation and climate of the best orbital analogues for the present interglacial in SW Europe. Clim. Dyn. 51, 667-686 (2018).

22. van der Hammen, T., Wijmstra, T. A. & Zagwijn, W. H. in The Late Cenozoic Glacial Ages (ed. Turekian, K. K.) 391-424 (Yale University press, 1971).

23. Polunin, O. & Walters, M. A Guide to the Vegetation of Britain and Europe. (Oxford University Press, 1985).

24. Gouveia, C., Trigo, R. M., DaCamara, C. C., Libonati, R. & Pereira, J. M. C. The North Atlantic Oscillation and European vegetation dynamics. Int. J. Climatol. 28, 1835-1847 (2008).

25. Gallardo-Lancho, J. F. Distribution of chestnut (Castanea sativa Mill.) forests in Spain: possible ecological criteria for quality and management (focusing on timber coppices). For. Snow Landsc. Res. 76, 477-481 (2001).

26. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. AA 428, 261-285 (2004).

27. Huybers, P. Early Pleistocene glacial cycles and the integrated summer insolation forcing. Science, https://doi.org/10.1126/science.1125249 (2006).

28. Luthi, D. et al. High-resolution carbon dioxide concentration record 650,000-800,000[thinsp]years before present. Nature 453, 379-382, http://www.nature.com/nature/journal/v453/n7193/suppinfo/nature06949_S1.html (2008).
29. Elderfield, H. et al. Evolution of ocean temperature and ice volume through the mid-Pleistocene climate transition. *Science* **337**, 704–709 (2012).

30. Hodell, D. A. et al. A reference time scale for site U1385 (Shackleton Site) on the Iberian Margin. *Glob. Planet. Change* **133**, 49–64 (2015).

31. Sánchez Goñi, M. F. et al. Pronounced northward shift of the westerlies during MIS 17 leading to the strong 100-kyr ice age cycles. *Earth Planet. Sci. Lett.* **511**, 117–129 (2019).

32. Sánchez Goñi, M. F. et al. Tropically-driven climate shifts in southwestern Europe during MIS 19, a low eccentricity interglacial. *Earth Planet. Sci. Lett.* **448**, 81–93 (2016).

33. Wright, A. K. & Flower, B. P. Surface and deep ocean circulation in the subpolar North Atlantic during the mid-Pleistocene revolution. *Paleoceanography* **17**, 20–21–20–16 (2002).

34. Davis, B. A. S. & Brewer, S. Orbital forcing and role of the latitudinal insolation/temperature gradient. *Clim. Dyn.* **32**, 143–165 (2009).

35. Lisiecki, L. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O records. *Paleoceanography* **20**, PA1003 (2005).

36. Mudelsee, M. *Climate Time Series Analysis: Classical Statistical and Bootstrap Methods*. First edition edn, (Springer-Verlag, 2010).

37. Mudelsee, M. *Climate Time Series Analysis. Classical Statistical and Bootstrap Methods*. Second edition edn, (Springer Cham, 2014).

38. Mudelsee, M. Estimating Pearson's correlation coefficient with bootstrap confidence interval from serially dependent time series. *Math. Geol.* **35**, 651–665 (2003).

39. 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).

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

41. Sun, Y. et al. Diverse manifestations of the mid-Pleistocene climate transition. *Nat. Commun.* **10**, 352 (2019).

42. Wu, B. & Wang, J. Winter arctic oscillation, Siberian High and East Asian winter monsoon. *Geophys. Res. Lett.* **29**, 31–34 (2002).

43. Chen, F. et al. East Asian summer monsoon precipitation variability since the last deglaciation. *Sci. Rep.* **5**, 11186 (2015).

44. Johnson, R. H. & Toy, M. D. in *The Global Monsoon System: Research and Forecast* (ed C. P. Chang et al.) 3–12 (World Scientific Publishing Co., 2017).

45. Chiang, J. C. H., Kong, W., Wu, C. H. & Battisti, D. S. Origins of East Asian summer monsoon seasonality. *J. Clim.* **33**, 7945–7965 (2020).

46. Yihui, D. & Chan, J. C. L. The East Asian summer monsoon: an overview. *Meteorol. Atmos. Phys.* **89**, 117–142 (2005).

47. Shackleton, N. J., Hall, M. A. & Vincent, E. Phase relationships between millennial scale events 64,000–24,000 years ago. *Paleoceanography* **15**, 565–569 (2000).

48. Ríos, A. F., Pérez, F. F. & Fraga, F. Water masses in the upper and middle North Atlantic Ocean east of the Azores. *Deep Sea Res. Part A. Oceanographic Res. Pap.* **39**, 645–658 (1992).

49. Hurrell, J. W. Decadal trends in the North Atlantic Oscillation: regional temperatures and precipitation. *Science* **269**, 676–679 (1995).

50. Trigo, R. M. et al. North Atlantic oscillation influence on precipitation, river flow and water resources in the Iberian Peninsula. *Int. J. Climatol.* **24**, 925–944 (2004).

51. Fiúza, A. F. D. G., Macedo, M. E. D. & Guerreiro, M. R. Climatological space and time variation of the Portuguese coastal upwelling. *Oceanologica Acta* **5**, 31–40 (1982).

52. Peinado Lorca, M. & Martínez-Parras, J. M. in *La vegetación de España* (eds. M. Peinado Lorca & S. Rivas Martínez) 163–196 (Universidad de Alcalá de Henares, 1987).

53. Loidi, J. in *La vegetación en España* (eds. Peinado Lorca, M. & Rivas Martínez, S.) 49–75 (Universidad de Alcalá de Henares, 1987).

54. Stow, D. A. V., Hernández-Molina, F. J., Alvarez Zarikian, C. A. & and, t. E. S. in *Proceedings of the IODP* (Proceedings IODP, **339**, Tokyo (Integrated Ocean Drilling Program Management International, Inc.), 2013).

55. Clopper, C. J. & Pearson, E. S. The use of confidence or fiducial limits illustrated in the case of the binomial. *Biometrika* **26**, 404–413 (1934).

56. Scherer, R. PropClS: Various Confidence Interval Methods for Proportions, R package version 0.3-0. https://CRAN.R-project.org/package=PropClS. (2018).

57. Juggins, S. Package "rijoa"-Analysis of Quaternary Science Data, The Comprehensive R Archive Network (2009).

58. Martrat, B. et al. Four cycles of recurring deep and surface water destabilizations on the Iberian margin. *Science* **317**, 502–507 (2007).

59. Villanueva, J., Pelejero, C. & Grimalt, J. O. Clean-up procedures for the unbiased estimation of C37-C39 alkenones sea surface temperatures and terrigenous n-alkane inputs in paleoceanography. *J. Chromatogr.* **757**, 145–151 (1997).

60. Prahl, F. G. & Wakeham, S. G. Calibration of unsaturation patterns in long-chain ketone compositions for palaeotemperature assessment. *Nature* **330**, 367–369 (1987).

61. Müller, P. J., Kirst, G., Ruhland, G., von Storch, I. & Rosell-Melé, A. Calibration of the alkenone paleotemperature index U37K' based on core-tops from the eastern South Atlantic and the global ocean (60°N–60°S). *Geochim. Cosmochim. Acta* **62**, 1757–1772 (1998).

62. Grimalt, J. O., Calvo, E. & Pelejero, C. Sea surfaces paleo-temperatures errors in U<sup>k</sup><sub>37</sub> estimation due to alkenone measurements near the limit of detection. *Paleoceanography* **16**, 226–232 (2001).

63. Polanco-Martínez, J., Medina-Elizalde, M., Sánchez Goñi, M. F. & Mudelsee, M. BINCOR: An R package for estimating the correlation between two unevenly spaced time series. *R. J.* **11**, 170–184 (2019).

64. Roche, D. M., Dumas, C., Bügelmayer, M., Charbit, S. & Ritz, C. Adding a dynamical cryosphere to iLOVECLIM (version 1.0): coupling with the GRISLI ice-sheet model. *Geosci. Model Dev.* **7**, 1377–1394 (2014).

65. Goosse, H. et al. Description of the Earth system model of intermediate complexity LOVECLIM version 1.2. *Geosci. Model Dev.* **3**, 603–633 (2010).

66. Opsteegh, J. D., Haarsma, R. J., Selten, F. M. & Kattenberg, A. ECBILT: a dynamic alternative to mixed boundary conditions in ocean models. *Tellus A* **50**, 348–367 (1998).

67. Goosse, H. & Fichefet, T. Importance of ice-ocean interactions for the global ocean circulation: a model study. *J. Geophys. Res.* **104**, 23337–23355 (1999).

68. Brovkin, V., Ganopolski, A., Claussen, M., Kubatzki, C. & Petoukhov, V. Modelling climate response to historical land cover change. *Glob. Ecol. Biogeogr.* **8**, 509–517, http://www.jstor.org/stable/2997870 (1999).

69. Roche, D. M. & Caley, T. δ18O water isotope in the iLOVECLIM model (version 1.0) – Part 2: Evaluation of model results against observed δ18O in water samples. *Geosci. Model Dev.* **6**, 1493–1504 (2013).

70. Berger, A. Long-term variations of daily insolation and Quaternary climatic changes. *J. Atmos. Sci.* **35**, 2362–2367 (1978).

71. Loulergue, L. et al. Orbital and millennial-scale features of atmospheric CH4 over the past 800,000[thinsp]years. *Nature* **453**, 383–386, http://www.nature.com/nature/journal/v453/n7193/Suppinfo/nature06950_S1.html (2008).

72. Ganopolski, A. & Calov, R. The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles. *Climate* **7**, 1415–1425 (2011).

73. Caley, T., Roche, D. M. & Renssen, H. Orbital Asian summer monsoon dynamics revealed using an isotope-enabled global climate model. *Nat. Commun.* **5**, 5371 (2014).
74. Bintanja, R., van de Wal, R. S. W. & Oerlemans, J. Modelled atmospheric temperatures and global sea levels over the past million years. *Nature* **437**, 125–128 (2005).

75. Kim, S.-T. & O'Neil, J. R. Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. *Geochim. Cosmochim. Acta* **61**, 3461–3475 (1997).

76. Coplen, T. B., Kendall, C. & Hopple, J. Comparison of stable isotope reference samples. *Nature* **302**, 236–238 (1983).

77. IPSL: LUDUS Framework, https://forge.ipsl.jussieu.fr/ludus, last access: 12 April 2023.

78. Porter, S. C. & Zhisheng, A. Correlation between climate events in the North Atlantic and China during the last glaciation. *Nature* **375**, 305–308 (1995).

79. Lofverstrom, M. A dynamic link between high-intensity precipitation events in southwestern North America and Europe at the Last Glacial Maximum. *Earth Planet. Sci. Lett.* **534**, 116081 (2020).

80. Railsback, L. B., Gibbard, P. L., Head, M. J., Voarintsoa, N. R. G. & Toucanne, S. An optimized scheme of lettered marine isotope substages for the last 1.0 million years, and the climatostratigraphic nature of isotope stages and substages. *Quat. Sci. Rev.* **111**, 94–106 (2015).

## Acknowledgements

This research used samples collected during the Expedition no. 339 "Mediterranean Outflow" of the Integrated Ocean Drilling Program (IODP). M.F.S.G. acknowledges funding from the GPR Human Past (University of Bordeaux). A.B. thanks Deutsche Forschungsgemeinschaft (DFG), project BA 3809/8. C.Z. acknowledges funding from IODP France and J.M.P.-M. from the Junta de Castilla y León and the European Regional Development Fund (Grant CLU-2019-03). T.R. acknowledges funding from FCT through projects Hydroshift (PTDC/CTA-CLI/4297/2021), WarmWorld (PTDC/CTA-GEO/29897/2017), UIDB/04326/2020, UIDP/04326/2020, LA/P/0101/2020 and EMSO-PT (POCI-01-0145-FEDER-022157). We thank Vincent Hanquiez for drawing Fig. 1 and Ludovic Devaux for pollen sample preparation.

## Author contributions

M.F.S.G. designed the research, performed the pollen analysis, and wrote the manuscript, T.E. and J.M.P.M. performed the model simulations and the statistical analyses, respectively. T.R. performed the

geochemical analyses (alkenone and %C₃₇:₄). T.E., J.M.P.M., C.Z., T.R. and A.B. participated in the discussion and final writing of the manuscript.

## 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-38337-4.

**Correspondence** and requests for materials should be addressed to María Fernanda Sánchez Goñi.

**Peer review information** *Nature Communications* thanks Manfred Mudelsee, Youbin Sun and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**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