Latest Pliocene Northern Hemisphere glaciation amplified by intensified Atlantic meridional overturning circulation

The global climate has been dominated by glacial–interglacial variations since the intensification of Northern Hemisphere glaciation 2.7 million years ago. Although the Atlantic meridional overturning circulation has exerted strong influence on recent climatic changes, there is controversy over its influence on Northern Hemisphere glaciation because its deep limb, North Atlantic Deep Water, was thought to have weakened. Here we show that Northern Hemisphere glaciation was amplified by the intensified Atlantic meridional overturning circulation, based on multi-proxy records from the subpolar North Atlantic. We found that the Iceland–Scotland Overflow Water, contributing North Atlantic Deep Water, significantly increased after 2.7 million years ago and was actively maintained even in early stages of individual glacials, in contrast with late stages when it drastically decreased because of iceberg melting. Probably, the active Nordic Seas overturning during the early stages of glacials facilitated the efficient growth of ice sheets and amplified glacial oscillations. Active overturning in the Nordic seas and strong Atlantic circulation coincided with enhanced glaciation in the Northern Hemisphere after 2.7 Ma, according to a multi-proxy investigation of cores from the IODP site U1314.

N orthern Hemisphere glaciation gradually developed after its onset 3.6 million years ago (Ma) 1 and further intensified after 2.7 Ma 2 . The latter intensification of Northern Hemisphere glaciation (iNHG) was an irreversible climatic deterioration 3 , and the 41-kyr component of glacial-interglacial variations was amplified 2 . The primary cause of iNHG has been debated, with two leading hypotheses: severe cooling induced by decreased atmospheric CO 2 concentration 4-6 , and increased rain and snow fall over Greenland and Europe attributed to the closure of the Central American Seaway (CAS) 7,8 . Maintaining the 41-kyr glacial cycles requires additional oceanic amplifier mechanisms, such as the Atlantic meridional overturning circulation (AMOC) 3 and storage of carbon dioxide 9 . The meridional overturning circulation is driven by deep waters produced by sinking of cold, salty, dense waters at high latitudes in both hemispheres, and in this process, Atlantic surface currents transport heat and moisture northward 10 . Although the AMOC played a dominant role in recent glacial-interglacial and stadial-interstadial climate conditions 11,12 , studies based on benthic foraminiferal carbon isotopes (δ 13 C) 13,14 , a major proxy of deep-water formation, suggested that North Atlantic Deep Water (NADW), the deep component of the AMOC, weakened after iNHG.
The weakened NADW and AMOC after iNHG are an enigma because the density of surface water in the high-latitude North Atlantic may have been increased during iNHG by a significant decrease in sea surface temperature (SST) 15 and by removal of a large amount of moisture by the development of continental ice sheets 16 . In fact, general circulation models 4,5 have predicted that the AMOC intensified during iNHG, whether the atmospheric CO 2 4-6 or CAS closure hypothesis 7,8 was adopted. Because benthic foraminiferal δ 13 C values are restricted by uncertainties associated with biological modification, air-sea gas exchange, and remineralisation effects [17][18][19][20] , it is necessary to examine the activities of the deep water and AMOC based on other proxies.
In this paper, we present glacial-interglacial variability of the Nordic Seas overturning circulation, northern sector of the AMOC, during iNHG. Our multi-proxy reconstructions of surface-water and bottom-water environments indicate that the Nordic Seas overturning intensified after the Marine Isotope Stage (MIS) G4/G3 deglaciation. More importantly, our highresolution proxy data reveal that the overturning remained relatively strong even in early stages of individual glacials after iNHG. We propose that the active AMOC during early glacial stages facilitated the efficient growth of continental ice sheets by mitigating heat and utilizing moisture supplied by the North Atlantic Current (NAC), resulting in the amplification of Northern Hemisphere glaciation.

