Onset and termination of Heinrich Stadial 4 and the underlying climate dynamics

Heinrich Stadial 4 during the last glacial period was marked by severe cooling at northern high latitudes along with the attendant changes in Asian Monsoon (Chinese Stadial 4) and South American Monsoon (South American Stadial 4). Here we present improved constraints on timings of Heinrich/Chinese/South American Stadial 4 onset and termination at sub-centennial precision based on speleothem records. We show that their initial onsets were essentially synchronous (40.20 ± 0.08 thousand years ago) and led the Antarctic warming by ~300 years. The Heinrich/Chinese Stadial 4 termination commenced at 38.34 ± 0.07 thousand years ago following a centennial-scale reduction in the Amazon River runoff and a poleward shift of the Southern Westerly wind belt. These two precursor events may have contributed to a reduced Amazon Plume Region and an enhanced Agulhas salt/heat leakage that led to an abrupt resumption of the Atlantic Meridional Ocean Circulation eventually triggering the Heinrich/Chinese Stadial 4 termination. Severe cooling at high northern latitudes, which marked the onset of Heinrich Stadial 4, was synchronous with changes in tropical monsoon systems and preceded Antarctic warming by around 300 years, according to a compilation of speleothem and ice core records.

T he last glacial period was characterized by a series of global-scale abrupt millennial-climate events known as Dansgaard-Oeschger (DO) cold stadials and warm interstadials 1 . Heinrich events contain layers of ice-rafted debris indicative of massive iceberg discharges into the North Atlantic [2][3][4] . The cold stadials containing the Heinrich events are also referred to as Heinrich Stadial (HS) 5,6 . HSs were marked by widespread cooler and drier conditions across large regions of North America, Eurasia, and North Africa [7][8][9][10] . A southward displacement of the Intertropical Convergence Zone (ITCZ) and a weakening of Northern Hemisphere (NH) Asian summer monsoon (namely Chinese Stadials (CSs)) 11 was associated, whereas many regions in the Southern Hemisphere (SH), especially the Amazon Basin, became wetter due to associated intensification of the summer monsoons (hereafter, South American Stadials (SASs)), which are nearly anti-correlated to CSs [12][13][14][15][16][17][18] . Meanwhile, Antarctica experienced warmer temperatures consistent with a reduction of the northward heat transport from SH to NH via for example, the Atlantic Meridional Ocean Circulation (AMOC) 19 . These scenarios manifest an important aspect of the global thermal seesaw [20][21][22][23] .
In the past three decades, six HSs and corresponding hydroclimate anomalies CSs/SASs (1)(2)(3)(4)(5)(6) have been recognized during the last glacial period with a~7 ky (thousand years) quasiperiodic recurrence 2,11 . Among them, HS4 contains the ice-rafted debris originated from the Laurentide Ice Sheet, similar to HS0 (the Younger Dryas) 24,25 . A large number of observational and modeling studies attribute HSs to the weakening of the AMOC in response to the release of freshwater by melting icebergs [26][27][28] . Consistently, recent studies have shown that the abrupt onset of Greenland warming (cooling) led the corresponding Antarctic cooling (warming) onset by~200 ± 100 y, implying a north-south propagation of abrupt climatic signals 29 . Subsequent ice-core studies have also found that the meridional migration of the SH Westerly Wind (SWW) belt and the ITCZ, two key components of the global atmospheric circulation implicated in the DOrelated climate changes, were in phase with Greenland HSs, suggesting a fast NH to SH directionality of global atmospheric teleconnection 21,23,30,31 . On the other hand, on the basis of precisely dated speleothem δ 18 O records, combined with polar icecore records, the evolution of SH or tropical Pacific temperature/ hydroclimate likely played an active role in DO variability/ Greenland climate dynamics 23 , as hypothesized previously in a number of studies but were based mostly on numerical models or low-resolution data with loosely constrained chronologies [32][33][34][35][36] . These developments, therefore, warrant further investigation of the phase relationship of millennial-scale events between different climate systems during the glacial period and their underlying mechanisms 23,37 .
Here, we report a set of new-generation speleothem δ 18 O records of CS4/SAS4 from Asian Monsoon (AM) and South American Monsoon (SAM) domains. We compare these data with a suite of previously published speleothem and ice-core records, and provide a synthesis of HS4 and associated CS4/SAS4 across different regimes. By focusing on the relative timing and structure of HS4, CS4, and SAS4 in each region, we aim to identify the dynamics associated with their onset and termination. Our results shed light on the underlying dynamics of the stadials, particularly with respect to their terminations.
The speleothem chronologies are based on extensive U-Th dating (~130 dates) by an improved technique 47 . The dating work was performed at the Isotope Laboratory of Xi'an Jiaotong University, using Multi-collector Inductively Coupled Plasma Mass Spectrometry (MC-ICP-MS, see "Methods" section). Typical age uncertainties (2σ) vary between 60 and 150 y for most key intervals (Supplementary Data 1). The age models of speleothem δ 18 O records were constructed using StalAge 48 ( Supplementary  Fig. 2). This algorithim estimates the uncertainty window of interpolated ages iteratively by using a Monte Carlo simulation (N = 300) (see "Methods" section). A total of~3550 oxygenisotope (δ 18

