Climate Evolution Through the Onset and Intensification of Northern Hemisphere Glaciation

The Pliocene Epoch (∼5.3–2.6 million years ago, Ma) was characterized by a warmer than present climate with smaller Northern Hemisphere ice sheets, and offers an example of a climate system in long‐term equilibrium with current or predicted near‐future atmospheric CO2 concentrations (pCO2). A long‐term trend of ice‐sheet expansion led to more pronounced glacial (cold) stages by the end of the Pliocene (∼2.6 Ma), known as the “intensification of Northern Hemisphere Glaciation” (iNHG). We assessed the spatial and temporal variability of ocean temperatures and ice‐volume indicators through the late Pliocene and early Pleistocene (from 3.3 to 2.4 Ma) to determine the character of this climate transition. We identified asynchronous shifts in long‐term means and the pacing and amplitude of shorter‐term climate variability, between regions and between climate proxies. Early changes in Antarctic glaciation and Southern Hemisphere ocean properties occurred even during the mid‐Piacenzian warm period (∼3.264–3.025 Ma) which has been used as an analog for future warming. Increased climate variability subsequently developed alongside signatures of larger Northern Hemisphere ice sheets (iNHG). Yet, some regions of the ocean felt no impact of iNHG, particularly in lower latitudes. Our analysis has demonstrated the complex, non‐uniform and globally asynchronous nature of climate changes associated with the iNHG. Shifting ocean gateways and ocean circulation changes may have pre‐conditioned the later evolution of ice sheets with falling atmospheric pCO2. Further development of high‐resolution, multi‐proxy reconstructions of climate is required so that the full potential of the rich and detailed geological records can be realized.


Introduction
Over the last ∼50 million years (Myr), there has been a long-term cooling trend in the Earth's climate, culminating in the transition to an "icehouse" climate state of bipolar and high amplitude glaciations which began in the Pliocene Epoch (5.3-2.58 million years ago, Ma) and was fully established across the Pliocene-Pleistocene transition (∼2.58 Ma) (Figure 1) (Westerhold et al., 2020).The Late Pliocene thus offers us an important opportunity to study the characteristics and pacing of major global climate change, and allows us to consider whether regions or parts of the climate system are vulnerable or resilient to climate forcings and feedbacks.
Late  climate was characterized by a sustained period (∼200,000 years) of global warmth, the mid-Piacenzian warm period (mPWP,.The mPWP has been an important window for collecting information about Earth's climate response to atmospheric carbon dioxide (CO 2 ) concentrations comparable to today and our near-future (∼400 ppmv; Figure 1) (Gulev et al., 2021).The Late Pliocene is also marked by the development of more intense cold stages (glacials) and an increase in ice-rafted debris (IRD) reaching the North Atlantic and North Pacific Oceans (Figure 1) (Blake-Mizen et al., 2019;Jansen et al., 2000;Kleiven et al., 2002;Mudelsee & Raymo, 2005).These changes mark the final step in the overall trend toward a bipolar-icehouse climate state (Bailey et al., 2013;Westerhold et al., 2020), and have been variously referred to as the "onset" or the "intensification of Northern Hemisphere Glaciation" (oNHG or iNHG).However, the terms "oNHG" and "iNHG" have been inconsistently used in the literature.A wide window of time in which a gradual expansion of the ice sheets developed, starting as early as ∼3.6 Ma (e.g.Mudelsee & Raymo, 2005), has been referred to as iNHG (e.g.Bolton et al., 2010Bolton et al., , 2018;;Kleiven et al., 2002).In contrast, iNHG has also described a narrower time window aligned with the stepped increase in IRD recorded in the high northern latitudes at ∼2.7 Ma (e.g.Bailey et al., 2011;Maslin et al., 1998;Ravelo et al., 2004;Teschner et al., 2016).Others have referred to the changes at ∼2.7 Ma as the "onset of major Northern Hemisphere Glaciation" (oNHG), to recognize that whilst ice was present and developing during the Pliocene (cf.Mudelsee & Raymo, 2005) the development of pronounced glacial cycles and larger continental ice sheets was a later significant climate shift (e.g.Bailey et al., 2013;Haug et al., 2005).Here, we make a distinction between (a) the longer-term transition toward bipolar glaciations, potentially beginning as early as ∼3.6 Ma and including the mPWP (which we refer to as the onset of NHG, oNHG), and (b) the relatively abrupt and later expansion of Northern Hemisphere ice sheets which occurs from ∼2.7 Ma (which we refer to as the intensification of NHG, iNHG).
One hypothesis to explain the iNHG is that the climate system crossed a threshold whereby a fall in CO 2 allowed large-scale glaciation in the Northern Hemisphere (Figure 1) (De Conto et al., 2008;Lunt et al., 2008).However, both the oNHG and iNHG occurred in the context of long-term shifts to meridional and zonal temperature gradients (e.g.Fedorov et al., 2015;Kaboth-Bahr & Mudelsee, 2022) as well as tectonic changes which can influence heat and moisture transport through the climate system, with potential impacts on ice-sheet growth (e.g.De Vleeschouwer et al., 2022;Karas et al., 2009;Knies, Mattingsdal, et al., 2014;Sánchez-Montes et al., 2020).It is important to determine the causes of the oNHG and iNHG, because they have been shown to affect the global climate system response to changes in the Earth's orbit and axial tilt ("orbital forcing"; Meyers & Hinnov, 2010;Turner, 2014;Westerhold et al., 2020).Orbital forcing impacts the distribution of energy received by the Earth from the sun through time and space, and has been shown to pace climate over tens and hundreds of thousands of years (e.g.Hays et al., 1976).A growing signal of global ice volume variability which could be linked to orbital forcing has been proposed to represent the development of a "deterministic" climate system associated with iNHG, whereby orbital forcing plays an important role in determining the climate system change (Meyers & Hinnov, 2010).In contrast, before iNHG the climate system is argued to be less strongly linked to orbital forcing, with a more stochastic (random) nature (Meyers & Hinnov, 2010).However, judging climate system response to forcing using only globally-integrated indicators, such as global ice volume, limits our assessment of the regional or local processes and feedbacks which may be facilitating the changes to external forcing as well as the underlying growth of the continental ice sheets.
We show that different regions of the climate system changed at different times, with some changing before the ice sheets expanded.The development of larger ice sheets in the Northern Hemisphere then impacted ocean temperatures and circulation, but there were many regions where no impacts were felt.Our analysis highlights regional differences in the timing and amplitudes of change within a globally-significant climate transition as well as in response to the current atmospheric CO 2 concentrations in our climate system.to the Pleistocene Epoch (2.58-0.01Ma) and the three time intervals used for our statistical analysis (P-1 3.3-3.0Ma, P-2 3.0-2.7 Ma, and P-3 2.7-2.4Ma; see Section 2.3).The time intervals defined as the onset of Northern Hemisphere Glaciation (oNHG), mid-Piacenzian warm period (mPWP) and the intensification of major Northern Hemisphere Glaciation (iNHG) are marked on panel (b), alongside the M2 glacial stage which has been suggested as a possible harbinger of more extensive glaciation (Westerhold et al., 2020).(a) Reconstructed atmospheric pCO 2 concentrations (de la Vega et al., 2020) including the ∼280 ppmv threshold for extensive glaciation proposed by De Conto et al. (2008) which aligns with Pre-Industrial pCO 2 recorded in ice cores (Loulerge et al., 2008); (b) a global stack of benthic oxygen isotope ratios (Ahn et al., 2017) which reflect a combined signal of the deep-water temperature and ice volume (see Technical Box).Key marine isotope stages referred to in the text are marked, blue for the cold glacial stages, red for the warm interglacial stage KM5c; (c) North Atlantic ice-rafted debris (IRD) inputs at Sites U1307 (southern Greenland margin) and U1313 (southern North Atlantic) (Blake-Mizen et al., 2019;Lang et al., 2014), which may be used to indicate when ice sheets expanded and reached the sea.(d) Antarctic IRD inputs from Site U1361 (Patterson et al., 2014).MAR refers to the Mass Accumulation Rate of the sediment.
signal is recorded, including whether the signal is generated by photosynthesizers or grazers, the preferred water depth or season of production, and the preservation of the SST signal during transport through the water column and into the seafloor sediments.By understanding both the similarities and differences between SST proxy records, we may identify the processes influencing our temperature signals and gain a more detailed view of environmental change associated with the oNHG and iNHG.
The information from Q1 is then combined with evidence for changes in global ice volume, including IRD, and other records of climate change to address three key questions linked to climate forcing and responses associated with oNHG and iNHG: (Q2) What were the characteristics of the mPWP interval?The mPWP is the most recent interval of geological time where pCO 2 and global temperatures were sustained above Pre-Industrial (by ∼50-110 ppmv for pCO 2 , by ∼2.5-4.0°C for temperature) (de la Vega et al., 2020;Gulev et al., 2021).For this reason, the mPWP has long been a target for data synthesis and data-model comparison efforts (e.g., Dowsett et al., 2012Dowsett et al., , 2016;;Haywood et al., 2020;McClymont et al., 2020;Salzmann et al., 2013), because it offers an opportunity to compare a warmer-than-modern climate with both modern observations and near-future projections (Burke et al., 2018;Gulev et al., 2021;Tierney et al., 2020).However, climate models of two interglacials within the mPWP show the potential for pronounced regional and seasonal sensitivity to variations to Earth's orbital configurations (Prescott et al., 2014), highlighting the importance of considering regional expressions of mPWP climates to better understand forcings and feedbacks on orbital timescales.(Q3) What were the characteristics of the iNHG?The mPWP and iNHG mark a shift in the patterns of climate change over the last ∼50 Ma, including a stronger global influence of complex high-latitude climate dynamics (Meyers & Hinnov, 2010;Turner, 2014;Westerhold et al., 2020).However, these analyses are focused on the integrated records of global ice volume and deep water temperatures recovered from stacks of benthic foraminifera stable oxygen isotope ratios (δ 18 O benthic ; Figure 1).Assessing a combination of globally-distributed surface ocean temperature and individual δ 18 O benthic records offers opportunities to identify and explain long-term trends in mean climate state as well as climate variability (see Question 4) for specific regions and/or circulation systems.(Q4) Did the amplitude or pacing of climate variability change with the iNHG?The climate system during the Pliocene and early Pleistocene was regulated on orbital timescales (10 3 -10 5 kyr) (Lisiecki & Raymo, 2005;Shackleton et al., 1984).However, the dominance of orbital obliquity periods (linked to Earth's axial tilt, ∼41 kyr) and the absence of a strong orbital precession signal (∼19-23 kyr) is puzzling, as daily summer insolation intensity has strongly influenced changes in polar ice volume for the last ∼800 kyr and is paced by precession (Hays et al., 1976).By exploring regional variability in ocean temperature and individual δ 18 O benthic records we can explore whether global ice volume is key to the pacing of orbital-scale climate variability, and how that might have changed with the iNHG.
Here, we draw on the high temporal resolution, multi-proxy marine sediment data synthesis which was generated and evaluated by the PAGES-PlioVAR working group (Figure 2).In previous work we have focused on a single interglacial within the mPWP (KM5c, ∼3.205 Ma; Figure 1) where well-constrained pCO 2 reconstructions allowed a detailed examination of the ocean temperature response to pCO 2 (McClymont et al., 2020).Here we apply the same stratigraphic and data evaluation protocols to the multi-proxy data over a longer time interval, to address the four questions outlined above.
We recognize that aspects of our data evaluation and analyses warrant technical explanations and justifications for a broad geophysics audience.In Section 2 we therefore first outline our approaches to ensuring sites meet our high-quality age control (Section 2.1) and briefly introduce both the climate records we have used (Section 2.2) and the statistical methods we employed to analyze our data compilation (Section 2.3).We direct readers who are interested in the details of the principles which underpin the climate reconstruction techniques we used to the Technical Box.The climate proxies we employ are indirect measures of climate variables, and through multi-proxy analysis we endeavor to tease out the most important signatures of past climate change.In particular, we employ both organic and inorganic proxies of ocean temperature change, which have a range of biological and environmental controls over their signatures: in Section 3 we address Q1 outlined above, by using three regions where we have the highest density of multi-proxy data to explore the differences and similarities expressed by these different approaches.In Section 4 we present our results examining changes in ocean temperature, ice volume and sea level through the late Pliocene to early Pleistocene (3.3-2.4 Ma).We first assess long-term 10.1029/2022RG000793 5 of 43 changes (Section 4.1) and then examine glacial-interglacial variability (Section 4.2) using the "PlioVAR data sets" which met our data quality thresholds as outlined in Section 2.1.We then address Q2 (Section 4.3), Q3 (Section 4.4), and Q4 (Section 4.5), by synthesizing our findings in the context of published information on changes to ice sheets and regional climates.1), which is a ratio describing the relative distribution of specific organic molecules (the C 37 alkenones) synthesized by selected haptophyte algae (Marlowe et al., 1984;Prahl & Wakeham, 1987).Many Pliocene studies have calibrated  U K ′ 37 to mean annual SST using the linear calibrations of seafloor surface sediments (core-tops) (Müller et al., 1998) or laboratory cultures (Prahl & Wakeham, 1987), which are almost indistinguishable.However, for some regions, these calibrations do not align with recorded mean annual SST.Non-linear  U K ′ 37 -temperature relationships are observed for high (i.e., warmest)  U K ′ 37 values (Conte et al., 2006;Pelejero & Calvo, 2003;Tierney & Tingley, 2018).Seasonality also likely impacts  U K ′ 37 reconstructions in the high latitudes of the North Pacific and North Atlantic Oceans (Tierney & Tingley, 2018).In our previous work on the interglacial KM5c (Figure 1), we showed that there was data-model agreement for both seasonal and mean annual SST reconstructions in the high latitudes, regardless of the calibration (McClymont et al., 2020).In the low-latitudes, the application of the non-linear BAYSPLINE calibration (Tierney & Tingley, 2018) elevated tropical SSTs relative to the Müller et al. (1998) core-top calibration (McClymont et al., 2020).Here, we present  U K ′ 37 -temperatures calibrated using BAYSPLINE (Tierney & Tingley, 2018) to acknowledge the non-linearity of the  U K ′ 37 -SST relationship at high values which is not captured by the Müller et al. (1998) calibration.However, we note that in doing so, there is increased uncertainty on the reconstructed SSTs above 24°C (from ∼1.5 to ∼4.4°C; Tierney & Tingley, 2018).