Results
Reconstruction of ice sheet-AMOC interactions. We performed multi-proxy analyses of sediments recovered from the southern Gardar Drift, subpolar North Atlantic (Integrated Ocean Drilling Project [IODP] Site U1314) ( Fig. 1) (see "Methods" section). This site is well suited to reconstruct ice sheet-ocean interaction, including the AMOC, because of its location under the strong influence of icebergs 21,22 , the Irminger Current (north-western branch of the NAC), and the Iceland-Scotland Overflow Water (ISOW). Currently, two kinds of cold, dense waters pass southward across the Greenland-Scotland Ridge as overflows: the ISOW to the east and Denmark Strait Overflow Water (DSOW) to the west (Fig. 1). The individual overflows are roughly 3 Sverdrups (1 Sverdrup= 10 6 m 3 s −1 ) 23 , and volume transport by both overflows corresponds to about one-third of total NADW production 24 . After passing the Greenland-Scotland Ridge, the overflows descend sills and further double their volume transport by entraining surrounding water 25,26 . Afterward, the overflows form the Deep Western Boundary Current around Cape Farewell and eventually become lower NADW around the Labrador Basin. In addition, the less dense Labrador Sea Water (4-6 Sv 26 ) forms upper NADW mainly above 2 km depth. Because ISOW, including the entrained downslope water, is the major source of lower NADW, its record from Site U1314 provides valuable information on the behavior of the lower limb of the AMOC during iNHG.
The analyses performed were (i) decomposition of isothermal remanent magnetization (IRM) of sediments as an ISOW activity proxy, (ii) measurement of coarse fractions (>150 μm) of ice rafted debris (IRD) and mineral amounts in sediments as a proxy of arriving and melted icebergs, (iii) alkenone biomarkers as proxies of SST and productivity, and (iv) abundances of diatoms as proxies of surface water productivity and NAC activity (see "Methods" section). IRM decomposition is an effective method to reconstruct ISOW activity; the high-coercivity component (HCC) of IRM reflects the amount of basaltic sediments from Iceland transported by the ISOW flow, whereas the low-coercivity component (LCC) is attributed to acidic sediments of southern origin 27 . An alternative explanation for the IRM-HCC, that the increased import of Icelandic basalts is a result of freshly weathered and eroded material in the aftermath of shrinking ice sheets, is unlikely because glacial HCC greatly increased with the development of ice sheets covering regolith after iNHG, 2.7 Ma ( Fig. 2; see below discussion). Our ISOW and iceberg records seem different from interpretations by previous studies (e.g., ref. 28 ) because of the much higher resolution of our records. For example, a peak in our iceberg record represented by one analysis point (a composite of a few hundred years) (Fig. 2) is undetectable in low-resolution IRD records on millennial time scales. Our high-resolution records succeed in reconstructing ice sheet-ocean interactions that were previously undetectable. increased at IODP Site U1314 in the MIS G4/G3 deglaciation (Fig. 2), reflecting the increase in ISOW formation. In addition, our high-resolution IRM record reveals two new important facts. First, after the MIS G4/G3 deglaciation, the fraction of the HCC doubled not only in interglacials but also in early stages of individual glacials (gray areas in Fig. 2c, hereinafter referred to as 'early glacial stages'), maintaining an average value comparable to interglacial values before the MIS G4/G3 deglaciation (dashed line in Fig. 2c). This finding indicates that the formation of ISOW was relatively strong even in the early glacial stages. Second, sharp declines in the HCC occurred in the late stages of individual glacials (hereinafter 'late glacial stages'). These declines corresponded to IRD and quartz peaks (Fig. 2a, b), indicating that suppression of ISOW flow accompanied the melting of icebergs. The MIS G4, G2, and 104 glacials each experienced a single prominent iceberg event (Fig. 2), in contrast to the MIS 100 glacial when icebergs repeatedly arrived on millennial time scales 22 . The scale of the individual iceberg events decreased step by step after the largest one in MIS G4. The HCC, the ISOW proxy, recovered to more than half of the pre-event level within about 2000 years after each iceberg event and continuously increased over the next several thousand years in subsequent deglaciations (Fig. 3). During MIS 104/103, however, final ISOW recovery was delayed by about 4000 years compared with the other deglaciations (Fig. 3a).
Surface water environments. Our proxy records show that surface water environments in the subpolar North Atlantic were under strong influence of glacial-interglacial variations, and that SST increased during interglacials (Fig. 2h). However, during the MIS G4/G3 deglaciation, before the MIS G3 interglacial, significant warming occurred, deviating from the glacial-interglacial variations. This pre-MIS G3 interglacial warming was also confirmed by the exceptional increase in the abundance of a   temperate and warm-water diatom species, Thalassionema bacillare ( Fig. 2l). On the other hand, concentrations of alkenone, calcite (reflecting nannofossil abundance 22 ), and diatoms decreased during glacials ( Fig. 2i-k), indicating decreases in surface water productivity.