Results and discussion
Speleothem δ 18 O records. Recently, we used a correlation strategy that relies on using "breakpoints" instead of conventional "mid-points" for correlating the HS0 among the AM-AW (Asian Westerlies) and Greenland ice-core δ 18 O records 23 in the NH ( Supplementary Fig. 3a). In this study, we have adopted a similar strategy to determine the initial onset and termination of HS4/ CS4/SAS4 from precisely dated high-resolution speleothem δ 18 O records. This strategy employs identification of the breakpoints and their temporal uncertainties from the visual inspection of δ 18 O profiles (Method 1, see "Methods" section) together with more objective determinations of breakpoints by using the BREAKFIT algorithm 49 . Using this approach, we conducted BREAKFIT analysis using both individual records as well as their regional composites. To account for the dating uncertainties and time-series "noise" (i.e., differences in δ 18 O profiles from the same region), we repeated the analysis using the 2.5th, 50th, and 97.5th percentile of StalAge-derived age model ensembles for each record and their composite (Methods 2 and 3, see "Methods" section).
Based on the visual inspections alone, we identified the initial onset of CS4 in AM-AW records (i.e., the approximate δ 18 O minima) at 40.24 ± 0.12 ky BP (MWS-2), 40.21 ± 0.09 ky BP (ZJD171), 40.20 ± 0.06 ky BP (Wulu-30), and 40.14 ± 0.10 ky BP (So-1). The So-1 record is from Sofular cave in the AW region 50 ( Figs. 1 and 2). The error-weighted average age of CS4 initial onset from these records is 40.20 ± 0.05 (2σ) ky BP. The initial onset of SAS4 in the SAM records (i.e., the approximate δ 18  . For both CS4 and SAS4 records, the error-weighted average of their breakpoint ages is 40.20 ± 0.04 ky BP with a 2σ uncertainty of ±0.07 ky BP. We surmise that the error-weighted uncertainty and standard deviation reflect the age model uncertainty and the proxy time-series noise, therefore, by quadratically combining these errors, we derive the breakpoint age of 40.20 ± 0.08 ky BP as the anchor point of the initial onset of CS4/SAS4 for correlating hydroclimate records, especially those that have large absolute but small relative chronology uncertainties, such as ice-core records. This estimate is consistent with the timing and uncertainty of the corresponding "breakpoints" determined via Method 2 and Method 3 ( Fig. 1  In contrast to their synchronous onset, the timings of the initial termination of the monsoon anomalies between CS4 and SAS4 are apparently different. The timings (i.e., the δ 18 O maxima around the initial termination) are 38.  Tables 1-3 and Supplementary Data 1), although the large uncertainties preclude a definitive conclusion. Yet, more robust supportive evidence comes from the records in eastern Brazil. The LSF13 record shows a rather earlier SAS4 initial termination at 38.85 ± 0.11 in central-eastern Brazil. The records from Toca da Boa Vista/Toca da Barriguda caves, northeastern Brazil show the drying trend inferred by δ 18 O records started even much earlier (at~39.5 ky BP) and ended by a growth hiatus at~38.6 ky BP, implying an arid condition similar to pre-SAS4 and post-SAS4 (Figs. 1 and 3 and Supplementary Fig. 6). As such, the observed initial termination of SAS4 appears to lead the anchor point (38.34 ± 0.07 ky BP) by a few to several hundred years ( Fig. 3

and Supplementary Figs. 4-6).
Correlation of Greenland ice-core records to Asian monsoon records. The current chronologies of Greenland ice-core records largely rely on the annual band counting of the North Greenland Ice-core Project (NGRIP) ice-core (i.e., the GICC05 chronology). The absolute age uncertainty of GICC05 remains significant, ca. ±1500 (2σ) y at the HS4 interval 51,52 , which theoretically precludes correlations with other climate records at sub-centennial/ centennial precision. Although speleothem records provide much more precise age constraints on HS4, the transfer of speleothem chronologies to ice-core relies on a proper correlation strategy. Recent studies of HS0 and Greenland Stadial (GS) 20 (or "HS7") have established a synchronization strategy that utilizes breakpoints rather than mid-points, assuming a causal synchroneity between abrupt Greenland temperature change and the initial AM-AW response 23,53 . Crucially, the breakpoint strategy was validated via dating of speleothem HS0 records from the North  39 , Zhangjia (ZJD171) (this study) and Wulu (Wulu-30, this study) 40 caves in the AM domain. e-j Speleothem δ 18 O records from El-Condor (ELC-B) 12 , Paraiso (PAR27 and PAR15, this study) 14 , Lapa Sem Fim (LSF13, gray and orange curves: low and high-resolution data, this study) 41 , Toca da Boa Vista (TBV5 (red): this study; TBV40 (dark blue) 42 ), and Toca da Barriguda (TBR10-13: this study) 42,43 caves in the SAM domain. Cave locations are shown in Supplementary Fig. 1. Error bars show U-Th dates (2σ) for each record. The Hulu (MSL) chronology is based on the previous study 38 and the error bar depicts the typical uncertainty. Two vertical dashed lines depict the initial onset of CS4/SAS4 at 40.20 ± 0.08 ky BP and initial CS4 termination at 38.34 ± 0.07 ky BP, respectively. Vertical dark gray bars show the uncertainties. Blue diamonds show the initial onset and termination points determined directly based on the Method 1.

Fig. 2 Correlation between
Greenland ice-core and AM speleothem records. a, b Greenland NGRIP δ 18 O and Ca 2+ records on the GICC05 chronology 1 with a +150 y shift (red arrow) based on the "breakpoint" correlation with AM records 23 (also see text and Supplementary Fig. 3). c-f Speleothem δ 18 O records from Hulu, Mawmluh, Zhangjia, and Wulu caves in the AM domain ( Fig. 1). g, h So-1 δ 13 C and δ 18 O records from Sofular cave in Turkey, indicating the AW ecosystem and climate changes 50 . i CH 4 record from the Antarctic WDC ice-core on the WD2014 chronology 60 . The WD2014 chronology is shifted by −60 years (red arrow) based on the breakpoint correlation of CH 4 with AM records 23 . Error bars in g show the U-Th dates (2σ). The vertical dashed lines, bars and blue diamonds are as same as in Fig. 1.
Atlantic and AM-AW domains showing synchroneity within a few decades 23 .
On the basis of the HS0 synchroneity tested on decadal-scale between Greenland ice-core and AM-AW speleothem δ 18 O records 23 in the NH, we correlate the initial δ 18 O decrease in the NGRIP record around the HS4 onset to the 40.20 ± 0.08 ky BP anchor point from AM-AW and SAM records, and the initial δ 18 O jump at the HS4 termination in the NGRIP record to the anchor point at 38.34 ± 0.07 ky BP from AM-AW records by shifting the GICC05 chronology to the older side by 150 y, well within its quoted absolute age uncertainty of~1500 y (Fig. 2). This shift does not require a change of the relative age of the icecore chronology, which is consistent with the fact that the ice-core record has a much smaller uncertainty in its relative age (< 200 y across HS4) 54 . An apparent age difference of~250 ± 500 (2σ) y around HS4 between cave and ice-core was found between the NGRIP (GICC05 chronologies) and previous cave δ 18 O records (U-Th ages) 37 . The age difference reported here (150 ± 80 y) is comparable, but has much smaller uncertainty because of the substantial improvement in resolution and dating precision of the cave records used in this study. From the dynamic point of view, the synchroneity on a sub-centennial timescale between the initial onset/termination of HS4 in Greenland and CS4/SAS4 in mid-low latitude records hints towards a fast (decadal-scale) teleconnection via global atmospheric circulations (e.g., changes in the tropical Hadley circulation, meridional shifts of the ITCZ and the mid-latitude westerly wind belts in both hemispheres) 21,23,30,31,55,56 . In contrast, the longer (~500 y) and more gradual AM-AW onset (from~40.2 to 39.4 ky BP) and termination (from~38.3 to 38.0 ky BP) processes imply a slower (centennial-scale) oceanic reorganization in response to the abrupt North Atlantic climate change (e.g., changes in AMOC, subsequent South Ocean temperature and feedbacks) [20][21][22][23]29,57,58 (Figs. 2 and 3).
Records from the South American Monsoon domain. Over the tropical South Atlantic and in the SAM domain, SAS4, emerges as the most remarkable millennial-scale event, recording the highest amplitude of SST and monsoon precipitation variability during the last glacial 41,59 . Speleothem δ 18 O records from the SAM domain spanning SAS4 are characterized by gradual shifts (Fig. 3 and Supplementary Fig. 6). Comparably, the SAS4 onset was initiated at~40.19 ky BP, synchronous within uncertainty with CS4 ( Fig. 1 and Supplementary Fig. 4). The SAS4 termination excursion in the SAM domain manifests as a progressive weakening of the SAM or decreasing rainfall 12,14,41 and a northward shift of the ITCZ 42,43 . Notably, the initial termination started apparently earlier than that of HS4/CS4 by~250 y in Paraiso, Lapa Sem Fim, and El Condor records ( Fig. 3 and Supplementary Figs. 5 and 6). Consistently, the growth interval of speleothems from Toca da Boa Vista/Toca da Barriguda cave (~40.1-38.6 ky BP) in Northeast Brazil provides precise constraints on the timing of a pluvial period, associated with a southward shift of the ITCZ 42,43 ( Fig. 3 and Supplementary Fig. 6). The SAS4 structure in Northeast Brazil is characterized by an abrupt change (within 30 y) to a pluvial phase at~39.5 ky BP 41,42 , followed by a persistent increase in δ 18 O (reflecting drier conditions), and eventually reaching the value of the early stage of SAS4 (or the growth phase I between 40.1 and 39.6 ky BP) 42 at~38.6 ky BP. Intriguingly, growths of all stalagmites (e.g., TBV5, 14, 40 and 63 and TBR10-13) 42 ceased ( Fig. 3 and Supplementary Fig. 6) prior tõ 38.6 ky BP (corresponding to the start of abrupt δ 18 O increase in the LSF13 record), suggesting that the hydroclimate became drier than during growth phase I 42 , with an arid state comparable to pre-SAS4 and post-SAS4. The SAS4 pluvial phase, therefore, started approximately at the same time as counterparts in Greenland and AM-AW domains but terminated notably earlier.
Additionally, Northeast Brazil is a semi-arid area whose precipitation variability is very sensitive to tropical South Atlantic sea-surface temperature (SST), and the earlier rainfall termination thus implies an earlier cooling in the tropical Western Equatorial Atlantic 59 . Collectively, the above observations are critical because the aforementioned sites are located in areas sensitive for tracking rainfall in South America including the Amazon Basin and thus Amazon River runoff into the Atlantic Ocean, as well as the ITCZ and tropical southwestern Atlantic SST (Fig. 4 and Supplementary Fig. 1).
Phasing between Antarctica and Greenland. Previous comparisons between Antarctic and Greenland ice-core records demonstrate that abrupt Greenland warming (cooling) led the corresponding onset of Antarctic cooling (warming) bỹ 200 ± 100 y (2σ) 21,29 or 122 ± 24 y (2σ) 59 in average for DO events, including HS and associated Antarctic hydroclimate variations. In a previous study, we have confirmed the coincidence between NGRIP ice-core δ 18 O, AM-AW speleothem δ 18 O, and atmospheric CH 4 across HS0/CS0 to within ± 20-40 y 23 . Using the same correlation strategy, we identified two corresponding breakpoints (see "Methods" section) in the CH 4 record from the West Antarctic Ice Sheet Divide ice-core (WDC) on the WD2014 chronology ( Fig. 3 and Supplementary Fig. 7), and then matched them to the two anchor points at 40.20 ± 0.08 and 38.34 ± 0.07 ky BP by shifting the WD2014 chronology by 60 y to the younger side (Fig. 3). The high WDC accumulation rate provides highresolution gas CH 4 and ice δ 18 O data with small uncertainty on their age difference around HS4 period (±75 y) 60 . As such, the WDC ice δ 18 O chronology around HS4 is precise and accurate at the sub-centennial-scale, provided that the CH 4 correlation to AM-AW and the Greenland records is robust, as tested on the decadal-scale by the previous HS0/CS0 study 23 .
Two breakpoints are identified in the HS4 interval of the WDC δ 18 O record (see "Methods" section) at~39.87 and~38.66 ky BP ( Fig. 3 and Supplementary Fig. 7 and Supplementary Tables 1-3 and Supplementary Data 1). The breakpoint at 39.87 ky BP marks the warming onset associated with the HS4 onset in Antarctica consistent with other Antarctic records ( Supplementary Fig. 7), which lags the HS4/CS4 initial onset (at 40.20 ± 0.08 ky BP) and associated atmospheric CH 4 change by~300 y. This temporal relationship is similar to the interpolar phasing revealed previously for DO events, including HSs 21,29 . In contrast, the phasing relation surrounding the HS4 termination appears to be rather intricate. The breakpoint at~38.66 ky BP (highest δ 18 O Fig. 3 Comparison between polar ice-core and SAM speleothem records. a Same as a in Fig. 2. b-e Speleothem δ 18 O records from the SAM domain. b Speleothem records from EL Condor Cave (ELC-B, this study). c Speleothem records from Paraiso Cave. Blue and purple curves: PAR27 and PAR15 records (this study). Gray curve: PAR07 record 14 . d Record from Lapa Sem Fim Cave (LSF13, gray and orange curves: low and high-resolution data, this study). e Speleothem records from Toca da Boa Vista (TBV)/Toca da Barriguda (TBR) caves (red and dark blue: TBV5 (this study) and TBV40 42 ; purple: TBR10-13, this study). The red vertical arrow marks the speleothem growth hiatus. Error bars show the U-Th dates (2σ) of TBR10-13. f Antarctic deuterium-excess records (d ln , a proxy for vapor source condition) on the WD2014 chronology 60 . Jade curve: 5-core average d ln anomaly 21 , gray curve: WDC d ln record 30 . g, h CH 4 and δ 18 O records from the Antarctic WDC ice-core on WD2014 chronology 60 . WD2014 chronology was shifted by −60 years based on the breakpoint correlation of CH 4 with AM records 23 . The black error bar depicts the uncertainty (~±75 years) of the WDC gas-ice age difference 60 . i The δ 18 O record from the Antarctic EDML ice-core on WD2014 chronology) 60 with −60 year shift. j Antarctic temperature stack (ATS) record relative to the present day (ΔT) on the WC2015 chronology 63 . k Atmospheric CO 2 records (gray; orange on WD2014 chronology) 60,91 . l WDC non sea-salt Ca 2+ (nssCa 2+ ) record (WD2014 chronology) 92 . High-frequency variability (>1/300 y −1 ) removed with a low-pass Butterworth filter. Two vertical dashed lines are the same as in Fig. 2. Dashed blue lines and squares depict trends and breakpoints. The vertical gray arrow in f depicts the direction of poleward shift inferred by the d ln change. The vertical orange arrow indicates the main termination (drying) trend~250 y before 38.34 ky BP in SAM records, and the olive bar highlights changes/states of the other records during the same period.
value) or 38.64 ky BP (using the BREAKFIT algorithm) 49 in the WDC δ 18 O record occurs~300 y before the initial termination of HS4/CS4 at~38.34 ± 0.07 ky BP (Fig. 3). This phasing relationship can be tested by the WDC CH 4 record on the same chronology (WD2014), due to a small uncertainty (±75 y) in the WDC ice-gas age difference 60 and a strong correlation of CH 4 with AM-AW and Greenland δ 18 O records 23,44,60,61 , controlled by the extent of wetlands and thus CH 4 emissions 62 . Essentially, the CH 4 initial jump associated with the HS4/CS4 termination lags the cooling onset (or hydroclimate changeover) inferred by δ 18 O in the same ice-core by about 300 ± 75 y, suggesting an earlier change in the West Antarctic δ 18 O record than the GS4/ CS4 termination (Fig. 3

and Supplementary Figs. 3, 5, and 7 and Supplementary Tables 1-3 and Supplementary Data 1).
Notably, a spatial disparity emerges from the existing ice-core records from Antarctica. While the EDML, WDC, and DF records in Northwest Antarctica (related centrally to the Atlantic sector, especially EDML) 21 show a changeover (or cooling onset) hundreds of years earlier than the initial HS4/CS4 termination, TAL, and EDC records (related to the Pacific and the Indian Ocean sectors) 21 show a changeover a few hundred years after the initial HS4/CS4 termination (Supplementary Fig. 7). The spatially heterogeneous pattern is also shown in the 5-core average anomaly with two breakpoints, hundreds of years before and after the initial HS4/CS4 termination 21 , similar to the Antarctic temperature stack (ATS) record that has two breakpoints as well 63 (Supplementary Fig. 7). The breakpoint after the HS4/CS4 termination is consistent with the previous notion that the abrupt Greenland warmings led the corresponding onset of Antarctic coolings by~200 ± 100 y; however, the breakpoint before the Greenland termination, linked closely to the South Atlantic Ocean, occurred without any clear Greenland precursors (Supplementary Fig. 7). The underlying dynamics of the Antarctic heterogeneous pattern surrounding the HS4/CS4 termination remain unclear, possibly because of the special SWW pattern in the Atlantic sector 21 , or regional effects, such as wind-driven changes to sea ice, gyre circulation or Weddell Sea deep convection 64 . Additional modeling and empirical studies are critically needed to further explore the mechanism(s).
Climate dynamics. In the chronological framework established here for HS4/CS4/SAS4, the associated onset of the Antarctic warming was~300 y after the initial onset of HS4/CS4/SAS4 at 40.20 ± 0.08 ky BP, suggesting a similar interpolar average phasing as observed previously for the DO events of the last glacial period 21,29 . Moreover, our analysis also provides evidence that mid-latitude to low-latitude hydroclimatic changes occurred synchronously within decades with abrupt changes in the North Atlantic, suggesting a fast teleconnection via meridional migrations of the Hadley circulation, ITCZ, and Westerlies in both hemispheres 20,21,23,30,55,56 . Conversely, the fact that AM-AW, SAM, and tropical hydroclimates show considerably longer and more gradual changes compared with Greenland also highlights a southward teleconnection influence from the northern high latitudes to the mid-low latitudes, which is in line with the role of the AMOC and the subsequent impacts on the monsoon systems 23,56,57,65,66 , demonstrating one of the important aspects of the global thermal seesaw [20][21][22][23]29 .
The events surrounding the HS4/CS4 termination are rather distinctive. On the basis of the correlation 23 , the CS4 termination began at~38.34 ky BP, synchronous with the abrupt northern high-latitude temperature rise evident in Greenland ice-core records (Fig. 2). The more gradual change observed in the CS4 termination process suggests a southward propagation from northern high latitudes to mid-low latitudes, similar to the CS4 onset (Fig. 2). Furthermore, hydroclimate changes around the HS4/CS4 termination appear to be earlier in the SH, particularly in the SAM records as well as in Northwest Antarctic δ 18 O records ( Fig. 3 and Supplementary Figs. 6 and 7). This earlier changes in the SAM realm were likely linked to cooling in the Northwest Antarctica which is related closely to the Atlantic sector ( Supplementary Fig. 7), shedding light on the HS4 termination dynamics. The Amazon River, originating mainly from the Amazon Basin and fed by SAM rainfall, is by far the largest river by volume of water (∼6600 km 3 y −1 ) in the world, representing~17% of the global riverine discharge to the oceans. In the SAM domain, including the Amazon Basin, the SAS4 termination manifests as a persistent drying trend lasting~250 y, which occurs mostly prior to the initial HS4/CS4 termination at 38.34 ± 0.07 ky BP ( Fig. 3 and Supplementary Fig. 5), consistent with marine records 67 (Supplementary Fig. 1). This observation suggests that the Amazon River discharge into the Atlantic Ocean decreased progressively for a few hundred years before the termination of HS4/CS4. Given the temporal relationship, the reduced discharge could plausibly act as a precursor of the HS4/ SC4 termination. Lower freshwater input would increase the seasurface salinity in the so-called Amazon Plume Region (APR) [68][69][70] and subsequently across the North Atlantic via increased northward salt transport by a shallow AMOC 71 or an overturning circulation involving the North Atlantic Intermediate Water (NAIW) in the late HS [72][73][74][75] (Fig. 4). The increased AMOC could have induced positive feedbacks via enhanced poleward heat transport, creating a positive temperature anomaly in the North Atlantic Ocean with a spatial structure similar to the positive phase of the Atlantic multidecadal variability (AMV) mode 76 . The persistence of a positive AMV mode would favor the development of a "seesaw" pattern in the Atlantic SST across the thermal equator, resulting in warming in the north and cooling in the south, thus displacing the ITCZ northward, further weakening the SAM, and reducing Amazon discharge 76 . Taken together, the above processes could eventually contribute to the resumption of a deep/strong AMOC mode, or the HS4/CS4 termination.
In contrast to Antarctic ice-core δ 18 O records that show disparities in their changeover timing to cooling around the end of the HS4/CS4, Antarctic ice-core d ln (a logarithmic definition of deuterium-excess) records show a rather coherent decreasing trend of several hundred years towards the Greenland HS4 termination (Supplementary Fig. 7). This suggests a progressive southward (poleward) shift of the strengthened SWW 21,30 , or a positive Southern Annular Mode 77,78 that tends to propel a poleward shift of the subtropical front (STF) and a contraction of the SWW, thus weakening the winds at mid-latitudes (40°S to 50°S) [77][78][79] . As a result, the SWW off the southern coast of South Africa would shift poleward, which enhances the Agulhas Current and its salt/heat leakage to the Atlantic Ocean from the Indo-Pacific Ocean (the Agulhas Leakage). This process would further contribute to a gradual strengthening of the AMOC 35,80-84 (Fig. 4).
Modeling studies show that a precursory gradual shift in the AMOC system 85 may ultimately have reached a tipping point 86 , allowing for an abrupt resumption of the deep/strong AMOC mode 26,27 and leading to the HS termination in the North Atlantic realm. Recent modeling and empirical studies also suggest a gradual recovery of the AMOC during the final part of HS events 87,88 . It remains unclear, however, if and how these gradual changes in the North Atlantic realm are dynamically linked to SH changes observed here because the present paucity of records with sufficient resolution and dating precision/ accuracy precludes a precise comparison. We argue that the aforementioned SH changes were one of the important factors that contributed to the gradual recovery of the AMOC, which warrants further testing by model simulations.