TEX 86 temperature proxy: principles and interpretations.
Only one published record (Site U1463) had TEX 86 temperature data meeting our age resolution criteria for analysis (Figure 1, Table 1).However, we draw on lower resolution TEX 86 records as secondary evidence in evaluating potential controls over other ocean temperature proxies for 4 other sites (Section 3).
The TEX 86 index is based on the distribution of isoprenoidal glycerol dialkyl glycerol tetraethers (isoG-DGTs) in marine sediments (Schouten et al., 2002), and has been correlated to surface or subsurface temperatures, where subsurface can be tens or a few hundred meters below the sea surface.Both linear (Schouten et al., 2002;Tierney & Tingley, 2014, 2015) or non-linear (Kim et al., 2010) calibrations have been used (see Inglis & Tierney, 2020, for an extensive review).Both calibration approaches yield similar values within the temperature range ∼5-30°C.TEX 86 has been particularly useful where  U K ′ 37 values have reached their upper limit in warm waters (i.e.,  U K ′ 37 = 1) (L.Li et al., 2011;O'Brien et al., 2014;van der Weijst et al., 2022;Y. G. Zhang et al., 2014).As for  U K ′ 37 , knowledge of past seawater chemistry is not required to calculate temperature values.However, care is required to ensure that data are not biased by contributions from archaea other than marine Thaumarchaeota, including terrestrial (Weijers et al., 2006), methanogenic (Inglis et al., 2015) and methanotrophic inputs (Y.G. Zhang et al., 2011).It has also been shown that changes to thermocline depth can influence TEX 86 values (Ford et al., 2015) (Section 3).We have screened all of the TEX 86 data included here using established indices for non-Thaumarchaeota inputs (e.g.BIT index, Ring index, Methane Index; Hopmans et al., 2004;Y. G. Zhang et al., 2011Y. G. Zhang et al., , 2016)).We also evaluated the potential contribution from "deep-water" Thaumarchaeota (those living below the permanent pycnocline) using the GDGT-2/GDGT-3 ratios (following Rattanasriampaipong et al., 2022;K. W. Taylor et al., 2013).We then applied a spatially-varying linear Bayesian regression model (BAYSPAR, Tierney & Tingley, 2014) to calculate TEX 86 -temperature estimates assuming a surface ocean origin (0 m water depth) for the isoGDGTs.

Foraminifera Mg/Ca temperature proxy: principles and interpretations.
The Mg/Ca-temperature proxy is underpinned by the observation that Ca is preferentially substituted by Mg in the calcite crystal lattice of foraminifera shells at higher temperatures (e.g.Anand et al., 2003;Dekens et al., 2002;Lea et al., 1999;Nürnberg et al., 2000).The Mg/Ca-temperature relationship is exponential, such that the proxy is more sensitive at higher temperatures.Planktonic foraminifera species commonly used for surface water temperature reconstruction include Globigerina bulloides, Globigerinoides ruber, and Trilobatus sacculifer (renamed from Globigerinoides sacculifer, Spezzaferri et al., 2015).However, depending on the oceanographic setting, there may be species-specific preferences in the preferred water depth of their habitat or the seasonality of their productivity, which can influence the recorded temperature signal.For example, the tropical and subtropical species G. ruber and T. sacculifer have photosynthetic symbionts and are confined to the photic zone, with T. sacculifer occupying a slightly deeper depth habitat than G. ruber (e.g.Bé, 1980;Fairbanks et al., 1982;Rebotim et al., 2017).Species-specific Mg/Ca-temperature calibrations also highlight species-specific differences in the substitution of Mg (Anand et al., 2003;Elderfield & Ganssen, 2000) and its distribution within foraminifera tests (Anand & Elderfield, 2005;Brown & Elderfield, 1996;Spero et al., 2015).Calcification temperature is the dominant control on planktic foraminiferal Mg/Ca (e.g.Anand et al., 2003;Dekens et al., 2002;Lea et al., 1999;Nürnberg et al., 2000), despite secondary controls over Mg incorporation into foraminifera tests.These secondary controls may be corrected as part of Mg/Ca-temperature calculations, and include (a) salinity, (b) surface water pH, (c) the Mg/Ca ratio of seawater, and (d) partial dissolution of calcite at the seafloor by sediment porewaters (Dekens et al., 2002;Evans et al., 2016;Regenberg et al., 2006;Rosenthal et al., 2022).
Six foraminifera Mg/Ca data sets met our criteria for analysis, from G. bulloides, G. ruber, and T. sacculifer (Table 1).The PlioVAR working group has previously presented the originally-published SSTs from these data sets for the KM5c interval (McClymont et al., 2020).The KM5c synthesis thus included temperature values generated across a range of calibrations and corrections, which were then compared to new PlioVAR calculations using the Bayesian regression model "BAYMAG" (Tierney et al., 2019).
Here, we calculate ocean temperature from Mg/Ca data sets where continuous data had been generated from a single species across our target interval (3.3-2.4Ma); we use an independent approach that maximized the fit between available core-tops and modern ocean temperatures (see Section 3 for results), while also maintaining consistency in our approaches to considering the secondary impacts on Mg/Ca incorporation noted in the previous paragraph, between sites and among species.We first accounted for the potential impact of dissolution at sites where modern-day carbonate saturation state is lower than 21.3 μmol kg −1 (identified by Regenberg et al. (2014) as the critical value below which dissolution starts).To do this, we identified the modern-day saturation state using values stated in the original publications (Table 1), or using World Ocean Atlas data and the CO2calc software (Robbins et al., 2010).Second, where reductive cleaning was included in the preparation method, 10% was added to the original Mg/Ca values before calibration, to enable comparison with non-reductively cleaned samples (Barker et al., 2003;Khider et al., 2015;Rosenthal et al., 2004); Sites ODP 806 and 1143).Next, we converted all corrected Mg/Ca data to temperature using the Elderfield and Ganssen (2000) calibration (G.bulloides, Site U1313) and the Anand et al. (2003) calibrations (all other sites/species), which yielded the best fit to core-top values (Section 3).
Finally, we account for the effect of changes in Mg/Ca seawater on foraminiferal Mg/Ca.Shifts in Mg/Ca seawater operate over million-year timescales, and thus should not affect the relative variability which is our focus here.It has nevertheless been shown that there can be a non-linear impact of Mg/Ca seawater on the calculated absolute SSTs (and thus the ranges) which would impact our assessment of glacial-interglacial variability (Evans et al., 2016).Although the magnitude and method of Mg/Ca seawater correction remains uncertain (e.g.Rosenthal et al., 2022;White & Ravelo, 2020), comparison of Mg/Ca-based and clumped isotope-based temperatures in the Pliocene show that a modest correction is best supported (Meinicke et al., 2021).Furthermore, since the application of a variety of Mg/Ca-based calibrations gave temperature estimates consistent within the calibrations' uncertainty (<±1°C) (Rosenthal et al., 2022), the Mg/Ca seawater calibration choice has minimal impact on the timescales we are targeting here.We therefore used the Evans et al. (2016) Mg/Ca seawater -sensitive temperature calibration, in conjunction with the BAYMAG seawater reconstruction (Tierney et al., 2019) to determine the magnitude of the Mg/Ca seawater correction which we then applied to our records.T. sacculifer Mg/Ca values were corrected upward by 10.3% to align with G. ruber before calculating the Mg/Ca seawater correction, as per Evans et al. (2016).The Mg/Ca seawater correction (about +0.5°C) was then added to each calculated sample temperature.Overall, our approach yields more reasonable core-top SSTs than when only using the Evans et al. (2016) calibration, the latter giving unrealistically cold core-top temperatures if interpreted as reflecting SST.A summary of the original published Mg/Ca data and the impact of our PlioVAR corrections is shown in Figure S2 in Supporting Information S1.

Foraminifera δ 18 O records for surface and deep ocean properties.
Seven of the PlioVAR synthesis sites provided planktonic foraminiferal δ 18 O records (δ 18 O planktonic ; Table 1).These δ 18 O planktonic records provide information about surface and near-surface temperatures and the δ 18 O composition of the seawater in which those species were living (Epstein et al., 1953).Surface seawater δ 18 O reflects the local precipitation-evaporation balance, local freshwater inputs like river runoff, and changes in global ice volume (Craig & Gordon, 1965;Shackleton & Opdyke, 1973).As for foraminifera Mg/Ca, biological and environmental factors influence foraminifera δ 18 O, including seasonality and depth habitat of the foraminifera, seawater pH, dissolution and post-depositional diagenesis (Pearson, 2012;Raymo et al., 2018).
In principle, the deep-sea δ 18 O benthic record is dependent on the δ 18 O composition of seawater and deepsea temperature (Shackleton & Opdyke, 1973).However, in contrast to the surface (planktonic) signals, the δ 18 O benthic record is considered to have a strong influence of global ice volume and a more minor influence of local temperature and salinity changes, as deep waters are expected to be more uniform through time and space (Ahn et al., 2017;Lisiecki & Raymo, 2005;Rohling, 2013;Waelbroeck et al., 2002).When quantifying past ice volume change, the temperature component of the δ 18 O benthic record is assumed to have the largest uncertainty (e.g.Evans et al., 2016;Raymo et al., 2018).Here, we compare δ 18 O benthic

Materials and Methods
To address the questions posed in Section 1, the PAGES-PlioVAR group considered the window of interest for this study to span from 3.3 to 2.4 Ma.This definition was used to ensure our target window included the M2 glacial stage and the mPWP which immediately followed it, as well as the onset and intensification of glacial-interglacial cycles in δ 18 O benthic records which align with iNHG (Figure 1) (McClymont et al., 2017).
The majority of the sites we examined were recovered by the Deep Sea Drilling Project (DSDP), the Ocean Drilling Program (ODP), the Integrated Ocean Drilling Program and the International Ocean Discovery Program (both IODP) (Figure 2, Table 1).The specific location at which each sediment sequence was recovered is termed the Site, and it is IODP nomenclature to name the sequence by a combination of the drilling program and the site number.We use this nomenclature in our text, tables and graphics (e.g.ODP Site 907 refers to Site 907 recovered by the ODP).

Age Models
Marine sediment sequences can be placed on an age scale (a "chronology") using a variety of methods.These include identifying events of known age, such as reversals in Earth's magnetic field and the first-appearance or extinction events in the fossil record (e.g., Gradstein et al., 2012).However, the intervals between these events can be hundreds of thousands of years apart.To explore climate changes on orbital timescales (see Section 1), we take advantage of the global impact of ice-sheet growth on the stable oxygen isotope ratio of seawater.The preferential storage of the 16 O isotope in continental ice during glacial stages leads to an increase in seawater δ 18 O, which is then recorded in δ 18 O benthic (see also the Technical Box).Although individual δ 18 O benthic records will have some local influences, the approach of "stacking" (or averaging) across multiple sites has been shown to increase the signal-to-noise ratio of δ 18 O benthic , and in turn provide a reference framework against which other sites can be aligned and then dated (Ahn et al., 2017;Lisiecki & Raymo, 2005).Regional differences in the expression of δ 18 O benthic and its relationship to orbital forcing have been determined (e.g.Caballero-Gill et al., 2019;Tian et al., 2002;Wilkens et al., 2017) which may hamper the accuracy of δ 18 O benthic -derived stratigraphy.However, for records from 22 sites which meet the PlioVAR criteria (Section 2.1), which span a range of regions and water depths to explore local/regional and global expressions of climate change.We did not include δ 18 O benthic records where more than one species had been used to generate the time series, so that we avoided any potential artificial introduction of shifts in the data when switching between species.Data are from Cibicides spp., Cibicides wuellerstorfi, Uvigerina spp., and Cibicides mundulus (Table 2).
To better understand the link between changes in ocean temperature and global ice volume, we examine indirect evidence for sea-level variability in the mPWP and early Pleistocene.Global sea-level variability is a measure of the magnitude and frequency of past ice volume fluctuations; direct records (e.g.paleo shorelines, speleothems, shallow marine sequences) are, however, spatially and temporally sparse (e.g.Miller et al., 2012 and references therein;Grant et al., 2019;Rovere et al., 2020).Age uncertainties associated with discontinuous records also makes it challenging to attribute sea-level maxima to single interglacials (e.g.Dumitru et al., 2019).To systematically assess sea-level variability through the mPWP and iNHG, we draw on the continuous but indirect sea-level estimates based on geochemical proxies from marine sediments (Miller et al., 2020;Rohling et al., 2021) and inverse modeling of δ 18 O benthic (Berends et al., 2021).Uncertainties come from the relative influence of temperature and diagenesis on δ 18 O noted above, as well as the unconstrained δ 18 O composition of polar ice sheets (Gasson et al., 2014) and varying ocean water masses.Here, we analyze three sea-level data sets generated using δ 18 O benthic , which have been calibrated to sea level by accounting for the influence of bottom water temperature variability through time (Miller et al. (2020), hereafter "Miller2020"), by accounting for non-linear changes to ice sheet δ 18 O during ice-sheet growth (Rohling et al. (2014), hereafter "Rohling2014"), or by undertaking an inverse modeling of temperature and ice volume contributions to δ 18 O (Berends et al., 2021, hereafter "Berends2021") Where the age model is noted as "PlioVAR" the record has been updated after the original publication.Basin refers to North Atlantic (natl), South Atlantic (satl), Pacific (pac), Indian (ind) and Southern Oceans (so).Table 1 Continued our interval of study, independent high-resolution δ 18 O benthic stratigraphies from the equatorial Atlantic confirm the LR04 age assignments (Wilkens et al., 2017).
Following protocols established by PAGES-PlioVAR (McClymont et al., 2017, 2020) (Gradstein et al., 2012).The age controls for all 32 sites which met these thresholds were reviewed, and the age models for several sites were revised compared to their published chronologies (Table 1 and Figure S1 in Supporting Information S1).Two sites were nevertheless excluded because although they met the age model requirements outlined here, the δ 18 O benthic data were unable to discern glacial-interglacial cycling as identified in the global stacks (n = 2; Figure S1 in Supporting Information S1).Fifty-three climate records met our criteria outlined above (the "PlioVAR data sets"), which include 5 climate proxies across 32 sites (Figure 1, Table 1).

