Increased abyssal ocean density stratification across the Middle Pleistocene Transition

The cause of the fundamental reorganisation of the climate system during the Middle Pleistocene Transition (MPT; ∼ 1.25–0.65 million years ago (Ma)), when glacial cycles intensified and lengthened from 41 thousand year (kyr) to quasi-100 kyr periodicity, remains one of the enduring unsolved questions in palaeoclimate science. Increasing attention has focused on the role of the marine carbon cycle, with enhanced carbon storage in the deep sea linked to changes in deep-water temperature and salinity that increase abyssal ocean density stratification. Here we report new high-resolution reconstructions of deep-water temperature and the isotopic composition of seawater (δ¹⁸O_seawater), a proxy for ice volume and salinity, from North Atlantic Site U1385, using paired benthic foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ and oxygen isotope (δ¹⁸O) analyses. Together with published records, we present basin-scale and globally distributed compilations of deep-water temperature and δ¹⁸O_seawater to assess changes in abyssal ocean stratification for the past 1.5 million years (Myr). Across the MPT, interbasinal gradients suggest North Atlantic deep-water became colder while Pacific deep-water became more saline during glacial periods after ∼ 930 thousand years ago (ka). Deep-ocean δ¹⁸O_seawater increased in both the Atlantic and Pacific across the MPT indicating increased continental ice volume. But, the increase was greater in the Pacific than the Atlantic, which we suggest reflects increased salinity of Southern Component Water (SCW). We propose the following sequence of changes during the MPT as a working hypothesis: (1) freshwater input to the marginal seas around Antarctica was reduced beginning at 930 ka by decreased melting of the Antarctic Ice Sheet and/or marine ice shelves and increased sea ice formation in the Southern Ocean; (2) the salinity and density of SCW increased, resulting in enhanced abyssal ocean density stratification and rendering the deep ocean a more effective carbon trap; (3) together with increased export production of organic matter and reduced deep-to-surface water exchange in the Southern Ocean, carbon storage in the deep ocean increased, lowering glacial atmospheric pCO₂; (4) lower pCO₂ permitted the growth of larger continental ice sheets, which reached a critical size and lengthened the glacial cycles. Our hypothesis supports an important role for abyssal ocean density stratification in the MPT, and requires further testing with additional benthic $\mathrm{Mg}/\mathrm{Ca}$-δ¹⁸O records, numerical model simulations, and forthcoming atmospheric pCO₂ and mean ocean temperature results from the Beyond EPICA–Oldest Ice core from Antarctica.

1 Introduction

The Middle Pleistocene Transition (MPT), which occurred between ∼ 1.25 and 0.65 million years ago (Ma), marks a major reorganisation of the global climate system and remains a central topic of palaeoclimate research since its first recognition over 50 years ago (Shackleton and Opdyke, 1976). Benthic foraminiferal oxygen isotope (δ¹⁸O) records are marked by a shift from dominant 41 thousand years (kyr) to quasi-100 kyr glacial-interglacial cycles across the MPT. There was no commensurate change in orbital forcing (Clark et al., 2006; Pisias and Moore, 1981), suggesting processes internal to the climate system were responsible. The exact cause(s) of the MPT remain(s) an open question (see reviews by Berends et al., 2021a; Clark et al., 2006; Herbert, 2023; and McClymont et al., 2013). Uncertainties include whether the transition was gradual or abrupt (Sosdian and Rosenthal, 2009; Elderfield et al., 2012; Ford et al., 2016; Legrain et al., 2023) and whether it was driven primarily by Northern or Southern Hemisphere processes (An et al., 2024; Basak et al., 2026; Chalk et al., 2017; Farmer et al., 2019; Hines et al., 2024; Pérez-Montero et al., 2026; Starr et al., 2021; Williams et al., 2024; Wirths et al., 2025; Yehudai et al., 2021). Most hypotheses invoke a long-term global cooling, commonly attributed to a reduction in atmospheric carbon dioxide (pCO₂) (Berends et al., 2021a; Clark et al., 2006; Willeit et al., 2019). Yet, despite extensive work, the relative roles and interactions of ice-sheet dynamics and carbon-cycle feedbacks, and their timing, remain poorly constrained.

Proposed explanations for the MPT can be broadly grouped into three categories: (i) internal ice-sheet dynamics, most notably the regolith hypothesis invoking progressive changes in basal conditions leading to the exposure of high-friction crystalline bedrock (Clark et al., 2006; Clark and Pollard, 1998); (ii) changes in ice-sheet geometry and extent, particularly involving North American and/or Antarctic ice sheets (An et al., 2024; Bintanja and Van De Wal, 2008; Elderfield et al., 2012; Raymo et al., 2006; Wirths et al., 2025); and (iii) feedbacks associated with long-term changes in the global carbon cycle, including deep-sea carbon storage and atmospheric pCO₂ variability (Berger et al., 1999; Chalk et al., 2017; Farmer et al., 2019; Lear et al., 2016; Qin et al., 2022; Shackleton, 2000; Thomas et al., 2022; Willeit et al., 2019). Increasing attention has been paid to the potential role of marine carbon cycle processes (Herbert, 2023), which are closely linked to changes in deep ocean circulation during the MPT (Hall et al., 2001; Hasenfratz et al., 2019; Kim et al., 2021; Li et al., 2026; Martínez-Garcia et al., 2010; Pena and Goldstein, 2014; Qin et al., 2022; Raymo et al., 1997). A “thermohaline circulation crisis” (Pena and Goldstein, 2014) (THC) was proposed for the deep Atlantic across the MPT, but the magnitude and extent of such an event have been questioned recently (Basak et al., 2026; Hines et al., 2024; Williams et al., 2024). Furthermore, recent findings from Allan Hills blue ice area of Antarctica show no significant decline in mean atmospheric pCO₂ over the MPT (Marks-Peterson et al., 2026), which does not necessarily preclude changes in the glacial values of atmospheric pCO₂ because the mean values are likely weighted towards interglacial periods. More recently, modelling simulations suggest that declining atmospheric pCO₂ and sea-level change may have driven a threshold response of the Antarctic Ice Sheet, resulting in expanded ice shelves and increasingly strong coupling between Antarctic ice-sheet variability, Northern Hemisphere ice volume, and other climate components following the MPT (Yun and Timmermann, 2026). Forthcoming atmospheric pCO₂ measurements from the Beyond-EPICA ice core are highly anticipated and will likely resolve the timing and magnitude of atmospheric pCO₂ changes across the MPT.

Whilst numerous studies implicate enhanced carbon storage in the deep sea and lowered glacial atmospheric pCO₂ across the MPT (Billups et al., 2018; Chalk et al., 2017; Farmer et al., 2019; Gildor and Tziperman, 2001; Hodell et al., 2003; Jaccard et al., 2013; Lear et al., 2016; Sigman et al., 2010; Starr et al., 2021; Thomas et al., 2022), no direct comparisons of the temperature and salinity controls on the density of deep water masses have been made between the North Atlantic and Pacific Ocean basins across the transition. This limits assessment of large-scale changes in abyssal density structure, and their potential role in carbon storage and atmospheric pCO₂ change (Adkins, 2013; Broecker and Peng, 1998, 1982). Here, we combine new and previously published records of $\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperature and oxygen isotopic composition of seawater (δ¹⁸O_seawater), a proxy for both ice volume and salinity, to evaluate how abyssal ocean density structure changed across the MPT.

Benthic δ¹⁸O records are invaluable tools for reconstructing glacial-interglacial variability of the Quaternary (Shackleton, 1967; Shackleton and Opdyke, 1973). However, interpretation of benthic δ¹⁸O records is complicated because the signal is dependent on both temperature and the δ¹⁸O_seawater, the latter varying with global ice volume and regional hydrography including salinity changes at sites of deep-water formation (Waelbroeck et al., 2002). Deconvolution of benthic δ¹⁸O into its temperature and δ¹⁸O_seawater components has been accomplished by paired measurement of $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O of calcite (Elderfield et al., 2010, 2012; Ford et al., 2016; Ford and Raymo, 2020; Hasenfratz et al., 2017, 2019; Sosdian and Rosenthal, 2009), with estimates of propagated uncertainty (Thirumalai et al., 2016). Other indirect approaches include inverse modelling that combines benthic δ¹⁸O data with ice sheet models of varying complexity (Berends et al., 2021b; Bintanja and Van De Wal, 2008) and an iterative process-based approach that assesses changes in relationships between the oxygen isotopic composition of ice (δ¹⁸O_ice), sea level, temperature, and δ¹⁸O_seawater (Rohling et al., 2021). More recently, changes in global mean sea surface temperatures (ΔGMSST) combined with proxy-based deep ocean temperature estimates, have been used to infer changes in mean ocean temperature (ΔMOT), enabling separation of temperature and ice volume components from a global benthic δ¹⁸O stack (Clark et al., 2024, 2025a). These different approaches can yield widely divergent interpretations of Quaternary temperature and ice volume evolution. For example, during Marine Isotope Stage (MIS) 22 (∼ 900 ka), Rohling et al. (2021) infer a substantial increase in ice volume with limited deep-water cooling, consistent with earlier $\mathrm{Mg}/\mathrm{Ca}$-based reconstructions (Elderfield et al., 2012; Ford and Raymo, 2020) (Fig. S1 in the Supplement). In contrast, the reconstruction of Clark et al. (2025a) shows pronounced deep-water cooling with relatively little change in ice volume across the MPT.

