Coupled atmosphere-ice-ocean dynamics during Heinrich Stadial 2

Our understanding of climate dynamics during millennial-scale events is incomplete, partially due to the lack of their precise phase analyses under various boundary conditions. Here we present nine speleothem oxygen-isotope records from mid-to-low-latitude monsoon regimes with sub-centennial age precision and multi-annual resolution, spanning the Heinrich Stadial 2 (HS2) — a millennial-scale event that occurred at the Last Glacial Maximum. Our data suggests that the Greenland and Antarctic ice-core chronologies require +320- and +400-year adjustments, respectively, supported by extant volcanic evidence and radiocarbon ages. Our chronological framework shows a synchronous HS2 onset globally. Our records precisely characterize a centennial-scale abrupt “tropical atmospheric seesaw” superimposed on the conventional “bipolar seesaw” at the beginning of HS2, implying a unique response/feedback from low-latitude hydroclimate. Together with our observation of an early South American monsoon shift at the HS2 termination, we suggest a more active role of low-latitude hydroclimate dynamics underlying millennial events than previously thought.

Our understanding of climate dynamics during millennial-scale events is incomplete, partially due to the lack of their precise phase analyses under various boundary conditions.Here we present nine speleothem oxygenisotope records from mid-to-low-latitude monsoon regimes with subcentennial age precision and multi-annual resolution, spanning the Heinrich Stadial 2 (HS2)a millennial-scale event that occurred at the Last Glacial Maximum.Our data suggests that the Greenland and Antarctic ice-core chronologies require +320-and +400-year adjustments, respectively, supported by extant volcanic evidence and radiocarbon ages.Our chronological framework shows a synchronous HS2 onset globally.Our records precisely characterize a centennial-scale abrupt "tropical atmospheric seesaw" superimposed on the conventional "bipolar seesaw" at the beginning of HS2, implying a unique response/feedback from low-latitude hydroclimate.Together with our observation of an early South American monsoon shift at the HS2 termination, we suggest a more active role of low-latitude hydroclimate dynamics underlying millennial events than previously thought.
The last glacial period was characterized by a recurrence of millennialscale cold-dry stadials and warm-humid interstadials, recognized from the Greenland ice-core records, known as Dansgaard-Oeschger oscillations 1,2 .Some Greenland Stadials have been associated with major iceberg discharge episodes and concomitant ice-rafted debris deposition in the central North Atlantic 2,3 .Previously marine studies have identified and characterized six such episodes known as Heinrich events during the last 60 thousand years 3,4 .The marine records show dramatic responses to Heinrich events for extended periods of time, often towards the end of the corresponding Greenland Stadial (e.g., ref. 5).We refer to these periods of extreme conditions as Heinrich Stadials (HSs), and note that they are not defined by the duration of the Greenland Stadial counterpart.Mid-to-low-latitude monsoon proxy records suggest a major weakening of the Asian summer monsoon (ASM) during the HSs 6,7 [hereafter Asian Heinrich Period (AHP)] and a strengthening of the South American summer monsoon (SASM) 8 [hereafter South American Heinrich Period (SAHP)].The weakening and strengthening of the two monsoon systems, respectively, are associated with the southward displacement of the Intertropical Convergence Zone (ITCZ) 6 .An abrupt positive (negative) isotope excursion during the AHP (SAHP) is well-expressed in the speleothem oxygenisotope (δ 18 O) records from the interhemispheric monsoon domains 8,9 (Supplementary Fig. 1c).It is noteworthy that Heinrich events have left a relatively small imprint in the Greenland ice-core δ 18 O records 10 .
Marine proxy records and numerical simulation experiments suggest that the Atlantic meridional overturning circulation (AMOC) weakened or even shut down during the HSs 11 .Iceberg discharge and associated freshwater flux almost led to a stop of the deep-water formation in the North Atlantic, which in turn weakened the AMOC 11 .Hence, the northward oceanic heat transport was reduced, resulting in Greenland cooling and corresponding Antarctica warming with a time lag of 122 ± 24 years 12,13 ; this phenomenon is known as bipolar or thermal seesaw 14,15 .It is also suggested that during Greenland's cooling and warming transitions, the global atmospheric teleconnections are characterized by a signal propagation from the northern high-latitudes to the tropics, and further, to the southern high-latitudes, which were synchronous within sub-centennial uncertainty [16][17][18][19] .Moreover, observational data points to an active role of hydroclimatic variations in the tropics and Southern Hemisphere around AHP and SAHP and other millennial-scale events 17,[20][21][22] .Recent research has highlighted the important role of the decreased Amazon River runoff and associated sea-surface salinity anomalies in resuming AMOC 6 , which in turn initiated the termination of AHP4.The noteworthy point is that the termination of SAHP4 [~38.62 ky BP (thousand years before present, where the present is 1950 CE)] predated the termination of AHP4 (~38.34 ky BP) and the associated Greenland warming transition (~38.34 ky BP) by hundreds of years 6 .Previous studies 6,17,20,22 point to a south-to-north signal propagation around the terminations of AHP and SAHP and other millennial-scale events suggesting a slow ocean response.Furthermore, during the early Last Glacial Maximum (LGM) (26 to 22 ky BP), the Greenland ice-core records have large age uncertainty (±500-±800 years) 2 and high-resolution absolutely-dated proxy records from low latitudes are still lacking 18 .These shortcomings largely impede the precise analyses of the climate signal propagation phase (north to south or vice-versa) inter-hemispherically.Consequently, it is yet unresolved whether the observed phasing remained consistent during the LGM 12,16,18 .Notably, LGM was an important period since it had a distinct boundary condition when global sea level and greenhouse gas (CO 2 and CH 4 ) concentrations were at their minima of the last glacial period (Supplementary Fig. 1) [23][24][25] , which may have influenced the hemispheric coupling.AHP2, SAHP2 and HS2 were typical millennial-scale events during the LGM (Supplementary Fig. 1) 3,26 , and it remains an important task to rigorously test the climate signal propagation and phase relationship between the events.Moreover, it remains poorly known how the global atmospheric circulation was tele-connected during HSs onsets, when there is a little-to-no imprint in Greenland ice-core δ 18 O.Therefore, an integrated understanding of hemispheric climate coupling including both high and low latitudes is still lacking.
In this study, we report nine speleothem δ 18 O records from ASM and SASM domains spanning 27 to 22 ky BP, including HS2 3,26 .We present speleothem δ 18 O records from the Indian summer monsoon (ISM) domain at a resolution of ~4 years.The age-model precision of our ISM δ 18 O records is <50 years and each record was obtained by combining annual lamina counting and 230 Th dating results (see Methods).The chronological precision and resolution of our ISM records facilitate the quantification of the chronological biases in Greenland ice-core records via comparison with the Greenland dust flux record (as reflected by ice-core [Ca 2+ ]).The robustness of the ISM speleothem-based Greenland chronology is further supported by extant volcanic evidence and radiocarbon dates.Speleothem δ 18 O records from the SASM domain are characterized by sub-centennial age-model uncertainties (<80 years) at the key intervals.Overall, we use speleothem δ 18 O and bipolar ice-core records to understand the causal link between the highand low-latitude climate systems in both hemispheres and to establish an interpretive framework of the rapid climate shifts during HS2/AHP2/ SAHP2, their characterization, and manifestation.
Speleothem chronologies are based on 127 230 Th dates (Supplementary Data 1).The age models of speleothem Cherrapunji-2 and Cherrapunji-2017-1 were obtained by annual lamina counting using confocal microscopy (Supplementary Fig. 3) in combination with 230 Th dates (Supplementary Fig. 4a, b, see Methods).The age models of speleothem YX-51, PX-07, MAG, BTV-4C, and NAR-C were obtained using the StalAge algorithm 38 (Supplementary Fig. 5, see Methods).Different age-modeling schemes yielded essentially identical age models (Supplementary Fig. 5).A total of ~2810 stable oxygen (δ 18 O) isotope data were obtained (Supplementary Data 2).The spatial resolution of the measurements varies between 0.05 and 1 mm (see Methods).The comparison between the δ 18 O records from the same and different caves in the same climatic region (Supplementary Fig. 6) suggests that the speleothem δ 18 O records broadly replicate, although there are minor differences in their absolute values (Supplementary Note 1.4).