Marine Proxies for Ocean Temperature, Ice Volume and Sea Level
A wide range of principles underpins the methods we employ here to reconstruct climate variables.The technical details and the considerations we made in our selections are provided in the accompanying Technical Box.In brief, we synthesized SST reconstructions generated by three proxies (  U K ′ 37 , TEX 86, and Mg/Ca).The  U K ′ 37 index is a ratio of organic molecules (alkenones) synthesized by selected haptophyte algae, which photosynthesize in the surface ocean (Marlowe et al., 1984;Prahl & Wakeham, 1987).The TEX 86 index is a ratio of organic molecules (isoprenoidal glycerol dialkyl glycerol tetraethers, isoGDGTs) synthesized by selected marine Archaea, in particular the ammonia oxidizing Thaumarchaeota (Schouten et al., 2002(Schouten et al., , 2013)).SSTs are also reconstructed from the Mg/Ca ratio in the shells of planktonic foraminifera: single-celled protists found at a range of depths in the upper parts of the ocean water column.As each of these three SST proxies has different biological sources and environmental influences, there is the potential to extract a rich suite of information by drawing on the similarities and differences of the reconstructed SSTs (discussed further in Section 3).We also synthesized foraminiferal δ 18 O records from species living in the upper water column (i.e., near the surface, δ 18 O planktonic ) and at the sea floor (δ 18 O benthic ).Foraminifera δ 18 O record both the influences of the temperature and δ 18 O composition of the seawater in which they live (Epstein et al., 1953).The seawater δ 18 O is also influenced by temperature and salinity, especially in the surface ocean, whereas in the deep sea an increasing influence of global ice volume is recognized (Shackleton & Opdyke, 1973) but not always easy to quantify (Evans et al., 2016;Raymo et al., 2018).The properties of the intermediate and deep waters which bathe the sites thus reflect surface-ocean conditions in the regions where they were formed through density-driven convection; in our compilation these include the high latitude oceans of both hemispheres and the Mediterranean Sea (e.g.Hodell & Venz-Curtis, 2006;Khélifi et al., 2014).As a result, multi-proxy analyses from a single site have the potential to record asynchronous changes between the surface and deep ocean.We also analyzed the results of several studies which reconstructed the sea-level variability contributions to δ 18 O benthic by accounting for differences in bottom water temperature variability through time (Miller et al., 2020), by accounting for non-linear changes to ice sheet δ 18 O during ice-sheet growth (Rohling et al., 2014), and by undertaking an inverse modeling of temperature and ice volume contributions to δ 18 O benthic (Berends et al., 2021).

Statistical Analysis
We adopted several statistical approaches to assess the character of the mPWP and to investigate how the climate system evolved across the iNHG.Previous analysis of the 5-0 Ma LR04 δ 18 O benthic stack using a Bayesian Change Point algorithm identified a change point at 2.73 Ma (±0.1 Ma), closely aligning with iNHG, which was marked by an increase in the obliquity signal (Ruggieri, 2013).To test whether regime boundaries also existed in the PlioVAR data sets (e.g. a shift in the mean or variance of the climate data), we applied the same Bayesian Change Point algorithm to the 53 individual proxy datasets and the three sea-level records which had high-resolution data for the whole 3.3-2.4Ma interval (Table 1).This algorithm calculates the probability density of the input time series, and uses a regression model and the combined application of forward recursion and Bayes rule to compute change point locations and their probability.As input, the program requires the number of change points to be found, the minimum distance between change points, and three hyperparameters that characterize the prior distribution of regression coefficients and residual variance (Ruggieri, 2013).
The hyperparameters are used to scale the prior distributions of regression coefficients (based on a multivariate normal distribution) and residual variance (scaled-inverse χ 2 distribution).Analyses using this program were performed using the same hyperparameters as used by Ruggieri (2013) for the Plio-Pleistocene LR04 δ 18 O benthic stack.
As we sought to identify whether iNHG represented a single shift in climate regime, as shown by Ruggieri (2013), we set the input to find only one viable change point within the analysis interval.This change point represents the mid-point of the climate shift that is, the point at which the preceding climate is different to the one that follows, rather than identifying the onset of a potential transition (Figure S3 in Supporting Information S1).To avoid single large glacial or interglacial peaks or data outliers artificially being indicated as change points, and accounting for the temporal resolution of our data sets, we fixed the minimum distance of detection to be equivalent to ∼10,000 years, regardless of the sedimentation rates at each site.Furthermore, some discontinuous records (8 out of 56) were interpolated using a cubic spline to their median time step prior to analysis.
We also assessed and compared the amplitude of glacial-interglacial variability spanning three periods for the δ 18 O benthic and SST data.We refer to these as intervals P-1, P-2, and P-3 throughout the manuscript (Figure 1).Interval P-1 broadly aligns with the mPWP (3.3-3.0Ma).Interval P-2 (3.0-2.7 Ma) represents the transition toward Interval P-3 (2.7-2.4 Ma), which is marked by the onset of consistently elevated IRD deposition in the high northern latitudes (Figure 1, Section 1).By comparing the amplitudes of glacial-interglacial variability we reduce the reliance on peak interglacial correlation and influence of individual age model uncertainty which might have impacted our comparison of globally distributed records with varying temporal resolution.For each record in the three intervals (P-1, P-2, P-3), we also generated the probability distribution of glacial-interglacial amplitudes for selected δ 18 O benthic and SST records (after Grant & Naish, 2021).To capture glacial-interglacial variability, the ranges of amplitude (maximum-minimum values) were calculated within a centered 41 kyr moving window, selected as the primary periodicity modulating glacial cyclicity, using R package Astrochron (Meyers et al., 2021).A minimum sampling of <10 kyr is required to resolve peak glacial-interglacial amplitudes.Data sets which met this requirement (n = 49) were linearly interpolated to the mean sampling resolution for each site (1-5 kyr), which was also used for the window time-step (i.e., the time between centered windows).
Finally, we assess the amplitude of glacial-interglacial sea-level cycles for the same three time intervals following Grant and Naish (2021).The full amplitude of sea-level change from colder-than-present glacial states to warmer interglacials provides insight on ice-sheet sensitivity, the magnitude of glacial-interglacial variability, and also allowed us to evaluate how changes in the cryosphere during the iNHG compared to the evolution of SSTs discussed in this manuscript.As we did for the δ 18 O benthic and SST records, the glacial-interglacial amplitudes were generated as the maximum range within a 41 kyr moving-window, for intervals P-1, P-2, and P-3, using a mean sampling step of the record (2 kyr).The amplitudes derived from this approach are on a floating scale (i.e., they sit unanchored to a Holocene reference).
3. Are There Proxy-Specific Differences in the Records of SST Change Through Time? (Q1) Previous global-scale syntheses of Pliocene climate have used reconstructed SSTs as key variables, especially for data-model comparison (e.g., Dowsett et al., 2012;Haywood et al., 2020;McClymont et al., 2020;Z. Zhang et al., 2021).This previous work has allowed quantification of climate response to pCO 2 forcing (Martinez-Boti et al., 2015), and enabled comparison to other geological time intervals where forcings and feedbacks may have been different (e.g.Pre-Industrial, Last Glacial Maximum, Eocene) (Fischer et al., 2018;Gulev et al., 2021).Although there is often good correspondence between different SST proxies, differences in absolute values and trends can arise which may vary by location but also through time (Lawrence & Woodard, 2017;McClymont et al., 2020).In turn, it can become challenging to calculate site-specific and global-scale temperature anomalies, and to undertake data-model comparison (e.g.Inglis et al., 2020;McClymont et al., 2020).Since each proxy system has a different biological origin with different biological and environmental controls on how the temperature signal is recorded (Section 2), we propose that there is the potential to understand the evolution of Pliocene and early Pleistocene environments more thoroughly by adopting a multi-proxy approach.
In this section, we examine SST records from three regional settings for which the most oNHG and iNHG temperature data have been generated: (a) the North Atlantic Ocean; (b) low latitude warm pools; (c) low latitude upwelling regions (Figure 2).The different oceanographic settings allow us to explore likely controls over the proxy signals.Since we have already accounted for a modest Mg/Ca seawater correction (Technical Box), offsets between foraminifera Mg/Ca reconstructions and other proxies need to be explained by other factors, which could include for example, the preferred season of growth, or the water depth where the source organisms were living (Bova et al., 2021).We sought to identify whether overlaps or offsets existed between data sets which might be consistent with our understanding of regional processes and biological controls.This was achieved by: (a) examining multi-proxy data from individual sites, where more than one proxy has been analyzed from the same sediment sequence; (b) examining SST data at a regional scale (Figure 3); and (c) comparing the PlioVAR data sets to published lower-resolution data.To facilitate comparison between sites where there may be different temporal resolutions of data, we plot PlioVAR SSTs as 20 kyr means centered on each glacial and interglacial maximum/minimum as determined by δ 18 O benthic (Figure 3).SST maxima or minima do not always align with the peak glacial or interglacial states as determined by the marine isotope stages (MIS) of the δ 18 O benthic record (e.g.Capron et al., 2014).However, our approach here is to use δ 18 O benthic as a simple framework for identifying mean SSTs within broad ∼20 kyr windows, to assess whether any proxy-specific bias occurs across iNHG.

A Comparison of Multi-Proxy Temperature Records From Single Sites
An important limitation to our analysis is that only 2 of the 22 PlioVAR sites examined here have both high temporal resolution and multi-proxy SST data recovered from the same location to enable direct proxy-proxy comparison across the whole interval of study (3.3-2.4Ma): at the North Atlantic site IODP U1313, and the South China Sea site ODP 1143 (Table 1; Figure 2).
Mg/Ca in G. bulloides consistently yields 1-2°C higher SSTs than  U K ′ 37 at the North Atlantic site (IODP Site U1313), excluding MIS 100 (∼2.52 Ma) where Mg/Ca G. bulloides temperatures are within error of, but slightly cooler than,  U K ′ 37 (Figure 3a).At this latitude (41°N),  U K ′ 37 is calibrated to mean annual SST (Müller et al., 1998;Tierney & Tingley, 2018), suggesting that the G. bulloides temperatures may be reflecting warmer seasonal SSTs.We did not analyze the G. ruber Mg/Ca data from IODP Site U1313 (Bolton et al., 2010(Bolton et al., , 2018) ) due to discontinuous sampling across our study interval, but the reconstructed G. ruber SSTs are warmer than from G. bulloides, consistent with modern summer production in G. ruber and spring production in G. bulloides, as also observed in other North Atlantic records (Friedrich et al., 2013;Hennissen et al., 2014;Robinson et al., 2008).The enhanced cooling in G. bulloides Mg/Ca during glacial stages (Figure 3a) thus indicates either enhanced spring cooling associated with iNHG, or a shift to deep (colder) growth temperatures, given that G. bulloides occupies a mixed-layer habitat (the upper ∼60 m) (e.g.Schiebel et al., 1997) compared to the upper 10 m for  U K ′ 37 (Müller et al., 1998).
In contrast, in the low latitude warm pools (Figure 2),  U K ′ 37 SSTs are within proxy calibration uncertainty of the Mg/Ca temperature estimates, as shown in the South China Sea (ODP Site 1143) with Mg/Ca in G. ruber (Figure 3), but also when comparing a lower time resolution  U K ′ 37 record to the high-resolution T. sacculifer Mg/Ca estimates from the Caribbean Sea (ODP Site 999, Badger et al., 2013;Figure S4 in Supporting Information S1).However,  U K ′ 37 temperatures are consistently higher than Mg/Ca in the South China Sea.Strong seasonal variations in SST and mixed-layer depth are recorded today in the South China Sea, linked to the East Asian monsoon system (Twigt et al., 2007).A summer signal from G. ruber (Lin & Hsieh, 2007) is unable to explain the cooler Mg/Ca temperatures we reconstruct relative to the (mean annual)  U K ′ 37 signal; rather, the high winter fluxes of G. ruber observed today (Lin et al., 2011) may also have imparted a relatively cool signal during the late Pliocene and early Pleistocene.However, the reconstructed SSTs at ODP Site 1143 sit in the warmest part of the  U K ′ 37 calibration, where errors approach 4.4°C (Tierney & Tingley, 2018), meaning that in effect there is agreement, within error, between the two proxies in these high resolution records.