Here we present a new 1.5-million-year-long record of North Atlantic deep-water temperature, using $\mathrm{Mg}/\mathrm{Ca}$ measured in calcite of two infaunal benthic foraminifera (Uvigerina peregrina and Globobulimina affinis) from sediment cores recovered at International Ocean Discovery Program (IODP) Site U1385 (37°34.3^′ N, 10°7.6^′ W, 2578 metres below sea level [m b.s.l.]) from the Iberian Margin (Hodell et al., 2015, 2023a) (Figs. 1, S2 and S3; Table 1). The $\mathrm{Mg}/\mathrm{Ca}$-derived temperature is paired with benthic δ¹⁸O to isolate the δ¹⁸O_seawater signal (Sect. 2.3; Figs. 2 and S4). Our record provides an important complement to existing North Atlantic data from Deep Sea Drilling Project (DSDP) Site 607 (Sosdian and Rosenthal, 2009), where $\mathrm{Mg}/\mathrm{Ca}$-derived temperature estimates were based partly on Cibicidoides wuellerstorfi. These results were called into question because this epibenthic taxon can be affected by changes in carbonate ion concentrations of bottom waters (Yu and Broecker, 2010) (Sect. 3.1; Supplement). This criticism was subsequently refuted (Sosdian and Rosenthal, 2010) prompting additional $\mathrm{Mg}/\mathrm{Ca}$ measurements of infaunal benthic taxa (Uvigerina spp.) at Site 607 (Ford et al., 2016).