Speleothem records in the ASM domain
The ASM is a vast climate system, composed mainly of the ISM and EASM subsystems, and closely coupled with the Asian westerly system 7,[39][40][41] .The timing of the AHP2 hydroclimatic transitions in the ASM domain is tightly constrained by the Cherrapunji δ 18 O record at high temporal precision (uncertainty <50 years) and resolution (~4 years) as well as its clear structure across AHP2 (Fig. 2a).work), and lasted for 76 ± 5 years (lamina counting uncertainty) (Supplementary Fig. 3a) before reaching the maximum value of AHP2 (Fig. 2a).A subsequent abrupt rebound commenced at 24.33 ± 0.02 ky BP and lasted 139 ± 2 years with a ~1‰ δ 18 O decrease (Supplementary Fig. 3b), which, together with the abrupt onset, manifests as a positive δ 18 O excursion during AHP2 (referred to as stage III) (Fig. 2a).Following the rebound a relatively stable period (stage II) occurred between 24.19 ± 0.02 ky BP and 23.87 ± 0.02 ky BP (Fig. 2a).The termination (stage I) was characterized by a gradual δ 18 O decrease of ~1.2‰, which was initiated at 23.87 ± 0.02 ky BP and lasted for ~350 years (Fig. 2a).In order to objectively select the speleothem records in the ASM domain, we used the Speleothem Isotopes Synthesis and Analysis (SISAL) database 42,43 to choose the speleothem δ 18 O records that have either a temporal resolution better than 40 years or a sub-centennial age uncertainty, and some high-quality ASM speleothem δ 18 O records that are not in the database were also used in this study (Supplementary Fig. 7).The comparison between the Cherrapunji δ 18 O record and other ASM domain speleothem records shows that their overall patterns and transitional timings across AHP2 are coherent (Fig. 2b and Supplementary Fig. 7).Notably, however, the stage III excursion is only evident in the records in the ISM domain and the transitional zone between the ISM and EASM domains, but is absent in the EASM domain (Fig. 2b and Supplementary Fig. 7), as discussed in following sections.In this study, significant changes in the isotopic profiles were calculated using the "Ramp-fitting" 44 and "BREAKFIT" methods 45 (see Methods; Supplementary Fig. 8 and Supplementary Table 1).However, these two methods do not consider the age model uncertainty, thus in our analyses, we quadratically combined the change point uncertainty and the age model uncertainty and obtained a combined uncertainty (see Methods; Supplementary Table 2).The timing of the significant changes in the isotopic file of our raw dataset is further supported by the "Trend-fitting" function analyses (see Methods; Supplementary Fig. 9).We also performed these calculations to identify key change points in bipolar ice-core records (Supplementary Figs. 9, 10 and Supplementary Table 2).

