Late Miocene climate cooling and intensification of southeast Asian winter monsoon

The late Miocene offers the opportunity to assess the sensitivity of the Earth’s climate to orbital forcing and to changing boundary conditions, such as ice volume and greenhouse gas concentrations, on a warmer-than-modern Earth. Here we investigate the relationships between low- and high-latitude climate variability in an extended succession from the subtropical northwestern Pacific Ocean. Our high-resolution benthic isotope record in combination with paired mixed layer isotope and Mg/Ca-derived temperature data reveal that a long-term cooling trend was synchronous with intensification of the Asian winter monsoon and strengthening of the biological pump from ~7 Ma until ~5.5 Ma. The climate shift occurred at the end of a global δ13C decrease, suggesting that changes in the carbon cycle involving the terrestrial and deep ocean carbon reservoirs were instrumental in driving late Miocene climate cooling. The inception of cooler climate conditions culminated with ephemeral Northern Hemisphere glaciations between 6.0 and 5.5 Ma.

T he late Miocene (~11.6 to 5.3 Ma) stands out as a period of exceptional interest within the long-term Cenozoic cooling trend toward icehouse conditions, as it represents a geologically recent interval of relative global warmth that was marked by profound environmental change in both terrestrial and marine ecosystems (e.g., refs. 1,2 ). This interval provides a unique opportunity to document climate-carbon cycle dynamics on a warmer-than-modern Earth and, thus, to help guide models and constrain predictions of climate change and sensitivity. The detailed sequence of climate events and the range of natural climate variability through the late Miocene remain, however, poorly understood, mainly due to the scarcity of continuous, high-resolution climate archives. Most available records are adequate for characterizing long-term trends or mean states, but do not capture short-term climate events and orbital-scale phase relationships required to assess, for example, changes in ice volume, monsoon intensity, and carbon fluxes.
Multiproxy temperature reconstructions indicated that a reduced sea surface temperature (SST) zonal gradient generally prevailed in the tropical Pacific Ocean during the late Miocene, in contrast to the sharper gradient that developed during the late Pliocene to Pleistocene 3,4 . This mean state, typically referred to as "permanent El Niño-like conditions" or "El Padre 5,6 , exerts a fundamental impact on regional and global climate because it is dynamically linked to the weakening of the Hadley and Walker circulation and the state of upper-ocean stratification 4,7,8 . However, it is difficult to reconcile the late Miocene warmth with inferred low atmospheric pCO 2 levels, close to preindustrial values (e.g., ref. 9 ). This apparent decoupling between climate warmth and atmospheric pCO 2 variations has prompted intense debate about the dynamics of warm climates and the role of pCO 2 as driver of climate variations under different background states (e.g., refs. 10,11 ). The widely held view of sustained equable warmth through the late Miocene was recently challenged by SST reconstructions, which revealed that a prolonged global cooling spell occurred between~7 and~5.5 Ma 2,12 . During this period, SST dipped below early Pliocene values in the Mediterranean, the Pacific, Atlantic and Indian Oceans and the latitudinal thermal gradient intensified 2,12 . Land records also indicate substantial aridification and vegetation changes between~8 and 5.5 Ma (see compilation in ref. 13 and references therein) with a reversal of this long-term trend 5.3 Myr ago, at the beginning of the Pliocene period 14 . However, the drivers of these major changes are still enigmatic. It remains unclear, in particular, whether climate cooling occurred as a response to changing climate boundary conditions, such as ice volume, pCO 2 , tectonic setting and ocean-atmosphere circulation and to what extent these changes were coupled to Northern or Southern Hemisphere climate dynamics.
In this study, we reconstruct the detailed evolution of deep and near surface water masses (using stable isotopes and mixed layer temperatures) at Ocean Drilling Program (ODP) Site 1146 in the South China Sea (Fig. 1) and we investigate relationships to highlatitude climate variability, global ocean circulation changes, and radiative forcing over the interval 9-5 Ma. This work extends published high-resolution benthic isotope and upper-ocean temperature time series from the same site [15][16][17][18] . During the late Miocene, the location and water depth of Site 1146 were approximately similar to that of today (19°27.40′ N, 116°16.37′ E, water depth: 2092 m) and the connection between the South China Sea and the western Pacific Ocean remained fully open 19,20 . The benthic isotope signal at Site 1146 is, therefore, representative of Pacific intermediate/deep water masses originating at higher latitudes. Site 1146 is located at the northwestern edge of the western Pacific warm pool (WPWP) and is also ideal to constrain meridional variations in the extent of the WPWP and to monitor changes in southeast Asian monsoon climate. The extended, carbonate and clay-rich succession recovered at this site 20 , thus, provides an outstanding archive of subtropical climate variations, allowing new insights into the dynamics and forcing processes of late Miocene climate evolution.

