Main

Precise records of climate variability from the past provide insights into the mechanisms of climate dynamics, which in turn helps to better assess ongoing and projected climate trends. The Last Interglacial (LIG; ~129,000–116,000 years ago) is the most recent warmer-than-present period in Earth history, when global sea surface temperatures (SSTs) were ~0.5–1 °C warmer than today1,2,3 and sea level 1.2 (ref. 4)–9 m (ref. 5) above the present-day level. Consequently, the LIG serves as a key interval for reconstructing climate variations that may be partially representative of the climate of the near future6,7. While the warmer mean SSTs during the LIG are reasonably well constrained, little is known about the short-term dynamics during this period, in particular the decadal to centennial rates of SST change. As future rapid changes will result in increased adaptational pressure to marine ecosystems and human societies8, data-based observations of the dynamic behaviour of SST changes from the most recent warmer-than-present analogue to future warming will provide critical information for assessing the oceans’ response. So far, reconstruction of continuous short-term climate variations from the LIG has been restricted by the limited availability of temporally well-constrained continuous archives that preserve undisturbed signals of climate variability in sufficient resolution and by the coarse temporal resolution provided by conventional analyses.

The Mediterranean Sea is one of the most sensitive regions to climate change, where the anthropogenic greenhouse gas forcing causes SST to increase at a faster rate than the global average9,10,11. Understanding the decadal climate variability from the LIG provides an opportunity to relate current and projected climatic trends to natural climate variability from a warmer world. The undisturbed sapropels deposited during the LIG in the eastern Mediterranean present an excellent opportunity for examining SST variability during this period at ultra-high temporal resolution. Sapropels are organic-rich sediment layers formed under anoxic conditions that developed in the eastern Mediterranean Basin due to reduced ventilation induced by monsoon intensification and subsequent freshwater flooding of the basin12. The LIG sapropel, known as sapropel S5, is typically only a few decimetres thick, precluding the extraction of a decadal climate record. Here we examined an exceptionally thick S5 sapropel layer recovered in the sediment core M51/3-SL104 (34.8247° N, 27.2925° E; water depth 2,155 m) from the Pliny Trench region of the eastern Mediterranean (Fig. 1a). The layer comprises a ~87-cm-thick section, where micrometre-scale annually formed laminations (varves) are preserved13,14 (see Supplementary Information for the detailed stratigraphy). Deposition of the varves provides a rare opportunity to generate a precise relative age model by counting annual layers in the varved part of the sapropel. Consequently, this section has the potential to provide a long and continuous subdecadal resolution climate record with a precise age model on the relative timescale.

Fig. 1: Location and stratigraphy of sapropel S5 and an exemplary MSI-based determination of SST via analysis of alkenones.
figure 1

a, Map of the Mediterranean Sea; the map was generated by The Generic Mapping Tools software. Core M51/3-SL104 is indicated with a black circle. Discussed cores from the eastern Mediterranean with preserved sapropel S5 and alkenone-based SST reconstruction (LC2121, ODP Site 97128, ODP Site 96728, KS 20528) are indicated with coloured circles. b, Sapropel S5 from core M51/3-SL104; sediment depth is presented in centimetres below seafloor (cmbsf). c–e, An illustrative example of MSI on a ~5 cm embedded subsample. c, Image of a slice on which MSI was performed, red rectangle indicates the area analysed. d, Spatial distribution of \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\)values obtained by MSI. Note that the scale is presented as x- and y-position spot, where each spot accounts for 200 µm. e, \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) and calculated SST time series for the selected sample; data points represent SST derived from \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) values calculated from the sum of alkenone intensities in one horizontal layer (only layers with >10 spots with successful detection of alkenones were considered). Depth is presented as y-position spots where each spot is 200 µm.

Full size image

Here we reconstruct molecular proxy-based SST at annual to subdecadal resolution by applying mass spectrometry imaging (MSI). MSI on sediments has been recently introduced to palaeoclimatology; this technique allows for the detection and visualization of biomarkers on intact sediment core sections at micrometre-scale resolution15,16,17 (Fig. 1c–e). Among other compounds, MSI allows for the measurement of long-chain alkenone abundances18, which provides access to the well-established \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) palaeotemperature proxy19,20, one of the most prominent approaches in SST reconstruction. We present alkenone-based SST at a spatial resolution of 200 µm (Fig. 2 and Extended Data Fig. 1), resulting in an annually to subdecadally resolved continuous millennial record indicative of climate variability during a past period with warmer-than-present climate. Moreover, we supplement this data with 200-µm-resolution data of isorenieratene, an aromatic carotenoid biomarker of anaerobic, phototrophic green sulfur bacteria that indicates photic zone euxinia21 (Extended Data Fig. 2), which provides important information on the stratification of the water column. Our goal was to examine the variability of SST and to constrain the rate of SST change on decadal to centennial timescales in the warmer-than-present LIG in the eastern Mediterranean.