Figure 1 Meridional profiles of deep-water δ¹⁸O_seawater showing the site locations. (a) δ¹⁸O_seawater at 2578 m b.s.l., corresponding to the depth of IODP Site U1385 (this study; yellow diamond) in the Northeast Atlantic, along with the geographic location of other sites referenced in the text: North Atlantic DSDP Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) (pink triangle); Southwestern Pacific ODP Site 1123 (Elderfield et al., 2012) (purple square); North Pacific ODP Site 1208 (Ford and Raymo, 2020) (lilac circle); and the Southern Ocean ODP Site 1094 (Hasenfratz et al., 2019) (red inverted triangle) located in the Atlantic sector of the Antarctic Zone. Meridional transects shown in panels (b) and (c) are indicated by white dashed lines. (b, c) Atlantic and Pacific meridional profiles of deep-water δ¹⁸O_seawater (‰, SMOW) showing site locations. Maps were created using GLODAPv2.2021 and GLODAPv2.2023 and GEOSECS in Ocean Data View (Schlitzer, 2019; http://odv.awi.de/, last access: 7 July 2026). Data for panel (c) are from (LeGrande and Schmidt, 2006).

Table 1 Site location details and modern-day deep-water properties.

^* Note: All salinity values reported follow the Practical Salinity Scale (PSS-78) and are unitless.
Unless otherwise stated, temperature and salinity values are from GLODAPv2.2021 and GLODAPv2.2023 and GEOSECS accessed via Ocean Data View (Schlitzer, 2019; http://odv.awi.de/, last access: 7 July 2026), while δ¹⁸O_seawater data for Sites 1123 and 1208 are from LeGrande and Schmidt (2006). Coordinates and water depths are taken from the original site publications cited in the first column.

Figure 2 Iberian Margin deep-water $\mathrm{Mg}/\mathrm{Ca}$-derived temperature and δ¹⁸O_seawater estimated using PSU Solver. (a) Composite benthic δ¹⁸O record, primarily Uvigerina peregrina, with Cibicidoides wuellerstorfi and Globobulimina affinis values adjusted to Uvigerina peregrina by +0.64 ‰ and −0.30 ‰ respectively (black line). Numbers denote Marine Isotope Stage (MIS) interglacials (odd, above) and glacials (even, below). (b) $\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperature (light green line and markers) with 1σ and 2σ uncertainty envelopes (dark and light grey, respectively) incorporating both analytical and calibration uncertainties. (c) Deconvolved δ¹⁸O_seawater (dark green line and markers; note the inverted axis) with uncertainties as in (b). All panels show Site U1385 (this study) together with records from Birner et al. (2016) (pink lines with crosses during MIS 41–37), MD01-2444 (Skinner and Elderfield, 2007) (light purple line) and MD99-2334 (Skinner et al., 2003) (dark purple line with crosses).

We combine (stack) North Atlantic records from Sites U1385 and 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) and compare them with equivalent Pacific records (Sites 1123 and 1208; Elderfield et al., 2012; Ford and Raymo, 2020) to reconstruct the basinal histories of temperature and salinity change (Sect. 2.5–2.6). The Atlantic and Pacific records are then combined with Site 1094 from the South Atlantic sector of the Southern Ocean (Hasenfratz et al., 2019) (Figs. 1, S2 and S3; Table 1) to produce globally distributed stacks. Although relatively few long reconstructions based on deconvolved benthic $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O records are available for stacking over the past 1.5 million years (Myr), some distinct differences emerge in the patterns of inferred deep-ocean temperature and salinity change between the North Atlantic and Pacific basins.

All study sites are deeper than 2500 m and are considered representative of the deep ocean (Figs. 1, S2 and S3; Table 1). Today, Sites U1385 (2578 m b.s.l.) and 607 (3427 m b.s.l.; Ford et al., 2016; Sosdian and Rosenthal, 2009) are bathed by North Atlantic Deep Water (NADW) and exhibit similar seawater physical properties (Table 1). Although both sites lie within the deep North Atlantic, they are influenced by different mixtures of deep-water masses. In the modern ocean, Sites U1385 and 607 occupy different parts of Lower North Atlantic Deep Water and receive varying proportions of Labrador Sea Water, Iceland–Scotland Overflow Water, Denmark Strait Overflow Water and Antarctic Bottom Water (Liu and Tanhua, 2021). The relative contributions of these water masses to each site also likely varied through glacial–interglacial cycles. Modern hydrographic observations nevertheless place both sites within the Lower North Atlantic Deep Water neutral density range (γⁿ ≈ 27.95–28.10 kg m⁻³). The North Atlantic stack is therefore intended to represent the average properties of deep North Atlantic waters rather than a single homogeneous water mass. Pacific Sites 1123 (3290 m b.s.l.; Elderfield et al., 2012) and 1208 (3346 m b.s.l.; Ford and Raymo, 2020) are located within Lower Circumpolar Deep Water (LCDW) and Pacific Deep Water (PDW) (Elderfield et al., 2012; Ford and Raymo, 2020; Insua et al., 2014), respectively, and share broadly similar physical properties originating from the Southern Ocean (Talley, 2013). The deep western boundary current (DWBC) of the Southwest Pacific is the main pathway through which most of the cold deep water formed around Antarctica enters the Pacific Ocean (Chandler et al., 2024). The deep water gradually mixes as it flows northward modifying its properties, leading to the formation of Pacific Deep Water (PDW). Given the locations of Site 1123 within the DWBC and 1208 in the PDW, it is reasonable to assume they broadly reflect conditions in the deep Pacific.

Deep Pacific sites are ∼ 1.2 °C colder and ∼ 0.3 lower in salinity (PSS-78) than those in the North Atlantic, reflecting the formation of relatively warm and salty NADW (Figs. S2 and S3; Table 1). Correspondingly, the δ¹⁸O_seawater values of the Atlantic sites are on average higher (∼ 0.25 ‰) higher than those of the deep Pacific (Fig. 1; Table 1) (Frew et al., 2000). The Atlantic and Pacific records are merged with Site 1094 (2807 m b.s.l.; Hasenfratz et al., 2019) from the South Atlantic sector of the Southern Ocean, which is currently bathed by LCDW.

2 Materials and Methods

2.1 Site details and chronology

Four holes drilled at Site U1385 were spliced to produce a continuous composite section extending to 166.5 m composite depth (mcd), corresponding to a basal age of ∼ 1.45 Ma during Marine Isotope Stage (MIS) 47 (Hodell et al., 2015, 2023a). The Site U1385 age model was established by alignment of the benthic oxygen isotope (δ¹⁸O) record to the LR04 reference stack (Hodell et al., 2023a; Lisiecki and Raymo, 2005), with chronostratigraphic ages assigned by linear interpolation between age-depth control points.

All previously published records used to compile the ocean stacks are presented on their original published age models. Most are based on alignment of benthic foraminiferal δ¹⁸O to LR04, including DSDP Site 607 combined with Chain 82-24-23PC (Ford et al., 2016; Sosdian and Rosenthal, 2009), and Ocean Drilling Program (ODP) Site 1123 (Elderfield et al., 2012). ODP Site 1208 (Ford and Raymo, 2020) was aligned to the Prob-stack using the HMM-Stack MATLAB code (Ahn et al., 2017; Ford and Raymo, 2020; Lin et al., 2014; Butcher et al., 2017), which is also based on the LR04 age model. The ODP Site 1094 age model (Hasenfratz et al., 2019) was derived by graphical alignment of benthic δ¹⁸O to LR04, with the original chronology of Jaccard et al. (2013) retained where benthic foraminifera were sparse (e.g. MIS 4).

The Iberian Margin piston cores MD99-2334 and MD01-2444 were independently aligned to Greenland ice-core records through synchronisation of Dansgaard–Oeschger events recorded in planktonic δ¹⁸O (Shackleton, 2000; Skinner et al., 2003; Skinner and Elderfield, 2007). Owing to the close agreement in both structure and absolute benthic δ¹⁸O values with Site U1385 over the interval of overlap, no further adjustments to these published age models were applied.

2.2 Stable isotopes

Oxygen and carbon isotopes were measured previously using Cibicidoidies wuellerstorfi from the >212 µm sediment coarse fraction (Hodell et al., 2023a, b). Additional analyses were made using Uvigerina peregrina and Globobulimina affinis selected from the 212–355 µm coarse fraction. All benthic oxygen isotope data were corrected to Uvigerina using an offset of 0.64 ‰ for C. wuellerstorfi and −0.3 ‰ for G. affinis. Tests were crushed between glass slides to facilitate cleaning. Isotopic measurements were made at the Godwin Laboratory for Palaeoclimate Research, Department of Earth Sciences, University of Cambridge following previously described methods (Hodell et al., 2015). Instrument precision was better than ±0.08 ‰ (1σ) for δ¹⁸O and ±0.06 ‰ (1σ) for δ¹³C.

Species-specific benthic δ¹⁸O offsets are commonly assumed to be constant when correcting between co-occurring taxa, although temporal variations have been reported (Hoogakker et al., 2010). Where same-species paired δ¹⁸O measurements were unavailable C. wuellerstorfi δ¹⁸O values were supplemented. Comparison of paired δ¹⁸O measurements from C. wuellerstorfi and G. affinis (n= 3580) at Site U1385 indicates that the corrected values agree closely with a mean residual difference of only −0.02 ‰ ± 0.23 ‰ (1σ). Restricting the comparison to 98 available paired samples used for the $\mathrm{Mg}/\mathrm{Ca}$ analyses yields a similarly small residual difference of −0.04 ‰ ± 0.26 ‰ (1σ). Together with the high sedimentation rate at Site U1385 (∼ 11 cm kyr⁻¹; Hodell et al., 2015) and interpolation of the reconstructions to 3 kyr resolution (Sect. 2.6), this suggests that any temporal offset associated with differing microhabitats is likely to be negligible relative to the temporal resolution of the study.

2.3 Foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ analyses

Benthic $\mathrm{Mg}/\mathrm{Ca}$, deep-water temperature and δ¹⁸O_seawater records from Site U1385 extend back to ∼ 1.45 Ma, spanning MIS 1–47. $\mathrm{Mg}/\mathrm{Ca}$ measurements on Uvigerina peregrina and Globobulimina affinis were conducted at ∼ 20 cm intervals, corresponding to a mean temporal resolution of ∼ 2000 years (yr), except for two intervals sampled at higher resolution: ∼ 4 cm spacing between ∼ 980–860 thousand years ago (ka) over MIS 27–21, and ∼ 5 cm spacing between ∼ 1300–1240 ka (MIS 40–37; Birner et al., 2016) yielding a temporal resolution of ∼ 500 years.

Three noticeable breaks (> 15 kyr) occur in the $\mathrm{Mg}/\mathrm{Ca}$-derived temperature and δ¹⁸O_seawater records because of low foraminiferal abundance between MIS 6–5 (∼ 132–116 ka), and MIS 16–15 (∼ 630–603 ka). In addition, a near 30 kyr hiatus removed the interval at Site 339-U1385 between MIS 12–11 (∼ 430–400 ka) at Termination V (Hodell et al., 2015, 2023a), which has now been recovered in companion Site 397-U1385 (Hodell et al., 2026). To produce a near-continuous Iberian Margin record spanning the past 1.5 Myr, Site U1385 records were supplemented through the last glacial cycle with $\mathrm{Mg}/\mathrm{Ca}$-derived temperature and δ¹⁸O_seawater data measured on Globobulimina affinis at nearby piston cores MD01-2444 covering 50–35 ka during the early glacial interval (Skinner and Elderfield, 2007), and MD99-2334 spanning 35–10 ka during the late glaciation (Skinner et al., 2003) (Fig. 2; Table 1). Following Skinner et al. (2007) a constant temperature offset of +0.6 °C was applied to $\mathrm{Mg}/\mathrm{Ca}$-derived temperatures from the deeper Core MD99-2334 to account for cooler modern deep-water temperature at that site location (Hodell et al., 2014; Skinner et al., 2003, 2007; Skinner and Elderfield, 2007; Table 1).

After data quality screening (Sect. 2.4), we report $\mathrm{Mg}/\mathrm{Ca}$ results for 864 samples of Uvigerina peregrina and 202 samples of Globobulimina affinis. All samples were oxidatively cleaned following the procedure of Barker et al. (2003). Trace elements were measured at the Godwin Laboratory for Palaeoclimate Research (Cambridge, UK) using inductively coupled plasma-optical emission spectroscopy (ICP-OES). Most analyses were performed on a Varian VISTA instrument (Birner et al., 2016) with the remaining ∼ 180 samples measured using an Agilent 5100 ICP-OES. Both instruments used the intensity-ratio calibration method of de Villiers et al. (2002). Interlaboratory comparison studies have validated the use of laboratory standards employed as reference materials for foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ analyses (Greaves et al., 2008). Repeated measurements of standards showed instrument precision for $\mathrm{Mg}/\mathrm{Ca}$ determinations to be better than 0.5 % relative standard deviation (r.s.d.) (Greaves, 2008; Greaves et al., 2008; de Villiers et al., 2002). Based on long-term standard reproducibility and 47 sets of replicate analyses (including four triplicates) measured on multiple aliquots of the same sample, the pooled standard deviation is ±0.09 mmol mol⁻¹ (1σ), corresponding to a mean precision of 8.2 % r.s.d. (see PANGAEA Data Repository, Thomas et al., 2026).

$\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperatures from Uvigerina peregrina and Globobulimina affinis show generally good agreement in both variability and mean values within uncertainty (1.4 ± 1.5 and 1.6 ± 1.0 °C, respectively) (Fig. S4a). Both species record similar mean δ¹⁸O_seawater values (∼  0.6 ‰ ± 0.4 ‰; Fig. S4b). This agreement is consistent with the application of species-specific $\mathrm{Mg}/\mathrm{Ca}$–temperature calibrations applied in this study (Elderfield et al., 2006, 2010, 2012; Weldeab et al., 2016) and justifies combining the two datasets. The individual Uvigerina peregrina and Globobulimina affinis time series were used as inputs for the PSU Solver calibration procedures described in Sect. 2.5 and combined into single records of deep-water temperature and δ¹⁸O_seawater (Fig. 2).

2.4 Data quality evaluation

$\mathrm{Mg}/\mathrm{Ca}$ ratios derived from foraminiferal calcite tests can be used to estimate palaeotemperatures provided that stringent cleaning procedures are followed (Barker et al., 2003; Larsson and Jung, 2025). To evaluate the efficiency of these cleaning procedures and to assess the effects of authigenic precipitation, we analysed element $/$ Ca ratios of co-measured trace elements – primarily Fe and Mn, but also Al, Ba, K, Na, Si, Ti and Zn (Barker et al., 2003; Elderfield et al., 2012; Greaves, 2008; Hasenfratz et al., 2017). Contamination by silicates and clay minerals can lead to significant uncertainty in measured $\mathrm{Mg}/\mathrm{Ca}$ values. Therefore, the removal of contaminant silicates is a critical step in the cleaning protocol, and its success can be evaluated through downcore covariation of $\mathrm{Mg}/\mathrm{Ca}$ and $\mathrm{Fe}/\mathrm{Ca}$ (or $\mathrm{Al}/\mathrm{Ca}$) records (Barker et al., 2003). $\mathrm{Mg}/\mathrm{Ca}$ records from the two foraminiferal species were assessed separately, as Globobulimina affinis is known to exhibit elevated $\mathrm{Mg}/\mathrm{Ca}$ ratios relative to other benthic foraminiferal species (Skinner et al., 2003; Weldeab et al., 2016). At Site U1385, five Uvigerina peregrina and three Globobulimina affinis samples were excluded from the $\mathrm{Mg}/\mathrm{Ca}$ dataset following identification of anomalously high co-occurring $\mathrm{Fe}/\mathrm{Ca}$ and $\mathrm{Mg}/\mathrm{Ca}$ values (Barker et al., 2003). Sample exclusion was based on co-occurring Fe/Ca–$\mathrm{Mg}/\mathrm{Ca}$ anomalies rather than a fixed $\mathrm{Fe}/\mathrm{Mg}$ threshold (Fig. S5). The $\mathrm{Mg}/\mathrm{Ca}$ dataset presented here includes Uvigerina peregrina samples with $\mathrm{Fe}/\mathrm{Ca}$ and $\mathrm{Fe}/\mathrm{Mg}$ values of up to 0.45 mmol mol⁻¹ and 0.33 mol mol⁻¹, respectively (and Globobulimina affinis values up to 0.79 mmol mol⁻¹ and 0.28 mol mol⁻¹), which exceed the typical $\mathrm{Fe}/\mathrm{Ca}$ (< 0.1 mmol mol⁻¹) and $\mathrm{Fe}/\mathrm{Mg}$ (< 0.03 mol mol⁻¹) values reported by Barker et al. (2003) and Elderfield et al. (2012).

Consistent with observations of Uvigerina spp. at Site 1123 (Elderfield et al., 2012), $\mathrm{Mn}/\mathrm{Ca}$ values for Uvigerina peregrina at Site U1385 show a slight downcore increase, with lowest glacial values rising from ∼ 0.012 mmol mol⁻¹ at 60 ka to ∼ 0.045 mmol mol⁻¹ at ∼ 1415 ka (Fig. S6). Although the maximum $\mathrm{Mn}/\mathrm{Ca}$ value for Uvigerina peregrina at Site U1385 (0.28 mmol mol⁻¹) is lower than the 0.36 mmol mol⁻¹ reported for Site 1123 (Elderfield et al., 2012), $\mathrm{Mn}/\mathrm{Ca}$ values exceeding ∼ 0.1 mmol mol⁻¹ may suggest the presence of authigenic Mn–Fe-oxide coatings on some foraminiferal tests (Elderfield et al., 2012; Hasenfratz et al., 2017, 2019). The incorporation of Mn into calcite tests has been shown to be species-specific (van Dijk et al., 2025). In Uvigerina peregrina, $\mathrm{Mn}/\mathrm{Ca}$ values range from 0.01 to 0.28 mmol mol⁻¹, while in Globobulimina affinis, values are generally higher, ranging from 0.04 to 1.88 mmol mol⁻¹, with two anomalous specimens exceeding this upper value. Despite these elevated $\mathrm{Mn}/\mathrm{Ca}$ values, the absence of a positive correlation between $\mathrm{Mn}/\mathrm{Ca}$ and $\mathrm{Mg}/\mathrm{Ca}$ indicates that the presence of Mn–Fe oxides does not significantly influence $\mathrm{Mg}/\mathrm{Ca}$ variability. We therefore interpret the higher $\mathrm{Mn}/\mathrm{Ca}$ values as reflecting variable contributions from authigenic Mn–Fe-oxide coatings (Skinner et al., 2003), or from species-specific Mn incorporation during biomineralization (van Dijk et al., 2025).

Diagenetic Mn–Fe-oxide coatings commonly form on foraminiferal shells in deep-sea sediment cores used in palaeoceanographic reconstructions. In oxygenated porewater dissolved manganese is precipitated onto carbonate tests as Mn⁴⁺ oxides, later becoming remobilized as Mn²⁺ and forming Mn-rich oxide coatings when organic matter undergoes anaerobic decomposition deeper in the sediment (Barker et al., 2003; Boiteau et al., 2012). These Mn-rich coatings can substantially impact bulk $\mathrm{Mg}/\mathrm{Ca}$ ratios and introduce bias into deep-water temperature reconstructions (Hasenfratz et al., 2017). Site specific studies reveal relatively consistent $\mathrm{Mg}/\mathrm{Mn}$ ratios in foraminiferal coatings, ranging globally between 0.17–0.32 mol mol⁻¹ [see Hasenfratz et al. (2017) Table 1 and references therein]. Regionally the upper end of this range – between 0.26 to 0.32 mol mol⁻¹ – is most applicable to the Atlantic Ocean, although individual sites may deviate notably from basin-wide averages (Hasenfratz et al., 2017). These $\mathrm{Mg}/\mathrm{Mn}$ ratios can be used to correct $\mathrm{Mg}/\mathrm{Ca}$ ratios associated with the calcite lattice using the following equation (Hasenfratz et al., 2017):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{Mg}/\mathrm{Ca}}_{\mathrm{corrected}}& ={\mathrm{Mg}/\mathrm{Ca}}_{\mathrm{measured}}\\ & -\left({\mathrm{Mn}/\mathrm{Ca}}_{\mathrm{measured}}\times {\mathrm{Mg}/\mathrm{Mn}}_{\mathrm{coating}}\right)\end{array}\end{array}$$