Greenland ice-core [Ca 2+ ] records
Previous studies have shown that the Taklimakan and Gobi deserts, hereafter referred to as the Asian dust source regions (Supplementary Fig. 11), were the main contributors of mineral dust transported to Greenland during 27-23 ky BP 46,47 (Supplementary Note 1.6.1).In modern times, the main dust emission season in the Asian dust source region is boreal spring 48 .Dry conditions and strong winds in the dust source areas favor the entrainment of dust into the upper atmosphere, which can be long-distance delivered to the downwind North Pacific and Greenland via westerly jet 48,49 .During boreal spring, the upper westerly jet is strong and its core axis is located south of the Tibetan Plateau, while its meridional range is the largest 50 (Supplementary Fig. 11a).This scenario is distinct from boreal summer, when the upper westerly jet weakens and its core axis migrates to the north of the Tibetan Plateau, allowing the low-level monsoonal moisture to reach inland areas 51 (Supplementary Fig. 11b).Paleoclimatic studies further showed that during stadials, the cooling in the Northern Hemisphere (NH) weakens the seasonal contrast, extending the dust season from spring to summer and suggesting that the jet axis maintained above the southern part of the Tibetan Plateau for a longer time 41,52,53 , which enabled more frequent dust emission and transport to Greenland.Based on these studies, we suggest that the variability of the Greenland ice-core mineral dust-derived [Ca 2+ ] (a dust proxy 54 , Supplementary Fig. 12) at decadal-to-millennial timescales during 27-23 ky BP primarily reflects the latitudinal position and intensity of the NH westerly winds, as well as the hydroclimate conditions in Asian dust source regions, consistent with the previous interpretation 55 .Such an argument gains support from the similarities between Greenland [Ca 2+ ] and an Asian westerly domain speleothem record from Turkey (So-1 record from Sofular Cave) for AHP2 (Supplementary Fig. 12).During the AHP2 onset, the So-1 δ 13 C record exhibits an abrupt positive transition, indicating a change to the colder/drier condition 56 , in line with a southward shift of the Westerlies.In the following discussion, Greenland ice-core [Ca 2+ ] is qualitatively used as a proxy for mid-latitude westerly winds in the NH and hydroclimate conditions in the Asian dust source regions, which is different from the Greenland ice-core δ 18 O records that mainly reflect NH high-latitude changes in temperature and precipitation seasonality.
Correlation between Greenland ice-core [Ca 2+ ] and ISM δ 18 O records Greenland ice-core δ 18 O records do not exhibit a distinct pattern during HS2 (Fig. 3a).One of several possibilities is that due to an AMOC slowdown, the cooling signal (lower δ 18 O) was obfuscated by the absence of winter precipitation in response to sea-ice expansion, as suggested by transient experiments for HS1 57 .The blurred signal precludes a direct correlation between Greenland ice-core δ 18 O records and the more precisely-dated ASM speleothem δ 18 O records.Alternatively, the Greenland ice-core [Ca 2+ ] record, as a dust proxy dominated primarily by Asian westerly winds that are strengthened rather than muted by sea-ice expansion 58 , exhibits a clear structure across Greenland Stadial 3 (Fig. 3b).Crucially, the westerly winds are closely coupled to the ASM circulation and tracked by speleothem δ 18 O 40,59 on seasonal-toorbital timescales 7,40,41,50,51 .We rely on this dynamic link to correlate Greenland ice-core [Ca 2+ ] to ASM speleothem δ 18 O records, assuming no time offset from lagged system responses at the decadal timescale.Indeed, Lidar data analyses and model simulations indicate that the transportation of modern dust is a fast process 60 .Modern Greenland snow-pit sample analyses and observations of dust-storm activities in the Asian dust source regions have further shown that no time offset exists between dust emission in the source region and dust deposition in Greenland at the interannual timescale (Supplementary Note 1.6.2).
A resemblance exists between Greenland ice-core [Ca 2+ ] and Cherrapunji δ 18 O records, which allows tuning the ice-core [Ca 2+ ] timeseries [on the widely used Greenland Ice Core Chronology 2005 (GICC05)] to the 230 Th-dated Cherrapunji δ 18 O record using prominent peaks and abrupt changes as tie points (Fig. 3 and Supplementary Table 3).These tie points have passed sensitivity tests (see Methods) and have been verified by the "Trend-fitting" method (Supplementary Fig. 9a-c), and are hence regarded as statistically robust.Furthermore, five of the tie points in the detrended z-score transformed Cherrapunji δ 18 O profile occurred at intervals with low z-score values (<−1) (Supplementary Fig. 13c), which we interpret to result from a higher frequency of extreme pluvials.Another tie point is at the start of AHP2, from which the Cherrapunji δ 18 O values increase abruptly to the high zscore values (>4) (Supplementary Fig. 13c), which we interpret to result from a higher frequency of extreme droughts.These lines of evidence show that the six tie points correlate with large changes in Asian hydroclimate.Based on these tie points, the result suggests a shift of the GICC05 timescale by +320 years (Fig. 3), which is well within its quoted age uncertainty (the maximum counting error, a measure of the potential bias around this time interval is 600-800 years 2 , Supplementary Fig. 14) and consistent with previous studies that also suggested that the GICC05 chronology needs to be shifted by several hundred years towards older ages 61,62 .Crucially, the +320-year shift is based on six tie points, much more precise than previous studies, which are merely based on one tie point with a larger age uncertainty (Supplementary Fig. 14).Moreover, although other millennial-scale ASM events (e.g., HS4/AHP4 6 and Younger Dryas 17 ) apparently correlate with the ice-core [Ca 2+ ] records (Supplementary Fig. 1), the dynamical mechanism underlying the correlation was not well explored in the previous studies.More importantly, no annually laminated speleothem record covers the full duration of 27-23 ky BP as of now, making it hard to precisely constrain the relative chronology, i.e., event durations.In this regard, considering that the GICC05 timescale was confirmed to be precise within 20-40 years of uncertainty during Greenland Stadial 1 (~12.9-11.7 ky BP) 17 , it is thus highly possible that ice layers were undercounted during Greenland Stadial 2 (~23-15 ky BP).
Applying the +320-year shift, the Greenland ice-core chronology has an uncertainty of 90 years (2σ), which comes from the age model errors of the Cherrapunji record plus a few other possible factors (Supplementary Note 1.7).We consider this uncertainty as a conservative estimate.Although ice-core [Ca 2+ ] and speleothem δ 18 O records show an impressive match with each other, differences exist.For example, the Cherrapunji δ 18 O record exhibits more gradual transitions compared with ice-core [Ca 2+ ] records, and the "double-spiked" Dansgaard-Oeschger (DO)-2 structure in the Cherrapunji δ 18 O record is relatively small (Fig. 3 and Supplementary Note 1.6.3).Clarifying these differences warrant further studies.Nevertheless, we argue that the abrupt shift of the large-scale atmospheric circulation (as reflected by the six tie points) was synchronous within multi-decadal uncertainty, especially considering that the ASM and the Asian westerlies are typically tightly-coupled 7,40,41,50,51 .This is also consistent with our proposition that the GICC05 timescale needs to be shifted as a whole between 27-23 ky BP without being stretched or compressed, which takes advantage of the annual layer-counted chronology of GICC05 and Cherrapunji record, ensuring the robustness of their direct comparisons.

Correlations to Antarctic ice-core records
The highest-resolution CH 4 24 and ice δ 18 O records 13 available around Antarctic Isotope Maximum 2 were derived from the West Antarctic Ice Sheet (WAIS) Divide Ice-core (WDC) on a recent chronology (WD2014 63 ).Remarkably, a unique volcanic triplet spike was identified in both Greenland and Antarctic ice-cores 12 , possibly associate with the Oruanui eruption from the Taupo volcano (Supplementary Fig. 2) 12 .After shifting the GICC05 chronology by +320 years, the volcanic triplet spike occurred at 24,939 ± 90 y BP in Greenland (Fig. 4), providing a robust tie point that requires the WD2014 ice chronology (which shows the volcanic triplet spike at ~24,539 y BP) to be shifted by +400 years (320 years according to the GICC05 shift in this study and 80 years based on a shift suggested previously 12 ) (Fig. 4).Moreover, after the +400-year shift, the tephra layer of the Oruanui eruption in WDC would occur at 25,718 y BP 64 , consistent within uncertainty with the calibrated radiocarbon-age of 25,675 ± 90 (1σ) calendar y BP (using the SHCal20 calibration curve 65 ) 66 .Furthermore, the WD2014 gas chronology is also shifted by +400 years, particularly considering the small uncertainty of the WDC ice-gas age difference (Δage) at this time (~110 years) 63,67 .Ultimately, the adjusted chronology is linked to our speleothem chronology, and thus, should have a similar precision of ~90 years (2σ, ice chronology) provided that the volcanic triplet is correctly aligned.WD2014 is an annual layer-counted chronology between 31.2-0ky BP 63 , which does not allow being stretched or compressed beyond the expected error, consistent with the shift of WD2014 as a whole between 26-23 ky BP.Similarly, the WDC ice layers were possibly undercounted during the late LGM.As a result, this tuning allows a precise comparison between bipolar ice-core records and mid-tolow-latitude speleothem records on a common chronology.We acknowledge the potential existence of other matching strategies and the risk of circular reasoning.However, because abrupt transitions of the large-scale atmospheric circulations as well as speleothem and radiocarbon chronologies are aligned so well and are connected (Figs. 3, 4), our tuning process must have high reliability within our stated uncertainty (~90 years) and constructs the hitherto strongest constraint on bipolar ice-core chronologies between 27-23 ky BP.
Although bipolar ice-core chronologies have been improved, the analyses of the phase relationship between Greenland Interstadial 2 warming transition and Antarctic Isotope Maximum 2 cooling transition cannot yet be resolved (Fig. 4).This is because the small signal-tonoise ratio in the Antarctic ice-core δ 18 O record 13 impedes a robust identification of change points (see Methods).Nevertheless, the precise chronologies provide important constraints to resolve the problem as soon as Antarctic ice-core δ 18 O records with higher resolution are available.2): NGRIP, GISP2, GRIP, and NEEM records (Supplementary Note 1.5).c Cherrapunji speleothem δ 18 O record (this study).Greenland ice-core records are shown on the GICC05 chronology (Supplementary Note 1.5), which is shifted by +320 years via tuning the Greenland [Ca 2+ ] time-series to the Cherrapunji δ 18 O record (see main text).Vertical red dashed lines depict the tie points between the ice-core and speleothem records (Supplementary Table 3).The leftmost two dashed lines correspond to Greenland Interstadial (GI)-2.1 and GI-2.2 peaks in Greenland 2 and their counterparts in the Asian summer monsoon regime (Dansgaard/Oeschger (DO)-2.1 and DO-2.2 peaks).Vertical gray bars depict the Greenland Stadial (GS)-3 dust peaks 94 .Error bars show uncertainty (2σ) of the 230 Th chronology of the Cherrapunji record (green) in comparison to the possible age bias of the GICC05 timescale 2 .