Conclusions
The temporal constraints of CS4/SAS4 datasets presented here are characterized mostly by sub-centennial age precision, which together with a correlation strategy based on breakpoints, provides an improved chronological and interpretive framework to explain the dynamics of global ocean-atmosphere teleconnections associated with the signal propagations during HS4. The initial onset of this stadial occurred at~40.20 ± 0.08 ky BP in the North Atlantic, synchronous within sub-centennial uncertainty with CS4/SAS4, implying a fast atmospheric teleconnection. In contrast, the longer and more gradual CS4/SAS4 onset suggests an oceanic reorganization in response to abrupt North Atlantic changes. The initial Antarctic warming appears to lag the initial HS4/CS4/SAS4 onset by~300 y. These observations manifest a southward climatic teleconnection via both atmospheric and oceanic processes, from the northern high latitudes to the midlow latitudes and finally the southern high latitudes. The initial termination of HS4/CS4 occurred at~38.34 ± 0.07 ky BP and the dynamic relationship between HS4 and CS4 termination appears similar to their onset. Hundreds of years prior to the HS4/CS4 termination at~38.34 ± 0.07 ky BP, however, the Amazon River discharge decreased and the SWW and STF shifted poleward. These precursor events, through positive feedbacks, gradually strengthened the AMOC and eventually culminated in deep/ strong AMOC resumption and the abrupt HS4 termination in Greenland. In this regard, the initial trigger might reside in SH, suggesting a northward propagation from SH to the North Atlantic-AM-AW regions.