Fig. 2: Comparison of alkenone-based SST, sediment lightness13 and Fe/Al ratio13 from core M51/3-SL104 with discussed proxies from nearby core LC2121.
figure 2

a–f, Ultra-high-resolution alkenone-based SST record (blue line, the black line is a 50-point running average; a), pixel-scale sediment lightness13 (note axis break; b) and Fe/Al ratio13 indicative of redox conditions (note axis break; c) from core M51/3-SL104 on the age scale26, compared with δ18O from planktic foraminifera (G. ruber; d), alkenone-based SST (e) and isorenieratene (f) from core LC2121 plotted on the depth scale. We refer to Extended Data Fig. 1, which additionally shows uncertainties for the SST record. Elemental ratios sensitive to redox conditions and primary productivity (see also Extended Data Fig. 6) support more reducing conditions and higher primary productivity in darker layers. Please note that M51/3-SL104 sapropel plotted on the age scale is only tentatively comparable to the LC21 depth scale as constant sedimentation rates are assumed for LC21. Grey area indicates the period of sapropel deposition, where darker grey indicates the non-varved part and lighter grey the varved part of the M51/3-SL104 sapropel. VPDB, Vienna PeeDee Belemnite.

Full size image

Eastern Mediterranean SST during the LIG

Sapropel formation is forced by the Northern Hemisphere summer insolation maximum22, and consequently sapropel S5 was deposited during the warmest part of the LIG in the Northern Hemisphere23. Precisely establishing the duration of the sapropel S5 deposition is challenging because the determination of the absolute age of the sapropel is infeasible, but it is considered to have been deposited over a period of ~4,800–6,800 years24,25,26. We applied an age model based on microscopic annual layer counting for the varved part, which covers 61% of sapropel duration (~2,900 years)26, while the non-varved part is considered to cover ~1,900 years (see Methods for further details). The absolute chronology is less certain and based on the analogy to the youngest S1 sapropel, where a lag of 3,000 years exists between the midpoint and the correlated maxima in insolation22,27 (see Methods). The down-core resolution of our dataset (200 µm) translates into a temporal resolution of ~4 years for the non-varved part of the sapropel (~126.3–124.4 thousand years ago; ka) and ~1 year for the varved part of the sapropel (124.4–121.5 ka). The detection of alkenones was unsuccessful at 6 intervals of 10–25 years duration between 125,685–125,658, 124,722–124,699 years bp in the non-varved part and between 122,546–122,521, 122,500–122,515, 122,335–122,223 and 122,267–122,253 years bp in the varved part.

The reconstructed SST range of 15.2 ± 1.9 to 23.8 ± 1.9 °C (for uncertainty determination see Methods) is consistent with the results from lower-resolution alkenone-based SST reconstructions from this study (Extended Data Figs. 3 and 4), and previous S5 sapropel records from the eastern Mediterranean21,28 (Fig. 2 and Extended Data Fig. 5). The temporal resolution provided by MSI, however, allows for an in-depth examination of the mode and amplitude of Mediterranean SST variability and its drivers during the studied part of the LIG. In terms of multicentennial trends, our record is consistent with previous low-resolution alkenone SST records from the eastern Mediterranean (Fig. 2 and Extended Data Fig. 5), suggesting a remarkable SST increase with the onset of sapropel deposition. Our record indicates that this warming was more abrupt than suggested by previous studies: the most prominent SST increase by ~3 °C to temperatures of ~23 ± 1.5 °C unfolded in just ~100 years between 125.67 and 125.57 ka (Fig. 2) and subsequently SST remained relatively high for the rest of sapropel S5 deposition. At site LC2121, SST increase during this early part of sapropel deposition is congruent with a decreasing trend in planktic δ18O, indicating an increase of freshwater flooding into the basin. This link between progressive freshwater flooding into the basin21,29 and the increase of SSTs can be explained by a reduction of deep mixing of the heat gained at the surface into deeper waters due to increased stratification. Accordingly, SST in the early part of the sapropel deposition was forced by the change in stratification of the water column, which allowed for large SST variations controlled by changes in the depth of the mixed layer.

The abrupt SST increase ~125.67 ka was followed by a pronounced increase in the concentrations of isorenieratene (Fig. 2 and Extended Data Fig. 2), reflecting the development of euxinic conditions in the photic zone of the stratified water column21. After the establishment of euxinic conditions until the onset of the varve formations, reconstructed SSTs were consistently higher than present-day SSTs within the calibration uncertainty level, with the exception of a few decades. We also observe prominent centennial SST variability in this part of the record, which was not revealed in low-resolution SST records21,28.