Regional-Scale Comparisons of Ocean Temperature Change
As an alternative to multi-proxy comparisons within individual sites, we compared the regional-scale glacial and interglacial patterns of SSTs from the PlioVAR data sets (Figure 3) to published lower-resolution data sets.In the North Atlantic Ocean, there is greater variability in  U K ′ 37 records between sites, than between the co-registered Mg/Ca and  U K ′ 37 data from IODP Site U1313 (Figure 3a).This variability demonstrates the sensitivity of these sites to changes in the position of the subtropical and subpolar gyres (e.g.IODP Site U1313, ODP Site 982) and the extent of Arctic surface water masses (e.g.ODP Site 907).With only one site containing both Mg/Ca and  U K ′ 37 data we are unable to assess whether there are any systematic differences in how SST is recorded by these proxies (e.g.season or depth of production).Several low-resolution multi-proxy time series are available from the North Atlantic Ocean, which did not meet the PlioVAR thresholds for analysis (Section 2.1; Figure 2), but they offer some insights into regional similarities and differences between our SST proxies.
For example, a reconstruction centered on MIS M2 (spanning MIS MG2-M1, ∼3.34-3.24Ma) showed comparable  U K ′ 37 trends among three sites along the path of the North Atlantic Current (spanning∼40 to ∼53°N) (De Schepper et al., 2013).SSTs from  U K ′ 37 were consistently ∼2°C warmer than Mg/Ca-temperatures from G. inflata (De Schepper et al., 2013) consistent with the interpretation that G. inflata reflect deeper (∼200-300 m) and thus cooler temperatures (De Schepper et al., 2013).Low-resolution  U K ′ 37 and TEX 86 temperatures ∼53°N (DSDP Site 610) had no consistent offset through time, although their SSTs were within analytical and calibration errors (Naafs et al., 2020).In contrast, the  U K ′ 37 and TEX 86 temperatures were consistently warmer than G. bulloides Mg/ Ca data from DSDP Site 610, where both proxies spanned the iNHG (Hennissen et al., 2014;Naafs et al., 2020).Cooler SSTs from G. bulloides might also reflect a spring or mixed-layer temperature signature here, as considered above for IODP Site U1313 (Friedrich et al., 2013;Hennissen et al., 2014;Robinson et al., 2008).Our regional-scale proxy comparison for the North Atlantic Ocean supports published interpretations of  U K ′ 37 and TEX 86 as proxies for mean annual or summer SSTs, compared to summer (G.ruber), spring or mixed layer (G.bulloides) or thermocline (G.inflata) proxies from foraminiferal Mg/Ca.
In the low latitude warm pools no long-term trends in glacial or interglacial SSTs were determined, regardless of the proxy (Figure 3b).However, both  U K ′ 37 and Mg/Ca SSTs in the South China Sea (ODP Site 1143) consistently exceed the Mg/Ca SSTs from all other warm pool sites, despite sitting more distal to the heart of the west Pacific warm pool (as recorded by proxy data from ODP Site 806, Figure 3b).The apparently warm Mg/Ca temperatures in the South China Sea may be explained by changes in salinity: a climate model which reproduces the lower zonal and meridional SST gradients shown in Pliocene proxy data (Burls et al., 2017) also generates saltier conditions at ODP Site 806 (+0.5 psu) and much saltier conditions at ODP Site 1143 (+2.5 psu).Given the salinity effect on Mg/Ca (+5.4% per psu) (Gray & Evans, 2019), increased Pliocene salinity would cause the Mg/Ca temperatures in the South China Sea to appear ∼1.3°C warmer than at ODP Site 806; without this effect, mPWP temperatures at the two sites are within ∼0.3°C (Figure 3b).A difficulty with this explanation is that the low-resolution TEX 86 temperatures from ODP Site 1143 are even warmer than both  U K ′ 37 and T. sacculifer Mg/Ca temperatures in the Pliocene (O'Brien et al., 2014), but cannot be explained by a response to elevated salinity alone.In contrast, Pliocene SSTs in the Caribbean Sea (ODP Site 999) are ∼2.5°Ccooler than the other warm pool sites (Figure 3b).Paleo-salinity changes are predicted to have been negligible (Burls et al., 2017;Rosenbloom et al., 2013), but seafloor dissolution may have been different (Groeneveld et al., 2014;Haug & Tiedemann, 1998); more corrosive deep waters in the Pliocene could have biased Mg/Ca-based temperatures to artificially cold values, which could be tested in the future at these sites by generating proxy data for bottom water carbonate ion saturation in tandem with Mg/Ca analysis (Rosenthal et al., 2022).
Although long-term evolution of the major upwelling systems have been intensively studied through the Pliocene and Pleistocene (e.g., Dekens et al., 2007;Etourneau et al., 2009;Rosell-Melé et al., 2014), the PlioVAR data sets are solely based on  U K ′ 37 analyses (Figure 3; Table 1).Yet, modern studies and multi-proxy analysis indicates the potential for seasonality or water depth to influence both the organic and inorganic proxies that we used (e.g.Hertzberg et al., 2016;Ho et al., 2011;Leduc et al., 2014;Lopes dos Santos et al., 2010;McClymont et al., 2012;Petrick et al., 2018;White & Ravelo, 2020).Additional low-resolution multi-proxy data are available at two upwelling sites from the south-east Atlantic, ODP Site 1087 (TEX 86 ) and ODP Site 1082 (G.bulloides Mg/Ca), where  U K ′ 37 consistently records the warmest SSTs.The relatively cold TEX 86 temperatures have been attributed to sub-surface production (Petrick et al., 2018;Seki et al., 2012), as is also observed in both late Pleistocene and modern upwelling systems (e.g.Hertzberg et al., 2016;Lopes dos Santos et al., 2010;McClymont et al., 2012) and discussed in detail for the Pliocene (Ford et al., 2015;White & Ravelo, 2020).Pliocene and Pleistocene SSTs from Mg/Ca in G. bulloides are consistently colder than  U K ′ 37 SSTs at ODP Site 1082, regardless of whether a Mg/ Ca seawater correction is applied (up to 10°C, Leduc et al., 2014; ∼1.5°C with Mg/Ca seawater correction, this study, Figure S2 in Supporting Information S1).The  U K ′ 37 -Mg/Ca offset was attributed to G. bulloides recording a deeper, winter signal during upwelling intensification at ODP Site 1082 (Leduc et al., 2014).In turn, if  U K ′ 37 is recording warm-season SSTs (Leduc et al., 2014), this may partially explain some of the observed data-model offsets in the mPWP, given that the model results were of mean annual SST (McClymont et al., 2020).Care is thus needed at upwelling sites to select the most appropriate SST reconstructions (e.g whether annual or seasonal) and interpret their results in relation to both complex local oceanographic change and the imprint of global climate transitions.

A Comparison of Ocean Temperature Proxy Records: Summary and Recommendations
Our synthesis and proxy comparison of regional SST data sets spanning the mPWP and iNHG demonstrates that for most sites,  U K ′ 37 , TEX 86 and Mg/Ca in surface-dwelling foraminifera can provide robust reconstructions of past mean annual SST.We have also shown that under certain circumstances, local or regional influences over the season or depth of production by the source organisms can influence the temperature signals, which can be identified through multi-proxy analysis from individual sites.
Our analysis highlighted that care should be taken to ensure that G. bulloides Mg/Ca data is reflecting a surface-dwelling or annual signal when deciding whether to include it into global-scale syntheses of mean annual SSTs.This is because G. bulloides-derived SSTs may be biased toward spring or subsurface temperatures, but with a relationship that may not be constant through time, for example, if there have been orbital or longer-term changes in seasonality, mixed-layer depth, or upwelling location/intensity.TEX 86 data may also show subsurface bias, and should be excluded from global SST syntheses where there is evidence for an increase in upwelling/ subsurface export, either directly (i.e., high GDGT-2/GDGT-3 ratios; van der Weijst et al., 2022) or indirectly (e.g.where microfossil assemblages highlight changing upwelling intensity).Incorporating seasonal or subsurface signals within global SST syntheses would impact data-model comparison at the affected sites, and would also impact the amplitude of changes in global SST used in assessments of climate sensitivity.
Our analysis has highlighted the value of multi-proxy temperature reconstructions of the same samples from the same core, since this approach offers the best opportunity to explore the processes which could influence individual proxies, including changing mixed layer depth, seasonality, or upwelling intensity through time.A multi-proxy approach also allows the assessment of contributions by non-thermal factors to all of our proxy meth-ods.Unfortunately, such an approach has rarely been applied at the sites we have examined here: the existing resolution of single-site multi-proxy data prevented us from achieving our goal of evaluating in detail the processes influencing differences and similarities in temperature records generated within and between sites.However, as the majority of our records do record mean annual SSTs, with some caveats, we utilize all of the SST proxies we have discussed here (  U K ′ 37 , TEX 86 , Mg/Ca) during our analysis and discussion in Section 4.

Assessing Trends and Variability in Ocean Temperatures, Ice Volume, and Sea Level Associated Across oNHG and iNHG
The transition from the late Pliocene to the early Pleistocene features ice-sheet expansion and a fall in pCO 2 (Figure 1; Section 1).An analysis of the LR04 δ 18 O benthic stack (Lisiecki & Raymo, 2005) showed that there was a shift away from a relatively noisy relationship between δ 18 O benthic and orbital forcing (a "stochastic" climate) ∼2.8 Ma (Meyers & Hinnov, 2010).After 2.8 Ma, a closer relationship between δ 18 O benthic and orbital forcing was identified (as a more predictable, or "deterministic" climate), and attributed to an increasing influence of larger ice sheets over climate feedbacks (Meyers & Hinnov, 2010).A change point in the same δ 18 O benthic stack at 2.73 Ma (±0.1 Ma) was driven by an increase in the orbital obliquity signal (Ruggieri, 2013), and also highlights a shift in the climate forcing-response relationship.As the oceans are a critical part of heat and moisture transport through the climate system, interacting with the ice sheets in a variety of ways, we seek to identify their response to, or influence over, the iNHG.
The focus of this Section is to present our analysis of the high-resolution, well-dated records of SST, ice volume, and sea level in our PAGES-PlioVAR synthesis, which we refer to here as the "PlioVAR data sets" (Figure 2, Table 1).Our aim is to constrain the timings of events associated with the iNHG, and to understand the potential interactions between climate variability and longer-term climate evolution.As outlined in Section 1, we make a distinction between the oNHG (a broad window of time with evidence for slow ice-sheet growth) and iNHG (a more focused interval with a rapid increase IRD to the high latitude northern oceans).We first assess the statistical signatures of long-term trends (Section 4.1) and shorter-term variability with 41 kyr-obliquity frequency (Section 4.2) in the PlioVAR data sets.We then draw on complementary published data to evaluate and explore the processes driving the patterns of climate change we observe, by targeting the three research questions we posed in Section 1: Q2 What were the characteristics of the mPWP (∼3.025-3.264Ma) (Section 4.3); Q3 What were the characteristics of the iNHG?(Section 4.4); Q4 Did the amplitude or pacing of climate variability change with iNHG?(Section 4.5).Alongside our analysis of the PlioVAR data sets, we draw on additional marine sediment evidence to illustrate environmental changes which occurred onshore, particularly ice-sheet expansion (using IRD) (Andrews, 2000), and atmospheric processes via aeolian deposition of terrigenous materials which can reflect wind intensity, atmospheric circulation shifts, and sediment availability for transport for example, from aridity (Rea, 1993).

PlioVAR Data-Sets: Long-Term Trends
To test whether climate change during the iNHG had a globally synchronous onset, we employed change point analysis on the SST, δ 18 O benthic , δ 18 O planktonic and sea-level data outlined in Section 2 (Figure 4, Table 3).The change point analysis serves to test whether there have been statistically significant regime shifts, by identifying the mean timing of the point which divides each record into two distinct intervals according to their mean state or orbital-scale variability (for the methodology see Section 2.3).Since δ 18 O benthic , δ 18 O planktonic and sea-level data should include a signature of global ice volume change, our null hypothesis was that a common change point would detail the onset of globally significant ice-sheet growth consistent with the iNHG.In contrast, regional or record-specific change points would reveal additional influences over δ 18 O including temperature and salinity, which may be most strongly expressed in δ 18 O planktonic data (Technical Box).Our SST records then give a contrasting analysis of surface ocean properties, providing information about ocean circulation and temperature changes which might be linked to other forcings (e.g.pCO 2 or gateway changes).
When considering the whole period of interest (3.3-2.4Ma), a wide temporal range of change points is detected (Figure 4, Table 3).This is not a reflection of our narrow time window of interest (3.3-2.4 Ma), since analysis of SST and δ 18 O planktonic records from some of the same sites using a different change point approach and a broader time window also identified a wide spread of change points with no strong regional patterns (Kaboth-Bahr &

Table 3
Change Point Analysis for the 3.3-2.4 Ma Time Interval as Shown in Figure 4 Mudelsee, 2022).One  U K ′ 37 (DSDP Site 593), one TEX 86 (IODP Site U1463) and three δ 18 O planktonic (ODP Sites 1143, 1148, and 214) records returned no change points.No change points were detected in the three sea-level time series (Table 3).
Our change point analysis challenges our long-held view that δ 18 O benthic records represent a globally synchronous signal dominated by ice volume.Instead, local differences in bottom-water temperatures or seawater δ 18 O may be more important in the Pliocene and early Pleistocene records (Mudelsee & Raymo, 2005).We found a broad spread in δ 18 O benthic change points, covering the full range from 3.29 to 2.42 Ma (Table 3, Figures 4a and 4e).This is consistent with a previous analysis of individual δ 18 O benthic records using a statistical "ramp" to describe the start, end and amplitude of events linked to iNHG, which identified a broad window of changes to δ 18 O benthic between 3.6 and 2.4 Ma (Mudelsee & Raymo, 2005).The earliest δ 18 O benthic change points in the PlioVAR data sets occurred between ∼3.3 and 3.0 Ma (our Interval P-1, Figure 1, Section 2.3), indicating that the climate system was already showing signs of change during the mPWP.The earliest mPWP δ 18 O benthic change points developed from ∼3.3 Ma, in the Southern Hemisphere, but these have low probabilities and suggest that a regime shift was unlikely to have occurred, or was weakly expressed in the data (Figures 4e and 4f).From ∼3.2 Ma, δ 18 O benthic change points with higher probabilities occurred in the northern low latitudes (<30°N) (Figures 4a,  4e, and 4f).Change points during the mPWP were also recorded in SSTs, first from ∼3.3 Ma at sites influenced by the extent of subpolar surface water masses and upwelling systems and the Mediterranean Sea, but also with relatively low probability (Figures 4a, 4c, and 4d).From ∼3.2 Ma, higher probability SST change points occurred in the mid-and high-latitudes of the Northern Hemisphere (Figures 4a, 4c, and 4d).After the mPWP, only two change points occurred between 3.0 and 2.8 Ma, both at ∼2.9 Ma (Figures 4c and 4e) in the Mediterranean Sea (δ 18 O benthic , SST).
From 2.8 Ma onwards a wide latitudinal spread of δ 18 O benthic change points occurred, first indicated by records influenced by mid-latitude climate change in both hemispheres.A cluster of change points between 2.8 and 2.7 Ma occurred in three sites which are bathed by the intermediate waters (∼1,000-1,500 m water depth) formed in the surface ocean of the Subantarctic region (South-east Atlantic, ODP Site 1082, and South-west Pacific, DSDP Sites 593 and 594) (Figure 4e).These change points coincided with a δ 18 O planktonic change point in the Subantarctic Atlantic (ODP Site 704), and occurred within a broader interval of southern mid-and high-latitude cooling as recorded by SSTs in the South-west Pacific and South-east Atlantic (Benguela) upwelling system at ∼2.8-2.6 Ma (McClymont et al., 2016;Petrick et al., 2018).In the mid-latitude North Atlantic (IODP Site U1313) there are synchronous and high probability change points in G. bulloides Mg/Ca and δ 18 O benthic at 2.77 Ma (noting that the δ 18 O planktonic change point comes later, at 2.56 Ma).
Regional differences in the timing and character of climate evolution between low-and high-latitude signals across the iNHG have previously been identified (Kaboth-Bahr & Mudelsee, 2022;Ravelo et al., 2004).Our results confirm comparatively late change points in sites from the low-latitudes, given a final cluster of change points in δ 18 O benthic from the South China Sea and Indian Ocean from 2.5 to 2. 4 Ma (ODP Sites 1143, 1148, andODP758), but acknowledging that for three of these sites the change point probabilities are very low (Figures 4a, 4d and 4e).A final cluster of low-latitude SST change points occurred between 2.6 and 2.4 Ma: all the change points in low-latitude Pacific and Indian Ocean Mg/Ca records occurred during this time interval (but with low probability), alongside high probability change points in  U K ′ 37 -derived SSTs (Figure 4).The timing of these change points is broadly aligned with the pronounced glacial-interglacial cycles of MIS 96-100 (Figure 4).

PlioVAR Data-Sets: Orbital-Scale (Glacial-Interglacial) Variability
In the absence of synchronous change points (Section 4.1), we analyzed three broad time windows to compare orbital-scale SST and δ 18 O benthic variability at the 41 kyr-obliquity frequency (Figure 1, Section 2.3): the sustained warmth of the mPWP (Interval P-1, 3.25-3.0Ma), the transition toward intensified glaciations in the late Pliocene (Interval P-2, 3.0-2.7 Ma) and the main interval of intensified glaciations spanning the latest Pliocene to early Pleistocene (Interval P-3, 2.7-2.4Ma) (Figure 5).Although the MIS M2 glacial stage is marked by cooling outside of the mPWP SST range at some sites (Table 4, Figure S6 in Supporting Information S1), its inclusion in the P-1 Interval does not impact the broad patterns in the data.To test whether iNHG was associated with climate cooling, we calculated the mean SSTs for each of the three time windows: in 19 ocean temperature records (17 sites), the mPWP (Interval P-1) was warmer than Interval P-3 (Table 5), but only 7 records (6 sites) showed cool- Site locations are shown in Figure 2. "Na" means data was analyzed but no statistical result was given.Blank spaces mean either data is not available or it is insufficient for analysis (cf.Table 1).
The probability that the change point is statistically significant is given in parentheses after each change point date.Table 3 Continued ing which exceeded calibration error, and three sites recorded warming (Table 5).The largest cooling generally occurred in the mid-latitude sites, whereas small SST changes occurred in the low latitudes.
We explored the probability distribution of SST variability by latitude using "violin" plots for SST (Figure 5) and δ 18 O benthic (Figure 6).In each plot, the mean variability is delineated by the solid line joining the sites, whereas the amplitude of the variability is shown by the vertical extent of data per site, and the probability of each amplitude is shown by the width of the "violin."For SST variability, a broad increase in the mean and the amplitude of SST variability with increasing latitude was identified in all three time intervals (Figure 5), noting that the maximum latitude recorded for the Northern Hemisphere (ODP 907, 69°N) is greater than in the south (DSDP Site 594, 46°S).In contrast to this broad latitudinal pattern, the largest increase in mean variability through time occurred in the mid-latitude North Atlantic (IODP Site U1313, from Mg/Ca), where the range of variability was also large (as seen in both Mg/Ca in G. bulloides and  U K ′ 37 ).Comparable latitudinal patterns of variability have been observed through the transition from the Last Glacial Maximum (∼27-19 ka) into the Holocene interglacial (e.g.Rehfeld et al., 2018) and across multiple late Pleistocene glacial-interglacial cycles (Rohling et al., 2012).These patterns reflect higher sensitivity to global radiative forcing at high latitudes, alongside interactions between ice albedo (mid and high latitudes) and potentially hydrological feedbacks in the low latitudes, which can enhance mid-latitude SST variability (Rohling et al., 2012).
The Southern Hemisphere sites displayed small increases in the mean and range of orbital-scale SST variability after the mPWP.First, variability increased in the mid-latitude sites as early as the transition from Interval P-1 to P-2 (∼3.0 Ma), especially in sites where the reconstructed SSTs are influenced by the position of the Subantarctic Note.SSTs were reconstructed using the  U K ′ 37 index, TEX 86 index, or Mg/Ca in three planktonic foraminifera species (G.bulloides, (MgCa-b), T. sacculifer (MgCa-s), and G. ruber (MgCa-r)).For each time interval, the mean, standard deviation (2σ) and the number of data points used (n) are shown.The KM5c-M2 anomaly is generated by subtracting the M2 mean from the KM5c mean, so that positive values equal a warming trend.The uncertainty in the KM5c-M2 SST difference is calculated as the sum of the uncertainties for each interval, calculated by dividing the 2σ range by the square root of the number of data points used in that time interval.

Table 4
Differences in Mean Sea-Surface Temperatures (SSTs) Between the MIS M2 Glacial  and the KM5c  front in both the and Atlantic Oceans (e.g.∼41°S, Southwest Pacific, DSDP Site 593) (e.g.McClymont et al., 2016).A second increase in SST variability occurred in the equatorial Pacific sites from Interval P-2 to P-3 (Figure 5; ODP 806 and 846).
A more complex picture of SST variability emerged temporally and spatially in the Northern Hemisphere, which may be linked to a long-term expansion of the subpolar gyres (Martinez-Garcia et al., 2010;Sánchez-Montes et al., 2020) and intensification of subtropical circulation (Fedorov et al., 2015).Between Interval P-1 and P-2 (∼3.0 Ma) there were increases in the mid-and high-latitude SST variability in several North Atlantic sites and in the Mediterranean Sea (Figure 5).Subsequently, moving from Interval P-2 to P-3 (∼2.7 Ma) there was a large increase in SST variability in the mid-latitude Atlantic (IODP Site U1313), and in the North-east Pacific (Gulf of Alaska, IODP Site U1417; California Current, ODP Site 1012).
The same statistical analysis was applied to the δ 18 O benthic data (Figure 6).A systematic increase in both the amplitude and mean value of orbital-scale variability occurred from P-1 to P-3; Figure 6.This is clearly demonstrated in the LR04 δ 18 O benthic stack (Lisiecki & Raymo, 2005), which showed a small increase in mean variability at  Note.The P3-P1 anomaly is generated by subtracting the P-1 mean from the P3 mean, so that negative values equal early Pleistocene cooling relative to the Pliocene.The uncertainty in the P3-P1 SST difference is calculated as the sum of the uncertainties for each interval, calculated by dividing the 2σ range by the square root of the number of data points used in that time interval.Sites where the P3-P1 difference is less than the uncertainty are italicized.If the M2 is excluded from the P-3 interval (i.e., 3.25-3.0Ma) the P-3 mean is changed by ≤0.3°C, except at ODP Site 907 (0.9°C difference) and ODP Site 982 (0.4°C difference).

Table 5
Quantifying Although no change points were detected in the three sea-level records (Section 4.1), a consistent signal of increasing orbital-scale variability moving from Interval P-1 to Interval P-3 was recorded (Figure 7).The highest Figure 6.Violin plots (symmetrical vertical histograms) of glacial-interglacial δ 18 O benthic amplitudes for each site, calculated as the maximum range within a moving 41-kyr window along the duration of each period, plotted using the same approach as for Figure 5. Colors refer to the water depth at the core site, except for the global LR04 stack (Lisiecki & Raymo, 2005) which is shown in black on the right.

What Were the Characteristics of the Mid-Piacenzian Warm Period (mPWP, ∼3.264-3.025 Ma)? (Q2)
We have identified several SST and δ 18 O benthic change points which occurred within the mPWP, showing that this period included long-term changes in climate regime despite relative stability in the pCO 2 record (Figure 1).The earliest change points spanned both the MIS M2 glacial stage, which immediately preceded the mPWP, and the earliest interglacial of the mPWP (KM5c), that is, events were developing throughout our Interval P-1 (3.3-3.0Ma).In this section, we compare our change point results with published records of ice-sheet size and configuration and complementary records of several key climate systems, to explore the nature of climate change occurring during both the M2 glacial stage and the mPWP.
The PlioVAR data sets reveal a wide range of SST responses to the MIS M2 glacial stage and the transition into interglacial KM5c (∼3.25 Ma, Figures 1 and 5).The MIS M2 has previously been identified by a short-lived but pronounced increase in δ 18 O benthic , argued to reflect the onset of more extensive glaciation in the Southern Hemisphere (e.g., McKay et al., 2012;Naish et al., 2009).There is also evidence for some ice-sheet growth in Greenland, Iceland and Svalbard during MIS M2 (De Schepper et al., 2014;Knies, Cabedo-Sanz, et al., 2014).However, despite MIS M2 being described as a possible "harbinger" of the onset of Pleistocene glacial cycles (Westerhold et al., 2020), there is a lack of widespread evidence for a large Northern Hemisphere ice-sheet expansion comparable to those of the Pleistocene (Kirby et al., 2020;Tan et al., 2017).Relative to the interglacial MIS KM5c, the PlioVAR data sets quantified the largest M2 cooling signals in SSTs from the mid-and high-latitudes of both hemispheres (2-4°C) (Table 4).Minimal cooling (<1°C) occurred in low-latitude SSTs (Table 4).This explains why the broad latitudinal increase in orbital-scale variability during the Interval P-1 was not altered if we included M2 in our analysis (Figure 5, Figure S5 in Supporting Information S1).Since only low probability change points in SST, δ 18 O benthic and δ 18 O planktonic occurred in the early part of Interval P-1, coinciding with MIS M2 and KM5c (3.3-3.2Ma), we conclude that the changes in climate regime were weakly expressed.In contrast, from ∼3.15 Ma change points with higher probabilities emerged, which were most significant in Northern Hemisphere SSTs and two low-latitude δ 18 O benthic records at ∼3.1 Ma, during the second half of the mPWP.Violin plots (symmetrical vertical histograms) of glacial-interglacial sea-level amplitudes for each site, calculated as the maximum range within a moving 41-kyr window along the duration of each period, plotted using the same approach as for Figure 5.The records shown were generated using different approaches to extract sea-level change from the global δ 18 O benthic stack (Berends et al., 2021;Miller et al., 2020;Rohling et al., 2014) as outlined in Section 2.5.

10.1029/2022RG000793 of 43
Long-term shifts in regional circulation systems through the mPWP may be reflected in the wide latitudinal distribution of SST change points, which have then been captured locally with apparent asynchroneity in response to migrations of the positions of surface ocean fronts or upwelling systems through time.The SST change points which occur during the mPWP, and which precede δ 18 O benthic change points, suggest a rapid surface ocean response to forcing, either through direct radiative forcing or via locally sensitive shifts to surface ocean circulation systems.Early SST changes preceding δ 18 O benthic have also been observed in the mid-and high-latitude North Atlantic region over mid and late Pleistocene glacial-interglacial cycles in response to shifts in the positions of the Subpolar or Arctic Fronts (e.g.Alonso-Garcia et al., 2011;Barker et al., 2015;Hernández-Almeida et al., 2012, 2013;Mokeddem et al., 2014;Wright & Flower, 2002).In contrast, the later development of δ 18 O benthic change points will include the additional signal of the slower expansion of continental ice (e.g.Alonso-Garcia et al., 2011;Elderfield et al., 2012).
Although the western Mediterranean Sea SST and δ 18 O benthic change points had low probability (ODP Site 978, Figure 4), they aligned with reconstructed increases in west Mediterranean Sea surface salinity (∼2 psu) and bottom water salinity (∼1 psu; also ODP Site 978), which are proposed to have contributed to a strengthening of Mediterranean Outflow Water into the Atlantic Ocean by 3.3 Ma (Khélifi et al., 2014).The Mediterranean change points at ∼3.3 Ma are also consistent with more restricted Indonesian Throughflow after ∼3.5 Ma (Auer et al., 2019;De Vleeschouwer et al., 2022;Karas et al., 2009Karas et al., , 2011)), which impacts Mediterranean SSTs and surface salinities by reducing African monsoon rainfall (Sarnthein et al., 2018).Increasing aridity in Eastern and Central Asia from ∼3. 6 Ma (Lu & Guo, 2014) or earlier (W.Zhang et al., 2018) have also been linked to a reduction in atmospheric water vapor as a result of a cooling climate.
Further evidence for longer-term cooling preceding, then continuing through, the mPWP comes from both highand low-latitude regions.Expanding polar water mass extent in the Subantarctic South Atlantic from ∼3.5 Ma (Martinez-Garcia et al., 2010) aligns with cooling and sea ice expansion within the mPWP in the Ross Sea (Riesselman & Dunbar, 2013).Enhanced dust deposition to the Subantarctic South Atlantic (ODP Site 1090) reflects intensification or equatorward displacement of southern hemisphere wind systems, likely causing cooling in the Benguela upwelling system (Southeast Atlantic: ODP Sites 1081, 1082, 1084, and 1087) (Martinez-Garcia et al., 2011;Petrick et al., 2018;Rosell-Melé et al., 2014).A more restricted Indonesian Throughflow would also have reduced heat transport through the Indian Ocean (Auer et al., 2019;De Vleeschouwer et al., 2022;Karas et al., 2009Karas et al., , 2011)), and led to cooling in the southeast Atlantic Ocean (Karas et al., 2011;Rosell-Melé et al., 2014).
In the North Atlantic (ODP Site 982), SSTs cooled from ∼3.6 Ma as Arctic sea-ice expansion was recorded at DSDP Site 610 (Karas et al., 2020 and references therein), and continued across the mPWP (ODP Site 982) (Lawrence et al., 2009).The slow growth in global ice volume starting as early as ∼3.6 Ma (Mudelsee & Raymo, 2005) further highlights shifting high latitude climates before and through the mPWP.In several low-latitude sites, the M2 glacial stage marks the inflection point between a cooling trend starting from ∼3.5 Ma and a subsequent warming through the mPWP: in the Atlantic (ODP Site 662), Pacific (ODP Site ODP846) (Herbert et al., 2010) and Indian (ODP Site 214) Oceans (Karas et al., 2009).The mPWP thus straddles an interval of long-term regional climate changes, which particularly connect the low latitudes and Southern Hemisphere sites.
Our analysis supports previous work indicating that the oNHG was a long-term transition, involving complex interactions between different climate system components and regions, potentially beginning before the M2 and extending through the mPWP (as early as ∼3.6 Ma;Kleiven et al., 2002;Meyers & Hinnov, 2010;Mudelsee & Raymo, 2005).Orbital-scale pCO 2 reconstructions are lacking before M2, but available data provide values which were either comparable to, or slightly lower than, the mPWP (de la Vega et al., 2020;Pagani et al., 2010;Seki et al., 2010), indicating that pCO 2 forcing is unlikely to explain the complex climate trends observed here.Herbert et al. (2010) proposed that the long-term trends in low-latitude SSTs (which spanned ∼300-500 kyr) could reflect a response to long-wavelength orbital forcing, favoring glacial-stage cooling during eccentricity minima.However, this mechanism does not explain signals of cooling from ∼3.5 Ma, since they develop during an eccentricity maximum.
Alternatively, it has been proposed that pre-mPWP shifts in heat transport through the Southern Hemisphere could have resulted from a reduction in the Indonesian Throughflow (Auer et al., 2019;De Vleeschouwer et al., 2022).As the Indonesian Throughflow is influenced by regional sea level, both tectonic changes to the gateway config-10.1029/2022RG000793 of 43 uration (Auer et al., 2019) and a fall in sea level in response to continental ice sheet growth (De Vleeschouwer et al., 2022) have been proposed to explain the reduction in the throughflow.In the high northern latitudes, the only likely major sources of IRD before ∼2.7 Ma were Svalbard (Knies, Cabedo-Sanz, et al., 2014) and northeast Greenland (Bachem et al., 2017;Jansen et al., 2000).Hence, it is likely that a significant component of the increasing δ 18 O benthic across the wider iNHG interval (from ∼3. 6 Ma;Mudelsee & Raymo, 2005) includes Antarctic ice sheet expansion.IRD deposition in Southern Ocean cores also indicate that extensive Antarctic glaciation developed prior to the mPWP (Figure 1) (e.g.Cook et al., 2013;M. A. Hansen et al., 2015;Passchier, 2011;Patterson et al., 2014).Although warming and retreat of the West Antarctic ice sheet recorded in the Amundsen Sea starts from ∼4.2 Ma and extends through the mPWP, more extensive West Antarctic glaciation from MIS M2 ∼3.3 Ma is indicated by ice-proximal cores in the Ross Sea (McKay et al., 2012;Riesselman & Dunbar, 2013).The relatively early onset of Southern Hemisphere δ 18 O benthic change points through the mPWP may reflect ice sheet impacts on Southern Ocean circulation.However, the complexity of the Antarctic signals, showing both advance and retreat of different parts of the ice sheet through time, may have influenced the observed asynchroneity of the δ 18 O benthic change points.We return to these complexities in Section 4.5, where we assess the evidence for changing orbital-scale variability in our PlioVAR data sets and evidence for ice-volume change both during and after the mPWP.