Results
Late Miocene astronomically tuned chronology. The 1146 benthic foraminiferal stable isotope records based on Cibicidoides wuellerstorfi and/or Cibicidoides mundulus (5 cm sample spacing along a composite sequence or splice from Holes 1146A and 1146C) were tuned to an eccentricity-tilt (ET) composite target generated from the La2004 orbital solution 21 Fig. 1 Location of ODP Site 1146 within a slope basin at the northern margin of the South China Sea. Satellite image and bathymetry from ref. 75 . Paleogeographic reconstruction at 10 Ma (simplified from ref. 19 ) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03950-1 mean sedimentation rate of~3 cm kyr −1 with a maximum of 5 cm kyr −1 and a minimum of 1 cm kyr −1 and a mean temporal resolution of~2 kyr over the interval 9-5 Ma ( Supplementary  Fig. 4B). We note that the short eccentricity (100 kyr) period is prominent in the untuned and tuned benthic δ 18 O records from 9.0 to 7.9 Ma and that the low amplitude of short eccentricity and high amplitude of obliquity between~7.7 and 7.2 Ma are clearly reflected in the benthic δ 18 O series, which exhibits pronounced 41 kyr variability over this interval ( Fig. 2c; Supplementary  Figs. 1-4). The interval 6-5 Ma includes prominent benthic δ 18 O maxima, identified as T8, TG4, TG12, TG14, TG20, and TG22. These globally traceable δ 18 O enrichments [22][23][24][25][26] provide additional stratigraphic control. Superimposed on higher frequency variations (mainly 41 kyr), the untuned benthic and planktic δ 13 C series display low-frequency oscillations that broadly relate to thẽ 400 kyr long eccentricity cycle ( Fig. 2a; Supplementary Fig. 1). Comparison of benthic and planktic δ 18 O and δ 13 C records plotted in the depth and time domains shows that the original spectral characteristics are preserved following the tuning procedure.
Temporal trends in benthic and planktic stable isotopes.  Fig. 2C). Between 7.3 and 6.9 Ma, the Site 1146 high-resolution benthic and planktic δ 18 O records reveal a series of previously unrecognized short-term climate events (Fig. 2b, c). The benthic δ 18 O curve resolves a~80 kyr long positive excursion (~0.3‰ amplitude) centered at 7.2 Ma followed by a rebound before a stepwise increase at 7.1-7.0 Ma (0.2‰ mean increase) ( Fig. 2c; Supplementary Fig. 2B). The benthic δ 18 O shift was coupled to a stepwise increase in planktic δ 18 O between 7.2 and 7.1 Ma (Fig. 2b, c; Supplementary Fig. 2A, B), which marked the onset of a long-term trend of substantially heavier values (0.3‰ mean  Benthic and planktic δ 13 C exhibit consistent long-term (400 kyr) and short-term (41 kyr) variability from~9 to~5 Ma Coherence (k) between benthic and planktic δ 18 O remains overall lower than between benthic and planktic δ 13 C over the  TG12   TG14  TG20  TG22   T8 TG4   T8  TG4   TG12  TG14   TG20   TG22   T8  TG4   TG12  TG14  TG20  TG22   23 (Fig. 3b). This transient warming is interrupted by a pronounced cooling step of~2°C between 7.1 and 6.9 Ma, previously documented at this site by a low-resolution study 34 . This cooling step marks the onset of a long-term trend of cooler temperatures, lasting until~5.5 Ma, which coincides with a longterm increase in the mean and amplitude variability of seawater

ARTICLE
In this study, we discuss relative changes in mixed layer temperatures rather than absolute values, since the history of Miocene seawater Mg/Ca composition is still poorly constrained (Supplementary Note 3). Furthermore, our interpretations are based on relatively short-term changes in mixed layer temperatures, which are not affected by the long-term variability of seawater Mg/Ca concentration with changes in the order of millions of years due to the long residence time of these elements in the ocean 29 .

Discussion
The planktic and benthic δ 18 O signals at Site 1146 differ markedly in their long-and short-term trends between 9 and 5 Ma, pointing to a decoupling of regional hydrology and the evolution of the Antarctic ice sheet, which formed the main component of the cryosphere during the middle and late Miocene (e.g., refs. 35,36 Supplementary Fig. 9A, C), as indicated by a previous low-resolution study 34 . These trends signal a change in the amount and/or δ 18 O composition of precipitation and runoff, likely associated with changes in the provenance and/or seasonality of precipitation toward a more pronounced monsoonal seasonality and a more temperaturecontrolled seasonality of rainwater δ 18 O (i.e., δ 18 O depleted winter precipitation 37 ). We attribute these hydrological changes in the northern South China Sea after~7 Ma to cooling and drying of the Asian landmass and a related southward shift of the average summer position of the Intertropical Convergence Zone (ITCZ), resulting in decreased influence of tropical convection and intensified dry winter monsoon over southeast Asia. Drying and cooling on the Asian continent at~7 Ma are supported by independent lines of evidence including enhanced dust accumulation rates in northern China 38 , vegetation change in central China 39 and an increase in the mean grain size of the terrigenous  40 . In addition, the predominance of a mollusc group preferring cold-arid conditions in the loess and paleosol layers of central China between 7.1 and 5.5 Ma is indicative of a dominant winter monsoon regime over this period 41 . In contrast to these major hydrological changes in the Northern Hemisphere, mean benthic δ 18 O suggests only a relatively modest, stepwise glacial expansion of the Antarctic ice sheet and/or deep water cooling at~7 Ma (Figs. 2c and 4f). However, the intensification of the southeast Asian winter monsoon after~7 Ma was associated with a long-term trend toward heavier benthic δ 18 O maxima, which culminated in the most intense maxima (TG22, 20, 14, and 12 between 5.8 and 5.5 Ma) within the entire late Miocene, before reversing in the   (Fig. 5e, f). Mixed layer temperatures additionally show concurrent sharp decreases of 2-3°C during these events (Fig. 5a), implying massive Northern Hemisphere cooling down to subtropical latitudes. The occurrence of ice rafted debris in North Pacific 43 and North Atlantic 44 sediment cores further indicates Northern Hemisphere ice buildups between 6 and 5 Ma. Expansion of Arctic sea ice during these intense cold spells would have increased the positive albedo feedback, amplifying cooling and favoring ice growth. Together these lines of evidence support the development of ephemeral Northern Hemisphere ice sheets (e.g., Greenland, Alaska, Labrador) between 6.0 and 5.5 Ma that were highly susceptible to insolation forcing. The Site 1146 records additionally reveal that climate cooling and intensification of the winter monsoon at~7 Ma coincided with the final stage of a long-term, global benthic, and planktic δ 13 C decline 27,28 (LMCIS, Figs. 2a and 4c, d). This major shift of 1‰, which started close to 7.8 Ma, has been interpreted as a global decrease in the δ 13 C of the dissolved inorganic carbon pool, although its causes remain debated (e.g., refs. 22,45,46 ). A long-held view among contending hypotheses is that this global δ 13 C decrease was linked to the late Miocene spreading of C4 grasslands, which are better adapted to low pCO 2 and to reduced seasonal precipitation. This large-scale expansion is thought to have resulted in a transfer of 13 C from the marine to the terrestrial carbon pool [47][48][49] . A decrease in atmospheric pCO 2 , linked, for instance, to long-term changes in the oceanic and/or terrestrial carbon inventories, could explain climate cooling after 7 Ma associated with equatorward migration of the ITCZ and contraction of the WPWP.
The gradient between benthic and planktic δ 13 C additionally provides insights into changes in atmospheric pCO 2 , since it is influenced by two main factors: the sequestration efficiency of the biological pump and equilibration processes between the upper ocean and the atmosphere (Supplementary Note 4; Supplementary Fig. 10). The equilibration time for δ 13 C in the mixed surface layer of the ocean exhibits a linear correlation to the ratio of dissolved inorganic carbon to pCO 2 , which leads to slow equilibration and elevated δ 13 C in the mixed layer of the ocean with respect to the atmosphere under low pCO 2 . Recent model simulations showed that accelerated equilibration under elevated atmospheric pCO 2 decreases the isotopic disequilibrium, leads to lower upper ocean δ 13 C and, thus, decreases the gradient between the δ 13 C of surface and deep water masses 50 . Consequently, the vertical δ 13 C gradient in the ocean exhibits a gentler slope under high atmospheric pCO 2 and steepens during intervals of declining pCO 2 .
Steepening of the gradient between planktic and benthic δ 13 C after~7 Ma at Site 1146, when mixed layer temperatures also declined (Fig. 3a, b), suggests that pCO 2 levels decreased, eventually reaching levels that enabled the formation of transient Northern Hemisphere ice sheets between 6.0 and 5.5 Ma. This steeper gradient also denotes a prolonged interval of substantially enhanced marine productivity and accumulation rates of biogenic components ("biogenic bloom" originally described in ref. 51 ) at numerous locations in the Pacific, Indian and Atlantic Oceans (e.g., ref. 52 and references therein). In the eastern equatorial Pacific Ocean, opal, and carbonate deposition reached a maximum between 7.0 and 6.4 Ma during the peak of the biogenic bloom in the region 46 . Thus, a plausible scenario is that changes in nutrient supply and/or pathways stimulated marine productivity after~7 Ma. Steepening of the equator to pole temperature gradient associated with global cooling after~7 Ma (Fig. 3a, b; ref. 2

) promoted intensification of the Hadley and
Walker circulation with repercussions on the wind-driven circulation and precipitation patterns (e.g., ref. 53 ). The strengthening of winds may have in turn fostered upwelling and ocean fertilization, helping to drive intense biogenic blooms through the Pacific Ocean, which enhanced carbon storage and decreased pCO 2 in the ocean in a positive feedback loop.
Previous work showed that the amplitude of the LMCIS differs in ocean basins (e.g., ref. 54 ). In particular, a comparison of benthic δ 13 C profiles indicates that the gradient between the Pacific and Atlantic Oceans intensified during the final stage of the LMCIS (Fig. 6b; ref. 54 ). The steeper inter-basinal gradient after~7 Ma cannot be explained by increased production and southward advection of North Atlantic Deep Water, as this relatively warm and/or fresh (lighter δ 18 O) and 13 C-enriched water mass appears not to have spread into the South Atlantic and Southern Ocean, which remained influenced by colder, denser (heavier δ 18 O) and δ 13 C depleted water masses through the late Miocene (Fig. 6a, b; ref. 54 ). Alternatively, the steeper inter-basinal δ 13 C gradient after~7 Ma may be driven by increased export of nutrient enriched waters with a lower preformed δ 13 C from the Southern Ocean into the Pacific Ocean (e.g., ref. 54 ) and/or to enhanced primary productivity and nutrient regeneration in the low-latitude Pacific Ocean.
Comparison of benthic δ 13 C profiles from Site U1338 in the abyssal equatorial Pacific Ocean 55 and the shallower Site 1146 in the northwestern subtropical Pacific Ocean (Fig. 6b) shows that the composition of Pacific water masses changed after 7.2 Ma. The convergence of the δ 13 C records after 7.2 Ma indicates expansion of a δ 13 C depleted central Pacific deep water mass into shallower depths during the peak of the biogenic bloom. If primarily driven by increased productivity and nutrient regeneration in the Pacific and Indian Oceans, the expansion of a 12 C-enriched deep water mass after 7.2 Ma also implies increased carbon storage in the deep Pacific Ocean. The global efficiency of the biological pump reflects a balance between high-and low-latitude regions with different  [15][16][17] and this work) with orbital parameters (eccentricity and obliquity from ref. 21 ) reveals similar sequence of climatic events during three Miocene cooling episodes with strikingly similar background orbital configuration. Blue shading marks cooling episodes following an extended period of high-amplitude variability in obliquity congruent with low variability in short eccentricity (gray shading). Light orange shading marks transient warming episodes coincident with high-amplitude variability in short eccentricity sequestration efficiency 56 . Thus, enhanced productivity and organic matter export in the tropical and subtropical ocean may increase global sequestration efficiency and lower atmospheric pCO 2 , even when deep water formation occurs in high-latitude areas with an inefficient biological pump.
The integrated benthic isotope data in Site 1146 provide the first continuous, highly resolved time series from a single site that span the last 16.4 Myr (Fig. 7a). These extended records track the transition from a warmer middle Miocene climate phase with a reduced and highly dynamic Antarctic ice cover (until~14 Ma) to an increasingly colder mode with more permanent and stable ice sheets in the late Miocene 17 . These records additionally allow us to evaluate the long-term relationship between radiative forcing and the response of the ocean/ climate that is imprinted on the benthic δ 18 O signal. For instance, the 41 kyr obliquity cycle is especially prominent in the benthic δ 18 O series between 7.7 and 7.2 Ma (Fig. 2c), during a configuration of the Earth's orbit, when high-amplitude variability in obliquity is congruent with extremely lowamplitude variability in short eccentricity ( Supplementary  Figs. 3, 4A, E). The onset of the~80 kyr long positive excursion in benthic δ 18 O centered at 7.2 Ma notably coincides with minima in obliquity (41 kyr) and in eccentricity (100 kyr, 400 kyr, and 2.4 Myr amplitude modulation) ( Fig. 7b; Supplementary Fig. 3). At obliquity and eccentricity minima, lower summer insolation at high latitudes inhibits melting of snow and ice. This conjunction of climatic forcing factors likely fostered a sustained cold phase in the high latitudes that lasted through two consecutive obliquity cycles, resulting in an extended benthic δ 18 O positive excursion. Renewed high-amplitude variations in eccentricity and precession together with maximum amplitude variability in obliquity probably drove the successive rebounds between 7.2 and 7.0 Ma. This sequence of climatic events, as well as their background orbital configuration were strikingly similar during two previous Miocene cooling episodes: the middle Miocene climate transition at~13.9 Ma, which resulted in major expansion of the East Antarctic ice sheet and the less pronounced late Miocene cooling step at~9.0 Ma (Fig. 7b). In all three instances, the 41 kyr cycle initially stands out in the benthic δ 18 O signal during a protracted period of high-amplitude variability in obliquity, congruent with low variability in short eccentricity. A marked enrichment in benthic δ 18 O (0.2-0.3‰), indicative of ice growth and/or deep water cooling toward the end of this interval, coincides with prolonged minima in eccentricity, lasting~100-200 kyr. Subsequent rebounds at peak insolation, linked to changes in eccentricity cadence (from 400 to 100 kyr variability), indicate episodes of transient ice sheet disintegration and deep water warming. This unusual orbital congruence appears propitious to high-latitude cooling in the Northern and Southern Hemispheres, although boundary conditions differed markedly during these three intervals of climate change. The middle Miocene cooling step occurred in a much warmer climate phase, characterized by substantially lighter mean benthic δ 18 O (Fig. 7a, b). At this time, the less extensive ice cover over Antarctica was likely more dynamic and highly responsive to Southern Hemisphere summer insolation 57,58 , in contrast to the more expanded Antarctic ice sheet during the late Miocene. This long-term perspective illustrates the nonlinear response of the ocean/climate system to orbital forcing and the role of internal feedback processes including ice sheet hysteresis, latitudinal temperature gradients, ocean circulation and CO 2 exchange between terrestrial, atmospheric and oceanic reservoirs.
Arguably, the uncertainty of the CO 2 forcing during the Miocene remains a major challenge for defining the characteristics and dynamics of warmer climate states. Although, current pCO 2 reconstructions show no significant change through the late Miocene, with levels staying close to or slightly above preindustrial levels, uncertainties in excess of 200 p.p.m. (see compilations in refs. 9,59,60 ) preclude assessment of variability and sensitivity to CO 2 forcing within the critical preindustrial to modern range. To test the sensitivity of outputs to pCO 2 uncertainties, some simulations of late Miocene climate using coupled atmosphere-ocean circulation models have applied atmospheric pCO 2 concentrations in the preindustrial range (~280 p.p.m.), as well as more elevated levels of 400-450 p.p.m. (e.g., refs. 60,61 ). These studies indicated major changes in vegetation distribution 60 and sea ice cover 61 in the Northern Hemisphere under these different pCO 2 states. In particular, forest areas decreased and the albedo of the Eurasian and North American landmasses increased under lower pCO 2 , due to the markedly lower (by 4-10°C) mean air temperatures and reduced precipitation during boreal winter 60 . These findings are in agreement with a previous modeling study 62 , which found that vegetation changes were more important than paleogeography in determining late Miocene climate. Simulated mean summer SST and sea ice concentrations in the Arctic Ocean 61 also showed that the region is highly sensitive to relatively small changes in pCO 2 , as a year-round sea ice cover prevails in the central Arctic Ocean at preindustrial levels, whereas summer conditions are ice-free at concentrations of 450 p.p.m. This difference in seasonal ice cover is critical because it implies very different feedbacks in terms of albedo and heat exchange with far-reaching repercussions for global climate 61 . Recent model simulations of Pliocene warmth additionally highlighted the importance of feedbacks associated with cloud albedo and ocean mixing in driving changes in meridional and zonal temperature gradients, despite relatively modest changes in pCO 2 63-66 . Data from this study support that subtropical climate cooling and intensification of the southeast Asian winter monsoon after~7 Ma were synchronous with decreasing pCO 2 (Figs. 3a and 4b) within a global context of steepening meridional thermal gradients 2 . We speculate that this late Miocene climate shift was associated with a relatively small decline in pCO 2 , which was amplified by a conjunction of positive feedbacks. Variations in Northern Hemisphere sea ice cover and vegetation in concert with changes in ocean-atmosphere circulation were likely instrumental in driving late Miocene climate, as illustrated by recent modeling simulations of late Miocene climate [60][61][62] . The dynamic behavior of the ocean-climate system between 9 and 5 Ma suggests a tight coupling between carbon cycle variations and low-latitude climate evolution. In particular, our results show that changes in Antarctic ice volume were not the primary driver of late Miocene climate development and that lowlatitude processes, including monsoonal wind forcing of upperocean circulation and productivity had a strong influence on climate-carbon cycle dynamics. Inception of colder climate conditions at~7 Ma during the final stage of the LMCIS coincided with intensification of the Asian winter monsoon and strengthening of the Pacific Ocean's biological pump, which persisted until~5.5 Ma. This suggests that changes in the global carbon cycle involved transfer of terrestrial carbon in a cooling, drying climate, as well as fluctuations in the carbon storage capacity of the deep ocean and the sedimentary carbon sink. Ephemeral Northern Hemisphere glaciations between 6.0 and 5.5 Ma additionally indicate that atmospheric pCO 2 levels hovered close to and occasionally reached the threshold necessary for Northern Hemisphere ice sheet growth during this period.