With the onset of varve formation, sedimentation rates increase and the resulting shift from ~4-year to ~1-year resolution of our record results in an apparent increase in SST variability (Fig. 2). A similar increase in variability can be observed in the other micrometre- to millimetre-scale resolved proxies from the M51/3-SL104 sapropel (Fig. 2 and Extended Data Fig. 6) and should not be considered as an indicator of the abrupt change in naturally occurring SST variability. However, the varve formation is indicative of major reorganization in the eastern Mediterranean oceanographic conditions. The deposition of the diatom-opal-rich laminae responsible for varve formation in M51/3-SL104 S5 sapropel is best explained by the mass sinking of diatom mats due to the onset of autumn/early winter mixing of the water column in the Mediterranean, indicating the breakdown of stratification in colder months14. By contrast, stratification probably persisted throughout the entire year during the early stage of sapropel deposition. Moreover, the onset of varve formation at 124.5 ka coincides with a decrease in isorenieratene concentrations in Mediterranean S5 sapropels (Fig. 2 and Extended Data Fig. 2). A previous study29 demonstrated that this decrease of isorenieratene is coupled with the deepening of the summer mixed layer, which reached an approximately modern level of ~30 m below surface30. This was further associated with a decrease in total organic carbon and a return of salinity to pre-sapropel conditions31, indicating a lowered input of freshwater. The presence of isorenieratene and high total organic carbon in varved sediment still indicate euxinic conditions in the photic zone and higher-than-modern marine productivity, that is, conditions that would affect carbon export but not the mechanism of SST variability. On the other hand, the return of salinity to pre-sapropel values31, the similar to modern depth of the summer mixed layer29 and the breakdown of stratification in colder months14 (present-day mixed layer depths for the studied location during winter is not deeper than 130 m30) suggest that the mechanism controlling SST dynamics during deposition of the varved interval was not fundamentally different to today. Consequently, in particular the varved part of the sapropel provides detailed information on the regional SST dynamics and their relation to regional atmospheric circulation during the LIG, and is thus a useful reference point for comparison with current and projected SST change in that area.

Decadal SST variability during the LIG

The reconstructed high-resolution SST record allows for an examination of regional decadal to multidecadal scale climate variability during the LIG. The varved part of the sapropel S5 is particularly well-suited to evaluate potential periodicities, as its age–depth model is constrained by counting of varves. As the non-varved part does not possess a robust relative age model, we considered both parts of the sapropel separately during cyclicity analysis (Fig. 3).

Fig. 3: 3–2π MTM spectral analyses.
figure 3

a,b, Power spectra of the alkenone-based SST time series from the non-varved (a) and varved (b) portion of sapropel S5 presented on the periodicity domain, with 95% CLs calculated by different techniques indicated in the legend (see Methods for more information on how CLs were determined).

Full size image

The non-varved part of the sapropel exhibits a spectral peak that exceeds 95% confidence levels (CL) at 23 years for both the conventional and robust CL (Fig. 3a; see Methods for explanation of CL) and does not exhibit any significant spectral peak in the multidecadal band. However, the 23-year peak does not exceed the 95% CL of the power-law test (see Methods) and hence cannot be reliably differentiated from red noise. The power spectrum of the varved part of the record, on the contrary, points to the potential existence of multidecadal atmospheric oscillations as it exhibits decadal to centennial spectral peaks that exceed the 95% CL at periods of 20, 32, ~100, ~115 and ~200 years (Fig. 3b). However, only 20 years cyclicity exceeds 95% CL in the power-law test, suggesting that other periodicities cannot be considered as robust. The observed ~20-year cycles closely match 22-year solar cycles32. However, we consider those cyclicities with caution as we did not observe the more prominent 11-year solar cycles, and we also highlight that such oscillation might also be related to an internal mode of ocean dynamics33. Moreover, wavelet analyses show that 20-year cycles do not persist throughout the time series (Extended Data Fig. 7).

To explore the amplitude of decadal LIG Mediterranean SST change, we used an approach with 30-year-wide sliding windows (which is the typical time period considered in climate studies) and determined the average SST change between two windows in 10-year steps (red and blue bars in Fig. 4 and Extended Data Fig. 8). Hereafter we refer to this metric as the decadal rate of SST change (DRCSST). With an identical approach, we also calculated current and projected future DRCSST using annually resolved time series of the Extended Reconstructed Sea Surface Temperature dataset (ERSST)34 in combination with SST model projections from Community Earth System Model (CESM1.2) simulations under Representative Concentration Pathway (RCP) RCP4.5 and RCP8.5 scenarios; see Methods. To test the sensitivity of our results on the selection of the specific climate model, we performed the same rate of SST change analyses on a set of Coupled Model Intercomparison Project Phase 5 (CMIP5)35 and 6 (CMIP6)36 projections (see Methods and Extended Data Figs. 9 and 10).