What Were the Characteristics of the iNHG? (Q3)
In this Section, we draw on published records of changes in the continental ice sheets and their interactions with the oceans, to consider the broader context of the patterns of long-term changes we identified.We focus here on the long-term trends; in Section 4.5 we consider how orbital-scale climate variability may have changed with iNHG.First, we consider the evidence recorded in δ 18 O benthic and IRD, which suggest two intervals of ice-sheet growth: during the mPWP (i.e., before ∼2.7 Ma) and an intensification of major ice-sheet growth between ∼2.8 and 2.4 Ma.We then reflect on the evidence for changes occurring elsewhere within the wider climate system to consider whether mechanisms to explain these differences can be detected.

Growth of Continental Ice Prior to ∼2.7 Ma (oNHG)
We have identified a wide range of δ 18 O benthic change points, which spanned our whole interval of interest (3.3-2.4 Ma) and with a wide range of probabilities.This result aligns with previous analysis of individual δ 18 O benthic data using a different approach (Mudelsee & Raymo, 2005).However, this broad interval of change contrasts with the change point analysis of the global LR04 δ 18 O benthic stack which placed the iNHG at 2.73 ± 0.1 Ma (Ruggieri, 2013), and with our analysis of sea-level variability which indicated an interval of globally significant ice growth between Intervals P-2 to P-3 (i.e., around ∼2.7 Ma) (Figure 7).The temporal spread of δ 18 O benthic change points indicates that factors other than global ice volume have contributed to the site-specific δ 18 O benthic shifts.The properties of intermediate and deep waters are initially set by high-latitude surface ocean conditions, and we have identified several δ 18 O benthic during the mPWP which align with some early mid-and high-latitude SST change points (Figure 4).However, the early SST change points are diachronous and of varying probability.Therefore, although there may be an influence of high-latitude SST changes on some of the δ 18 O benthic signals (or vice versa), a direct connection is difficult to establish based on the data presented here.
Alternatively, the asynchronous advance and retreat of smaller ice sheets (Section 4.3) prior to ∼2.7 Ma might account for the difficulties experienced in isolating signals of ice volume response to orbital forcing using δ 18 O benthic (Meyers & Hinnov, 2010).Evidence for growth of relatively small ice sheets includes reduced terrestrial organic matter and increased meltwater inputs to the Gulf of Alaska from ∼3.1 Ma (IODP Site U1417), attributed to an expanding Cordilleran ice sheet which had not yet reached the sea to generate IRD (Sánchez- Montes et al., 2020Montes et al., , 2022)).In the deep North-west Pacific (ODP Site 1208), by extracting the temperature influence from the δ 18 O benthic signal, a gradual increase in global ice volume was determined between ∼3.15 and 2.73 Ma (Woodard et al., 2014).In the absence of extensive Northern Hemisphere Glaciation at this time, a growing Antarctic ice sheet was proposed to account for the >11 m sea-level fall by 2.73 Ma (Woodard et al., 2014).For the same time interval, long-term cooling of ∼2°C in North Atlantic Deep Water (ODP Site 607) in turn suggests that there was cooling in the high latitude (northern) surface ocean where North Atlantic Deep Water is formed (Sosdian & Rosenthal, 2009).North Atlantic SST change points through 3.15-2.99Ma provide support for decreasing mid-and high-latitude SSTs in this region (ODP Site 907, IODP Site U1313, ODP Site 982; 10.1029/2022RG000793 of 43 Figure 4).A challenge is to isolate whether these cooling patterns were conducive to ice-sheet growth, or whether they contributed a cooling signature to the individual records of δ 18 O benthic, potentially masking evidence for the expansion of relatively small and dynamic ice sheets (Mudelsee & Raymo, 2005).