Methods
Revision of shipboard sediment splice. Cores were sampled in~5 cm intervals (~2 kyr time resolution) from a composite sequence (shipboard splice) of Holes 1146A and 1146C (Cores 1146C-30X to 1146C-39X). After comparison of the shipboard natural gamma ray, color reflectance, magnetic susceptibility data, and overlapping benthic isotope records over the splice tie points, we made the following modification to the original shipboard splice: we defined a new tie point between Cores 1146C-38X and 1146A-39X (1146C-38X-4, 122 cm at 359.52 m below sea floor (mbsf) tie to 1146A-39X-2, 107 cm at 359.87 mbsf), based on the match of isotope data from Holes 1146A and C. This adjustment resulted in the addition of an 80 cm splice segment from Hole 1146A to the meter composite depth (mcd) scale.
Benthic and planktic foraminiferal stable isotopes. All samples were oven dried at 40°C and weighed before washing over a 63 µm sieve. Residues were oven dried at 40°C on a sheet of filter paper, then weighed and sieved into different size fractions. We measured δ 18 O and δ 13 C in the epifaunal benthic foraminifers C. wuellerstorfi and/or C. mundulus and on the mixed layer foraminifer G. sacculifer. Well-preserved tests were broken into large fragments, cleaned in alcohol in an ultrasonic bath, then dried at 40°C. Stable isotopes were measured with a Finnigan MAT 253 mass spectrometer at the Leibniz Laboratory, University of Kiel. The instrument is coupled on-line to a Carbo-Kiel Device (Type IV) for automated CO 2 preparation from carbonate samples for isotopic analysis. Samples were reacted by individual acid addition (99% H 3 PO 4 at 75°C). On the basis of the performance of international and lab-internal standard carbonates, the precision is better than ±0.09‰. Paired measurements in middle Miocene samples from ODP Sites 1146 and 1237 previously indicated no significant offset in δ 18 O and δ 13 C between C. wuellerstorfi and C. mundulus 16 . Results were calibrated using the National Institute of Standard and Technology (Gaithersburg, MD) carbonate isotope standard and NBS (National Bureau of Standard) 19 and NBS 20, and are reported on the Vienna PeeDee Belemnite scale.
Astronomically tuned chronology. The chronology is based on minimal tuning 67 of the benthic oxygen isotope record to a computed orbital solution 21 (Supplementary Note 1). We used an ET composite with equal weight of eccentricity and obliquity as tuning target and tuned δ 18 O minima to ET maxima ( Supplementary Fig. 4). This strategy is based on the notion that relatively warm summers at high obliquity promote ice sheet melting at high latitudes and that a low summer insolation gradient between low and high latitudes at high obliquity leads to a decrease in poleward moisture transport, inhibiting ice sheet buildup 68 . We did not adjust our tuning for possible phase lags between δ 18 O and insolation forcing, since the response time of smaller Miocene ice sheets is unknown. We did not include high-latitude summer insolation and/or precession (P) parameters in our tuning target to remain independent of possible changes between dominant northern (ET−P) or southern (ET+P) hemisphere precessional insolation forcing. We note that the obliquity cycle (41 kyr) and the eccentricity cycles (100 and 400 kyr) are also prominently imprinted in the benthic δ 13 C record ( Supplementary Fig. 6A). The detection of these astronomical frequencies in the δ 13 C record supports the age model, based on independent tuning of the δ 18 O series. Age correlation points are given in Supplementary Fig. 4D and Supplementary Table 1.
Mixed layer temperature reconstructions. Ma/Ca measurements were performed on~30 well-preserved specimens of G. sacculifer from the 250-315 μm fraction (average weights of~400 μg). Sample spacing is 20 cm from Sample 1146C-30X-1, 33 cm (298.43 revised mcd) to Sample 1146A-39X-4, 88 cm (387.28 revised mcd), corresponding to a mean temporal resolution of~7 kyr between 8.11 and 5.07 Ma. Sample spacing was decreased to 5 cm, corresponding to a temporal resolution of~2 kyr, over the prominent benthic δ 18 O maxima TG4, TG12, TG14, and TG22 between 6.0 and 5 Ma. Tests of G. sacculifer were crushed between glass plates and cleaned with methanol and reductive (hydrazine) and oxidative (hydrogen peroxide) steps, as detailed in refs. 69,70 . Samples were then leached and diluted with nitric acid and analyzed with a SPECTRO Ciros CCD SOP ICP-OES at the ICP-MS Laboratory of the Institute of Geosciences, University of Kiel. We monitored Fe/Ca, Al/Ca, and Mn/Ca ratios to assess the efficacy of the cleaning procedure. These trace-element/Ca ratios do not show significant correlation with Mg/Ca. Mixed layer temperatures were estimated from Mg/Ca ratios using the exponential equation for G. sacculifer with an assumed A constant of 0.09 71 . Fifty-four duplicate measurements show a mean Mg/Ca standard deviation of 0.13 mmol mol −1 . The error associated with mixed layer temperature estimates was defined by propagating the errors introduced by the Mg/Ca measurements and the Mg/Ca temperature calibration assuming no covariance among these errors, following ref. 72 .
Time frequency analysis. Orbital tuning, bandpass-filtering, calculation of sedimentation rates, and Blackman-Tukey cross-spectral analyses are performed with AnalySeries 2.0.4.2 73 . Cross-wavelet analysis is performed following ref. 74 . We utilize wavelet coherence to quantify coherence and phase relationships between the benthic δ 13 C and δ 18 O time series in frequency space. We employed the Morlet wavelet and a Monte Carlo count of 300 to establish the statistical significance level (5%). Prior to analysis, both time series were interpolated to 2 kyr time steps and standardized (zero mean and unit standard deviation). Software is available at http://www.pol.ac.uk/home/research/waveletcoherence/.
Data availability. Data are archived at the Data Publisher for Earth and Environmental Science (https://doi.pangaea.de/10.1594/PANGAEA.887393).