Fig. 4: LIG, present and future Mediterranean SST, their decadal SST rate of change and centennial SST trends.
figure 4

a,b, Observational and projected SST of the eastern Mediterranean for the period from 1860 to 2099 (splice of ERSST34 observational data and CESM1.2 model projections under the RCP8.5 (a) and RCP4.5 (b) scenarios; model projections are highlighted by grey background). c, Reconstructed alkenone-based SST of the eastern Mediterranean during the period of sapropel S5 formation (grey area marks non-varved part of the sapropel; see Extended Data Fig. 1 for data uncertainty). Red and orange dashed lines indicate the maximal projected SST for the twenty-first century under RCP8.5 (25.1 °C) and RCP4.5 (23.4 °C) scenarios, respectively, while the dashed grey line indicates the average ERSST34 for the period from 1860 to 2020 (20.5 °C). d,e, Decadal rate of SST change for the period from 1860 to 2099 under the RCP8.5 (d) and RCP4.5 (e) scenarios. f, Decadal rate of SST change for the period of sapropel S5 deposition (see Extended Data Fig. 8 for data uncertainty). Red dashed line indicates the maximal projected decadal rate of SST change in the twenty-first century under the RCP8.5 scenario (0.48 °C per decade), and orange dashed line indicates 0.40 °C per decade wich is the same values of maximal projected decadal rate of SST change in the twenty-first century under the RCP4.5 scenario and maximal decadal SST rate of change from 1860 to 2020. g,h, Observational and projected centennial SST changes for the period from 1860 to 2099 under the RCP8.5 (g) and RCP4.5 (h) scenarios. i, Reconstructed centennial SST changes for the period of sapropel S5 deposition (see Extended Data Fig. 8 for data uncertainty). Red and orange dashed lines indicate the maximal centennial SST change in the twenty-first century under the RCP8.5 (4.3 °C per century) and RCP4.5 scenarios (2.8 °C per century), respectively, while the dashed grey line indicates SST warming of the last century (1920–2019; 0.68 °C per century). Warming and cooling trends are indicated with red and blue colours, respectively.

Full size image

We observe that the amplitude of DRCSST in sapropel S5 is generally larger in the non-varved part than in the varved portion. The highest positive DRCSST value of 1.1 °C in the non-varved portion is almost three times higher than the highest positive value of ~0.4 °C in the varved part. Conversely, the maximal negative DRCSST value in the varved part of −0.65 °C is similar to the corresponding value in the non-varved portion of −0.55 °C per decade. We assign the higher positive amplitudes in the non-varved portion to the highly stratified water column and reduced deep mixing of the heat gained.

In comparison, during the past 160 years, the eastern Mediterranean DRCSST values obtained from the ERSST have generally not exceeded warming or cooling above ~0.25 °C per decade, with the exception of the warming of ~0.4 °C per decade in the last DRCSST in our evaluation (SST difference between 1991–2020 and 1981–2010). DRCSST in the non-varved portion of the sapropel repeatedly exceeded 0.4 °C per decade. By contrast, in the varved portion, DRCSST did not exceed this threshold. While some DRCSST values during deposition of the varved portion of the sapropel exceeded the maximal cooling of −0.25 °C per decade of the instrumental era, the majority of cooling decades (269 out of 293) are within the modern envelope. This suggests that the highest amplitude of decadal variability in the instrumental era is similar to the highest amplitude of decadal SST variability from the varved sapropel portion, where the most recent and the highest instrumental DRCSST of 0.4 °C per decade is in the range of the highest values recorded for the varved sapropel portion (Fig. 4 and Extended Data Fig. 8).

According to CESM1.2 model projections, in the twenty-first century maximal DRCSST in the Mediterranean is not expected to exceed 0.4 °C per decade in the medium-to-low-emissions scenario RCP4.5 (0.30 °C per decade is obtained as an average maximum value by CMIP5 projections; Extended Data Fig. 9) and 0.48 °C per decade in the high-emissions scenario RCP8.5 (Fig. 4). The average value of the set of CMIP6 projections suggests a slightly higher maximal DRCSST of 0.46 °C per decade under the medium-to-low emissions Shared Socioeconomic Pathway (SSP)-245 scenario, whereas under the high-emissions SSP-585 scenario the maximal DRCSST reaches 0.50 °C per decade (Extended Data Fig. 10). This indicates that the highest DRCSST in the medium-to-low emissions scenario is in the similar range as maximal DRCSST recorded in the varved part of the sapropel. On the other hand, projected DRCSST in the high-emissions scenario exceeds DRCSST from the varved part of the sapropel, but is still lower than the highest DRCSST from the non-varved portion of the sapropel. This underscores that increased ocean stratification can serve as an additional amplifier of SST change.