Intensification of Major Continental Ice Sheet Growth ∼2.7-2.4 Ma (iNHG)
The δ 18 O benthic records reveal both an increased amplitude of orbital scale variability from ∼2.7 Ma (P-3, Figure 6) and a clustering of change points ∼2.8-2.7 Ma (Figure 4) which suggest a growing global impact of increasing ice volume in δ 18 O benthic .The timing of these events broadly aligns with the iNHG (Section 1; Haug et al., 2005).Although the change points in δ 18 O benthic are clustered, they span several glacial-interglacial cycles rather than a relatively abrupt and synchronous event (Figure 4).
Asynchroneity in ice-sheet expansion from ∼2.7 Ma is also evidenced by two separate increases in Northern Hemisphere IRD starting from ∼2.7 Ma (Figure 1).It is overly simplistic to relate IRD abundance to glacial extent (Andrews, 2000), and we recognize that the early stages of ice-sheet expansion occur inland without any commensurate IRD deposition.However, abundant and widespread IRD deposition in high-latitude marine settings provides a powerful marker of continental-scale glaciation that has extended to sea level (Kleiven et al., 2002;Shackleton et al., 1984).The first dramatic increase in IRD abundance occurred during MIS G6 (∼2.72 Ma; Figure 1) in sediments from all sectors of the Nordic Seas (Jansen et al., 2000;Jansen & Sjøholm, 1991;Knies et al., 2009), the subarctic northwest Pacific (Bailey et al., 2011) and as far south as ∼53˚N in the subpolar North Atlantic (Bailey et al., 2013;Blake-Mizen et al., 2019;Kleiven et al., 2002).The IRD increase starting from MIS G6 is supported by additional sedimentary evidence demonstrating that it reflected a synchronous expansion of glaciation in Greenland, Scandinavia, the British Isles, and parts of North America, and the advance of ice caps to their iceberg-calving margins (e.g.Chow et al., 1996;Kleiven et al., 2002;Lein et al., 2022;Thierens et al., 2012;Vanneste et al., 1995).A large increase in aeolian dust deposition in the mid-latitude North Atlantic during cold stages from ∼2.72 Ma may also reflect strengthening of glaciogenic dust sources (Naafs et al., 2012) or a shift in regional wind fields and vegetation biomes (Lang et al., 2014) in North America following ice-sheet expansion.
The second large increase in IRD occurs from MIS G4 (∼2.63 Ma, Figure 1) in sites recording ice expansion beyond the circum-Nordic Sea landmasses, suggesting a delayed expansion of ice sheets on North America (as originally proposed by Maslin et al., 1998).The onset of abundant IRD deposition in the Gulf of Alaska ∼2.63 Ma (Sánchez-Montes et al., 2020) lies within the dating error of glacio-fluvial gravels marking the existence of the early and most extensive Cordilleran ice sheet (∼2.64 +0.20 / −0.18 Ma) (Hidy et al., 2013).In the North Atlantic, the first abundant IRD deposition in the center of the Last Glacial Maximum IRD belt (∼50°N at IODP Site U1308) occurs during MIS G4 (Bailey et al., 2010) and later on the southernmost fringe (∼41°N at IODP Site U1313) during MIS 100 at ∼2.52 Ma (Bolton et al., 2010;Lang et al., 2014).
The diachronous expansion of IRD deposition in the subpolar North Atlantic between 2.72 and 2.52 Ma is not a product of iceberg survivability, because significant glacial excursions of Subarctic Front surface waters into the mid-latitude North Atlantic occurred earlier, during cold stages from MIS 104 (Bolton et al., 2018;Hennissen et al., 2014) and perhaps as early as MIS G6 (Naafs et al., 2010).Iceberg calving models for the Last Glacial Maximum suggest that the source of IRD to this region of the North Atlantic is the Gulf of St. Lawrence in midlatitude North America (Bigg & Wadley, 2001).The relatively late arrival of IRD deposition at IODP Site U1313 therefore likely reflects the first time during the iNHG that the proto-Laurentide Ice Sheet extended into the mid-latitudes of North America, consistent with other evidence for the timing of the onset of midlatitude glaciation here (Balco & Rovey, 2010;Shakun et al., 2016).

Climate Changes Accompanying the Intensification of Major Continental Ice Sheet Growth ∼2.7-2.4 Ma (iNHG)
The range of change points and the diverse signatures of orbital-scale variability in δ 18 O benthic (Figure 6) are a reminder that factors other than ice volume are also required to explain the differences we observed between sites, which could include temperature or salinity-driven changes in seawater δ 18 O in the regions of deep water formation.Only one SST record has a change point between 2.9 and 2.6 Ma (Figure 4; G. bulloides Mg/Ca from North Atlantic site IODP U1313).Yet, SST and dust records from both hemispheres detail signatures of temperature changes and intensification of westerly wind systems in the mid-and high-latitudes during this iNHG interval.The reconstructed changes to North Pacific surface ocean conditions were conducive to preconditioning Northern Hemisphere ice-sheet growth, and impacted the configuration of deep ocean circulation.Warmer SSTs from ∼2.7 Ma onwards (North Pacific sites ODP 882 and IODP U1417) would have provided a source of moisture to enhance the growth of the North American ice sheets, alongside the cooler winter SSTs and associated sea-ice expansion (Haug et al., 2005;Sánchez-Montes et al., 2020).These patterns reflect the development of the halocline in the North Pacific, whereby a strong vertical gradient in salinity leads to isolation of deep waters from the surface, encouraging elevated summer warming but also potentially enabling sequestration of atmospheric CO 2 into the deep ocean (Haug et al., 2005).Model-data comparison suggests that as this strong halocline developed, there was a cessation of the North Pacific Deep Water formation which occurred during the Pliocene (Burls et al., 2017;Ford et al., 2022).
Through MIS G6 (∼2.72 Ma), an increase in deep ocean connectivity between the Pacific and Atlantic Oceans is evidenced by the convergence of North Atlantic and North Pacific δ 18 O benthic values (Woodard et al., 2014).The pronounced cooling of North Atlantic Deep Water ∼2.7 Ma (Site 607) marked the end of a long-term trend started in the mPWP (Sosdian & Rosenthal, 2009;Woodard et al., 2014).This deep water cooling is likely linked to the acceleration of surface ocean cooling in the North Atlantic at ∼2.7 Ma (ODP Site 982, Lawrence et al., 2009) and an expansion of sea ice, which reached its maximum winter extent by ∼2.6 Ma in the Fram Strait (Clotten et al., 2018;Knies, Cabedo-Sanz, et al., 2014).
While the Northern Hemisphere showed a general pattern of cooling over the ∼2.7 Ma iNHG interval, the Subantarctic Atlantic was characterized by an interval of relative warmth (ODP Site 1090, Martinez-Garcia et al., 2010).The intermediate waters which are formed in the Subantarctic surface ocean were also impacted.Change points in δ 18 O benthic from sites bathed by Antarctic Intermediate Water at ∼2.7 Ma (Figure 4) reflect both a shift in the mean and a range of variability (Figure 6) in the Subantarctic surface ocean.In the South-west Pacific (DSDP Site 593), reconstructed Antarctic Intermediate Water temperatures confirm that a pronounced cooling occurred in the zone of intermediate water formation during MIS G6, set within an extended interval of surface ocean cooling in the mid-latitude South-west Pacific ∼2.8-2.6 Ma (McClymont et al., 2016).
As the mid-and high latitudes cooled and climate became more variable, there was an intensification of wind-driven dust deposition in the mid-latitudes of both hemispheres and across multiple ocean basins from MIS G6 at ∼2.73 Ma (Abell et al., 2021;Lang et al., 2014;Naafs et al., 2012).Both the marine sediment data and climate models suggest that globally synchronous equatorward shifts in the westerlies occurred in response to stronger atmospheric temperature gradients as the high latitudes cooled and ice volume increased (Abell et al., 2021;X. Li et al., 2015).In North America, MIS G6 marks a shift from a wetter-than-modern continental climate during the Pliocene to a more arid climate for the Pleistocene (Lang et al., 2014).A long-term increase in aridity in the lower latitudes has also been recorded from ∼2.7 Ma and explained as a response to a globally cooler climate with reduced evaporation supplying monsoon rains (e.g.X. Li et al., 2015).In central Africa, a shift toward more arid conditions was followed by an increase in dust fluxes associated with the oNHG, suggested to reflect strengthening winds in response to the increasing latitudinal temperature gradients driven by ice-sheet growth (Crocker et al., 2022;de Menocal, 2004).Increasing aridity in North-east and South-east Africa also developed with iNHG (Liddy et al., 2016;A. K. Taylor et al., 2021), linked to long-term cooling in east Indian Ocean SSTs (A.K. Taylor et al., 2021).
Although we identified a wide range of climate changes aligned with the iNHG, we do not observe regime shifts (change points) or changes to variability in the low-latitude regions ∼2.7 Ma, consistent with a previous proposal that the iNHG is a transition which largely reflects high-latitude climate change (Ravelo et al., 2004).Even some mid-latitude sites show little long-term SST response to the iNHG ∼2.7 Ma, for example, MIS G6 is relatively cool at North Atlantic Site IODP U1313, but otherwise forms part of a gradual long-term cooling which began ∼3.1 Ma (Naafs et al., 2012).In the Southern Hemisphere, long-term stability or gradual cooling is observed in SSTs at several sites: in the Benguela upwelling system (Petrick et al., 2018;Rosell-Melé et al., 2014), the subantarctic South Atlantic (Martinez-Garcia et al., 2010), and equatorial Pacific (Lawrence et al., 2006).An impact of sea-level change (and orbital forcing) on the Indonesian Throughflow was demonstrated by an intensification of the Leeuwin Current between 2.9 and 2. 6 Ma,but this was strongly influenced by ice sheet expansion (De Vleeschouwer et al., 2022) and so has a strong link to high-latitude changes.
The stronger signals of ice-sheet expansion and mid-and high-latitude climate changes with iNHG are consistent with a fall in pCO 2 increasing the potential for ice sheets to expand into the mid-latitudes.A puzzle is to explain 10.1029/2022RG000793 of 43 why pCO 2 fell below mPWP values during MIS G10 (∼2.8 Ma) whereas IRD increases are detected from MIS G6 (∼2.72 Ma) (Figure 1).One possibility is that with the time resolution and spatial distribution of our existing records, we have not yet fully identified the expression of climate changes and ice-sheet evolution across these relatively short glacial intervals, including the pCO 2 reconstructions.Alternatively, the second lowering of pCO 2 during MIS G6 may have been required to facilitate the development of major ice-sheet growth (Figure 1), or a temperature threshold for positive ice-mass balance may finally have been crossed.Although much of our focus has been on the Northern Hemisphere ice sheets, several climate patterns associated with iNHG also suggest a climate system response to changes in the Antarctic ice sheet, including polar water expansion in both hemispheres (Haug et al., 2005;Martinez-Garcia et al., 2010) and shifting Southern Hemisphere circulation patterns (e.g.De Vleeschouwer et al., 2022), which might also have been connected to the monsoons, Mediterranean Outflow, and North Atlantic circulation (Sarnthein et al., 2018).In turn, it has been suggested that the increased Mediterranean Outflow from ∼2.9 to 2.7 Ma might have exerted a negative impact on Northern Hemisphere ice sheet growth, delaying a response to North Atlantic cooling (Kaboth-Bahr et al., 2021).
Finally, our analysis identified a last clustering of change points focused largely on the low latitudes (<20°N/S) in both SST and δ 18 O benthic data, which occurred after ∼2. 6 Ma and seem to be linked to the development of the first large glacial cycles (MIS 100-96; Figure 4e).The drivers of these final change points are unclear.Several are clustered around the Indo-Pacific warm pool (e.g.DSDP Site 214, ODP Site 806, ODP Site 1143, ODP Site 1148) or reflect changing zonal gradients in the Pacific Ocean (e.g. between ODP Site 806 and ODP Site 846).We acknowledge that caution is required since the Mg/Ca change points have very low probability (Figure 4), but these low-latitude signals may be hinting at a late response to iNHG or even an independent evolution of the Walker Circulation system.Previously, Kaboth-Bahr and Mudelsee (2022) identified high spatial and temporal variability in change points spanning the Pliocene-Pleistocene from sites influenced by Walker Circulation, and suggested that different ocean basins affected by this complex system evolved asynchronously from as early as 3.5 Ma.Our analysis supports this previous work, and shows that many of the high-latitude changes associated with iNHG occurred much earlier (before ∼3.0 Ma), whereas the low latitudes tend to show a later long-term response aligned with the development of pronounced glacial maxima from MIS 100.