To maintain consistency, we applied a $\mathrm{Mg}/\mathrm{Mn}$_coating correction following Eq. (1), with the $\mathrm{Mg}/\mathrm{Mn}$_coating term fixed at 0.32 mol mol⁻¹, and applied uniformly to all $\mathrm{Mg}/\mathrm{Ca}$ data (Hasenfratz et al., 2017; de Lange et al., 1992) (Fig. S7). This value represents the maximum adjustment in the direction of lower temperatures and yields the least discrepancy between Site U1385 temperatures and those from the Pacific and South Atlantic. The correction has minimal influence on glacial $\mathrm{Mg}/\mathrm{Ca}$ values.

In total, $\mathrm{Mg}/\mathrm{Ca}$ was initially determined on 911 Uvigerina peregrina and 264 Globobulimina affinis samples. After quality control, we retain $\mathrm{Mg}/\mathrm{Ca}$ results for 864 samples of Uvigerina peregrina (including replicates) and 202 samples of Globobulimina affinis measurements, having excluded 42 and 59 samples, respectively, based on co-measured trace element indicators of potential contamination. Samples were rejected according to the following criteria:

• Clay contamination: Low [Ca] (< 10 ppm) coupled with high $\mathrm{Na}/\mathrm{Ca}$ (> 5.5 mmol mol⁻¹), and $\mathrm{K}/\mathrm{Ca}$ (> 0.4 mmol mol⁻¹), almost always associated with elevated $\mathrm{Mg}/\mathrm{Ca}$.

• Additional indicators of clay influence: $\mathrm{Al}/\mathrm{Ca}$ > 0.4, $\mathrm{Si}/\mathrm{Ca}$ > 0.4, $\mathrm{Ti}/\mathrm{Ca}$ > 0.005 mmol mol⁻¹, or very low $\mathrm{Mn}/\mathrm{Fe}$ (< 0.1 mol mol⁻¹) due to elevated $\mathrm{Fe}/\mathrm{Ca}$.

• Analytical blank effects: Low [Ca] combined with anomalously high $\mathrm{Mg}/\mathrm{Ca}$ values, suggesting contamination during measurement.

2.5 PSU Solver: Calibration equations for $\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperature and δ¹⁸O_seawater

Deep-water temperature and δ¹⁸O_seawater estimates were calculated for Site U1385 over the past 1.5 Ma using an updated Python version of Paleo-Seawater Uncertainty Solver (PSU Solver) originally developed in MATLAB (Thirumalai et al., 2016). Uncertainties associated with $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration slopes and intercepts were propagated explicitly within a Monte Carlo framework by sampling calibration slopes and intercepts within their reported uncertainties during each iteration (n= 1000). Input data comprise foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ and benthic δ¹⁸O measurements. Standard laboratory uncertainties of 0.08 ‰ for δ¹⁸O and 0.05 mmol mol⁻¹ for $\mathrm{Mg}/\mathrm{Ca}$ were applied as inputs to the PSU Solver (Thirumalai et al., 2016). Species-specific $\mathrm{Mg}/\mathrm{Ca}$–temperature calibrations were employed within the PSU Solver framework (Table 2).

For Uvigerina peregrina, we used the linear $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration of Elderfield et al. (2012):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{Mg}/\mathrm{Ca}\left[\mathrm{mmol}{\mathrm{mol}}^{-\mathrm{1}}\right]& =\mathrm{1.0}\\ & +\left(\mathrm{0.1}\pm \mathrm{0.013}\times T\left[\mathit{°}\mathrm{C}\right]\right)\end{array}\end{array}$$

Core-top $\mathrm{Mg}/\mathrm{Ca}$-derived temperatures (∼ 2.5 °C; Figs. 2 and S4) are consistent with modern deep-water temperatures (∼ 2–3 °C) at the Iberian Margin (Hodell et al., 2014; Skinner et al., 2003), supporting the applicability of the Elderfield et al. (2012) $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration at this site.

For Globobulimina affinis, we applied the calibration of Weldeab et al. (2016):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{Mg}/\mathrm{Ca}\left[\mathrm{mmol}{\mathrm{mol}}^{-\mathrm{1}}\right]=\mathrm{2.22}& \pm \mathrm{0.19}\\ +(\mathrm{0.36}& \pm \mathrm{0.02}\times T\left[\mathit{°}\mathrm{C}\right])\end{array}\end{array}$$

To reconstruct δ¹⁸O_seawater we used the oxygen-isotope palaeothermometry equation (Elderfield et al., 2010, 2012):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}T\left[\mathit{°}\mathrm{C}\right]=\mathrm{16.9}& -\mathrm{4.0}\times ({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{calcite}}\left[\mathrm{‰}\right]\\ & -{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{seawater}}\left[\mathrm{‰}\right]+\mathrm{0.27})\end{array}\end{array}$$

Age model uncertainty was set at ±2 kyr, which represents both the optimal resolution for combining the records of differing resolution and the mean resolution of the five records.

Table 2 Species-specific $\mathrm{Mg}/\mathrm{Ca}$–temperature calibrations.

Abbreviations: LDEO, Lamont-Doherty Earth Observatory (Columbia University, New York, USA); ETHZ, Geological Institute, ETH Zurich (Zurich, Switzerland).

2.6 Stacking temperature and δ¹⁸O_seawater records

The rationale for stacking records of temperature and δ¹⁸O_seawater is the same as that given for oxygen isotope records (Imbrie et al., 1984; Lisiecki and Raymo, 2005); that is, to improve the signal $/$ noise ratio and reduce the influence of local variability. Although stacked records provide more representative regional and/or global mean signals of temperature and ice volume they also tend to dampen the amplitude of the signal owing to misalignments and smoothing. Stacked deep-water temperature and δ¹⁸O_seawater records were constructed for the North Atlantic using Sites U1385 (this study; Birner et al., 2016), MD01-2444 (Skinner and Elderfield, 2007), MD99-2334 (Skinner et al., 2003) and Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009). The deep Pacific stack was derived using Sites 1123 (Elderfield et al., 2012) and 1208 (Ford and Raymo, 2020) (Figs. S8 and S9). Globally distributed mean deep ocean temperature (MDOT) and mean deep ocean δ¹⁸O_seawater stacks (Figs. 3a and 4a respectively) incorporate these records together with ODP Site 1094 located in the South Atlantic sector of the Southern Ocean (Hasenfratz et al., 2019) (Fig. S10).