Discussion
Our finding that the ISOW strengthened after the MIS G4/G3 deglaciation is well consistent with the radiogenic neodymium isotope variability of bottom-water masses in the northeast Atlantic, which indicated that advection of more radiogenic overflow from the Nordic Seas was enhanced after 2.7 Ma 29 . This finding does not completely contradict the interpretation of benthic foraminiferal δ 13 C records; that is, decrease in NADW formation after iNHG 13,14 . The IRM-HCC record shows that the reduced ISOW formation with large-scale iceberg melting in the late glacial stages became more severe from MIS G2 to MIS 100 (arrow in Fig. 2c). Therefore, it is suggested that low-resolution δ 13 C records reflect the reduced ISOW formation in the late glacial stages (Fig. 2d). Nevertheless, the most significant fact in our findings is that ISOW formation remained relatively strong during the early glacial stages, as opposed to the late glacial stages when it significantly decreased. It is reasonable to infer that the increased ISOW formation after the MIS G4/G3 deglaciation represents intensified overturning north of Iceland. However, the NAC is considered to not have reached the Nordic Seas during glacials after iNHG, because increases in surface alkenone productivity with SST declines in the mid-latitude North Atlantic (IODP Site U1313) were interpreted as southward movement of cold and productive Arctic waters 30 . Our research site is more suited to monitoring the NAC movement because of its location, just south of the Arctic Front and under the influence of the northern NAC ( Fig. 1). Nonetheless, surface water productivity clearly decreased during the glacials ( Fig. 2h-k). Thus, we exclude the possibility of the glacial southward movement of Arctic waters and the NAC. If the NAC did not reach the Nordic Seas during the glacials, the close relationship between the reduced ISOW formation and iceberg melting would not be detected in our results. A late Pliocene-earliest Pleistocene SST record from IODP Site U1313 suggested that the NAC did not move southward but only modestly reduced poleward heat transport, even during the glacials 31 .
Additional evidence for the intensified Nordic Seas overturning after the MIS G4/G3 deglaciation were obtained from SST and bottom water temperature. The pre-MIS G3 interglacial warming recorded in our research site was also confirmed in the midlatitude North Atlantic 30 and was most pronounced in the South Atlantic 32 . Simultaneously, diatoms in the subpolar North Atlantic achieved their highest productivity, based on the drastic increase in abundance of the temperate and warm-water species Thalassionema bacillare (Fig. 2k, l). The current sporadic occurrence of T. bacillare in the Norwegian Sea is considered to reflect transport from farther south by the NAC 33 . Therefore, the enhanced Nordic Seas overturning allowed the NAC to transport heat that was previously stored in southern seas into the subpolar North Atlantic. In addition, the increase in the ISOW production in the MIS G4/G3 deglaciation occurred synchronously with major cooling of deep water in the North Atlantic 34 (Fig. 2f), consistent with their relationship during the Holocene and lastglacial interstadials 12 .
As we have discussed, the increased ISOW formation after the MIS G4/G3 deglaciation indicates the intensification of the Nordic Seas overturning circulation (northern sector of the AMOC). During the Holocene, however, activity of the lower limb of the AMOC remained relatively stable, as upstream tributaries of NADW (ISOW, DSOW, and Labrador Sea Water) fluctuated to compensate for each other on a millennial time scale [35][36][37][38] . This is because inflow and freshening of Arctic sea ice reduced convection in the Nordic Seas, whereas in the more southerly Labrador Sea, surface cooling, more important than sea ice during the interglacial, increased convection 38 . In addition, antiphase behavior between ISOW and DSOW may have been generated as a result of atmospheric and hydrographic reorganizations in the Nordic Seas driven by the decrease in Northern Hemisphere summer insolation 37 .
However, although the upstream tributaries of NADW were dominated by decadal-scale, centennial-scale, and millennialscale fluctuations with different mechanisms 37 , antiphase variability was found only in the millennial-scale ones [35][36][37][38] . Rather, ISOW and DSOW showed in-phase variability on decadal-tocentennial scales (see Fig. 5 in ref. 37 ). In addition, there is no convincing evidence of antiphase variation in instrumental records on shorter time scales 37 . On a larger, orbital time scale, Holocene-like antiphase variation was unlikely to occur during glacials, given that the Labrador Sea became covered by a large amount of sea ice and icebergs, and the atmospheric and hydrographic responses to insolation were different from those in the Holocene.
Our IRM-HCC record shows that significant development of ISOW occurred not only during interglacials but also during glacials after the MIS G4/G3 deglaciation, which was also coincident with the strengthening of the surface NAC. Therefore, it is highly likely that the changes in the IRM-HCC reflect the orbitalscale or even larger time scale development of NADW, rather than ISOW alone. Of course, it is still possible that Holocene-like antiphase variations of the upstream tributaries of NADW occurred more or less during interglacials around iNHG. Nonetheless, it is important to note that the Nordic Seas overturning was clearly strengthened after the MIS G4/G3 deglaciation, regardless of the presence or absence of Holocene-like antiphase variations.
The intensified Nordic Seas overturning after the MIS G4/G3 deglaciation was consistent with the final phase of the gradual development of the AMOC, which was brought about by tectonic events. In the Nordic Seas, current water movement was not attained until at least the opening of the Arctic-Atlantic gateway by late Miocene-early Pliocene subsidences of the Hovgaard Ridge to the north 39 and the Greenland-Scotland Ridge to the south 13,40 . However, these subsidences were prerequisites rather than triggers because they greatly preceded the latest Pliocene intensification of the northern sector of the AMOC. Both the Bering and Canadian Arctic Archipelago straits were temporarily closed during the mid-Piacenzian (ca. 3.3-3.0 Ma) 41 . Model simulations suggested that the closure of the straits inhibited freshwater transport from the Pacific to the Arctic Ocean and from the Arctic Ocean to the Labrador Sea, thereby increasing deep-water formation between the Labrador Sea and the subpolar North Atlantic, and consequently strengthening the AMOC at latitudes 40-60°N 42,43 . At the same time, however, increased freshwater (particularly sea ice) transport from the Arctic through the Fram Strait freshened the Nordic Seas and greatly decreased deep-water formation 43 . Thus, other mechanisms were still needed for further development of the AMOC, including the Nordic Seas.
The closure of the CAS had various impacts on the AMOC by strengthening the Gulf Stream-NAC system 7 . Although the NAC delivered heat to high latitudes and increased ablation of ice sheets 4,5 , its northernmost extent may have increased salinity by supplying large amounts of moisture to ice sheets. This increased moisture supply may have facilitated freshwater delivery to the Arctic Ocean 44 , and sea ice expanded to its modern winter maximum extent during iNHG 45 . Unlike the modern Arctic Ocean, sea-ice formation produced large amounts of dense brines in the Eurasian shelf regions, but these brines did not contribute to the deep-water formation in the Nordic Seas and iNHG 46 .
Instead, sea ice expanding to the Nordic Seas counterbalanced heat transported by the NAC and promoted cooling of the highlatitude North Atlantic through ice-albedo feedback 45 .
The CAS was initially considered to have mostly closed by about 3.0 Ma 47,48 or 4.2 Ma 49 . However, much earlier closure in the middle Miocene was suggested based on a finding of detrital zircons that were interpreted to have been transported by a river from the Panama Arc to middle Miocene sediments of the lower Magdalena Basin, northern Andes 50 . This much earlier closure of the CAS was supported by some palaeoceanographic and biological studies [51][52][53] . According to the earlier closure theory of the CAS, deep to middle water connection between the Pacific and the Atlantic through the CAS ceased during the Pliocene, but shallow water exchange may have continued and played a role in regulating the AMOC strength in stages. Recently, however, the earlier closure theory 50 was refuted based on the existence of alternative source areas of the detrital zircons other than the Panama Arc and resolution of sedimentological misunderstanding 54 . Although the controversy about the timing of the final closure of the CAS has not yet been resolved 55,56 , a study based on palaeontological and molecular records suggested that this closure occurred around 2.8 Ma 54 , shortly before the intensification of the northern sector of the AMOC, which was recorded in our proxy records.
A final trigger for the AMOC intensification including the Nordic Seas sector was probably global cooling induced by a low atmospheric CO 2 concentration. The atmospheric CO 2 concentration gradually declined throughout the Cenozoic because of decreased volcanic release and enhanced chemical weathering of continental rocks in uplifted regions 57 . Furthermore, latest Pliocene changes in oceanic carbon storage reduced the atmospheric CO 2 concentration to a low threshold at which orbital forcing could induce glaciation in the Northern Hemisphere 6 . The lowered atmospheric CO 2 with the above ice-albedo feedback caused a drastic SST decline at 2.7 Ma in the high-latitude North Atlantic 15 (Fig. 2g). Then, the convection of cold, densified water was facilitated in the Nordic Seas, enhancing the AMOC. Once the AMOC had intensified, it would have served as an amplifier of glaciations. The relatively strong AMOC in early glacial stages indicates that the heat transport of the NAC was mitigated by  c Age model construction 27 . A normalized χ+ natural gamma radiation index record (light blue curve) 21,27 was tuned to the normalized LR04 δ 18 O stack record (black curve) 16 .
glacial coldness and expanding sea ice. Furthermore, the moisture supply by the NAC would likely have been actively maintained. Therefore, the vigorous AMOC, including the Nordic Seas sector, during early glacial stages facilitated efficient growth of continental ice sheets and amplified 41-kyr glacial oscillations after iNHG.