Did the Amplitude or Pacing of Climate Variability Change With iNHG? (Q4)
In this Section we investigate climate changes which occurred on orbital timescales (tens to hundreds of thousands of years), to consider whether iNHG impacted the relationships between climate forcings and feedbacks for different components of the climate system.We have identified that the mid-and high-latitudes had higher SST variability than the low-latitudes, which increased in association with the iNHG (Figure 5).Here, we explore the pacing of this variability, because it has been proposed that with the growth of larger Northern Hemisphere ice sheets, two impacts on the climate response to orbital forcing were experienced.First, that larger ice sheets were slower to respond to orbital forcing (Lawrence et al., 2009;Meyers & Hinnov, 2010), and second, that the ice sheets enhanced climate feedbacks related to changing ice albedo, meridional temperature gradients and winds (e.g.Herbert et al., 2015;Martinez-Boti et al., 2015).In particular we consider the relative importance of climate cycles related to obliquity pacing (∼41 kyr) and precessional pacing (∼19-23 kyr; see Section 1) by comparing our results to the published results of site-or proxy-specific time-series analysis.
In contrast, precession signals dominate SST and δ 18 O benthic records from the Mediterranean, reflecting a tight connection between local insolation and ventilation of bottom water masses (Herbert et al., 2015;Khélifi et al., 2014).Precession signals also dominate multiple low-latitude SST records (Herbert et al., 2010) and the wind-and river-driven delivery of African sediments to marine sediment sequences offshore, linked to changing African monsoon strength (Crocker et al., 2022;de Menocal, 2004).Significant precession signals are also found 10.1029/2022RG000793 of 43 in the high latitudes: in ice-marginal records from the marine-based East Antarctic ice sheet margin (Bertram et al., 2018;Cook et al., 2013;Patterson et al., 2014;Williams et al., 2010), and records of Antarctic Ice Sheet variability from the Scotia Sea (South Atlantic) (Reilly et al., 2021).Precession signals are also reflected strongly in the New Zealand sea level record, which is paced by Southern Hemisphere local insolation changes (Grant et al., 2019).The presence of both strong obliquity and precession signals indicates that even within the Antarctic Ice Sheet there were likely different regional catchment sensitivities to orbital forcing (Golledge et al., 2017), which might explain the challenges of isolating a strong link between orbital forcing and the global δ 18 O benthic stack prior to iNHG (Meyers & Hinnov, 2010).
From ∼2.7 to 2.6 Ma (MIS G6-G4; Interval P-3), an increasing influence of obliquity pacing appears in several components of the climate system, and argued to reflect strengthening feedbacks associated with the expansion of continental ice sheets into the mid-latitudes (Section 4.3) (e.g.Herbert et al., 2015).In turn, these stronger connections may have strengthened the relationship between the LR04 δ 18 O benthic stack and orbital forcing, which then persisted until ∼1.3 Ma (Meyers & Hinnov, 2010).For example, in the Mediterranean Sea, a shift away from precession-driven local insolation changes and the development of enhanced cooling during glacial stages with an obliquity period develops from ∼2.8 Ma (Herbert et al., 2015).Obliquity continues to dominate dust deposition to the North Atlantic, but the relative timing shifts: having led δ 18 O benthic in the mPWP and late Pliocene, dust varies in-phase with δ 18 O benthic by ∼2.8 Ma, suggesting a close link between dust supply and glaciation (Naafs et al., 2012).
However, we acknowledge that not all regions respond to the events marking the iNHG (MIS G6-G4) that is, the intensification of glacial cycles at 2.7-2.6 Ma (Section 4.4).Continued complexity in regional glacial-interglacial variability signals after 2.6 Ma is indicated by several records where precession, rather than obliquity, continues to dominate.These include Subantarctic Atlantic SSTs (ODP Site 1090; Martinez-Garcia et al., 2010), and meltwater discharges from a proto-Laurentide (North American) ice sheet into the Gulf of Mexico (Shakun et al., 2016).The precession-paced meltwater discharges are proposed to reflect a strong ice sheet mass balance response to local summer insolation (Shakun et al., 2016), yet the North American-sourced wind-driven dust deposition to the mid-latitude North Atlantic continued to be paced by obliquity and glaciation-driven sediment availability (Naafs et al., 2012).This apparent discrepancy in the observed pacing of ice-sheet dynamics during the iNHG might be explained by the fact that dust generation on North America during this time was dominated by non-glaciogenic processes (Lang et al., 2014).Alternatively, it may simply reflect that glaciogenic processes driving dust generation on North America during the iNHG operated at high latitudes under the influence of changes in axial tilt, whereas meltwater discharges from the mid-latitude margins were more strongly influenced by precession.This example demonstrates that complex and sometimes apparently contradictory evidence for climate pacing and forcing/response can be determined even in single regions, but further emphasizes the need for multi-proxy and multi-region assessments of climate evolution associated with the iNHG.In summary, we have found evidence for ice sheet growth, high latitude ocean changes, and low-high latitude teleconnections developing on both orbital and longer-term timescales through the mPWP and late Pliocene.The complex mixture of different orbital-scale variability being represented by different regions and climate variables, and their evolution through the mPWP and late Pliocene, may explain why the identification of common signals of climate and ice sheet evolution through this time window were not apparent in our change point analysis.The signals we observe also support the proposal that as well as this time interval being one with a stochastic character, the expansion of continental ice sheets with iNHG strengthened some low-high latitude climate feedbacks, and led to a more stable response to orbital forcing (Meyers & Hinnov, 2010).Yet even the Antarctic ice sheet showed variable regional responses to insolation, making the oNHG and iNHG challenging intervals to characterize, and perhaps an unexpectedly complex time window given the past research focus on the mPWP as an interval of equilibrium climate response to modern and near-future pCO 2 .
We observed that some regions of the climate system, particularly the low latitudes, showed little response to the onset of major ice-sheet expansion ∼2.7 Ma (iNHG).This pattern contrasts with the more pronounced low-latitude cooling signatures which observed during the Mid and Late Pleistocene glacial stages (e.g.Lawrence et al., 2006).As the impact of CO 2 forcing is highest in the low latitudes (Rohling et al., 2021), the relatively muted response to oNHG in these regions suggests that only a small decrease in pCO 2 occurred at this time.In contrast, the stronger expression of change in the mid-and high-latitudes is consistent with strengthening ice-albedo feedbacks which amplified the response to pCO 2 (Martinez-Boti et al., 2015;Rohling et al., 2012).
By combining these observations of various aspects of late Pliocene and early Pleistocene climate change, we have identified varying influences of ocean gateways, radiative forcing associated with pCO 2 , ocean circulation and ice-sheet feedbacks.The early changes observed in sites which have connections to the Indonesian Throughflow, in the absence of any clear changes to pCO 2 , suggest that changes to this ocean gateway may have had regional climate impacts prior to and during the mPWP.These events are accompanied by relatively localized and likely asynchronous patterns of growth in smaller ice masses, leaving a complex signature in both δ 18 O benthic and SST signals.It is uncertain whether the ice sheets in both hemispheres were undergoing a long-term expansion during this interval, whereby shifts in ocean gateways and ocean circulation may have pre-conditioned a later expansion (iNHG), so that a small change in forcing led to a large ice-sheet response.Alternatively, the ice masses may have remained relatively small, in which case a much larger forcing would have been required to cause the larger expansion of continental ice defined by iNHG.Ice-proximal marine sediment records may offer some additional clues (see Section 4.6.2),but where ice margins remained on land there is a need to generate terrestrial data sets which can detect evidence for changes to air temperature, precipitation or ice margin position.Although some data exist (e.g.Hidy et al., 2013) it is generally sparsely distributed.An increase in terrestrial data density would also be valuable for data-model comparisons in warmer climates of the past (e.g.Haywood et al., 2020).
The later, and more abrupt, iNHG around 2.8-2.4 Ma has a stronger signal of northern hemisphere ice-sheet expansion as well as impacts on SSTs.The preceding fall in pCO 2 and latitudinal expressions of SST change are consistent with an ice-sheet response to reduced radiative forcing.The delayed response of the ice-sheets to the initial fall in pCO 2 at 2.8 Ma suggests either a role for slow ice-sheet feedbacks during growth, a need for further pCO 2 decline to have an impact, or some combination of the two.The asynchronous expansion of ice masses around, and then beyond, the Nordic Seas region suggests that there were different regional responses to these scenarios.Our analysis shows that once these larger ice masses grew, their impact on other parts of the climate system also grew, with some low latitude regions showing a response depending on the nature of the teleconnections at play.

Reflections and Future Work
In this work, we synthesized a large set of SST data and assessed their long-term trends and orbital-scale variability.The development of similar multi-proxy, spatially-distributed temperature data sets is also required to generate robust estimates of global mean surface temperature (GMST), which can be informative in data-model comparisons and calculations of climate sensitivity to changing radiative forcing (e.g.Inglis et al., 2020;Martinez-Boti et al., 2015).GMST estimates from MIS KM5c were derived by calculating a globally weighted mean temperature anomaly (see McClymont et al., 2020).However, this approach is sensitive to sampling density and does not account for unevenly distributed proxy data sets.Future studies can assess the influence of this approach by resampling the modern temperature field but with the sampling density of the Pliocene and Pleistocene.This has been assessed during other time intervals (e.g, the early Eocene) and indicates that this method is able to estimate GMST within 1.5°C (Inglis et al., 2020).Although δ 18 O benthic values have also been used to estimate GMST at high-resolution during the Pliocene (e.g., J. Hansen et al., 2013), this approach is subject to several uncertainties (e.g., changes in ocean vertical stratification; discussed by Inglis et al., 2020).Furthermore, using δ 18 O benthic values in this way assumes that they represent a global climate signal.Our study shows a wide range in δ 18 O benthic change points starting from 3.3 Ma, suggesting that caution is required when using δ 18 O benthic values to infer GMST, especially between glacial-interglacial cycles.
A limitation of our work was that to understand ice-sheet evolution, we relied on the expected global signature of ice volume from δ 18 O benthic records, alongside published IRD evidence for ice sheets whose margins reach the sea.The diversity of the glacial-interglacial variability we identified between sites emphasizes the need to better understand both the global and local factors influencing the δ 18 O benthic signal.More studies of the temperature-controlled Mg/Ca signal in benthic foraminifera would allow us to start to evaluate the varying environmental contributions to δ 18 O benthic records.We identified some differences between the timing of changes recorded by δ 18 O benthic and IRD (e.g.longer-term shifts in δ 18 O benthic starting in the mid-Pliocene, a more narrowly focussed shift in Northern IRD ∼2.7-2.6 Ma).While several deep-sea ice-marginal records have provided continuous archives of ice-sheet variability, complexities exist in assessing the orbital pacing of ice sheets from these records.In part, this is because of the difficulty in separating out the specific processes controlling individual ice-sheet proxies (e.g.McKay et al., 2022).Detailed lithofacies analysis is required to ensure depositional processes are carefully considered, since otherwise there can be confusion even between studies undertaken at near-by sites (e.g.M. A. Hansen et al., 2015;Patterson et al., 2014).
For several regions of the ocean (e.g. the Indian Ocean, the Pacific sector of the Southern Ocean, the North Pacific Ocean), a paucity of published data also limits assessment of key systems for regulating regional and global climate, including monsoons and low-latitude hydroclimate, circulation systems connecting the low and high latitudes, and ocean-ice sheet interactions close to the ice sheets.Promising new high-resolution sequences recovered from the northern Indian Ocean by IODP will allow examination of the core monsoon region of the northern Bay of Bengal and Andaman Sea (Clemens et al., 2016).In the South-west Indian Ocean, IODP Expedition 361 has recovered sediment sequences offering new opportunities to reconstruct the transport of heat and salt from the low-to mid-latitudes and between the Indian and Atlantic Ocean, and to link them to onshore environmental changes in southern Africa (e.g. A. K. Taylor et al., 2021).A relative scarcity of data from the South Pacific sector of the Southern Ocean will soon be resolved as a result of IODP Expedition 383 (Lamy et al., 2021).Results from recent expeditions adjacent to the Antarctic Ice Sheet margin (IODP Expeditions 374, 379, 382) (Gohl et al., 2021;McKay et al., 2019;Weber et al., 2021), and scheduled expeditions adjacent to the Greenland margin (IODP Expedition 400) and close to Svalbard (Expedition 403) are also expected to shed new light on regional variability over high-latitude oceans and their proximal ice sheets.These expeditions will in future help to better inform the scientific community on the role of natural forcing and thresholds that may exist in warmer climates like the Pliocene Epoch compared to the relatively cool and more heavily glaciated Pleistocene Epoch.Reflecting on our analysis here, we also encourage the development of co-ordinated and strategic multi-proxy sampling strategies and generation of the highest quality age controls, to maximize the information which can be recovered and explored to better understand the iNHG and the warmth of the mPWP which immediately preceded it.
Finally, our focus on marine evidence to examine the iNHG neglects the rich archive of environmental information which can be returned from continental sequences, but whose age controls can be more challenging in the absence of time-continuous sedimentation or the global temporal framework of a δ 18 O benthic stack.Yet climate feedbacks linked to changes in vegetation and soils are important components of Earth system change (e.g.Feng et al., 2022).To maximize the potential for linking ocean, ice, and terrestrial information and gain a holistic understanding of the iNHG will require further analysis for continental environmental changes using both continental inputs to marine sediment sequences (cf., Crocker et al., 2022) and the generation and integration of continental sequences into a global framework (e.g.Feng et al., 2022).We also did not examine marine ecosystem responses to the environmental changes we identified, yet our interval of study includes evidence for diachronous extinction events (e.g.Brombacher et al., 2021;De Schepper & Head, 2008), and the mPWP offers an opportunity to consider how biota responded to a warmer climate with elevated pCO 2 (e.g.Todd et al., 2020).

