MD01-2444 core location and study interval

Core MD01-2444 (37°33.68′N; 10°08.53′W; 2637 m water depth; 27.45 m in length) was recovered near the location of core MD95-2042 (ref. 20) (37°48′N, 10°10′W; 3146 m water depth), in 2001, using the CALYPSO Giant Piston corer aboard the French research vessel Marion Dufresne II. This study focuses on the 22.50–19.50 m section in MD01-2444, (~140–108 ka), with a mean sedimentation rate of ~10 cm kyr−1. Samples for pollen and palaeoceanographic analyses were obtained from the same levels, typically at 1-cm intervals.

MD01-2444 stable isotopes

Measurements were carried out in the Godwin Laboratory for Palaeoclimate Research (University of Cambridge) using a VG PRISM, VG SIRA, or Thermo MAT253 mass spectrometer. For PRISM/SIRA measurements, samples were dried in an oven at 50 °C overnight prior to sealing the vials with a septum and screw cap. The samples were analysed using a Micromass Multicarb Sample Preparation System attached to a VG SIRA or PRISM Mass Spectrometer. Each run of 30 samples is accompanied by 10 reference carbonates and 2 control samples. For the Thermo MAT253 measurements, foraminifera are transferred into sample vials, crushed and then dried in an oven at 50 °C. The vials were loaded into the carousel and analysed using a Thermo Kiel device attached to a Thermo MAT253 Mass Spectrometer in dual inlet mode. The preparation system operates automatically analysing samples in sequence. One hundred percent orthophosphoric acid is dropped onto the evacuated vial and reacts with the calcium carbonate sample. The evolved carbon dioxide is cryogenically dried and then admitted to the dual inlet mass spectrometer for isotopic analysis by comparison with a reference gas. Each run of 30 samples was accompanied by 10 reference carbonates and 2 control samples. The results are reported with reference to the international standard Vienna PeeDee Belemnite (VPDB) and the precision is better than ±0.06‰ for 12C/13C and ±0.08‰ for 16O/18O.

For the planktonic-isotope record about 15–30 specimens of Globigerina bulloides were selected from the 250 to 300 μm size fraction. Replicates were carried out for a selection of intervals using ~20–30 specimens from the 300 to 350 μm size fraction, confirming no apparent offset or differences in down-core variability in the different size fractions. Where replicate analyses were performed, measurements have been averaged.

For the benthic δ18O record, several different species (Cibicidoides sp., Cibicidoides wuellerstorfi, Uvigerina peregrina, Hoeglundina elegans and Globobulimina affinis) from the >212 μm size fraction were analysed as no single species was present in sufficient numbers to generate a continuous 1-cm sampling resolution record. Where possible two or three separate analyses of different species were made in each sample; a correction factor was applied (Cibicidoides sp.: +0.57; U. peregrina and similar specimens: 0.0; G. affinis: −0.3; C. wuellerstorfi, +0.64; H. elegans: −0.67). These adjustments are optimized for this particular core63 in accordance with the long-standing convention by which U. peregrina is assumed to deposit oxygen in isotopic equilibrium. The average of all the corrected values at each level is shown in the figures. The epibenthic δ13C record was derived from measurements on monospecific samples (1–4 individuals) of C. wuellerstorfi selected from the >212 μm size fraction wherever they were present. These new measurements were combined with those previously reported by ref 25.

Results are reported in Supplementary Data 1.

MD01-2444 Mg/Ca measurements

Trace element measurements were conducted at the University of Cambridge via inductively coupled atomic emission spectrometry, using a Varian Vista machine64. New measurements have been combined with those previously reported25, with replicate analyses being averaged to yield a continuous record. All Mg/Ca analyses have been screened for dissolution artefacts by reference to shell weights and for contamination by reference to measured iron, aluminium and manganese concentrations. Only analyses with Fe/Ca and Mn/Ca ratios <0.5 mmol mol−1 were retained. Mg/Ca ratios in G. bulloides have been calibrated65,66 to temperatures, which yield absolute temperatures that correspond to the modern habitat ranges of these species. It must be stressed that any proxy record derived from a biological ‘proxy-carrier’ necessarily reflects a record of ‘habitat change’, and that this will reflect local climate change to the extent that the organism’s habitat has been coupled with the ‘synoptic’ local/regional climate. The ‘shallow water’ (here referred to as ‘surface water’, as distinct from ‘deep water’) temperature records derived here from G. bulloides Mg/Ca ratios thus represent a record of the impact of local climate change on the habitats of this planktonic foraminifer. At present, G. bulloides thrives near the surface under eutrophic (‘bloom’ or upwelling) conditions during spring and summer, or at frontal upwelling zones67.