Centennial trends of LIG SST versus projected future SST

Our study shows that the maximum amplitude of decadal SST variability in the eastern Mediterranean was higher in an early stage of sapropel formation of the LIG than during the instrumental era, but similar during deposition of the varved sapropel portion. However, we rarely observe persistent warming and cooling trends that last for more than several decades during the LIG, suggesting that decadal trends were cancelling each other on the centennial scale. Consequently, examination of the longer-term SST change from the period of sapropel deposition may add relevant information on the long-term trends of SST evolution during the LIG in comparison with the projected twenty-first century SST increase.

We determined centennial SST trends using a sliding window approach and calculating the linear regression in 100-year windows in 10-year steps (Fig. 4). Centennial warming and cooling trends had an average value of 0.7 °C and −0.7 °C per century, respectively (Fig. 4). The average centennial SST warming trend from the LIG slightly exceeds the centennial SST increase of 0.6 °C per century in the eastern Mediterranean from the past 100 years (1921–2020; Fig. 4). The maximum centennial SST increase of about 3.6 °C per century observed in the LIG SST exceeded the centennial warming from the last century by a factor of five. As this maximum centennial SST increase from 125.67 to 125.57 ka is forced by strong stratification, its comparison with current and projected climate trends has to be considered with caution. The maximum centennial SST increase of 1.8 °C per century in the varved sapropel portion suggests that during some periods of the LIG the amplitude of natural centennial SST variability in the eastern Mediterranean indeed exceeded the centennial SST increase in the last century (Fig. 4).

This 1.8 °C centennial SST warming is, however, lower than the highest projected centennial SST increase in the twenty-first century of 2.8 °C under the RCP4.5 scenario in the CESM1.2 projection (Fig. 4). This is in line with the average of the highest projected centennial SST increases from CMIP5 projections under the RCP4.5 scenario of 2.3 °C (Extended Data Fig. 9), as well as the highest projected centennial SST increases from all CMIP6 projections under the medium-to-low emissions SSP-245 scenario with an average of 3.1 °C per century. On the other hand, the projected SST increase in the high-emissions scenarios indicates that the centennial SST increase could outpace any analogous LIG SST increase already within the next few decades and be as high as 4.4 °C per century by the end of the twenty-first century (Fig. 4 and Extended Data Fig. 10). Such persistent warming would expose the eastern Mediterranean and its ecosystem to a climatic state whose consequences can no longer be assessed from analyses of past analogues such as the LIG. Consequently, limiting future greenhouse gas emissions to the medium-low scenario and below could perhaps keep the twenty-first century Mediterranean SST variability in the scale comparable to natural SST variability from the LIG and thus reduce the risk of major perturbations of eastern Mediterranean marine ecosystems.

Methods

Chronology

Sapropel S5 from sediment core M51/3-SL104 is a ~87-cm-thick sapropel layer (Supplementary Fig. 1; see Supplementary Information for more details on stratigraphy). Starting from 9 cm above the base to the top of the sapropel, a coupling of bright layers (formed by annual blooming of mat-forming diatoms) and dark layers (composed of terrigenous silt, clay and mixed assemblages of mostly solitary diatoms) is evident, indicating the formation of annual lamination (varves)13,14. Generating a precise age model for the whole period of sapropel formation is challenging due to large uncertainties in the determination of absolute dates. This makes a determination of the absolute timing, as well as the exact duration of the period of sapropel deposition, difficult. We adopted the age model introduced in ref. 26, which is based on the microscopic counting of 2,889 varves from the exact same core and further correlation of detected sapropel events to the same events in other sapropels. This correlation of the sapropel events relies on a previous study37, which established a spatially unambiguous, consistent and repeatable sequence of events in δ18O, δ13C and foraminifera abundance records, leading to the establishment of common tie points among the majority of S5 sapropels in the eastern Mediterranean. Following the same approach, ref. 26 confirmed these tie points in δ18O and δ13C values and foraminifera abundance of Globigerinoides ruber and Globigerina bulloides in M51/3 SL104 core and correlated M51/3-SL104 to ODP 971 A core as a reference core, concluding that the varved part of the sapropel spans 61% of the total duration of the sapropel deposition and the non-varved part 39% (that is, ~1,900 years). The absolute chronology is based on the lag of 3,000 years between the midpoint of the youngest S1 sapropel (8.5 ka) and the correlated maxima in insolation (11.5 ka)22,27, suggesting the onset of the sapropel S5 formation at ~126.4 ka26.