Prior to stacking, temperature and δ¹⁸O_seawater estimates were interpolated to a regular 3 kyr time step using MATLAB's interp function. Uncertainty envelopes were estimated using bootstrap resampling (n= 1000), assuming normally distributed errors defined by the reported standard deviations (1σ) (Figs. 3 and 4). Gaps in the original age models were identified, and interpolated values were excluded where age gaps exceeded ∼ 6 kyr. At each interpolated time step, stacked estimates were calculated using either arithmetic or volumetrically weighted averaging of available basin-scale records. Uncertainties in non-weighted stacks were estimated by propagating the combined calibration and analytical errors in quadrature and incorporating inter-record variance, such that the uncertainty of the mean reflects both measurement uncertainty and the spread between paired estimates. For weighted stacks, basin weights were defined using fixed deep ocean volume fractions following the volumetric approach of Lisiecki and Stern (2016) (see their Table S2) and renormalised to unity at each time step to reflect the ocean volume represented by the available records. This ensures that stacked estimates reflect the relative deep ocean volume of the contributing basins without introducing bias when one or more basins are absent. Only half of the stacked data points (n= 268) include contributions from all three basins (North Atlantic, Southern Ocean and Pacific), together representing 43.6 % of the global ocean volume (Table S1). The impact of data gaps in the Site 1094 record appears to exert little influence on the weighted MDOT and mean δ¹⁸O_seawater reconstructions (Fig. S11).

Figure 3 Interpolated and bootstrapped $\mathrm{Mg}/\mathrm{Ca}$-derived temperature records for all sites. (a) Non-weighted mean deep ocean temperature stack (MDOT; this study, orange line and markers), with MIS numbers denoting glacial (below) and interglacial (above) stages. $\mathrm{Mg}/\mathrm{Ca}$-derived temperature records for: (b) DSDP Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) (light blue); (c) IODP Site U1385 (this study, including additional Iberian Margin records (Birner et al., 2016; Skinner et al., 2003; Skinner and Elderfield, 2007) (light green); (d) ODP Site 1094 (Hasenfratz et al., 2019) (yellow); (e) ODP Site 1123 (Elderfield et al., 2012) (light pink); and (f) ODP Site 1208 (Ford and Raymo, 2020) (dark purple). $\mathrm{Mg}/\mathrm{Ca}$-derived temperature records were interpolated on a 3 kyr interval and subsequently bootstrapped. Shading represents 1σ and 2σ error margins (dark and light grey respectively), propagated from both analytical and calibration sources.

Figure 4 Interpolated and bootstrapped δ¹⁸O_seawater records for all sites. (a) Non-weighted mean δ¹⁸O_seawater stack (this study; blue line and markers) with MIS numbers denoting glacial stages (below) and interglacial stages (above). Interpolated and bootstrapped δ¹⁸O_seawater records for: (b) DSDP Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) (blue); (c) IODP Site U1385 (this study, including additional Iberian Margin records (Birner et al., 2016; Skinner et al., 2003; Skinner and Elderfield, 2007) (green); (d) ODP Site 1094 (Hasenfratz et al., 2019) (orange); (e) ODP Site 1123 (Elderfield et al., 2012) (pink); and (f) ODP Site 1208 (Ford and Raymo, 2020) (purple). All δ¹⁸O_seawater records were interpolated to a 3 kyr interval and subsequently bootstrapped. Shading denotes 1σ and 2σ error margins (dark and light grey respectively), propagated from both analytical and calibration uncertainties.

At each time step t, volumetrically weighted stacked estimates were calculated as

$$\begin{array}{}\text{(5)}& \overline{x}\left(t\right)={\sum }_{i=\mathrm{1}}^{N\left(t\right)}{w}_i\left(t\right)x_i\left(t\right),\end{array}$$

where x_i(t) is the basin-scale estimate (temperature or δ¹⁸O_seawater) and N(t) is the number of available basins. Weights were defined as

$$\begin{array}{}\text{(6)}& w_i\left(t\right)={\displaystyle \frac{V_i}{{\sum }_{j=\mathrm{1}}^{N\left(t\right)}{V}_j}}\end{array}$$

where V_i denotes the deep ocean volume fraction represented by basin i, and the denominator represents the total volume of all basins available at time t.

Total uncertainty was estimated by propagating intra-basin variance and incorporating an inter-basin variance term that captures the spread of basin-scale estimates about the weighted mean:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\sigma }}^{\mathrm{2}}\left(t\right)& ={\sum }_{i=\mathrm{1}}^{N\left(t\right)}{w}_i^{\mathrm{2}}\left(t\right){\mathit{\sigma }}_i^{\mathrm{2}}\left(t\right)\\ & +{\sum }_{i=\mathrm{1}}^{N\left(t\right)}{w}_i\left(t\right){\left(x_i\left(t\right)-\overline{x}\left(t\right)\right)}^{\mathrm{2}},\end{array}\end{array}$$

where σ_i(t) represents the reported 1σ uncertainty of basin i, and the second term represents inter-basin variability. Where only a single basin contributed to the stack, the uncertainty reduces to the reported basin 1σ value. This approach assumes that residual chronological misalignment between records is small relative to the interpolated 3 kyr resolution. Uncertainty envelopes incorporate analytical, calibration and age uncertainties propagated through a Monte Carlo framework, with basin-scale uncertainty additionally reflecting inter-site variability.

Change-point analysis was applied to the interpolated stacks (3 kyr resolution) using MATLAB's findchangepts function to identify statistically significant shifts in the mean state of the record. The analysis detects a primary breakpoint at ∼ 930 ka in both mean averaged and basin-volume weighted δ¹⁸O_seawater stacks, marking the most prominent transition in the time series.

3 Results

3.1 Comparison of Sites U1385 and 607 deep-water records

North Atlantic Sites U1385 and 607 temperature and δ¹⁸O_seawater records generally show similar values and trends over the past 1.5 Myr (Figs. 5b and 6b), but there are some differences in the magnitude of the $\mathrm{Mg}/\mathrm{Ca}$-derived temperature changes. Greater cooling is expressed at Site 607, where temperatures show decreasing trends of ∼ 1.2 °C between 1500–900 ka and ∼ 1.6 °C over the past 1.5 Myr. In contrast, Site U1385 shows no statistically significant long-term cooling prior to 900 ka, and a decreasing trend of ∼ 0.8 °C across the full record. At times, differences between the two temperature records (e.g. between 1180–1150 ka and ∼ 400–350 ka) correspond to $\mathrm{Mg}/\mathrm{Ca}$ measurements of Cibicidoides wuellerstorfi at Site 607, which is more prone to carbonate ion changes and dissolution than infaunal benthic taxa. However, the cooler temperatures at Site 607 cannot be attributed solely to a carbonate ion effect because Uvigerina spp. also give high temperatures between ∼ 1400–1300 ka (Ford et al., 2016; Sosdian and Rosenthal, 2009) and may instead reflect hydrographic differences between the eastern and western North Atlantic basins (Chalk et al., 2019). Despite the subtle differences in the records, the close similarity in both temperature and δ¹⁸O_seawater data from the two sites supports their integration into a North Atlantic stack.

Figure 5 Comparison of $\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperature from Site U1385 versus Sites 607 and 1123. (a) Prob-stack benthic δ¹⁸O (Ahn et al., 2017) (black line), with MIS numbers denoting glacials (below) and interglacials (above). $\mathrm{Mg}/\mathrm{Ca}$-derived deep-water temperatures from Site U1385 (grey) are compared with (b) Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) (blue), and (c) Site 1123 (Elderfield et al., 2012) (pink). Site 607 exhibits a gradual cooling through time (Ford et al., 2016; Sosdian and Rosenthal, 2009), whereas Site 1123 documents relatively stable, near-freezing temperatures during glacial maxima (Elderfield et al., 2012). Dark and light shading represent 1σ and 2σ uncertainties, respectively, estimated by propagating analytical and calibration uncertainties in quadrature, and accounting for inter-record variance. Records are interpolated to 3 kyr resolution and bootstrapped.

Figure 6 Comparison of benthic δ¹⁸O_seawater records from Site U1385 with Sites 607 and 1123. (a) Prob-stack benthic δ¹⁸O (Ahn et al., 2017) (black line), with MIS numbers denoting glacial (even) and interglacial (odd) stages. Deep-water δ¹⁸O_seawater from Site U1385 (dark grey) is compared with (b) Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009) (dark blue), and (c) Site 1123 (Elderfield et al., 2012) (dark pink). Records are interpolated at a 3 kyr interval and bootstrapped. Dark and light shading indicates 1σ and 2σ uncertainty envelopes, respectively. No clear long-term trend in δ¹⁸O_seawater (ice volume/salinity) is observed at Site 607 (Ford et al., 2016; Sosdian and Rosenthal, 2009), whereas Site 1123 records an abrupt increase in Antarctic ice volume (salinity) ∼ 900 ka (Elderfield et al., 2012).

3.2 Deep Atlantic-Pacific gradients in deep-water temperature and δ¹⁸O_seawater