The absolute uncertainty in the calculated temperature estimates (due to statistical calibration uncertainty and analytical reproducibility) is estimated at ~±0.7 °C. However, the ‘noise’ in the temperature time series, which can be approximated by the average standard deviation of paired adjacent measurements (which essentially provides an estimate of the degree of autocorrelation in the record), is closer to ±0.8 °C or about ±0.2 mmol mol−1. This estimate is obtained by assuming that, as their depth offset tends to zero, adjacent measurements will tend to represent replicates of a single mean value.

Estimates of seawater δ18O (δ18O sw ) have been derived by combining calcite δ18O measurements with Mg/Ca-derived calcification temperature estimates, both performed on the same planktonic species. We use the palaeotemperature equation68,69, whereby:

$${\mathrm{T}}_{{\mathrm{Mg}}/{\mathrm{Ca}}} = 16.9 - 4.38\left( {\delta _{\mathrm{c}} - \delta _{\mathrm{sw}}} \right) + 0.1\left( {\delta _{\mathrm{c}} - \delta _{\mathrm{sw}}} \right)^2.$$ ( 1)

In the above equation \(\delta _{\mathrm{c}}\) and \(\delta _{\mathrm{sw}}\) are the calcite and seawater δ18O, respectively, and are both referenced to the same primary standard (i.e. \(\delta _{\mathrm{c}}\) corrected by +0.2‰ to convert from VPDB to Vienna Standard Mean Ocean Water (VSMOW)).

Results are reported in Supplementary Data 1.

MD01-2444 pollen analysis

Subsamples of 3–7 g of sediment were prepared for pollen analysis using the standard hot acid digestion technique. Fine sieving, through a mesh of 10 μm or less, was not used as it has been found to result in a loss of pollen, particularly Gramineae. Residues were mounted in silicone oil for microscopic analysis at magnifications of 400, 630 and 1000 times on a Leica DM2000 light microscope. Nomenclature follows ref 70. Abundances are expressed as percentages of the main sum, which includes all pollen except Pinus, Pteridophyte spores and aquatics. Pinus is conventionally excluded from the main sum, as it is strongly over-represented in marine sediments because of its extensive dispersal ability and buoyancy71. Following the convention for marine pollen analyses, a minimum of 100 pollen grains, excluding Pinus, spores and aquatics, was counted in each sample. Pollen totals (including Pinus) ranged from 123 to 540 grains. Shown here are the summary pollen curves of Mediterranean sclerophylls (Olea, Pistacia, Phillyrea and evergreen Quercus) and temperate trees, which includes Mediterranean sclerophylls and Eurosiberian taxa (the latter includes deciduous trees/shrubs and Abies, and excludes Juniperus, Hippophae, Salix, Betula). Results are reported in Supplementary Data 1.

Pollen studies from continental shelf sequences suggest that palynomorph transport to these areas is controlled primarily by fluvial and secondarily by aeolian processes72. Studies on modern pollen deposition in fluvial systems, considering the transfer of pollen from vegetation to the channel, transport in the channel and deposition in coastal waters, indicate the rapid incorporation of pollen to marine sediments72. At the Portuguese Margin margin, aeolian pollen transport is limited by the direction of the prevailing offshore winds and pollen is mainly transported to the abyssal site by the sediments carried by the Tagus River73,74. Comparison of modern marine and terrestrial samples along western Iberia has shown that the marine pollen assemblages provide an integrated picture of the regional vegetation of the adjacent continent74. In MD01-2444, the coherence of the pollen and marine isotope signals in MIS 3 and MIS 6 (ref. 63) argue against reworking of non-contemporaneous material.

MD01-2444 XRF analysis

Archive halves of sections from Core MD01-2444 were analysed using an Avaatech XRF core scanner at the University of Cambridge55. The core surface was carefully scraped cleaned and covered with a 4-mm thin SPEXCertiPrep Ultralene foil to avoid contamination and minimize desiccation. XRF data were collected every 2.5 mm. The length and width of the irradiated surface was 2.5 and 12 mm, respectively, with a count time of 40 s. Results are reported in Supplementary Data 1.

Corchia Cave analytical methods