Because the major part of this age model is constrained by varve counting, the varved interval of the core provides a detailed account of the change in sedimentation rates and a precise relative timescale. The age model for the non-varved part is characterized by higher uncertainty. The obtained total duration of ~4,800 years is ~2,000 years shorter relative to suggestions by refs. 24,38, but within the dating uncertainties of the same study (see Supplementary Discussion). If the varved section of the sapropel is considered continuous, then the duration of the non-varved portion derived from the biostratigraphic correlation is unlikely to represent less time than estimated, but may potentially contain more time (such that the total duration is closer to the estimate by refs. 24,38.). The chronology of the non-varved section should therefore be considered with caution (see Supplementary Information and Supplementary Fig. 2 for alternative age model).

MSI

MSI is a laser-based technique that collects mass spectra of organic compounds at micrometre-scale spatial resolution and thus allows for the detection and visualization of biomarkers on intact sediment core sections at exceptionally high resolution15,16,17. Application of MSI for ultra-high-resolution biomarker-based SST reconstruction has been recently validated on sediments from the Santa Barbara Basin and its comparison with instrumental data39 (Supplementary Discussion). In our study, data were generated on the sediment core M51/3-SL104. Two U-channels were subsampled with four 25-cm-long LL-channels. Subsequently, X-ray photographs of LL-channels with clearly marked 5 cm subsamples were taken. Sample preparation was performed according to ref. 15. In short, samples were divided into 5-cm-long subsamples and freeze-dried. To stabilize subsamples, they were afterwards embedded in a mixture of 5% gelatin and 1% sodium carboxymethyl cellulose and frozen. Each embedded subsample was sectioned into 60-µm-thick slices using a cryomicrotome. The obtained slices were placed on indium-tin-oxide-coated glass slides and dried in a desiccator. MSI was performed using a 7 T solariX XR Fourier-transform ion cyclotron resonance mass spectrometer (FT-ICR-MS) coupled to a matrix-assisted laser desorption/ionization (MALDI) source equipped with a Smartbeam II laser (Bruker Daltonik). Analyses were performed in positive ionization mode. To enhance the sensitivity in the m/z range of targeted biomarkers, data were acquired in the continuous accumulation of selected ions (CASI) mode with m/z range at 554 ± 12. All data were acquired with a 50% data reduction in ftmsControl 2.1.0 (Bruker Daltonik). Spatial resolution was obtained by rastering the ionizing laser across the sample in a defined rectangular area at a 200 µm spot distance. Alkenones and isorenieratene were detected as Na+ adducts. Data were mass-calibrated to pyropheophorbide a in Data Analysis 4.4 (Bruker Daltonik). On average, ~13,000 spectra were obtained for each 5 cm subsample with a resolution of 200 μm. Afterwards, selected mass spectra in the region of interest were exported. The export routine generates a comma-separated file including spot position, m/z, intensity and signal-to-noise ratio (SNR) for every measured position. Exported data were further evaluated with MATLAB R2017a (MathWorks), where compounds of interest were identified according to a mass accuracy of ±0.003 Da. As data were already reduced during acquisition (50%), only the most prominent peaks were recorded. Consequently, a very low SNR threshold of 0.8 was applied.

Owing to the potential expansion of the sedimentary matrix due to sample preparation, the mapped areas were corrected for initial dimensions. The correction was done by visual identification of tie points in the X-ray photograph and imaging data, following the procedure from ref. 39. In short, three previously marked and clearly visible teach points were chosen on each ~5 cm section image and three corresponding teaching points on the X-ray photograph. The corresponding xy-spot coordinates of the teaching points on the measured ~5-cm slices were read off from Bruker Daltonics flexImaging software (version 4.1). Using the spot coordinates of the teaching points from the individual cryomicrotome slices as moving points and the pixel coordinates of the teaching points on the X-ray image as fixed reference points, the data maps of the molecular proxies were loaded as xyz coordinates and mapped onto the X-ray photograph using fitgeotrans in MATLAB.

The intensity of the two alkenone species relevant to the \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) proxy (C37:2 and C37:3) were recorded for each individual laser spot. Only spots in which both compounds were detected were further considered. Intensity values were then summed up over the depth corresponding to one 200 µm layer. If at least ten spots presenting both compounds were available for a single horizon, data quality criteria were satisfied18. \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) was calculated as follows20:

$${{U}}_{37}^{{{{\mathrm{K}}}}\prime } = {{{\mathrm{C}}}}_{37:2}/({{{\mathrm{C}}}}_{37:2} + {{{\mathrm{C}}}}_{37:3})$$
(1)

Examination of the sediment trap and core-top \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) index values at various sites of the Mediterranean Basin suggested that the spring bloom period does not remarkably imprint the temperatures recorded in the sediments, but that the sedimentary temperature estimates would rather reflect annually integrated SST, with a major influence of the autumnal post-bloom development of the coccolithophores in the euphotic zone40,41,42,43. Thus, and in line with our approximately annual resolution for the varved part of the record, data points were transformed into SST using the calibration from ref. 44 to obtain annual average SST:

$${{{\mathrm{SST}}}} = ({{U}}_{37}^{{{{\mathrm{K}}}}\prime } - 0.044)/0.033$$
(2)

Variability of \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) values in each horizontal layer was estimated for each 200 µm depth increment using a bootstrap approach45 (Extended Data Fig. 1a). To assess the uncertainty of SST estimates, we considered the analytical uncertainty related to the conversion of MSI data to quasi-gas chromatography data, which historically serve as the basis for \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) calibration to SST; this uncertainty ±0.012 \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) units is defined by the standard error of the linear regression between the \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) values from 12 samples from the same depth intervals obtained by both techniques (Extended Data Figs. 3 and 4). On top of this range, the uncertainty of the global calibration from ref. 44 (0.05 \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) units) was applied, resulting in a range of ±1.9 °C (Extended Data Fig. 1b).

Gas chromatography flame ionization detector

Ten samples in 1 cm and 2 samples in 2 cm resolution were analysed by gas chromatography flame ionization detector (GC-FID) to crosscheck FT-ICR-MS data with the conventional method. Around 0.5–0.6 g of freeze-dried sediment samples were extracted using the modified Bligh and Dyer technique46,47,48. For a 50% aliquot of the total lipid extract, aminopropyl solid-phase extraction was performed, where the ketone fraction containing the alkenones was eluted using dichlormethane/acetone (9:1). The extract was dried under a stream of nitrogen and redissolved in 100 µl n-hexane for GC-FID analysis. The alkenones were detected using a ThermoFinnigan Trace GC-FID equipped with a splitless injector and a Restek Rxi-5ms capillary column (30 m × 0.25 mm ID). The GC-FID was programmed to an initial temperature of 60 °C (1 min isothermal) to 150 °C at a rate of 10 °C min−1, and then to 310 °C at a rate of 4 °C min−1. CG-based \({{U}}_{37}^{{{{\mathrm{K}}}}\prime }\) and SST were calculated in the same way as for MSI data following refs. 20,44. Data are presented in Extended Data Fig. 3.

Spectral analysis

Prior to spectral analysis, all time series were detrended through the fitting and removing of a lowess smoother. This smoother only contained variability with wavelengths >600 years, hence not impacting the centennial and decadal oscillations that form the focus of this study. To detrend the data, the noLow function with the freely available R library astrochron49 was used. Subsequently, proxy time series were interpolated to equally spaced intervals of 4 years in the non-varved portion of the sapropel, and 0.8 years in the varved portion of the sapropel. These sampling intervals closely correspond to the median sampling intervals of the respective portion. All spectral analyses in this study were carried out using the multitaper method (MTM) with three Slepian 2π tapers50. Spectral analyses were conducted separately for the non-varved and varved parts of the sapropel because the age model of the non-varved part is considerably more uncertain than the age–depth model for the varved part.

Despite the well-constrained chronology, testing whether SST periodicity reflects a true climatic oscillation or just a transient feature on a ‘red noise’ background remains challenging51. Hence, we here present a critical evaluation of the spectral character of our SST time series, using multiple techniques to calculate CLs. These are estimated through the conventional and through the robust lag-1 autoregressive model (AR1)52. The former has been demonstrated to be less biased towards false positives than the robust AR1 model53, especially in the low frequency part of the spectrum. This is relevant because our focus lies with decadal and centennial climate oscillations and the centennial band corresponds to the low-frequency part of the spectrum. To further scrutinize the risk of false positives, we calculated noise levels that are the product of purely red-noise time series with the same power-law spectrum (β)54,55 as the examined time series. The key parameter β, determining the slope of the red-noise spectrum, is derived by fitting a power-law through the SST data MTM power spectra. Using a Monte Carlo approach, we then generated 1,000 purely random datasets with the same β as our data and calculated 95% CL from their spectral distributions.

Sliding-window analysis of rates of SST change