Whereas the benthic δ¹⁸O_calcite records are similar between the Atlantic and Pacific (Fig. 7a), the deconvolved temperature and δ¹⁸O_seawater records are not. Comparisons between Atlantic and Pacific deep-water temperature and δ¹⁸O_seawater stacks show greater inter-basin differences before ∼ 930 ka than afterwards. From 1500 to 930 ka (MIS 49 to 24), North Atlantic deep-water was on average ∼ 2.5 ± 1.5 °C warmer than the deep Pacific during both glacial and interglacial stages (Fig. 7b). After ∼ 930 ka, this temperature difference decreases to ∼ 1.1 ± 1.4 °C primarily reflecting cooling of glacial North Atlantic deep-water while deep Pacific temperatures show only minor changes. Prior to 930 ka, deep ocean glacial δ¹⁸O_seawater values in the Pacific were distinctly less than those of the Atlantic. Across the MPT, glacial deep ocean δ¹⁸O_seawater increased markedly in the Pacific and more modestly in the deep Atlantic (Fig. 7c). After 930 ka, glacial deep-water temperature and δ¹⁸O_seawater values from the Atlantic and Pacific were more similar.

Figure 7 Comparison of North Atlantic and Pacific Ocean deep-water temperature and benthic δ¹⁸O_seawater stacks. (a) Prob-stack benthic δ¹⁸O (Ahn et al., 2017) (black line), overlain with contributing benthic δ¹⁸O records for the Atlantic (this study; Ford et al., 2016; Sosdian and Rosenthal, 2009) (bright green) and Pacific (Elderfield et al., 2012; Ford and Raymo, 2020) (bright pink) stacks in (b) and (c). MIS numbers denote glacial (even) and interglacial (odd) stages. (b) Stacked deep-water temperatures for the Pacific (red) and North Atlantic (blue) with dark and light shading indicating 1σ and 2σ uncertainties respectively. (c) Stacked deep-water δ¹⁸O_seawater for the Pacific (purple) and North Atlantic (green) with uncertainties as in (b). Stacks and associated uncertainties represent mean of manually aligned, interpolated (3 kyr resolution) and bootstrapped records, with uncertainties reflecting both analytical error and the spread between paired estimates (see Sect. 2.6).

Deep-water temperatures at Site 1094 in the Atlantic sector of the Southern Ocean (Hasenfratz et al., 2019) are similar to those of the deep Pacific prior to ∼ 900 ka (Fig. S10), but are generally slightly colder thereafter consistent with the modern temperature difference between the two sites (Fig. S2; Table 1). No substantial change in glacial temperatures occurred at Site 1094 across the MPT (Fig. S10). The δ¹⁸O values of seawater at Site 1094 are more similar to those of the deep North Atlantic record. Prior to the MPT, a gradient in δ¹⁸O_seawater existed between Site 1094 and the deep Pacific, whereas the values converge after 930 ka.

All five records were stacked to produce global signals of deep ocean temperature (> 2500 m) and δ¹⁸O_seawater change. The non-weighted MDOT stack (Fig. S12c) shows a modest gradual decrease in glacial deep-water temperatures from ∼ 1050 to 925 ka primarily related to cooling in the deep North Atlantic. At ∼ 930 ka, the non-weighted mean δ¹⁸O_seawater stack (Fig. S12b) increases abruptly. This change is more pronounced in ocean volume-weighted records because of the proportionally greater influence of Pacific deep ocean temperature, which shows a warming from MIS 24 through 22 (Fig. 8). Mean glacial δ¹⁸O_seawater values in the ocean volume-weighted stack increase from 0.38 ‰ ± 0.08 ‰ prior to ∼ 930 ka to 0.80 ‰ ± 0.11 ‰ thereafter, representing a mean difference of ∼ 0.42 ‰ that is statistically significant (p< 0.001).

Figure 8 MDOT and mean δ¹⁸O_seawater stacks weighted by ocean-basin volume. (a) Basin-volume–weighted benthic δ¹⁸O stack compiled from Atlantic, Pacific and Southern Ocean records used in this study (pink) overlaid on Prob-stack benthic δ¹⁸O (Ahn et al., 2017) (black). MIS numbers denote glacial (even) and interglacial (odd) stages. (b) Globally distributed basin-volume–weighted mean deep-water δ¹⁸O_seawater stack (blue). Horizontal dashed lines indicate the glacial mean values for intervals before and after 930 ka. (c) Ocean volume–weighted MDOT stack (orange); the horizontal dashed lines indicate the total mean values for the intervals before and after 930 ka. The vertical dashed line marks 930 ka. Ocean-basin volume weighted stacks were constructed following the approach of Lisiecki and Stern (2016) (Sect. 2.6). Uncertainty envelopes incorporate analytical, calibration and age uncertainties propagated through a Monte Carlo framework, with basin-scale uncertainty additionally reflecting inter-record variability.

4 Discussion

4.1 Basinal differences in abyssal temperature and salinity

Although the deep Atlantic and Pacific stacks contain only two records each, we suggest they are broadly representative of deep-water mass properties in the two basins. Comparison of Atlantic and Pacific temperature and δ¹⁸O_seawater records indicate differences in the magnitude of temperature and δ¹⁸O_seawater changes across the MPT (Fig. 7b, c). The deep North Atlantic cooled during glacials between ∼ 1050 to 925 ka whereas the deep Pacific remained cold, both before and after the MPT. The magnitude of cooling in the deep North Atlantic is supported by a decrease (∼ 2 °C) in high-latitude North Atlantic sea surface temperature (SST) at Site 982 on the Rockall Plateau (58° N, 16° W; Lawrence et al., 2009) (Fig. 9b), and nearby Site 983 (McClymont et al., 2008), both located close to the source region of deep-water formation. Because deep water derived from the Southern Ocean is colder than that sourced from the North Atlantic, part of the cooling may also result from an expansion of southern sourced water into the deep North Atlantic during glacial stages. Cooling coincides with an abrupt increase in neodymium isotopes (ϵNd) at Site 607 that has been previously interpreted as indicating an increased contribution of Southern Component Water (SCW) in the deep Atlantic (Kim et al., 2021; Pena and Goldstein, 2014) (Fig. 9d), with similar deep circulation changes inferred for the Pacific (Li et al., 2026). An alternative interpretation is that neodymium isotope changes in the Atlantic reflect changes in the end-member Nd value of NCW alongside more modest circulation changes (Pöppelmeier et al., 2020, 2022; Williams et al., 2024; Zhao et al., 2019). Nevertheless, an increase in the δ¹³C gradient between the intermediate and deep North Atlantic also supports a greater proportion of SCW in the North Atlantic across the MPT (Fig. 9e; Lang et al., 2016; Lisiecki, 2014).

Figure 9 Comparison of deep-water temperature, salinity, ocean circulation and atmospheric pCO₂ records over the last 1.5 Ma. (a) Benthic δ¹⁸O Prob-stack (Ahn et al., 2017) (black line), MIS numbers denote glacial (below) and interglacial (above) stages. (b) Alkenone U${}_{\mathrm{37}}^{{\mathrm{k}}^{\prime }}$-derived, sea surface temperature (SST) record from ODP Site 982 (Lawrence et al., 2009) (green line with markers) overlaying the Atlantic stack for deep-water temperature (this study; blue line with markers). (c) Atlantic (green) and Pacific (purple) δ¹⁸O_seawater stacks (this study). (d) DSDP Site 607/V30-97 ϵNd proxy for North Atlantic Ocean circulation changes (Kim et al., 2021) (dark pink line with pink diamonds). (e) Carbon isotope gradient between regional δ¹³C stacks (Lisiecki, 2014) for the intermediate depth North Atlantic (black) and middle deep Atlantic (orange). Low δ¹³C reflects reduced ventilation in the deep Atlantic. (f) Atmospheric pCO₂ records (presented as in Thomas et al., 2022): Antarctic 800 kyr ice core data (Bereiter et al., 2015) (black line); Allan Hills Blue Ice (blue triangles (Yan et al., 2019) and blue squares (Higgins et al., 2015)); and δ¹¹B-based reconstructed pCO₂ [(Dyez et al., 2018) green line with diamonds; (Hönisch et al., 2009) solid green squares; (Chalk et al., 2017) early Pleistocene light green line/small circles, and late Pleistocene light green circles]. Chronology for Atlantic and Pacific stacks is available in Sect. 2.1, all other records are presented on their original timescales.

Cooling of Northern Component Water (NCW), assuming its salinity remained approximately constant, would increase its density and create a stratification inversion with deeper water masses in the absence of a density increase of SCW (Knorr et al., 2021). The δ¹⁸O_seawater in both the deep North Atlantic and Pacific increases across the MPT reflecting expansion of continental ice sheets, but the increase is considerably greater in the Pacific than the Atlantic (Fig. 7c). We attribute the greater increase in Pacific δ¹⁸O_seawater to an increase in the salinity and density of SCW, which may have shoaled the boundary between NCW and SCW in the deep Atlantic if the density contrast between the two water masses increased (Fig. 9c, e) (Song et al., 2025). Glacial-interglacial change in the global ocean density stratification in much of the deep ocean (> 2000 m) is controlled by surface freshwater forcing in the Southern Ocean (Sun et al., 2016). Increased salinity can be achieved by decreasing meltwater input to the marginal seas around Antarctica (e.g. Weddell and Ross Seas) either from the Antarctic ice sheet and/or grounded ice shelves. Sea ice expansion around Antarctica also raises salinity through brine rejection and expands the depth range of SCW (Ferrari et al., 2014), but sea ice formation does not significantly alter the δ¹⁸O_seawater of seawater (Kim and Timmermann, 2024; Toyota et al., 2013).