Methods
Materials. All analysis samples were collected from complete spliced sections of sediment cores recovered from IODP Site U1314. At this site, three holes (Holes A, B, and C) were drilled with an advanced piston corer; the core recovery was nearly 100%, with the quoted rate exceeding 100% at each hole influenced by core expansion 58 . The complete spliced section was obtained down to 281.00 meters composite depth (mcd), and the two deepest cores of Hole B were appended to the splice down to 299.00 mcd 58 . Note that the mcd units in this paper are shifted up by 0.90 m relative to the shipboard mcd units for depths greater than 260.00 mcd (see Appendix B of ref. 21 ). Only one lithologic unit was defined at IODP Site U1314, and it consisted of nannofossil-rich and clay-rich sediments 58 . An age model between 2.88 and 2.11 Ma for the core was obtained by tuning a hybrid environmental proxy record (a combination of magnetic susceptibility and natural gamma radiation) to the LR04 global δ 18 O stack record 16 by Hayashi et al. 21 , and was further developed by Sato et al. 27 (Fig. 4). The tuning of the age was very successful (Fig. 4c), but sedimentation rates were high, up to 31.2 cm/kyr, in the MIS G12 glacial, G8 glacial, and G5 interglacial (Fig. 4b). The high sedimentation rates during the MIS G12 glacial may reflect tuning uncertainty near the end of the age model, whereas those during the MIS G8 glacial and the MIS G5 interglacial probably represent uncertainties inherent in the tuning target (LR04 global δ 18 O stack) itself 16 . Excluding the uncertainties in these three periods, the sedimentation rates from MIS G13 to MIS 99 varied between 9.9 and 24.7 cm/kyr (average: 17.3 cm/kyr) without a consistent and systematic correlation with glacial-interglacial variability, implying that our palaeoenvironmental proxy records were not impacted by sediment dilution. The high sedimentation rates (uncertainties of the age model) in MIS G12, G8, and G5 have potential to interfere with data obtained by three quantitative proxies (IRD [ Fig. 2a], alkenone concentration [ Fig. 2i], and diatom abundance [ Fig. 2k]) in this study. However, their variation patterns are in good agreement with variation patterns of corresponding qualitative proxies (quartz weight percentage [ Fig. 2b], calcite weight percentage [ Fig. 2j], and relative abundance of diatom species [Fig. 2l]), indicating that they were virtually unaffected by the high sedimentation rates. In addition, IRM-HCC and alkenone SST are also qualitative proxies; thus, their interpretations are unaffected by the age-model uncertainties and changes in sedimentation rates.
Prior Decomposition of IRM gradient curves. The method for decomposition of IRM gradient curves was the same as that used in ref. 27 , whereas the end-member curves were chosen from the new data set, and the decomposition of IRM gradient curves was recalculated in this study. The IRM acquisition curves were measured for the samples at 2-cm intervals between 269. 13 22 ) were used for calculations, but the samples with the highest IRD fractions (267.95 and 267.93 mcd) were excluded from analysis because of their weak magnetic signals. The IRM acquisition curves were obtained by the application of stepwiseincreasing uniaxial fields to the samples with 29 steps from 1 mT to 1 T, using a MicroMag 2900 Alternating Gradient Magnetometer (Princeton Measurement Corporation). The IRM acquisition curves were normalized by the IRM intensity at 1 T, and the first derivatives of the acquisition curves were then calculated with respect to the log 10 field. The observed IRM gradient curves (f obs ) were decomposed into a linear combination of two end-member curves (f HCC and f LCC ), where f HCC and f LCC are the IRM gradient curves for the HCC and LCC, respectively. The Sratio was defined as a ratio of two IRM intensities at 1 T and −0.1 T and is shown in Fig. 5a. The samples with low S-ratio values (<0.455) and high S-ratio values (>0.835) were chosen for the HCC and LCC, respectively, and the end-member curves were determined by averaging the IRM gradient curves of selected samples (Fig. 6). The optimal mixing fractions (C HCC and C LCC ) were chosen, and the sum of squared residuals (S) was minimized in the ranges C HCC ≥ 0 and C LCC ≥ 0: where B i and f res are a magnetic field and residual curve, respectively. Note that analytical uncertainties for the mixing fractions cannot be calculated in this method. Alternatively, the integral of the residual (R) was calculated as an indicator of fitting error: The average of the R through the analytical period was 5.0%. The R value intermittently showed large values (Fig. 5d), although these data points were sharp pulses and had little influence on our interpretations. The IRM gradient curves represent the distribution of coercivity of magnetic minerals included in the experimental samples. Previous rock-magnetic studies based on the IRM acquisition curves, S-ratios, and low-temperature and high-temperature magnetometry (ref. 22 and references therein) revealed that the dominant magnetic mineral of the Site U1314 sediments was magnetite. The coercivity of a magnetic mineral varies with the mineral's properties, such as the chemical composition, grain size, grain shape, and stress. The various states of magnetite grains contained in the lithic particles form the coercivity spectra of basaltic and acidic sediments (Fig. 6), and the mixing of basaltic and acidic sediments transported from different provenances result in systematic change in coercivity spectra among these endmembers (Fig. 5). Note that it is difficult to apply this method to estimate contributions from sediments containing trace amounts of magnetic minerals. At Site U1314, noise of IRD was negligible in IRM because it contributed little to magnetic properties 22 . -99]) using the same core, we re-analyzed this record from the core because it was likely to include noise caused by gypsum (Figs. 7d, 8, 9a) that formed diagenetically after the deposition of IRD. The samples for IRD counts were freeze-dried, weighed, and then sonicated in ethanol for 2 min. We counted all non-biological coarse grains per 0.5 g collected on a 150-µm sieve with the use of a stereoscopic microscope (SMZ18, Nikon). Non-biological coarse grains included in the analysis samples were IRD, volcanic glass, pyrite, and gypsum ( Figs. 7 and 8). Large peaks of volcanic glass during the MIS 100 glacial (Fig. 7b) represent ash falls, but their synchronous occurrence with IRD suggests that some quantity of volcanic glass was also transported by icebergs. Diagenetic pyrite and gypsum occurred throughout the samples (Fig. 7c, d). To reconstruct a pure iceberg record, IRD counts were made by carefully excluding volcanic glass, pyrite, and gypsum grains.
Mineral analysis. To reveal changes in the amounts of minerals in the sediments, we performed X-ray diffraction (XRD) analysis at (1) 22 . However, the profile fitting for that data set possibly failed to be separated from gypsum (q.v., Fig. 10), and thus, we re-analyzed the data. Each bulk sample was ground and homogenized with an agate mortar. Simultaneously, powdered zincite (ZnO) was added as an internal standard in a fixed percentage (5 wt%) to the samples. All XRD measurements were performed on an X-ray diffractometer (RINT 2100 V, Rigaku), using CuKα radiation monochromatized by a curved graphite crystal in a step of 0.02°with a step-counting time of 2 s. The profile fitting of the XRD patterns was performed with the scientific graphical analysis program XRD MacDiff 59 . Relative amounts of gypsum and pyrite in the sediments were estimated based on the integrated peak-intensity ratios of the individual minerals to the internal standard zincite (5 wt%) (Fig. 9). We also determined relative amounts (wt%) of quartz and calcite using calibration curves measuring artificial standards comprising 5 wt% zincite plus a matrix consisting of specified weight percentages of quartz and calcite (Fig. 2b, j; Fig. 9c, d). Errors for all XRD data (miscounting errors) were within 1%. Gas chromatography was conducted using a gas chromatograph (GC-17A, Shimadzu) with on-column injection and a flame ionization detector. Samples were dissolved in hexane. Hydrogen was used as a carrier gas, and the flow velocity was maintained at 30 cm/s. The column used was a CP-Sil 5 CB ( 37 value corresponds to that of estimated SST better than ±0.3°C when the SST is calculated according to the relationship between U K0 37 and SST described below. In addition, irreversible adsorption of C37:3 on the chromatographic column is often a major source of error when the total amount of C37 alkenone injected into the system for analysis is less than 5 ng 61 , which leads to overestimation of SST. Therefore, the error should be minimized in this study because more than 5 ng was used in each injection except for a few samples, and excess negative error based on the amount of alkenones injected into the system was estimated using the reported relationship 61 between the injected amount of alkenones and the difference between analytical and true U K0 37 values. Alkenone concentration (the sum of two  C37 alkenones) is shown as values normalized by the total organic carbon (TOC) content in the sediment samples. TOC content was measured using the dry combustion method with an elemental analyzer (EuroEA3000, EuroVector). Prior to TOC analysis, about 5 mg of each dried sediment sample was weighed into a silver capsule, and a few drops of 6 N HCl were added to remove inorganic carbon; then, the sample was dried under vacuum. The dried silver capsule was wrapped with tin film and then installed into the elemental analyzer.
The SSTs were calculated using the equation of ref. 62 because their calibration was determined by laboratory culture experiments on Emiliania huxleyi collected in the eastern North Pacific: U K0 37 = 0.034 T (°C) + 0.039. However, E. huxleyi, which is a member of the class Prymnesiophyceae, is a main producer of alkenones in the modern ocean (e.g., ref. 63 ). The oldest record of carbonate microfossil coccoliths derived from E. huxleyi has been reported from mid-latitude North Atlantic sediments of the middle Pleistocene 64 , but alkenones have been found in much older sediments. This suggests that other extant and extinct members of the class Prymnesiophyceae presumably produced the alkenones. Although the alkenone producers have changed over time, previous work has shown that U K0 37 and global core-top calibration are applicable beyond the first occurrence of E. huxleyi 65,66 . In recent years, the alkenone thermometer has therefore been used to produce several long-term SST records (e.g., Refs. 15,67 ), suggesting that it provides reliable temperature estimates for at least the last 5 Ma 30 .
Diatom valve counts. We chose 66 samples from 299. .99 mcd to reconstruct the diatom productivity record corresponding to glacial-interglacial variations for 2.88-2.50 Ma (MIS G13-99). The samples were freeze-dried, weighed, and then sonicated in deionized water for approximately 0.2 s following the procedure of ref. 68 . We counted more than 300 valves of all diatom species collected on membrane filters (pore size: 0.45 μm) with the use of a scanning electron microscope (JSM-IT100, JEOL). Total valves per 10 µg sample were counted for 11  . In the counting of species in the genera Thalassiothrix and Thalassionema (including Thalassionema bacillare), a head or foot pole was identified as a half valve because very long valves of those species were usually broken and fragmented into small pieces in taphonomic and experimental processes, except for the two robust poles 68 . Abundances of all diatoms and Thalassionema bacillare (a temperate and warm-water species 33 ) per microgram were determined based on the count data. The relative valve abundance of Thalassionema bacillare was calculated based on the count data, except for the samples with very few valves.

Data availability
The data sets analyzed during this study are available from the PANGAEA data repository (https://doi.org/10.1594/PANGAEA.922223).