AHP2/SAHP2 onsets and the large excursion
A recent study 68 has established age models for the 92 published marine sediment records from the Atlantic Ocean via correlating to the Greenland ice-core GICC05 chronology.This, in turn, allows us to simply shift the chronologies of these marine records by +320 years across the HS2 period (Supplementary Fig. 15 and Supplementary Note 1.8) to compare the marine records with the precisely-dated speleothem records.
It has been known that global atmospheric shifts are synchronous within decades at Greenland cooling and warming transitions, with a signal propagation from north to south [16][17][18][19] .However, only limited data are available about the global atmospheric teleconnections at HSs onsets when there are relatively small changes in Greenland δ 18 O records (e.g., ref. 69).In our temporal framework, the AHP2 onset at 24.43 ± 0.05 ky BP (combined uncertainty) inferred from the Cherrapunji δ 18 O record is synchronous within sub-centennial uncertainties with changes in the SASM hydroclimate (Fig. 5j and Supplementary Table 2).The abrupt and anti-phase hydroclimatic changes in the monsoon domains of the two hemispheres reinforce the hypothesis (e.g., ref. 70) that iceberg discharge or freshwater forcing originating from the NH subtropics could push the tropical rain-belt far southward.A positive anomaly can also be observed in the Δε LAND record (Fig. 4d), which also supports a southward shift of the tropical rainfall and terrestrial oxygen production during the HS2 71,72 .This is also consistent with the observation that the weakening of the AMOC coincides with the waning of the ASM within the marine age uncertainties (Supplementary Fig. 15), similar to the results of water-hosing experiments 73 .Additionally, the large and abrupt transition provides a direct example of a fast and widespread monsoon climate swing that occurred on a decadal timescale, as constrained by lamina counting (i.e., onset in ~75 years, Supplementary Fig. 3a).Furthermore, we compared our records with proxy records from the Southern Hemisphere representing mid-to high-latitude atmospheric circulations and/or different moisture origins 16,19,74 (Supplementary Note 1.9), i.e., Antarctic ice-core d ln (the logarithmic definition of deuterium excess) and non-sea-salt soluble Ca 2+ ([nssCa 2+ ]) records (Fig. 5h, i).The results show synchronous changes within uncertainties (Fig. 5j and Supplementary Table 2), suggesting parallel changes in Southern Hemisphere high-to low-latitude atmospheric circulations and NH mid-latitude atmospheric circulations (Fig. 5b-i).This synchronicity hints towards a fast (decadal-scale) atmospheric teleconnection, involving meridional shifts of the ITCZ and the mid-latitude westerly wind belts in both hemispheres, changes in the tropical Hadley circulation as well as monsoon circulations, which is dynamically consistent with the interhemispheric atmospheric teleconnections in model simulations 75,76 .
After the abrupt AHP2/SAHP2 onsets, a rapid rebound is observed in Cherrapunji, PA-LA-1, MAG and PX-07 δ 18 O records (Fig. 5c-f).Notably, three of the four records are derived from tropical regions (Fig. 2f, g), and Cherrapunji Cave is close to the tropical regimes (Fig. 2f) with a major moisture trajectory from tropical oceans during summer monsoon months 36 .Further comparison shows that the rebound feature is absent in Southern Hemisphere mid-to highlatitude atmospheric circulation records (Fig. 5h, i), and the Southern Hemisphere subtropical record from Botuverá Cave (BTV-4C record) also does not exhibit a clear rebound structure (Fig. 5g).In the NH, Greenland [Ca 2+ ] records do not exhibit a rebound feature (Fig. 3b), and the same is true for speleothem δ 18 O records from the EASM domain (Fig. 5b and Supplementary Fig. 7a-e).Additionally, the NGRIP deuterium-excess record, a proxy for the past changes in evaporation conditions or shifts in moisture sources in the NH mid-to highlatitudes 69 , does not show a rebound feature (Fig. 5a).We, therefore, identify simultaneous changes in low-latitude regions without any fingerprint in mid-to high-latitude climates.The direction and speed of the rebound are nearly opposite to the onset transition, corresponding to a breakpoint during stage III (Fig. 5c-f).Change point detection results show that the breakpoints in our records (Cherrapunji, MAG and PX-07) are coherent within sub-centennial uncertainties (Fig. 5c, e, f and Supplementary Table 2), hinting at rapid atmospheric teleconnection within the tropics.The δ 18 O reversal in our monsoonal  13 and EDML 95 , respectively.Antarctic ice-core records are plotted on the WD2014 chronology 63 , which is shifted by +400 years via a correlation based on the volcanic triplet spike (depicted by the orange diamond and the vertical bar at ~24.94 ky BP).It is noted that the EDML δ 18 O record has been tuned to the WD2014 chronology in a previous study 16 .Red dots and associated error bars in (e) and (f) indicate the timing and combined uncertainties (2σ) of the change points in Greenland and Antarctic icecore δ 18 O records (Supplementary Table 2).The Brown error bar depicts the 14 C age 66 and uncertainty (1σ) of the Oruanui eruption from the Taupo volcano.The vertical bar at ~25.72 ky BP marks the tephra layer of the Oruanui eruption recorded in the WDC ice core 64 .The two vertical dashed lines depict the DO-2.1 and DO-2.2 peaks.Cave, ice-core and volcano locations are shown in Supplementary Fig. 2. Vertical bars are the same as in Fig. 2. AIM Antarctic isotope maximum 95 , ppb parts per billion, ISM Indian summer monsoon, SASM South American summer monsoon.
record is rapid, which occurred within decades as constrained by annual lamina data (Fig. 5k), manifesting the existence of a tipping point in the tropical climate system.We refer to this oscillation as a "tropical atmospheric seesaw".Oceanic processes would last longer than a few decades 17 .In this regard, the tropical atmospheric seesaw involves synchronous and opposite changes (positive δ 18 O excursion in Cherrapunji and PA-LA-1; negative δ 18 O excursion in MAG and PX-07) in monsoon domains of both hemispheres (Fig. 5c-f), possibly associate with the meridional shift of the ITCZ.It is noteworthy that oceanic changes would accompany atmospheric changes in a coupled climate system, and our results may help to differentiate the roles of oceanic processes and atmospheric processes mentioned in this study, which is vital for understanding climate dynamics 13 .Although Greenland ice-core [Ca 2+ ] and EASM speleothem δ 18 O records show a general coherent pattern with Cherrapunji δ 18 O record across much of 27-23 ky BP, the NH mid-latitude changes were decoupled with ISM domain during stage III (Figs. 2b, 3b and Supplementary Fig. 7a-e).The mid-latitude records lack the prominent Cave (Cherrapunji, this study), Larga Cave (PA-LA-1 96 ), Marota Cave (MAG, this study), Paixão Cave (PX-07, this study), and Botuverá Cave (BTV-4C, this study), respectively.h Antarctica 5-core averaged d ln anomaly 16 .i Antarctica WDC ice-core non-sea-salt Ca 2+ [nssCa 2+ ] record 74 and overlain by the thick line (high-frequency variability (>1 cycle per 300 year) removed by a low-pass Butterworth filter, same as in ref. 74).Antarctic ice-core records in (h) and (i) are on the WD2014 chronology and shifted by +400 years.The yellow dashed lines in (c-f) are generated via the BREAKFIT algorithm 45 .Black dots and associated error bars in (c-f) depict the timing of breakpoints and combined uncertainties (2σ) (Supplementary excursion, which is exhibited when comparing the YX-51 and Cherrapunji records (Fig. 2b) (Supplementary Note 1.11).After stage III, the NH mid-to-low-latitude atmospheric circulations coupled again (Figs.2b,  3b).This suggests that during AHP2 and SAHP2 onsets, the tropical atmospheric system amplified the disturbances from the north, and subsequently triggered the rebound after reaching a tipping point in the tropical hydroclimate.Moreover, the prominent excursion occurred during the warming phase of Antarctic Isotope Maximum 2 (Fig. 4f, g), while Greenland temperatures maintained low (Fig. 4e).Meantime, changes in low-latitude regimes suggest large-scale reorganizations of the atmospheric circulation, probably associate with meridional ITCZ shifts.In this context, the conventional bipolar seesaw concept 15 needs to be expanded in order to fully reflect the complexity of processes at play, including both high and low latitudes and various boundary conditions.A previous study using a heuristic model of tropical rainfall distribution hypothesized that the CH 4 overshoots within HS1, HS2, HS4, and HS5 correlate with the onset of Heinrich events 24 .Although subsequent studies confirmed the correlation during some HSs by speleothem records of high-resolution and high-precision 6 , the phase relationship between the ~40 ppb (parts per billion) CH 4 increase and Heinrich event 2 remains unclear (Supplementary Fig. 16).It is noteworthy that the previous relationship was derived by comparing the small WDC CH 4 peak (~15 ppb) with the GSIP2 δ 15 N record (both in the gas phase, Supplementary Fig. 16), assuming that the δ 15 N peak represents the DO-2 peak 24 .The ice-core δ 18 O record, however, shows two abrupt warmings and associated peaks 2 (Fig. 4e), unlike the single DO-2 peak as observed in previous low-resolution data.Besides, the GISP2 δ 15 N record bears very low resolution 2 (Supplementary Fig. 16).Therefore, the result would be different when matching the small WDC CH 4 peak with the DO-2.1 or DO-2.2 peak (Supplementary Fig. 16).Indeed, if based on the CH 4 record on the WD2014 chronology (i.e., without the +400-year shift proposed in this study), the abrupt CH 4 increase would correlate with the termination of AHP2 rather than its onset (Supplementary Fig. 16c, d), which is distinct from the previous assumption 24 .In our correlation framework, we confirm that the small CH 4 peak correlates with the DO-2.2 peak (Fig. 4a, b and Supplementary Fig. 16h, i).Moreover, the abrupt CH 4 increase (~40 ppb) coincides with stage III, the abrupt δ 18 O shift shared by our records from the ASM and SASM domains (Fig. 4a-c and Supplementary Fig. 16h-j).This common signal provides compelling evidence that enhanced rainfall in the Southern Hemisphere tropical wetlands indeed coincided with a significant CH 4 increase early in AHP2, favoring the previous hypothesis 24 .

AHP2/SAHP2 terminations
During AHP2/SAHP2 terminations, a notable trait is a variance in trends between ASM and SASM records across stages I and II.The SASM speleothem δ 18 O records exhibit a long-term trend across stage II showing a δ 18 O increase (drying) in South America (Fig. 6a-d).After the termination of SAHP2 (~23.85-24.05ky BP), all SASM speleothem δ 18 O records exhibit maximum aridity, comparable to or even drier than pre-SAHP2 conditions (Fig. 6a-d and Supplementary Fig. 17g-j).In contrast, the ASM was still weak during stages I and II compared to the pre-AHP2 conditions, and AHP2 termination occurred ~300 years after the SAHP2 termination (Fig. 6f, g and Supplementary Fig. 17).Other records from the SASM domain (Supplementary Fig. 18c-e), however, lack sufficient 230 Th dates control around the SAHP2 in order to evaluate the transitional timing at sub-centennial precision.Nevertheless, the long-term drying trend during the SAHP2 termination, as observed in previous records (Supplementary Fig. 18c-e), is overall consistent with our observation (Supplementary Fig. 18a, b, f, g).The long-term drying trend in the SASM domain was possibly associated with a cooling of the tropical western equatorial Atlantic 77 .Additionally, a few centennial-scale events are superimposed on the SASM long-term drying trend.Of note is a large dry event (inferred by a positive δ 18 O excursion) in stage II in three SASM speleothem δ 18 O records [(MAG and PX-07, this study, Fig. 6a, b) and arguably LSF-3 (Supplementary Fig. 18c)], in contrast to a rather stable climate in the ASM realm (Fig. 6e).Indeed, MAG, PX-07, and LSF-3 records are from caves located close to the tropical South Atlantic (Supplementary Fig. 2), and thus might be sensitive to changes in the surface thermohaline circulations nearby.Nevertheless, this dry event does not affect our conclusion of the much earlier termination in the SASM domain compared to the ASM domain, as demonstrated by various change point detection methods and sensitivity tests (see Methods, Supplementary Figs. 8, 9).
The observed long-term drying in the vast SASM domain would have caused a decrease in Amazon River discharge into the Atlantic Ocean 6 (Supplementary Fig. 19b).Reduced freshwater input may have induced a positive sea-surface salinity anomaly in the Amazon Plume Region 78,79 , ultimately advected to the deep-water formation areas in the North Atlantic 80 (Supplementary Fig. 19b) and contributing to the strengthening of the AMOC 81,82 .The sea-surface salinity reconstructed for the eastern subpolar North Atlantic (MD95-2002 Core) indeed exhibits an increasing trend during the AHP2 termination (Supplementary Figs. 2, 15e).A strengthened AMOC would induce positive feedback via transporting more saline water to the north, facilitating northward heat transport and a northward shift in the tropic rain-belts, intensified ASM 6,17,22 , and Greenland warming 14 , consistent with the proxy data for stage I (Fig. 4a, e).Therefore, the temporal phasing established here shows that the initiation for the long-term drying in the SASM domain occurred hundreds of years before the AHP2 termination (the start of stage I) (Fig. 6f and Supplementary Fig. 17), consistent with marine records (Supplementary Fig. 2) and shedding light on the termination dynamics.The lead of the SASM in the termination process is prominent, which is based on the assessment of the combined uncertainties of proxy records (Fig. 6f, g, Supplementary Fig. 17, and Supplementary Table 2).This indicates that long-term drying in the SASM domain and associated positive sea-surface salinity anomalies in the Amazon Plume Region may have acted as a precursor for the AHP2 termination, which reinforces the hypothesis 6 recently proposed for the AHP4 termination, albeit having a different boundary condition compared to AHP2 (Supplementary Fig. 1).On the other hand, in the current theoretical framework, the hydroclimatic events surrounding the earlier termination of SAHP2 may have been causally linked to the hydroclimatic swerve (or the end of the Antarctic Isotope Maximum 2 warming transition) observed in Antarctic ice-cores (particularly the EDML ice core from the Atlantic sector) (Fig. 4 and Supplementary Fig. 19) (e.g., refs.6,83), but the cause of the earlier Antarctic change itself requires further investigation.
In conclusion, our study provides a causal relationship between ASM speleothem δ 18 O records and Greenland dust records, improving the bipolar ice-core chronologies by nearly an order of magnitude.Our observations support the notion that rapid atmospheric changes in the tropics is not merely a response to forcing from the north, since the amplification effect in the tropical atmospheric system is vital for the distinct isotope excursion, which is absent in mid-to high-latitudes.In terms of AHP2/SAHP2, the signal of tropical changes could be transported to the high-latitudes of both hemispheres only when the oceanic mode is dominant, such as a resumption of the AMOC.The emerging picture is that both atmospheric and oceanic processes in the tropical regions and the Southern Hemisphere played active roles during millennial-scale event terminations, superimposed on the bipolar seesaw.In this regard, our findings extend the bipolar seesaw theory in both temporal (including LGM) and spatial domains (including mid-to-low latitudes), which is vital for an integrated understanding of hemispheric coupling during rapid climate changes.This is also important for model simulations seeking to capture the climate dynamics of abrupt climate changes like Heinrich events.

Paleoclimate records
We reconstructed and considerably improved nine speleothem δ 18 O records.The sample information is summarized in Table 1.Cherrapunji Cave (~1100 m above sea level) and Mawmluh Cave are located at Cherrapunji on the southern fringe of the Meghalayan Plateau in Northeast India (Fig. 2f).The mean annual air temperature recorded at the nearby meteorological station in Cherrapunji is 17.4 °C.The mean annual rainfall in this region is ~12,000 mm, 70% of which falls during the peak ISM months (June to September) 84 (Supplementary Fig. 7i).Dongqinghe Cave (~1600 m above sea level) is located at Guangyuan city, Sichuan province, Southwest China (Fig. 2f).The mean annual air temperature recorded at the nearby meteorological station is 16 °C.The mean annual rainfall in this region is ~1100 mm, 70% of which falls during the summer monsoon months (June to September) (Supplementary Fig. 7i).Other caves are also located in the monsoon domains of both hemispheres 8,[27][28][29][30] .

Proxy interpretation
The simulation work by ref. 36 examined the South Asian summer monsoon and showed that the interannual speleothem δ 18 O variability in India is strongly tied to the overall monsoon circulation.The intraseasonal to interannual rainfall variability over the Indian subcontinent exhibits a quasi-east-west precipitation dipole with anomalies of one sign over Northeast India and of an inverse sign over North, Northwest, and Central India 36,85 .The dynamical constraints 36 explained the situation that the oxygen-isotope composition of precipitation (δ 18 O p ) of summer monsoon rainfall at Northeast India, albeit being located at the opposite end of the precipitation dipole, reflects upstream changes in monsoon precipitation amount over North, Northwest, and Central India (~15-28°N and 70-84°E) (via the moisture source effect instead of the "classical" amount effect).Following these reasonings and given the resemblance between Cherrapunji δ 18 O and MWS-1 δ 18 O (from Mawmluh Cave, which is located in the same climate region as Cherrapunji Cave) records (Fig. 2f and Supplementary Fig. 6d), we interpret low and high δ 18 O values in the Cherrapunji record to reflect strong and weak large-scale monsoonal circulation, respectively.At millennial timescale, the interpretation of other speleothem δ 18 O series utilized herein are consistent with the intensity of large-scale monsoonal rainfall/circulation and the northsouth shifts of the ITCZ 7,8,[27][28][29][30] .Supplementary Note 1.3 gives a detailed description of the ASM domain speleothem δ 18 O interpretation.

Th dating method
Subsamples for 230 Th dating were drilled on the polished speleothem section using a 0.3-mm carbide dental drill.We used standard chemistry procedures 86 to separate U and Th.A triple-spike ( 229 Th-233 U-236 U) isotope dilution method was used to correct instrumental fractionation and to determine U-Th isotopic ratios and concentrations 87,88 . 230 Th dating was performed at Xi'an Jiaotong University, China, using a Thermo-Finnigan Neptune Plus multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS).U and Th isotopes were  2), red lines in (a) and (c-e) show the ramps and break lines generated by measured on a MasCom multiplier behind the retarding potential quadrupole in the peak-jumping mode using standard procedures 87 .Uncertainties in U and Th isotopic measurements were calculated offline at the 2σ level, including corrections for blanks, multiplier dark noise, abundance sensitivity, and contents of the same nuclides in the spike solution.The most recent values of the decay constants of 234 U 88 and 230 Th 88 and 238 U 89 were used.Corrected 230 Th ages assume an initial 230 Th/ 232 Th atomic ratio of (4.4 ± 2.2) × 10 −6 and those are the values for material at secular equilibrium with the bulk earth 232 Th/ 238 U value of 3.8 88 .The corrections for samples in this study are small because their uranium concentrations are high (~200-2500 ng/g) and detrital 232 Th is low (less than 4000 pg/g) (Supplementary Data 1).Moreover, some studies applied an initial 230 Th/ 232 Th ratio of (8.8 ± 8.8) × 10 −6 to provide a maximum uncertainty estimate for speleothem datasets all over the globe 18 .This value is too large for our clean samples, nevertheless, we recalculated our 230 Th ages using the 230 Th/ 232 Th ratio of (8.8 ± 8.8) × 10 −6 (Supplementary Data 1).The comparison shows that the recalculated ages are broadly similar to the original ages for our speleothems within uncertainty and do not show considerable changes (Supplementary Data 1).Therefore, in this study, a 230 Th/ 232 Th ratio of (4.4 ± 2.2) × 10 −6 is applied to the final age models.

Annual lamina counting
In order to construct a precise floating lamina chronology (counting chronology), the annual nature of the laminae needs to be established.This can be achieved by comparing the number of laminae counted between absolutely-dated ages.After polishing, samples Cherrapunji-2 and Cherrapunji-2017-1 have clear lamina bands observed using confocal laser fluorescence microscopy (CLFM) (Supplementary Fig. 3).In this study, we utilized a Nikon A1-1024 instrument from the State Key Laboratory for Manufacturing Systems Engineering, Xi'an Jiaotong University.A 40 mW, 488 nm blue laser line was used for excitation 90 , and fluorescence images were collected using an emission filter which allows light with wavelengths between 500 and 550 nm (visible, green) to pass through 90 .Lamina counting was conducted five times (Supplementary Note 1.2).The result shows that the laminae can be continuously identified, and the numbers of laminae between different depths are identical to the age difference of the respective 230 Th dates within analytical uncertainty (Supplementary Data 1).This agreement supports our interpretation that the laminae are annual 91 and allows constructing a precise age model 92 .

Age models
Age models for annually laminated speleothems were reconstructed using the floating lamina chronology anchored by 39 and 9 highprecision 230 Th dates for samples Cherrapunji-2 and Cherrapunji-2017-1, respectively (Supplementary Data 1 and Data 2) (Supplementary Fig. 4), using the least-squares fitting method to get the best fit of the 230 Th dates to the floating lamina chronology 92 .These age models have an estimated maximum age uncertainty of ±21 years (2σ, 99.5-329.5 mm in Cherrapunji-2), ±23 years (2σ, 331.5-437 mm in Cherrapunji-2), and ±43 years (2σ, 40.5-56 mm in Cherrapunji-2017-1) (Supplementary Data 2).The age model constructed using this method considered counting uncertainty and the uncertainty when anchoring the floating lamina chronology to the absolutely-dated 230 Th dates 92 .This method does not depend on the analytical uncertainties of 230 Th dates, thus the uncertainty of the age model is generally smaller than the dating uncertainty (Supplementary Data 2).We used StalAge algorithm 38 to construct the age models for YX-51, MAG, PX-07, BTV-4C and NAR-C (Supplementary Fig. 5).StalAge is particularly suited for speleothems creating objective age models based on two assumptions: (1) the age model is monotonic and (2) a straight line is fitted through all data or through as many data points as possible within error bars 38 .StalAge produces 300 realizations of age models by the Monte-Carlo simulation to account for the 95% confidence limits 38 .Major outliers are identified by disagreement with at least two data points, and minor outliers are screened if more than 80% of the simulated straight lines fail to have a positive slope.In our case, no major or minor outliers were detected because all ages in each age model increase monotonically within dating uncertainties (Supplementary Fig. 5).The uncertainty is large (more than 150 years) in the age model boundary of sample BTV-4C, YX-51 and NAR-C, due to limited dating controls.Nonetheless, the key change points (the onset and termination of AHP2/SAHP2) are robust, which is away from the age model boundaries and well constrained by several 230 Th dates (Supplementary Fig. 5).Moreover, age models of our speleothem records were also calculated using other age-modeling schemes (Supplementary Fig. 5).All these schemes yielded nearly identical results and the conclusions of this study are thus insensitive to the choice of the age model (Supplementary Fig. 5).

Change point determination
In order to objectively identify change points in various records, we used the Ramp-fitting 44 and BREAKFIT 45 algorithms.The choice of method is based on the shape of the time-series (Supplementary Table 1)."Ramp-fitting" finds two change points (t 1 and t 2 ), it estimates two constant levels of pretransition (t > t 2 ) and post-transition (t < t 1 ) and fits a ramp between them; Each data series was calculated for 3 × 10 6 steps using a Markov Chain Monte-Carlo (MCMC) sampler.These features make Ramp-fitting suitable for BTV-4C, Antarctica d ln and Greenland δ 18 O records, as well as the AHP2/SAHP2 termination process in Cherrapunji and MAG δ 18 O records.Both EDML and WDC δ 18 O records are characterized by an initial increase, a plateau, and subsequent decrease across the Antarctic Isotope Maximum 2 (Fig. 4), suitable for the Ramp-fitting routine in this case.Ramp-fitting also evaluated the influence of the addition of autocorrelated noise on the identified transitions, which in turn creates a more conservative estimate (with larger uncertainties) in the ramp parameters than other similar change point detection methods.BREAKFIT finds one change point and fits a linear slope on either side; these features make it suitable for the large excursion in Cherrapunji, MAG, PX-07 and PA-LA-1 δ 18 O records, as well as the determination of tie points between Cherrapunji δ 18 O and Greenland [Ca 2+ ] records.The BREAKFIT algorithm can provide statistical uncertainties of the timing of breakpoints using 2000 block bootstrap simulations, it also considers autocorrelation coefficients for the case of uneven time spacing 45 .In this study, we report the uncertainties of change points based on Rampfitting and BREAKFIT methods, and the analytical time intervals are shown in Supplementary Table 1.The main criteria to choose analytical time intervals are as follows: (1) the interval contains two prominent change points for Ramp-fitting and (2) the same time intervals are used for records from the same region if possible.To assess the sensitivity of the identified change points, we performed tests in which we randomly change the search time windows (Supplementary Figs. 8,  10).The results are regarded as robust only in the case where these techniques provide an unequivocal solution, i.e., the timing of the identified change points do not vary by more than 60 years when the width of the search time window changed (Supplementary Fig. 8).The results show that most change points are insensitive to the choice of search time windows, one exception is the Antarctic Isotope Maximum 2 cooling transition in Antarctic ice-core δ 18 O records where the identified change points do not meet the aforementioned criteria (Supplementary Fig. 10).
It is noteworthy that in most cases, our search time intervals meet the requirements of Ramp-fitting (contains two prominent change points) and BREAKFIT (contains one obvious change point), and no cropping is needed.However, the large dry event in the MAG record makes it hard to obtain the overall transitional trend during the SAHP2 termination and associated change points; in this scenario, we need to crop this event (Supplementary Fig. 8c).The cropped sections of the MAG raw data are replaced with constant values which equal the boundary values of the uncropped part of the record (Supplementary Fig. 8c).
In order to further check the robustness of the identified change points, we used another method by fitting "trends" to the whole timeseries which does not require setting the search time window in advance (namely "Trend-fitting" method) (Supplementary Code 1).This method does not calculate the uncertainties for change points.The Trend-fitting results are shown in Supplementary Fig. 9.In this study, we only display results that have passed the sensitivity test and are confirmed by "Trend-fitting" method, i.e., the change points picked by "Trend-fitting" method fall within the uncertainties determined via BREAKFIT and "Ramp-fitting" algorithms.In our case, the change points derived via different methods are comparable (Supplementary Fig. 9).The exceptions are the Antarctic Isotope Maximum 2 cooling transition in the WDC δ 18 O record and the Antarctic Isotope Maximum 2 transition in the EDML δ 18 O record (Supplementary Fig. 9), which are excluded in our analyses.

Combined uncertainty
In order to conduct a conservative estimate, we calculated the combined uncertainties for critical change points.Because Ramp-fitting and BREAKFIT merely calculate the change point uncertainty (Supplementary Table 2), we then quadratically combined the change point uncertainty and the age model uncertainty ( ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 + y 2 p ) to obtain the combined uncertainty (Supplementary Table 2).We contend that these uncertainties should be regarded as conservative estimates.

Additional test
Considering the age models could be asymmetrical and may influence the change points, we performed an additional test to check this situation.In this regard, we conducted the BREAKFIT and Ramp-fitting analyses for our speleothem records using the 2.5 th and 97.5 th percentile of age model ensemble of each record (Supplementary Table 4).We then calculated the average of the change points derived via these two age models and compared it with the change point derived via the median (50 th ) percentile of age models (Supplementary Table 4).The result shows that the asymmetry of age models does not significantly influence the change points of our speleothem records (Supplementary Table 4) as they are all constrained by several precise 230 Th dates, the conclusion of this study is therefore unaffected by this problem.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.© The Author(s) 2022 1

Fig. 1 |
Fig. 1 | Cherrapunji Cave speleothem δ 18 O records.a Speleothem δ 18 O records from Cherrapunji Cave (Cherrapunji-2 and Cherrapunji-2017-1, this study) were used to construct a composite record.Over their contemporary growth intervals, we exclusively use the Cherrapunji-2 record (Supplementary Note 1.1).Error bars show 230 Th dates with uncertainties (2σ) for each record (color coded).b Annual lamina thickness of speleothem Cherrapunji-2 and estimated annual growth rate of speleothem Cherrapunji-2017-1.The black box in (b) shows the period with an extremely slow growth rate in Cherrapunji-2017-1, and the corresponding δ 18 O data in (a) is depicted by dots.

2 Fig. 2 |
Fig. 2 | Speleothem δ 18 O records from Asian summer monsoon and South American summer monsoon domains and cave locations.a-e Speleothem δ 18 O records from Cherrapunji Cave (composite Cherrapunji record, this study) (a), Yongxing Cave (YX-51, black, this study) in comparison with the Cherrapunji record (light green) (b), Marota Cave (MAG, this study) (c), Paixão Cave (PX-07, this study) (d), and Botuverá Cave (BTV-4C, this study) (e).Note the inverted y-axis for (a) and (b).Error bars show 230 Th dates with uncertainties (2σ) for each record (color coded).Greenland Stadial (GS) and Greenland Interstadial (GI) 2 are shown at the top.It is noted that the timing of Greenland events are based on the GICC05 chronology (Supplementary Note 1.5), which is shifted by +320 years (see main text).AHP2 and SAHP2 durations are marked by double-sided arrows.The vertical bars from left to right during AHP2 depict stages I, II, and III.f Spatial pattern of June-July-August (JJA) minus December-January-February (DJF) precipitation amount-weighted δ 18 O (‰) in the Asian summer monsoon domain using the Isotope-incorporated Global Spectral Model version 2 (IsoGSM2 93 ) for 1979-2017.g is the same as in (f) but for DJF minus JJA in the South American summer monsoon domain.

Fig. 4 |
Fig. 4 | Comparison between Asian summer monsoon/South American summer monsoon and ice-core records.a Cherrapunji δ 18 O record (this study) from Cherrapunji Cave in the ISM domain, b Antarctica WDC ice-core CH 4 record 24 , c MAG δ 18 O record (this study) from Marota Cave in the SASM domain, d Antarctic ice-core composite Δε LAND proxy71 , e Greenland NEEM ice-core δ 18 O record (Supplementary Note 1.5) on the improved chronology (GICC05 timescale +320 years) (see main text), f, g are Antarctic ice-core δ 18 O records from WDC13 and EDML95 , respectively.Antarctic ice-core records are plotted on the WD2014 chronology63 , which is shifted by +400 years via a correlation based on the volcanic triplet spike (depicted by the orange diamond and the vertical bar at ~24.94 ky BP).It is noted that the EDML δ 18 O record has been tuned to the WD2014 chronology in a previous study16 .Red dots and associated error bars in (e) and (f) indicate the timing and combined uncertainties (2σ) of the change points in Greenland and Antarctic icecore δ18 O records (Supplementary Table2).The Brown error bar depicts the 14 C age66 and uncertainty (1σ) of the Oruanui eruption from the Taupo volcano.The vertical bar at ~25.72 ky BP marks the tephra layer of the Oruanui eruption recorded in the WDC ice core64 .The two vertical dashed lines depict the DO-2.1 and DO-2.2 peaks.Cave, ice-core and volcano locations are shown in Supplementary Fig.2.Vertical bars are the same as in Fig.2.AIM Antarctic isotope maximum95 , ppb parts per billion, ISM Indian summer monsoon, SASM South American summer monsoon.

Fig. 5 |
Fig.5| Global atmospheric teleconnections during stage III. a NGRIP deuteriumexcess record69 on the GICC05 chronology and shifted by +320 years.b-g speleothem δ18 O records from Yongxing Cave (YX-51, this study), Cherrapunji Cave (Cherrapunji, this study), Larga Cave (PA-LA-196 ), Marota Cave (MAG, this study), Paixão Cave (PX-07, this study), and Botuverá Cave (BTV-4C, this study), respectively.h Antarctica 5-core averaged d ln anomaly16 .i Antarctica WDC ice-core non-sea-salt Ca 2+ [nssCa 2+ ] record74 and overlain by the thick line (high-frequency variability (>1 cycle per 300 year) removed by a low-pass Butterworth filter, same as in ref.74).Antarctic ice-core records in (h) and (i) are on the WD2014 chronology and shifted by +400 years.The yellow dashed lines in (c-f) are generated via the BREAKFIT algorithm45 .Black dots and associated error bars in (c-f) depict the timing of breakpoints and combined uncertainties (2σ) (Supplementary Table2).Latitudes and longitudes for the various records in (a-g) and (i) are shown.Note the

Fig. 6 |
Fig. 6 | Speleothem records from South America and their phasing relationship with Asian summer monsoon records during AHP2 and SAHP2 terminations.a-d SASM speleothem δ 18 O records from Marota Cave (MAG, this study), Paixão Cave (PX-07, this study), Cueva del Diamante Cave (NAR-C, this study and ref. 8), and Botuverá Cave (BTV-4C, this study).e ASM speleothem δ 18 O record from Cherrapunji Cave (this study) (note the inverted y-axis).Cave locations are shown in Fig. 2f, g.Error bars in (c) depict 230 Th dates and 2σ errors of the NAR-C record.Red dots and associated error bars in (a) and (c-e) indicate timing and combined uncertainties (2σ) of the change points in speleothem δ 18 O records (Supplementary Table2), red lines in (a) and (c-e) show the ramps and break lines generated by

Table 1 |
Cave locations and climate domains