Penultimate deglacial warming across the Mediterranean Sea revealed by clumped isotopes in foraminifera

The variability of seawater temperature through time is a critical measure of climate change, yet its reconstruction remains problematic in many regions. Mg/Ca and oxygen isotope (δ18OC) measurements in foraminiferal carbonate shells can be combined to reconstruct seawater temperature and δ18O (δ18OSW). The latter is a measure of changes in local hydrology (e.g., precipitation/evaporation, freshwater inputs) and global ice volume. But diagenetic processes may affect foraminiferal Mg/Ca. This restricts its potential in many places, including the Mediterranean Sea, a strategic region for deciphering global climate and sea-level changes. High alkalinity/salinity conditions especially bias Mg/Ca temperatures in the eastern Mediterranean (eMed). Here we advance the understanding of both western Mediterranean (wMed) and eMed hydrographic variability through the penultimate glacial termination (TII) and last interglacial, by applying the clumped isotope (Δ47) paleothermometer to planktic foraminifera with a novel data-processing approach. Results suggest that North Atlantic cooling during Heinrich stadial 11 (HS11) affected surface-water temperatures much more in the wMed (during winter/spring) than in the eMed (during summer). The method’s paired Δ47 and δ18OC data also portray δ18OSW. These records reveal a clear HS11 freshwater signal, which attenuated toward the eMed, and also that last interglacial surface warming in the eMed was strongly amplified by water-column stratification during the deposition of the organic-rich (sapropel) interval known as S5.


Supplementary information
Data smoothing. The  47 -based time series is a mixed signal between the ambient seawater 25 temperature recorded in the foraminiferal calcite, instrument shot-noise limit, measurement 26 error, sample preparation artefacts, and other unidentified sources of uncertainty 1 . In order to 27 reduce the high variability that these inflict on the  47 -data and highlight the main trends, we 28 apply a Moving Gaussian Window (MGW) filter. We use this smoothing method because it 29 allocates weights for averaging the  47 -replicates based on the size of the Gaussian window.

30
This allows one of the key advancements of our data processing approach, i.e., it does not rely 31 on any a-priori knowledge to obtain the final  47 -temperature records, as opposed to the 32 traditional binning method. That is, for the traditional  47 -small method it is necessary to define 33 the boundaries for averaging  47 - 18 O C replicates (binning) from neighbour samples.

34
Unfortunately, that requires a-priori information that is not always available. Although the 35 simultaneously generated and higher-temporal-resolution  18 O C record can be used to define the 36 binning boundaries, this relies on the assumption that temperature and  18 O C co-vary, which is 37 not always the case. In the wMed, for instance, temperature and  18 O C records across HS11 are 38 markedly asynchronous 2,3 ( Fig. 4b-c). The approach we propose in this study only requires 39 defining the size of the MGW, which here is referred to in terms of the 99.7 % coverage (±3) 40 of a Gaussian distribution in kyr. In general, the larger the size of the window, the smaller the 41 estimated error of the reconstructed parameter, and the smoother the record. In our case we use a 42 ~5 kyr window (±1 = 0.8 kyr) because it provides a reasonable compromise between associated 43 uncertainties and the length of the climate events we are investigating. Additionally, the 44 spans the same time interval of our data (Supplementary Fig. 5) shows that it satisfies the goal of 46 reducing the high frequencies in the record. That is, the 5 kyr MGW should minimize the noise 47 and highlight the main trend of the  47 -data.