Methods
Paleoclimate records. We reconstructed or considerably improved ten speleothem δ 18 O records spanning HS4. The sample information and cave settings are as follows: MSL from Hulu Cave (32°30′N, 119°10′E) Fig. 1). Zhangjia cave is located near Guangyuan city, central-western China. The cave formed in limestones of the lower Triassic Feixianguan Formation. The columnarshaped calcite stalagmite ZJD171, 22 cm in height and 6 cm in width, was collected in the first cave chamber, at~800 m behind the cave entrance. According to instrumental records from Guangyuan meteorological station (~54 km southwest of the cave), the mean annual air temperature is 16.1°C and the mean annual precipitation is~950 mm,~75% of which occurs during summer (June-September) (1951-2019). Cave monitoring results from 2017 to 2020 show a stable cave temperature of~16.2°C and constant relative humidity of~100%.
U-Th dating method. Stalagmites were halved along the growth axis and polished. About 20-150 mg of powder was drilled near the central axis for each U-Th subsample. These subsamples were obtained by drilling the polished stalagmite section along the growth axis with a carbide dental burr. U-Th dating work was performed at the Isotope Laboratory, Xi'an Jiaotong University using multicollector inductively coupled plasma mass spectrometers (MC-ICP-MS) (Thermo-Finnigan Neptune-plus). We used standard chemistry procedures to separate U and Th for dating 4 . A triple-spike ( 229 Th-233 U-236 U) isotope dilution method was employed to correct for instrumental fractionation and determine U-Th isotopic ratios and concentrations. The instrumentation, standardization and half-lives are reported in Cheng et al. 47,89 . All U-Th isotopes were measured on a MasCom multiplier behind the retarding potential quadrupole in peak-jumping mode. We followed similar procedures of characterizing the multiplier as described in Cheng et al. 47 . Uncertainties in U-Th isotopic data 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. Corrected U-Th ages assume an initial 230 Th/ 232 Th atomic ratio of 4.4 ± 2.2 × 10 −6 , the values for a material at secular equilibrium with the bulk earth 232 Th/ 238 U value of 3.8. Most samples have high U/Th ratios and thus the corrections are negligible. A total of 130 U-Th dates were obtained from eight speleothem samples (Supplementary Data 1). All dates are in stratigraphic order within error ( Supplementary  Fig. 2).
Age models. We used StalAge software 48 to construct age models for MWS-2, ZJD171, Wulu-30, ELC-B, PAR27, PAR15, and LSF13 records reported in this study ( Supplementary Fig. 2). StalAge is particularly suited for speleothems creating objective age models based on two assumptions: (1) the age model is monotonic, and a straight line is fitted through all data or through as many data points as possible within error bars 48 . StalAge produces 300 realizations of age models by the Monte-Carlo simulation to account for the 95% confidence limits 48 . The 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 the ages in each age models increases monotonically within dating uncertainty. in the boundary of sample MWS-2, ELC-B due to limited dating controls. Nonetheless, the key breakpoints (the initial HS4/CS4/SAS4 onset and termination) in detection periods are robust, which is away from the boundary areas and well constrained by several U-Th dates ( Supplementary Fig. 2). The age model of MSL was constructed by Oxcal 90 , consistent with previous publication 37 . Age models of TBV5 and TBR10-13 were constructed by linear interpolation.
Oxygen isotope analysis. For each oxygen isotope measurement, ∼200 μg of powder was drilled from the sample. A total of~3550 oxygen isotope (δ 18  Breakpoint determination. We used three methods to identify the timings and age uncertainty of initial onset and termination of CS4/SAS4 from the AM-AW and SAM domains. Method 1: Through visual inspection, we picked the minima (maxima) in different records between 40.0 and 40.4 ky BP indicating the beginning of persistent positive (negative) excursions in the AM-AW (SAM) records marking the initial onset of CS4/SAS4. Similarly, we chose the maxima (minima) in the records between 39.0 to 38.2 ky BP indicating the beginning of persistent negative (positive) excursions in the AM-AW (SAM) records marking the initial termination of CS4/SAS4. The age uncertainties of these points were calculated by combing 2σ uncertainty and error-weighted mean uncertainty quadratically, which reflect the age model uncertainty and partially the proxy time-series "noise". Method 2: We used BREAKFIT 49 algorithm to identify the onset and termination of HS4/CS4/SAS4 in various speleothem and ice-core records. The algorithm employs a continuous function, consisting of two linear parts that are joined at the breakpoint. The break model is fitted to data using a weighted least-squares method with a brute-force search for the breakpoint. While BREAKFIT provides an objective estimate for the change points in a given dataset, the choice of "fit interval" is subjective and can influence the results. The main criteria to choose analytical time intervals are as follows: (1) the interval contains one breakpoint, and (2) use the same time intervals for records from the same region if possible. Though BREAKFIT algorithm can give statistical uncertainties in the timing of breakpoints using 2000 block bootstrap simulations, which preserved the distribution and serial dependence of the data over the length of a block (Supplementary Tables 1 and 2), it does not take into account the dating uncertainty associated with the record. Additionaly, the differences in δ 18 O records from each domain, referred here as the "noise", also contribute to the uncertainty. To account for these additional sources of uncertainties, we conducted the BREAKFIT analysis for our speleothem records using the mean (50th) as well as the 2.5th and 97.5th percentile of age model ensemble of each record. The final breakpoint and 2σ uncertainty of each record was determined by the mean and standard deviation of three BREAKFIT results. The uncertainty mainly encompasses the age model uncertainty. Method 3: To further assess both the dating uncertainty and the noise of records within each domain, we z-scored and composited the record from the same climate domain as follows. We combined the ZJD171, MWS-2, and Wulu30 records as the representative of the AM-AW domain, and the PAR15, PAR27, LSF13, and ELC-B records as the representative of the SAM domain. For each domain, we again conducted the BREAKFIT analysis on the composite record using the 2.5th, 5th, and 97.5th percentile of age model ensemble. The final breakpoint and σ uncertainty were determined by the mean and standard deviation of the three BREAKFIT results. The reported uncertainties incorporate uncertainties from the age model and proxy time-series noise. The selected time intervals for BREAKFIT analyses are listed in Supplementary Tables 1-3. Results and  additional discussion are presented in Supplementary Figs. 4 and 5 and Supplementary Tables 1-3. Generally, breakpoints determined via BREAKFIT algorithms agree well with visual inspections.
Correlation strategy. The direct comparison between AM-AW and North Atlantic records supports the correlation strategy via matching "breakpoints" 23,53 rather than their "mid-points" 13,29 (Supplementary Fig. 3).

Data availability
The new data used in this study are reported in Supplementary Data 1, and will be available at https://www.ncdc.noaa.gov/paleo/study/34255

Code availability
A version of the breakfit software (without timescale simulation) 49