Raymo et al. (2006) suggested the East Antarctic ice sheet was terrestrial-based with melting margins before the MPT and transformed to a marine-based ice sheet afterward, which would have decreased meltwater to marginal seas. Hemispheric ΔSST reconstructions exhibit little precession-band variability (Clark et al., 2025b) and have been interpreted as evidence against this mechanism. However, precession is often underestimated by spectral analysis in relatively low-resolution stacked records with chronologic error and noise (Liautaud, 2021; Liautaud et al., 2020). We therefore include the mechanism proposed by Raymo et al. (2006) as one of several possible explanations for the increase in salinity of Antarctic marginal seas. An unconformity in the Ross Sea Embayment that spans the MPT has been interpreted to represent widespread expansion of a marine-based ice sheet (McKay et al., 2012) into the Ross Sea where AABW is formed today. Adkins (2013) proposed that cooling of NADW during glacial stages of the late Pleistocene decreased melting rates of floating Antarctic ice shelves. Adkins (2013) suggested decreased meltwater input led to the formation of a cold and salty version of Antarctic Bottom Water, and this water mass filled the deep ocean basins. We observe cooling of the high-latitude and deep North Atlantic across the MPT at the same time as the inferred salinity increase of the deep Pacific, which may support the mechanism outlined by Adkins (2013).

The freshwater fluxes to marginal seas can also be decreased by exporting freshwater to lower latitude via iceberg discharge and sea ice (Ferrari et al., 2014; Jansen, 2017). Starr et al. (2021) reported an increase in IRD in the Subantarctic Zone across the MPT and suggested this reflected increased transport of meteoric freshwater away from Antarctica. This process could have contributed to increased salinity in source areas of SCW formation, especially in the Weddell Sea sector. At the same time, increased transport of icebergs and sea ice from Antarctica increased the freshwater flux to the subantarctic, strengthened the halocline and slowed deep-to-surface exchange that controls carbon dioxide release to the atmosphere (Hasenfratz et al., 2019).

4.2 Implications for carbon storage

Although we cannot pinpoint which of the many processes discussed above was responsible for increased salinification and densification of SCW – some or all may have contributed – increased stratification of the glacial deep ocean may have played a role in enhancing ocean carbon storage and lowering atmospheric pCO₂ concentrations (Adkins, 2013). Increased deep ocean stratification provides a physical mechanism by which sinking organic carbon can be remineralised and accumulate in the deep sea. In the process, oxygen is consumed, nutrients are released, and total dissolved inorganic carbon (DIC) increases, although carbonate ion decreases. Several studies have reported an increase in deep-sea nutrients and lowered carbonate ion and oxygen concentrations in the deep Atlantic (e.g. Fig. 9e; Farmer et al., 2019; Lear et al., 2016; Lisiecki, 2014; Thomas et al., 2022) and Pacific (Diz et al., 2020; Peng et al., 2026; Qin et al., 2022) across the MPT.

Increased stratification of the deep ocean by itself is insufficient to increase deep-ocean carbon storage. Instead, enhanced storage requires a more efficient biological pump, whereby a greater proportion of nutrients in the deep ocean are regenerated through remineralisation of sinking organic matter rather than supplied as preformed nutrients during deep-water formation (Marinov et al., 2006, 2008; Sigman et al., 2010). The biological pump exports organic matter from the surface ocean that is remineralised at depth. In the Southern Ocean, biological export production is thought to have increased owing to iron fertilization from dust over the MPT (Martínez-Garcia et al., 2011). The return of deep water towards the surface in the upwelling region of the Antarctic Circumpolar Current (ACC) ventilates the deep ocean through gas exchange in the Southern Ocean. This surface/deep exchange of water was reduced (“isolated”) during glacial periods because of a strengthened halocline and reduced upwelling (Sigman et al., 2021). At Site 1094 across the MPT, Hasenfratz et al. (2019) suggested a reduction in deep water supply to the surface and a concomitant freshening of surface waters that strengthened the halocline, thereby reducing carbon dioxide release to the atmosphere. The strength of the Westerlies and ACC are also important for the upwelling of deep-water masses from the ocean interior in the Antarctic Zone, although changes may be asynchronous in different sectors of the Southern Ocean (Lamy et al., 2024; Sun et al., 2016).

Most hypotheses and models to explain the MPT invoke a glacial atmospheric pCO₂ decline to promote the growth of larger ice sheets; indeed, it is difficult to explain the MPT without this radiative cooling especially in the tropics (Herbert, 2023). Some studies have indicated a modest drop (∼ 20 ppm) in concentrations of glacial atmospheric pCO₂ across the MPT (Chalk et al., 2017; Higgins et al., 2015; Hönisch et al., 2009; Yan et al., 2019) (Fig. 9f), whereas others have found no mean change (Marks-Peterson et al., 2026). Our findings are also broadly consistent with recent modelling simulations suggesting that declining atmospheric pCO₂ and sea-level change may have triggered a threshold response of the Antarctic Ice Sheet, resulting in enhanced coupling between Antarctic ice-sheet variability, atmospheric pCO₂ and Northern Hemisphere glaciation following the MPT (Yun and Timmermann, 2026).

4.3 Mean Deep Ocean Temperature

A critical question for the MPT is by how much and how fast did global deep ocean temperature and ice volume change? Global stacking of mean deep ocean temperature and δ¹⁸O_seawater records (derived from paired benthic $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O) can address this question and be compared to other more indirect approaches of estimating temperature and ice volume changes from benthic δ¹⁸O records. Although only 5 long records of tandem benthic δ¹⁸O–$\mathrm{Mg}/\mathrm{Ca}$ measurements exist spanning the last 1.5 Myr, we stacked the records first with equal weighting to reconstruct global records and then weighted by the volumes of different ocean basins (Lisiecki and Stern, 2016) (Sect. 2.6) to better estimate global averages. Whereas the unweighted glacial MDOT record shows a slight decrease after ∼ 900 ka, it is biased by the cooling of the deep Atlantic at Site 607 (Fig. S12c). Considering the much smaller volume of the deep North Atlantic compared to the deep Pacific, the volume-weighted MDOT stack is more representative and shows minimal change across the MPT (Fig. 8c).

We compare our stacked volume-weighted MDOT record with similar estimates derived by other methods. Rohling et al. (2021) (R21) calculated deep-sea temperatures (DST) using benthic δ¹⁸O and estimates of δ¹⁸O_seawater from sea level reconstructions. Our MDOT reconstruction agrees well with R21 before ∼ 1200 ka, although after ∼ 1200 ka R21 exhibits larger-amplitude changes (Fig. 10a). Neither record shows substantial cooling of deep ocean temperatures across the MPT. Chandler and Langebroek (2024) similarly used benthic δ¹⁸O, $\mathrm{Mg}/\mathrm{Ca}$ and sea level reconstructions to estimate changes in mean Circumpolar Deep Water temperature (MCDWT) for the past 800 ka, including Sites 1123 and 1094 but using slightly different calibrations than this study (Fig. 10b). For the last 800 kyr, our MDOT agrees well with MCDWT, indicating the slightly different approaches yield similar results. Lastly, we compare our MDOT estimates with changes in mean ocean temperature (MOT) reconstructed using noble gas ($\mathrm{Xe}/\mathrm{Kr}$) from the EPICA ice core (Grimmer et al., 2025) and shallow ice cores recovered from the Allan Hills blue ice area, Antarctica (Shackleton et al., 2026) (Fig. 10c). The latter likely records a weighted averaging of glacial and interglacial conditions. There is good agreement between the two, especially considering that our MDOT does not include the entire ocean recorded by noble gas MOT. Both records show little change in MDOT or MOT across the MPT (Shackleton et al., 2026).

Figure 10 Comparison of MDOT with similar mean ocean and deep-water temperature reconstructions. All panels show the change in volume-weighted MDOT (ΔMDOT) (this study; black line) relative to an average present day value (1.73 °C) versus: (a) the change in deep-sea temperature (ΔDST) estimated relative to the present (Rohling et al., 2021) (purple line); (b) the change in mean CDW temperature (ΔMCDWT) relative to present (Chandler and Langebroek, 2024) (dark orange circles); (c) the change in MOT (ΔMOT) derived from noble gas measurements in ice cores (Grimmer et al., 2025) (purple and pink circles) and, in Allan Hills blue ice (pink diamonds) (Shackleton et al., 2026), error bars on individual MOT estimates represent absolute uncertainties; and, (d) ΔMOT inferred from ΔGMSST and proxy-based deep-water temperature (Clark et al., 2025a) (red), with both records referenced to the Early Holocene (0 °C at 10 ka) and offset by −1.73 °C following (Clark et al., 2025a). Dark and light shading of the same colour represent 1σ and 2σ uncertainty envelopes, respectively. Uncertainties associated with ΔMOT are taken directly from the published datasets of Clark et al. (2025a). MIS numbers mark many of the glacial (even) and interglacial (odd) stages.