49
Outlier test. The outlier test is based on the box-plot approach, which is an intuitive and broadly 50 recommended method for data analysis, e.g., in medical research 4,5 . It has the advantage of 51 detecting outliers regardless of the distribution of the sample and for this reason we implemented 52 it in the approach. The box-plot uses the 25 th (Q1), 50 th (Q2), and 75 th (Q3) percentiles of the 53 data to characterize the sample (see ref. 4 for a detailed explanation). Data points ( 47 -replicates) 54 are considered outliers if they fall outside the Q3+W*(Q3-Q1) and Q1-W*(Q3-Q1) range. The 55 standard approach of Tukey 6 uses W = 1.5. In the case of Gaussian distributed data, that value 56 would identify  47 -replicates outside a ±2.698 range (99.3 % coverage) as outliers. We 57 consider this range too large to detect outliers in our dataset and therefore we also tried a 58 conservative definition that would consider outliers  47 -replicates falling outside a ±2 (~95.5 % 59 coverage) range in a normal distribution (W = 0.9826). Supplementary Figure 6 shows that the 60 final ODP975 record using W = 1.5 is not much different from those obtained using W = 0.9826 61 (main figures). The records only diverge where the density of the  47 -replicates is low at, e.g.

62
122-115 kyr in LC21. Yet, because (i) we focus our discussion on the interval between 140-122 63 kyr and (ii) the MGW filter is sensitive to outlying data, we favoured the conservative test, even 64 if this implies that a few "potentially good"  47 -replicates might be wrongly removed from the 65 dataset. 66 67 record, there is a slight effect on how smooth the confidence intervals appear. Fewer simulations 71 result in a slightly less smoothed record ( Supplementary Fig. 7). We selected N = 5,000 for our 72 study because it offers a good compromise between the smoothness of the confidence intervals 73 and the time that the simulations take.   (c) are homoscedastic (i.e., have constant variance) 8 . We checked the efficiency of the non-116 traditional  47 data-analysis approach by analysing the residuals between the measured  47 -117 (Supplementary Fig. 9-11). 119

(a)
We checked that the residuals were normally distributed by applying the so-called 120 'Lilliefors Test' 9 at the 0.05 significance level. The histogram and Quantile-Quantile plot of the 121 residuals are also shown, along with the residuals (Supplementary Fig. 9-11a-c), to visually 122 corroborate the results of the Lilliefors test. 123 We tested that the residuals were not auto-correlated by estimating their persistence time 124 () using a first-order autoregressive model (AR (1) We calculated S(a) for the  47 residuals 5,000 times using the unique ages (i.e., no replicated 132 ages), with the option of randomly drawing a  47 residuals value when a sample has  47 -133 replicates >1. We compare the results to that of a system that has no memory in it. For this, we 134 randomised (i.e., induced = 0yrs) each of the 5,000  47 residual realizations and calculated 135 S(a). Supplementary Figures 9-11d shows the 97.5 th , 50 th , 2.5 th percentiles for the 5,000 136 instances that S(a) was calculated before (blue) and after (red) residuals were randomised. For 137 all the cores the two systems have similar behaviours implying that , and in turn the 138 autocorrelation, in the residuals is very small. This is corroborated when looking at the kernel 139 density estimate of Supplementary Fig. 9-11e), a rough estimate of the memory in the 140 residuals, which is always <1000years.

(c)
The heteroscedasticity was evaluated by checking that the squared residuals 11 were not 142 autocorrelated, following the same procedure explained in (b  variation of S(a) -the least-squares estimator that minimizes a (a = exp(-1/-when a varies 293 from 0.8 to 1 and we use the residuals (blue) or the randomised residuals (red). e, kernel density 294 estimate of f-g are like d-e, respectively, but using the squared-residuals. 295 296 during the onset of S5 with respect to late-Holocene. Summer mixed layer value at present is 298 ~50m and it is expected to shoal during sapropel periods 15 . Insolation was obtained from ref. 16 299 using AnalySeries software. We averaged summer (June 21 st ) and winter insolation (December 300 21 st ) at 34°N during the late-Holocene (0-2 ka) and onset of the S5 (127-129 ka). We calculated 301 the resultant mixed layer heating for three scenarios: a) using late-Holocene insolation and 302 present mixed layer; b) onset of S5 insolation and present mixed layer; c) onset of S5 insolation 303 and shallower mixed layer (35 m). Attenuation (a) in all cases was varied between 0.4-0.6 (20 % 304 of that in Athens).