The speleothem time series used in this study is derived from four Corchia Cave stalagmites (CC1, CC5, CC7 and CC28) collected from a single chamber (44°02′N, 10°18′E) between 1999 and 2005. The interval reported spans 140 to 108 ka. The cave setting, the collection details and the general characteristics of each stalagmite have been described in detail elsewhere21,32,33,75,76. Except for two small sections of CC5, each segment under study was microsampled at 100-μm resolution using a Taig CNC micromilling lathe. A ~41-mm interval of rapid growth in CC5 during the early LIG was sampled at intervals of between 0.2 and 2.0 mm, while another CC5 section (~38 mm length) at ~121 ka was sampled at an average interval of 0.12 mm.

The δ18O and δ13C measurements were conducted on two identical continuous-flow mass spectrometers at the Scottish Universities Environment Research Centre (SUERC) at East Kilbride, UK (an Analytical Precision AP2003) and the School of Earth Sciences, University of Newcastle, Australia (a GV Instruments GV2003). In both cases, subsamples of 0.8 ± 0.1 mg were digested in 105% orthophosphoric acid at 70 °C and measurements made on the evolved CO 2 . Sample results were converted to the VPDB scale using the known values of Cararra Marble house standards (NEW1 at Newcastle and MAB2C at SUERC) previously calibrated to the VPDB scale using NBS19 and NBS18. Each batch of 132 vials included between 22 and 27 house standards ~evenly spaced through the sequence. At Newcastle, scale corrections (small, <0.10‰) for δ18O were applied where necessary based on values obtained for two measurements each of NBS19 and NBS18. At SUERC, scale corrections were unnecessary, based on results of regular measurement of samples against both NBS19 and NBS18 using dual-inlet mass spectrometry (VG Prism II). The repeatability precision (1σ) on the AP/GV2003 machines throughout the analyses was ≤0.10‰ for δ18O and ≤0.05 for δ13C. Results are reported in Supplementary Data 2.

U–Th age data for each speleothem are shown in Supplementary Data 3. The U–Th dating was conducted either on samples drilled using a 1-mm drill bit or on composites of leftover powders originally microsampled for stable isotope analysis. The latter samples represented between 0.1 and 0.5 mm of speleothem distance (i.e. between one and five 100-μm-thick samples), depending on the amount of powder remaining. Conservative estimates were made of the 100% sampling uncertainty (Supplementary Data 3). Chemical separation of U and Th, spike procedures and the method of isotopic measurements followed the protocol of refs. 77,78. All age calculations use the latest U and Th decay constants79.

A Bayesian Monte Carlo method ‘finite growth-rate age model’ was used to calculate depth-age models for each stalagmite and to generate 95% age-uncertainty estimates for each stable isotope sampling depth position21,32,80. The depth-age plots and the modelled 2σ uncertainties for each speleothem are shown in Supplementary Fig. 1. Two U–Th ages were immediately rejected from CC5 because of their outlying position in depth-age space (red symbols in Supplementary Fig. 1; red text in Supplementary Data 3). Based on the updated U–Th age calculations, the time periods covered by each speleothem considered in this paper are:

CC1: 141.0 ± 1.6 to 128.5 ± 1.1 ka.

CC5: 141.0 ± 1.5 to 126.0 ± 1.0 ka, and 124.4 ± 0.99 to 107.0 ± 0.7 ka.

CC7: 128.6 ± 1.4 to 126.4 ± 1.5 ka, and 125.6 ± 1.3 to 121.3 ± 1.6 ka.

CC28: 118.4 ± 1.3 to 107.0 ± 1.0 ka.

Corchia Cave composite (‘stacked’) time series

All but ~4 kyr of the 32-kyr interval is replicated and the isotopic patterns between coeval stalagmite sections display good consistency (not shown here). Accordingly, we constructed a single composite stable isotope record that attempts to maximize the level of detail by utilizing sections with the highest sampling resolution where possible, bearing in mind the challenges of splicing individual records into a single time series.

We used the stalagmite possessing the most complete record, CC5, as the tuning target. CC5 spans all but ~1.6 kyr of the 140–108 ka period. To produce the master time series, each individual speleothem depth series was first smoothed using a 5-point window. The stable isotopic data for CC1, CC7 and CC28 were then tuned to the CC5 depth scale using tie points common to overlapping growth segments. Tuning points were identified by eye using both δ18O and δ13C profiles. In difficult-to-match sections, the δ18O was used as the main tool for pattern matching, as this ratio is more robust than δ13C, which is more susceptible to local (including ‘in-cave’) influences.