The central estimates of both our MDOT and noble-gas MOT reconstructions (Shackleton et al., 2026) differ from the ΔGMSST-derived deep ocean temperature reconstruction of Clark et al. (2025a) prior to and across the MPT, particularly through MIS 22 (Fig. 10d). Although the uncertainty envelopes of these independent reconstructions overlap, they differ in the pattern of reconstructed deep-ocean temperature variability, including its amplitude and evolution through time. In particular, the Clark et al. (2025a) ΔGMSST-derived reconstruction exhibits relatively little glacial-interglacial variability in the earlier part of the record prior to MIS 26 (∼ 960 ka), reflecting differences in the approaches used to estimate MOT. The ΔGMSST estimates are converted to ΔMOT using a non-constant ratio of ΔMOT $/$ ΔGMSST over the MPT. For the period from 1.5 to ∼ 0.9 Ma, these estimates give warmer deep-water temperatures relative to our study, followed by a marked cooling in deep ocean temperature (> 2000 m) at ∼ 900 ka (Fig. 10d). The Clark et al. (2025a) compilation is strongly biased towards the deep North Atlantic (predominantly Site 607) with comparatively sparse Pacific deep ocean temperature (> 2000 m) data over this interval (Clark et al., 2025a; Lear et al., 2003). Site U1385 does not show cooling between 1.5 to ∼ 0.9 Ma, but the records from Site 607 show a ∼ 1.2 °C long term cooling trend, which is unlikely to represent a global mean change in deep-ocean temperature (Figs. 2b, 3b, c, 5b, S4a, and S8b, c; Sect. 3.1). Instead, our findings support previous interpretations for near-freezing glacial bottom water temperatures in the Pacific and the South Atlantic Southern Ocean throughout the period from 1.5 to ∼ 0.9 Ma, with no substantial change across the MPT (Elderfield et al., 2012; Ford and Raymo, 2020; Hasenfratz et al., 2019; Siddall et al., 2010) (Figs. 8c, S10c and S13c). These results appear to be consistent with recent ice core noble gas-based reconstructions of MOT (Shackleton et al., 2026) (Fig. 10c), but higher resolution MOT data are needed to capture the full amplitude of glacial-interglacial variations.

4.4 Mean δ¹⁸O_seawater and Ice Volume

Both weighted and unweighted mean ocean δ¹⁸O_seawater stacks were compiled and show a step-like increase across the MPT at 930 ka (Figs. 8b and S12b), supporting previous interpretations of increased continental ice volume during most glacial periods after ∼ 900 ka (Elderfield et al., 2012; Ford and Raymo, 2020). Our volume-weighted mean δ¹⁸O_seawater stack is consistent with process-based estimates of δ¹⁸O_seawater over the past 1.5 Myr (Rohling et al., 2021) (Fig. 11a) and with reconstructions of relative sea level (RSL) change for the past ∼ 500 ka (Grant et al., 2012, 2014) (Fig. 11b). However, our value for MIS 22 is considerably greater than R21. This is the result of a warming trend of MDOT from MIS 24 through MIS 22 that yields greater δ¹⁸O_seawater values when the signal is deconvolved from the benthic δ¹⁸O of calcite. While all records exhibit warming to varying degrees over this interval, the more prolonged warming in the Pacific records, particularly at Site 1208, has a greater influence on the basin-volume-weighted MDOT stack and therefore contributes to the discrepancy with R21.

Figure 11 Comparison of mean δ¹⁸O_seawater with similar process-based and relative sea level reconstructions. Volume-weighted mean δ¹⁸O_seawater stack shown in each panel (this study; black) is compared against: (a) process-based estimated ice-volume related deep-sea δ¹⁸O_seawater (Rohling et al., 2021) (light blue circles); (b) Red Sea relative sea level (RSL) reconstructions (Grant et al., 2012, 2014) (dark and light blue circles) scaled so that δ¹⁸O_seawater of 1.0 ‰ is equivalent to 10 m relative sea level; and, (c) δ¹⁸O_seawater derived from the deconvolution of the benthic δ¹⁸O Prob-stack using inferred changes in MOT (Clark et al., 2025a) (red). MIS numbers denote glacial (below) and interglacial (above) stages. Dark and light shading of the same colour represent 1σ and 2σ uncertainty envelopes, respectively. Uncertainties associated with δ¹⁸O_seawater are taken directly from the published datasets of Clark et al. (2025a).

Our results diverge from Clark et al. (2025a) in that our δ¹⁸O_seawater values are lower prior to the MPT, particularly in the interval older than 1200 kyr, suggesting less glacial ice volume before the MPT (Fig. 11c). Clark et al. (2025a) infer surprisingly large ice sheets before the MPT and relatively little change at ∼ 900 ka, which is contrary to our results and most previous studies (Bintanja and Van De Wal, 2008; Clark et al., 2006; Elderfield et al., 2012; Ford and Raymo, 2020; Rohling et al., 2021). The δ¹⁸O_seawater reconstruction of Clark et al. (2025a) incorporates the long-term benthic δ¹⁸O correction as originally proposed by Raymo et al. (2018) (0.083 ‰ Myr⁻¹), but our reconstruction does not. This correction amounts to a maximum adjustment of ∼ 0.12 ‰ at 1.5 Ma, which is small relative to the propagated uncertainties of the two δ¹⁸O_seawater reconstructions and would increase divergence between them. The method of reconstruction used by Clark et al. (2025a) also discounts the temperature and δ¹⁸O_seawater records from Pacific Sites 1123 and 1208 (Elderfield et al., 2012; Ford and Raymo, 2020) claiming they are unrepresentative of global values and reflect regional hydrographic changes. While we agree there is a hydrographic (salinity) component to the increase in the glacial deep Pacific δ¹⁸O_seawater signal that exceeds the Atlantic signal, there is also a global increase in δ¹⁸O_seawater reflecting greater ice volume across the MPT.

4.5 Limitations and Future Opportunities

The small number of long, continuous benthic $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O records available across the MPT leads to questions about how well they represent basin-wide or global averages. Low abundance of preferred benthic foraminiferal taxa in some of the records results in time gaps where fewer than 5 sites are available for stacking. The water depths of the studied sites are all below 2500 m b.s.l. and a greater depth range would provide a better comparison with noble gas-based MOT estimates across the MPT (Shackleton et al., 2026). The temporal resolution is relatively low (3 kyr) and may not reflect the full amplitude of changes and certainly do not capture millennial-scale variability. Direct assessment of test preservation (e.g. using the Foraminiferal Preservation Index; Poirier et al., 2021), together with complementary $\mathrm{B}/\mathrm{Ca}$ measurements on the same samples, would help evaluate the extent to which dissolution may have influenced $\mathrm{Mg}/\mathrm{Ca}$ palaeotemperature reconstructions from different benthic foraminiferal species. The propagation of error when using paired benthic $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O to calculate the δ¹⁸O_seawater is fairly large (±0.24 ‰ for U. peregrina and ±0.16 ‰ for G. affinis (1σ), with the dominant contribution arising from temperature uncertainty). Using benthic foraminifera taxa with a greater sensitivity of $\mathrm{Mg}/\mathrm{Ca}$ to temperature could potentially reduce the errors (Yang et al., 2025).

The recent recognition that the time rate of change of benthic δ¹⁸O is a proxy for Earth's Energy Imbalance (EEI) offers the opportunity to investigate the processes that modify the global energy budget on long timescales (Shackleton et al., 2023). With more records, the derivative of the deconvolved stacks would permit separation of the sensible and latent heat components and provide valuable insight into how the Earth's energy balance changed during the MPT.

All of the points discussed above underscore the importance of acquiring additional paired benthic $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O records to improve the records presented herein. Furthermore, the proposed causes for the changes in deep-ocean stratification and their impact on deep-ocean carbon storage should be tested with numerical models.

5 Conclusions

The causes of the MPT have been sought since it was first recognized in oxygen isotope records over 50 years ago (Pisias and Moore, 1981; Shackleton and Opdyke, 1976). Our results are consistent with the view that the MPT involved a fundamental change in the interactions among deep-ocean circulation, atmospheric pCO₂, and the cryosphere that permitted continental ice sheets to grow large enough to survive modest rises in summer insolation that would have resulted in deglaciation prior to the MPT (Ganopolski, 2024; Raymo et al., 1997; Willeit et al., 2019). Here we suggest that the MPT involved changes in the physical properties of the abyssal ocean, including a cooling of the deep North Atlantic and salinification of the deep water in the Pacific. This resulted in increased density stratification of the deep ocean in glacial periods after ∼ 900 ka.

Global ocean density stratification in much of the deep ocean (> 2000 m) is controlled by surface freshwater forcing in the Southern Ocean (Sun et al., 2016). We suggest that across the MPT, meltwater input to the marginal seas around Antarctica decreased and sea ice increased during glacial periods, resulting in formation of cold (near freezing), salty deep water. Increased deep ocean stratification provides a physical mechanism by which sinking organic carbon from the surface can accumulate in the deep sea, especially if accompanied by increased export production and surface stratification in the Southern Ocean. Many complementary studies support increased carbon storage in the deep Atlantic (e.g. Fig. 9e; Farmer et al., 2019; Lear et al., 2016; Lisiecki, 2014; Thomas et al., 2022) and Pacific (Diz et al., 2020; Peng et al., 2026; Qin et al., 2022) across the MPT, and some atmospheric pCO₂ proxy records show a moderate decline in glacial atmospheric concentrations (Chalk et al., 2017; Higgins et al., 2015; Hönisch et al., 2009; Yan et al., 2019). Forthcoming analyses of the newly recovered continuous Antarctic ice core spanning the MPT (Wolff et al., 2022) will provide critical information on changes in greenhouse gases and MOT, which are needed to test many of the inferences made herein.