Conclusions
By compiling, evaluating and exploring the high-resolution reconstructions of ocean temperatures, benthic and planktonic stable oxygen isotope signals (δ 18 O benthic , δ 18 O planktonic ), and sea level we have identified complex signals of ice sheet and climate evolution between 3.3 and 2.4 Ma associated with the late Pliocene and earliest Pleistocene "onset and intensification of Northern Hemisphere Glaciation" (oNHG and iNHG).Here, we directly respond to the questions we posed in Section 1.
Q1: Are there proxy-specific differences in the evolution of ocean temperatures?We did not identify any proxy-specific signatures where high-quality multi-proxy SST data were available.We confirmed that caution is required when considering whether data are recording mean annual SST: seasonality, depth of production, and how these two factors have evolved through time, as well as the possibility of non-thermal influences, are all important considerations when synthesizing and interpreting proxy data.The existing resolution of single-site multi-proxy data prevented us from fully addressing this question.
Q2: What were the characteristics the mPWP (∼3.264-3.025Ma)? Immediately preceding the mPWP was the MIS M2 glacial stage, argued to reflect a short-lived increase in Southern Hemisphere glaciation and some minor ice-sheet expansion in the north.Our analysis shows that the M2-KM5c transition was marked by a reduction in global ice volume, and mid-and high-latitude ocean warming of 1-2 and 2-4°C, respectively.Low latitude ocean temperature response to the M2-KM5c transition was comparatively muted and within proxy calibration uncertainties (<1°C).
The mPWP was marked by warmer-than-Pre-Industrial climates, consistent with reduced global ice volume and elevated CO 2 .It was also a time of climate evolution: sectors of the Antarctic ice sheet were expanding, and there is some evidence for growth in the smaller Northern Hemisphere ice sheets.There is a consistent pattern of increasing climate variability at the 41 kyr-obliquity scale moving from the low to the mid-and high-latitudes.Our analysis thus confirms that targeting a single interglacial for data-model comparison offers the best chance to minimize the influence of orbital-scale high-latitude ocean temperature variability for data-model comparison (Haywood et al., 2020;McClymont et al., 2020).In the absence of high-resolution pCO 2 reconstructions before the mPWP, it is difficult to assess the role played by radiative forcing in the long-term trends.The observed patterns are consistent with ocean circulation responses to gateway-driven changes in Indonesian Throughflow and long-term intensification of meridional temperature gradients.

Q3: What were the characteristics of the iNHG?
We have shown that the iNHG was not a globally synchronous event.We confirmed previous suggestions that long-term changes in climate and ice sheets began before and during the mPWP (a longer-term "oNHG").Asynchronous transitions in δ 18 O benthic records from individual sites confirm the importance of identifying and assessing the relative importance of local influences over δ 18 O benthic (e.g.temperature, salinity) which are in turn linked to high-latitude surface ocean properties, alongside the more integrated signature of global ice volume.
Our review of existing data confirms that the iNHG was complex in terms of ice-sheet advances/retreats, in both timing and magnitude.A global expression of ice volume change is not indicative of how dynamic the individual ice masses were across the complete oNHG and iNHG.The pattern of high-latitude cooling proceeding first, followed by an expansion of ice towards the mid-latitudes in the northern hemisphere, is consistent with intensification of ice-albedo feedbacks which amplified the decline in CO 2 and enabled ice sheets to expand and survive further to the south.The lag between the fall in CO 2 during MIS G10 and the subsequent onset of IRD and glacial-stage intensification (Figure 1) may reflect the influence of changes in North Atlantic circulation (e.g.Kaboth-Bahr et al., 2021) or the time required to grow and sustain a larger continental ice mass, particularly in North America.

Q4: Did the amplitude or pacing of climate variability change during iNHG?
Yes, but in a complex and non-uniform way.In almost all SST and δ 18 O benthic records we examined, we observed an increase in the mean or range of variability at the 41 kyr-obliquity frequency from the mPWP through to the early Pleistocene, but there was strong regional scale variability.The low-latitude SST records show minimal or no shifts, whereas larger increases in SST variability emerge in the mid-and high-latitude Northern Hemisphere oceans.There was diversity in the signals of δ 18 O benthic which may reflect evolving surface ocean conditions in the regions of intermediate and deep water formation, though no clear patterns could be isolated.In some proxy records, previous work has established switches from precession to obliquity in response to strengthened low-high latitude teleconnections as the Northern Hemisphere ice sheets expand, but not all records or regions express this pattern.
Our analysis has demonstrated the complex, non-uniform and globally asynchronous nature of climate changes associated with the oNHG and iNHG.The evidence available today indicates that shifting ocean gateways and ocean circulation changes may have pre-conditioned the later evolution of ice sheets with falling pCO 2 , but there remains complexity in confirming the mechanism(s) or sequence of events which influenced the timing of ice sheet expansion and its expression in marine records.We recommend that future work targets multi-proxy data sets with high time-resolution, which can integrate our understanding of the environmental processes controlling and responding to the changes in SST and ice volume outlined here, which will continue to be critical to develop our understanding of the causes and impacts of past warm climates and the transitions between different climate regimes.group on Pliocene Climate Variability over glacial-interglacial timescales (PlioVAR), and we thank all workshop participants and members of the PlioVAR steering committee for their discussions.Funding support has also been provided by the UK Natural Environment Research Council (PA, HLF [NE/N015045/1]), the Leverhulme Trust (Philip Leverhulme Prize, ELM), the Ministry of Science and Technology in Taiwan (SLH, Grants 110-2611-M-002-028-, 111-2611-M-002-027-), MICIN (MAG, the Grant PICTURE, PID2021-128322NB-I00), a Royal Society Dorothy Hodgkin Fellowship (GNI, DHF\R1\191178), the Research Council of Norway (BR, Grant 221712), the German Research Foundation (JG, project DePac-GR 3528/8-1), the Agencia Nacional de Investigación y Desarrollo, Anillo (CK, Grant ACT210046), the Universidad de Santiago de Chile (CK, Grant DICYT 092112KSSA), a UKRI Future Leaderships Fellowship (BH, MR/S034293/1), the USA National Science Foundation (JA, Grant NSF-OCE-PRF #2126500) and the Agence Nationale de la Recherche in France (CTB, ANR-16-CE01-0004) and the IODP in the UK (PA) and France (CTB).This research used samples and/ or data provided by the International Ocean Discovery Program (IODP), Ocean Drilling Program (ODP), and Deep Sea Drilling Project (DSDP).We also acknowledge PANGAEA and the NOAA Paleoclimatology databases for free access to published data.The work presented here resulted from the efforts of the late Antoni Rosell-Melé to generate a community effort to facilitate a better understanding of Pliocene and Pleistocene climate evolution: we thank him for all of his work and enthusiasm getting PlioVAR up and running, and for contributing to our discussions and assessments over the years that followed.

Figure 1 .
Figure 1.Key climate parameters across the Pliocene-Pleistocene transition.The uppermost panel delineates the transition from the Pliocene Epoch (3.6-2.58Ma) to the Pleistocene Epoch (2.58-0.01Ma) and the three time intervals used for our statistical analysis (P-1 3.3-3.0Ma, P-2 3.0-2.7 Ma, and P-3 2.7-2.4Ma; see Section 2.3).The time intervals defined as the onset of Northern Hemisphere Glaciation (oNHG), mid-Piacenzian warm period (mPWP) and the intensification of major Northern Hemisphere Glaciation (iNHG) are marked on panel (b), alongside the M2 glacial stage which has been suggested as a possible harbinger of more extensive glaciation (Westerhold et al., 2020).(a) Reconstructed atmospheric pCO 2 concentrations (de la Vega et al., 2020) including the ∼280 ppmv threshold for extensive glaciation proposed by De Conto et al. (2008) which aligns with Pre-Industrial pCO 2 recorded in ice cores(Loulerge et al., 2008); (b) a global stack of benthic oxygen isotope ratios(Ahn et al., 2017) which reflect a combined signal of the deep-water temperature and ice volume (see Technical Box).Key marine isotope stages referred to in the text are marked, blue for the cold glacial stages, red for the warm interglacial stage KM5c; (c) North Atlantic ice-rafted debris (IRD) inputs at Sites U1307 (southern Greenland margin) and U1313 (southern North Atlantic)(Blake-Mizen et al., 2019;Lang et al., 2014), which may be used to indicate when ice sheets expanded and reached the sea.(d) Antarctic IRD inputs from Site U1361(Patterson et al., 2014).MAR refers to the Mass Accumulation Rate of the sediment.

Figure 2 .
Figure 2. The PlioVAR-synthesized records of climate change spanning the 3.3-2.4Ma interval, superimposed upon present-day mean annual sea-surface temperatures (SSTs).Four climate proxies for ocean temperature were assessed:  U K ′ 37 , TEX 86 , Mg/Ca in planktonic foraminifera, and δ 18 O in planktonic foraminifera.Benthic δ 18 O records were analyzed to assess changes in deep water temperature and global ice volume.The three regions which are discussed when comparing SST proxies (Section 3) are annotated, noting that upwelling sites are found in the Atlantic, Indian and Pacific Oceans.Not shown are the three sea-level records, since they assess global stacks of available δ 18 O benthic .Sites which did not meet the PlioVAR thresholds for age control or temporal resolution of data are indicated in gray.

Figure 3 .
Figure 3. Regional comparisons of glacial and interglacial means for selected marine isotope stages (MIS) spanning the late Pliocene and early Pleistocene.Only the MIS with the highest frequency of sea-surface temperature (SST) data and where the glacial or interglacial duration was ≥20 kyr were selected.Data are presented as 20 kyr means centered on the glacial or interglacial maximum as defined by Lisiecki and Raymo (2005), with the 2σ range of SSTs within that 20 kyr interval shown by the vertical error bars.For data sources see Table 1, for site locations see Figure 2. (a) δ 18 O benthic stack (Lisiecki & Raymo, 2005) and the MIS analyzed in panels (a-c) where even numbers (blue shading) highlight glacial stages, and odd numbers (orange shading) highlight interglacial (warm) stages.(b) North Atlantic sites.Note that the Site SPP (Mediterranean Sea) lies outside the North Atlantic but its SSTs are likely to be closely connected to changes there; (c) Low-latitude sites which are unaffected by upwelling ("warm pools"); (d) Low-latitude upwelling sites; (d) Temperatures are reconstructed using  U K ′ 37 (10 sites, squares) and foraminifera Mg/Ca (5 sites, triangles).

Figure 4 .
Figure 4. Change point analysis of the four high time-resolution ocean proxies spanning the 3.3-2.4Ma interval.(a) Summary map of change point results mapped by proxy.Symbol size represents the probability (large = high probability); (b) benthic δ 18 O stack (Ahn et al., 2017) with key marine isotope stages (MIS) indicated as per Figure 2; (c) Change points in sea-surface temperature (SST) records plotted by time and latitude; (d) Change point probabilities in the SST records plotted by time; (e) Change points in benthic and planktonic δ 18 O plotted by time and latitude; (f) Change point probabilities in benthic and planktonic δ 18 O records plotted by time.
Globally integrated time series of sea-level using benthic δ

Figure 5 .
Figure 5. Violin plots (symmetrical vertical histograms) of glacial-interglacial sea-surface temperature (SST) amplitudes for each site, calculated as the maximum range within a moving 41-kyr window along the duration of each period.SSTs recorded by  U K ′ 37 (u on the x-axis), Mg/Ca in planktonic foraminifera (m), and TEX 86 (t).Sites are ordered by latitude, from northernmost (left) to southernmost (right).Site locations are shown in Figure 2 and include Atlantic, Indian, and Pacific Ocean records.(a) a comparison of the evolution of mean SST variability across the three time intervals, extracted from panels (b-d); (b-d) frequency distributions of SST amplitudes within each time interval (colors).Box-whisker plots are superimposed on the frequency distribution data: means are joined between sites by the solid line.The 16th to 84th percentiles are indicated by the white vertical box, and the full range of data is shown by the shading for each site.For the frequency distributions which include the marine isotope stages M2 stage in panel (a) see Figure S5 in Supporting Information S1.
the Mean Sea-Surface Temperature (SST) Differences Between the Three Time Intervals of Interest, Using SSTs Reconstructed by the ′ 37 Index, TEX 86 Index, or Mg/Ca in Three Planktonic Foraminifera (G.bulloides, (MgCa-b), T. sacculifer (MgCa-s), and G. ruber (MgCa-r)) and a larger shift from ∼2.7 Ma.However, the δ 18 O benthic stack masks the diverse patterns of in individual δ 18 O benthic records, which are not easily grouped by latitude or water depth (Figure6).For example, at 2,082 m water depth in the South Atlantic, ODP Site 1088 consistently showed lowest mean variability in δ 18 O benthic despite recording a wide range.In contrast, the deep ODP Site 849 (3,296 m water depth, equatorial Pacific) showed consistently low variability in both the mean and range throughout all of the time intervals we studied.

Figure 7 .
Figure7.Violin plots (symmetrical vertical histograms) of glacial-interglacial sea-level amplitudes for each site, calculated as the maximum range within a moving 41-kyr window along the duration of each period, plotted using the same approach as for Figure5.The records shown were generated using different approaches to extract sea-level change from the global δ 18 O benthic stack(Berends et al., 2021;Miller et al., 2020;Rohling et al., 2014) as outlined in Section 2.5.

Table 2
Sites Used for δ 18 O benthic Probability Distribution Analysis