For the long period of overlap between CC5 and CC28, we used the stable isotope data of the latter due to its superior time resolution. Data from this speleothem commence from a point as close as possible to its base where a clear isotopic matching with CC5 is evident. For the periods of overlap between CC7 and CC5, CC7 stable isotope data is only used across the hiatus in CC5 and where the time resolution is superior and there is confidence in the splice with CC5.

To develop the stack depth scale, the depth scale of CC5 was reset to an arbitrary zero corresponding to the termination of the longest period of continuous growth in this stalagmite (124.4–87.8 ka32), which occurs at an absolute depth position of 14.4 mm from the tip of CC5. To accommodate the hiatus in CC5 growth, a depth increment equivalent to the segment of the CC7 section was inserted into the stack depth scale. The pre-hiatus depth of CC5 was adjusted accordingly.

The results of the cross-tuning are shown in Fig. 2a, b and reveal good replication between the overlapping portions of the speleothems for both δ18O and δ13C. For CC7, we added 0.3‰ to the δ18O data to bring its profile into a more convincing alignment with CC5, although it appears this adjustment is somewhat inconsistent with the older segment. The need to apply an offset, and the inconsistency of such an offset over several thousand years, is not entirely surprising: the cave chamber lies at great depth below the ground surface (at least 400 m), and this plus the recharge elevation range of the structurally complex carbonate geology that comprises the aquifer31 would conspire to permit recharge waters from markedly different altitudes (and thus with different isotopic composition31) to reach the chamber at different times, according to the evolution of the flow path for individual drip points. The final stable isotope series on the stack depth scale is shown in Fig. 2b and reported in Supplementary Data 2.

The advantage of a single-stacked-record approach is that it allows the U–Th age data from all speleothems to be combined into a single depth-age model, thereby reducing overall age-model uncertainties for the entire 140–108 ka interval21. Accordingly, we transposed the original U–Th age-depth positions from the individual speleothems onto the stack depth scale and derived a new stack depth-age model using the Monte Carlo technique described previously (Fig. 2c, d). In implementing this model, we increased the sample-depth error (by ±1 mm) of all non-CC5 ages to account for uncertainty in the tuning procedure. The composite depth-age model reveals a further age outlier not detected in the original single-speleothem depth-age model, giving a total of three rejected outliers from a total of 90 U–Th age determinations. These outliers are highlighted by red symbols in Fig. 2c and have been excluded from the depth-age model.

Finally, a useful test of the robustness of the chronology of a stacked record is how closely the model stack age for each isotope data point agrees with its original model age. Supplementary Fig. 2 shows the results of this comparison. Only five very brief intervals display stack model ages outside the 95% uncertainties of the original model age, equivalent to ~2% of the entire record. In each case, the model-age difference is small (within 99% model-age uncertainties) and in no instance do these intervals correspond to the inferred aridity events identified in Fig. 3.

MD01-2444 age model

The age model of the LIG section of MD01-2444 is based on aligning its temperate tree pollen curve to the δ18O speleothem record from Corchia Cave (Supplementary Fig. 3), on the premise that, under the influence of mid-latitude westerlies and cyclogenesis over the Mediterranean, rainfall amount in southern Europe exerts a dominant control over both the composition of vegetation and isotopic signature of speleothems21,31,34. Supplementary Data 4 provide the alignment tie points, ages and associated uncertainties. It is important to emphasize that the comparisons of proxy records from the same samples in the MD01-2444 sediment sequence are independent of the choice of age model.

Detrending the MD01-2444 and Corchia Cave records

To remove the low-frequency component from the temperate tree pollen, planktonic foraminifera and speleothem-isotope records, a 6-kyr Gaussian interpolation was subtracted from the data sets on their original time step. The residuals were then smoothed using a 600-year window Gaussian interpolation, significantly reducing the amplitude of periodicities below (shorter than) 400 years (Fig. 4).

Comparison of intra-interglacial aridity events and ocean changes

To identify interglacial samples with negative temperate tree pollen and positive δ18O planktonic values, the detrended data on their original time step (common to both data sets) were examined (Fig. 4). To avoid selecting minor variations that may be analytical noise or artefacts of the detrending method, temperate tree pollen and δ18O plantonic values in the same sample had to reach or exceed thresholds of –1.13% and +0.04 ‰, respectively. These thresholds are equivalent to one-third the mean absolute deviation of the respective data sets.

Alignment of ODP984 to MD01-2444