A sliding-window approach is adopted to calculate the rate of SST change between two subsequent analysis windows. Sliding windows are 30 years wide and slide shift with increments of 10 years on the age model from ref. 26. The settings of this sliding window approach are motivated by classical reference periods in climate science, as to allow for a direct comparison between past climate reconstructions and future climate projections. Within each analysis window, the average SST is calculated by taking the non-weighted mean of all data points within the window. The rate of SST change is subsequently calculated by taking the difference between the average temperatures of two subsequent (that is, 10 years apart) analysis windows.

To assess the uncertainty on the average SSTs within each window, we calculated the standard error of the mean within each analysis window. Uncertainties are then propagated by means of a Monte Carlo approach that estimates the uncertainties on the calculated rates of SST change: the Monte Carlo procedure resamples the SST temperature within each window, assuming a normal distribution described by the average SST and the standard error of the mean SST. After resampling, the procedure re-calculates all rates of SST change between subsequent windows. The Monte Carlo procedure was run 1,000 times to obtain a robust handle on rate of SST change uncertainty. This approach provides an individual measure of uncertainty for each rate of SST change calculation. For the sake of clarity, individual uncertainty assessments are available in Extended Data Fig. 8.

With an identical approach, we calculated decadal rates of SST change between 1860 and 2099, where for the period between 1860 and 2020 we used the ERSST dataset in annual resolution for the coordinates of our studied core34, and for the period between 2021–2099 we used annually resolved SST model projections from CESM1.2 simulations under RCP4.5 and RCP8.5 scenarios for the same coordinates. To calculate the offset/bias of the modelled SST and the ERSST time series we compared ERSST and CESM1.2 historical simulation from 1860 to 2000. We observed that the decadal SST amplitude is in the same range, however, ERSST exhibits values systematically lower by 0.4 °C on average. We have accounted for this offset when combining ERSST and CESM1.2 future projections by adding 0.4 °C to the SST values obtained by CESM1.2. Considering the close match between the SST amplitude of the CESM1.2 historical simulations and ERSST data, we rely on the CESM1.2 model projections for the visualization of the future climate trends in Fig. 4. However, to test the sensitivity of our results on the selection of the specific climate model and to use a larger dataset to determine the statistics of SST rate of change, we also performed the same rate of SST change analyses on the set of 31 CMIP535 and 22 CMIP636 ensembles (Extended Data Figs. 9 and 10).

Centennial SST trends were also calculated by means of a sliding window approach, with 100-year-wide windows and 10-year steps. Within each analysis window, a linear regression through all data points within the window forms the basis for calculating the SST rate of change on centennial timescales. The sliding window analyses were carried out on the R platform56 for statistical computing.

Climate model projections

Transient climate simulations were performed with version 1.2.2 of the Community Earth System Model (CESM1.2), which is a fully coupled global climate model57. The coupled components include the Community Atmosphere Model (CAM5), the Community Land Model (CLM4.0), the Parallel Ocean Program (POP2) and the sea-ice model Community Ice CodE (CICE4). The ocean and sea-ice models share the same horizontal grid with a nominal 1° resolution. The ocean/sea-ice model grid resolution varies and is higher around Greenland, with the North Pole displaced, as well as around the Equator. The ocean model grid has 60 levels in the vertical. The atmosphere model has a finite volume dynamical core with a uniform horizontal resolution of 1.25° × 0.9° which is shared with the land model grid. The atmosphere model runs with 30 levels in the vertical. Two twenty-first century scenarios were run with the model, the medium-low emissions scenario RCP4.5 and the high-emissions scenario RCP8.558. For further information on the model and experimental setup the reader is referred to ref. 59. The modelled surface temperatures shown and analysed in this study were taken as annual means from the grid box that includes the location of core M51/3-SL104 in the eastern Mediterranean.

To test whether our results from the CESM1.2 model projections are robust, we also analysed the set of 31 coupled atmosphere–ocean model runs from the CMIP535, as well as 22 model runs from the CMIP636 (see Supplementary Table 1 for information on the models used). For CMIP5 models, we use the SST fields of the historical simulations (hist) covering 1850–2005 and the future projections (rcp45) covering 2006–2099. The historical simulations are forced by reconstructed natural and anthropogenic radiative forcing from solar variations, greenhouse gas concentrations, and volcanic and anthropogenic aerosols. The twenty-first century scenarios were run under the medium-low emissions scenario RCP4.5. For the CMIP6 simulations, we use the SST fields of the historical simulations covering 1850–2014 and the future projections under the medium-to-low emissions scenario (SSP-245) and the high-emissions scenario (SSP-585) covering 2015–2099. The modelled surface temperatures were taken as annual means from the grid box that includes the location of core M51/3-SL104 in the eastern Mediterranean.