The LIG record9 of ODP984 south of Iceland (61.25°N, 24.04°W; 1648 m water depth) is here aligned to that of MD01-2444 (Supplementary Fig. 4) on the basis that iceberg discharges, expansion of cold surface-water and attendant SST changes are quasi-synchronous (within the sampling resolution) between the two sites. Large changes in XRF Zr/Sr and IRD and also SST and percentages of Neogloboquadrina pachyderma, sinistral coiling (s.), considered a cold indicator species9, provide an unambiguous stratigraphic correlation of HS11.1–HS11.3 and C24. Smaller decreases in SSTs also constrain the correlation of C28 and C25. In between these events, the alignment of the two records is based on correlations of their respective δ18O planktonic records. Supplementary Data 4 provide the alignment tie points, ages and propagated uncertainties.

Alignment of MD03-2664 to ODP984

The LIG record10,38,39 of MD03-2664 from the Eirik Drift, south of Greenland (57°26.34′N, 48°36.35′W; 3442 m water depth) is here aligned to that of ODP984 (ref. 9) (Supplementary Fig. 5), on the basis that iceberg discharges, expansion of cold surface water and attendant SST changes are quasi-synchronous (within the sampling resolution) between the two sites. Large changes in IRD, SSTs and percentages of N. pachyderma (s.), provide an unambiguous stratigraphic correlation of HS11.1, HS11.2 and C25 and C24. Although the deglacial section of MD03-2664 is condensed, the end of HS11 and start of the LIG is clearly recorded in the IRD and SST records, which then allows the correlation of C28. The alignment of the intra-interglacial sections is based on correlations of changes in Mg/Ca SSTs in MD03-2664 and in the percentages of subpolar N. incompta percentages, representing a warm indicator species in ODP984 (ref. 9). Supplementary Data 4 provide the alignment tie points, ages and propagated uncertainties.

Numerical experiments performed with LOVECLIM

North Atlantic meltwater experiments were performed with the Earth System model LOVECLIM44, which comprises a spectral T21 quasi-geostrophic atmospheric model and an ocean general circulation model coupled to a thermodynamic-dynamic sea-ice model, with a horizontal resolution of 3° longitude by 3° latitude and 20 vertical levels. The experiments were performed under varying LIG boundary conditions, that is, appropriate orbital parameters81, Northern Hemisphere ice-sheet extent and albedo82, atmospheric CO 2 content83 and constant CH 4 (700 ppb) and N 2 O (310 ppb) values. Two experiments were performed with meltwater added in the northern North Atlantic, to the area south of Greenland (55.5°N–66.3°N, 58.5°W–31.5°W) for 400 years. For the first experiment, 0.05 Sv (103 m3 s−1) was added between 126.8 and 126.4 ka (FWF005), and for the second one, 0.03 Sv was added between 124 and 123.6 ka (FWF003). Both experiments run for another 400 years without any forcing and were compared to an LIG control experiment. Following the freshwater addition, the AMOC weakens from ~23 to 15 Sv (~35%) at 126.7 ka and from 24 to 19 Sv (~20%) at 123.9 ka. We note that LOVECLIM’s North Atlantic and European climate response to an AMOC weakening is within the range of coupled climate models84.

Numerical experiments performed with CESM

To further assess the link between changes in North Atlantic SST and precipitation over Europe, we perform additional experiments with an independent model, the NCAR CESM. The atmospheric component of the CESM, namely the Community Atmospheric Model version 4 (ref. 85) is forced with monthly SST and sea-ice conditions from LOVECLIM. The CAM4 configuration consists of a spatial resolution of 1.9° latitude by 2.5° longitude and 26 vertical levels in a hybrid sigma-pressure coordinate. Four experiments are performed with an atmospheric CO 2 concentration fixed at 284.7 ppmv and integrated as follows: (i and ii) two control run experiments forced with a repeating 12-month climatology of SST and sea-ice concentration from the LOVECLIM LIG control run at 126.8 and 124 ka; (iii) an experiment forced with a climatology of both North Atlantic SST and sea-ice concentration derived from the LOVECLIM FWF005 simulation when the AMOC is reduced by 35% (Supplementary Fig. 6b, c); and (iv) an experiment forced with a climatology of both North Atlantic SST and sea-ice concentration derived from the LOVECLIM FWF003 simulation when the AMOC is reduced by 20% (Supplementary Fig. 6e, f).

Outside the forcing domain, a 10° latitude/longitude band was used to linearly damp the North Atlantic SST and sea-ice conditions to the control run climatology. Unless otherwise specified, SST and sea-ice concentration were the same as for the control run. The control runs were integrated for 270 and 100 years and the perturbation experiments were integrated for 90 and 80 years, respectively.

Analysis of the variance of LIG and Holocene records