Tipping point-induced abrupt shifts in East Asian hydroclimate since the Last Glacial Maximum

Domain: earth_science
Published: 2025-01-08
Editorial outcome: minor_revision

Abstract
MultipletippingpointsintheEarthsystemcouldbetriggeredwhenglobalwarmingexceedsspecificthresholds.However,thedegreeoftheirimpactontheEastAsianhydroclimateremainsuncertainduetothelackofquantitativerainfallrecords.HerewepresentanensemblereconstructionofEastAsiansummermonsoon(EASM)rainfallsincetheLastGlacialMaximum(LGM)usingninestatisticalandmachinelearningmethodsbasedonmulti-proxyrecordsfromamaarlakeinsouthernChina.OurresultsdefinefivetippingpointsintheEASMrainfallsincetheLGM,whicharecharacterizedbyabruptandirreversibleregimeshiftswithamedianamplitudeof387±73mm(24±5%).Combinedwithmulti-modelsimulationsandexistingrecords,weattributethesetippingpointstocascadesofabruptshiftsintheAtlanticmeridionaloverturningcirculation(AMOC)andSaharanvegetation.OurfindingsunderscorethenonlinearbehavioroftheEASManditscouplingwithothertippingelements.

Introduction
Thetippingpointisacriticalthresholdbeyondwhichasmallperturbationcanqualitativelytransformthestateorevolutionofasystem1.Aclimatetippingpointhappenswhenaminorvariationinforcingtriggersaprofoundnonlinearresponseinthedynamicprocesseswithintheclimatesystem,leadingtoanabruptandirreversibleregimeshift2.Therehasbeenincreasingawarenessandconcernthatglobalwarmingcouldtrigger tipping points in the climate system, introducing far-reaching influence on natural ecosystems and human societies3, which call for urgent political and economic actions to reduce greenhouse gas emissions and keep atmospheric CO2 concentration below a safe level to mitigate the risk4. Several components of the Earth system have been identified as tipping elements that may reach or pass a tipping point when global warming exceeds a critical threshold such as 1.5°C above preindustrial temperature1,3. The Earth system is a complex dynamical system with different components highly interconnected and interdependent, thus tipping in one element can lead to cascading changes in another5–7, and even could trigger a planetary-scale tipping point8,9.

The global hydrological cycle, particularly the global monsoon system, has been suggested to be vulnerable to climate tipping points and may experience abrupt shifts in the future10. However, the direction and amplitude of abrupt shifts are unclear, mainly because the available observational records are too short to provide useful information11, and rainfall response is usually more difficult to constrain than temperature12. Paleoclimate records can provide the only long-term context for evaluating the impact of climate tipping points11. However, the precision, resolution, and spatial coverage of paleoclimate records are insufficient and need to be improved11. In the EASM domain, our understanding of paleomonsoon evolution mainly depends on the precisely-dated speleothem δ18O records13. However, the interpretation of speleothem δ18O records in East Asia remains debated, excluding it as a local rainfall proxy13–16. For instance, a speleothem trace-element record from Haozhu Cave in central-eastern China suggested wetter conditions during North Atlantic cold stadials when speleothem δ18O records indicated a weaker summer monsoon15. The situation is further complicated by inconsistent results of other proxy records in the same region17–21. Besides, the nonlinear behavior of the EASM is understudied, making it unclear whether changes in the EASM rainfall are gradual or abrupt. The potential for abrupt shifts in the EASM rainfall under excessive global warming were largely overlooked in both assessment reports and academic literature1,3,22,23.

Absent from these debates is an independent quantitative reconstruction of past rainfall changes in the EASM realm, particularly in southern China. Southern China is the core region of EASM influence where annual rainfall exceeds 800 mm and rainfall variability is the largest (Fig. 1).

In this study, we aim to provide a robust reconstruction of the EASM rainfall in southern China and use it to detect tipping points in the EASM rainfall since the LGM and assess the cascading impacts of climate tipping points on the EASM rainfall. The first goal is achieved by applying nine statistical and machine learning methods to retrieve rainfall signals from multi-proxy records of a maar lake in southern China. Our analysis indicates that the rainfall variability at the lake location is representative of rainfall variability in southern China (Supplementary Fig. 1). The second goal is realized by comparing a large number of proxy records with multi-model simulations. We focus on two well-known tipping elements, the AMOC and Saharan vegetation cover1,3. The AMOC has a strong impact on global climate by conveying surface warm water northward and deep cold water southward, redistributing the solar energy received by the Earth24. Through its impact on the interhemispheric thermal gradient, the AMOC can alter the position and intensity of low-latitude rainfall systems, including the intertropical convergence zone (ITCZ) and the EASM25,26. In situ observations indicate that the AMOC has weakened during recent decades27, and the sea surface temperature (SST) fingerprint suggests that this weakening trend extends back to the mid-20th century28. Climate models further project the AMOC to weaken under future anthropogenic warming29. Recent studies even warn that the AMOC may collapse in the coming decades30,31. However, the climatic impact of AMOC collapse in the context of anthropogenic warming remains unclear. The last deglaciation is characterized by a large increase in the atmospheric CO2 concentration and abrupt shifts in the AMOC32, providing an ideal analog for examining its future influence. The Sahara is currently the largest hot desert and source of dust emission on Earth, which plays a crucial role in regulating global climate and biophysical feedbacks33. Paleoenvironmental records indicate the greening of the Sahara during the African Humid Period (AHP) when the African summer monsoon intensified and rainfall increased in response to higher summer insolation34. Conversely, the termination of the AHP and collapse of Saharan vegetation during the mid-to late Holocene are abrupt events occurring faster than orbital forcing10,34,35, which is a paradigm for going through climate tipping points1. Recent studies suggest that the changes in Saharan vegetation cover and dust emissions can not only influence local climate in northern Africa36 but also the distant climates ranging from the tropics37,38 to the Arctic39. The Asian summer monsoon could be affected by Saharan vegetation changes through large-scale atmospheric circulation such as the upper-level Rossby wave train and the Pacific Walker circulation40,41. However, the timing, amplitude, and structure of EASM rainfall changes in response to the end of the Green Sahara and AHP are not well defined due to the lack of high-resolution quantitative rainfall records. Given that the Sahara/Sahel vegetation and the West African monsoon have been identified as a tipping element that may experience abrupt shifts under future anthropogenic warming1,3, investigating the effect of Saharan vegetation change during the past will have practical implications for better understanding of vegetation feedback on future climate.

Results
EnsemblereconstructionoftheEASMrainfall
Quantitative paleoclimate reconstructions are closely related to climate proxy and reconstruction methodology42,43. We, therefore, employ a multi-proxy, multi-method probabilistic framework to generate a large ensemble for robust reconstruction of the EASM rainfall. The reconstruction target is anomalies in the annual rainfall (Pann) relative to 1961–1990 CE. The period 1953–2005 CE was chosen for calibration because it has the most complete climate data and overlaps with the proxy data. The period before 1953 CE was not used for calibration due to significant uncertainty in the climate data. However, gridded rainfall data from the Global Precipitation Climatology Centre (GPCC)44 was utilized as an additional independent validation of our reconstructions. Because annual rainfall in East Asia is dominated by summer monsoon rainfall and their changes are highly correlated with similar amplitudes (Supplementary Fig. 2), we only reconstruct the annual rainfall. To preserve the full amplitude of rainfall variability and avoid spurious signals, both rainfall observations and proxy records were decadally smoothed with a cubic-smoothing spline45 to have roughly the same frequency (Supplementary Fig. 3), following the suggestion of previous studies46–48. Smoothing can also enhance the signal-to-noise ratio and stabilize the relationship between climate and proxies, reducing the impact of age model uncertainty (Supplementary Fig. 4). From a sedimentological perspective, smoothing is necessary because climate signal is naturally smoothed during the deposition of lake sediments, and lake records are unable to capture the year-to-year climate variability as precisely as tree-ring records47,48.

In November 2016, a sediment core (QTY-2-1) was recovered from the center of the Qingtongyang Maar Lake (110°10′E, 20°34′N, 131 m above sea level) in the Leizhou Peninsula, southern China. This lake is a closed-basin crater lake with a very small catchment and lack of river input, making it a natural gauge for past rainfall amount change. The core was split and subject to high-resolution color extraction and X-ray fluorescence (XRF) element scanning (Methods). The age model of the core is based on linear interpolation and extrapolation of the median values of sixteen AMS (accelerator mass spectrometry) 14C dates with two outliers excluded, which have been calibrated to calendar years before present (BP), where present refers to 1950 CE (Methods, Supplementary Table 1, and Supplementary Fig. 5). The age model over the past 150 years is constructed by the piecewise constant rate of supply (CRS) model49 of 210Pb dates constrained by the 137Cs peak at 1963 CE and two post-bomb 14C dates at 1955 CE and 2014 CE (Supplementary Fig. 6). The full ranges of these post-bomb 14C dates are only 2 to 3 years (Supplementary Table 1), ensuring a precise age model and allowing for meaningful calibration with rainfall observations on the decadal timescale. Even when assuming a larger age uncertainty by treating the full age range as one standard error and applying Bayesian age modeling to assess its impact, the rainfall-proxy relationship remains robust, regardless of changes to the age model (Supplementary Fig. 4).

Four climate-proxy records with nearly annual to decadal time-resolution were used as predictors for quantitative rainfall reconstructions, including the color b* and the concentrations of iron (Fe), manganese (Mn), and calcium (Ca) (Fig. 2). Color b* is a proxy for the concentration of goethite and hematite in lake sediment50. Fe and Mn concentrations are sensitive proxies for redox conditions in the lake at the time of sediment deposition51. Ca concentration indicates the content of carbonate in lake sediment52. Therefore, they are excellent indicators of past rainfall changes (Supplementary Discussion). Both correlation and partial-correlation analyses indicate a significant correlation with rainfall (P < 0.02) but a poor correlation with temperature during the calibration period 1953–2005 CE when accounting for the degrees of freedom (Supplementary Table 2), demonstrating that rainfall is the dominant controlling factor of changes in these proxy records and thus can be faithfully reconstructed. We used the expressed population signal (EPS)53 and mean inter-series correlation (Rbar)53 to evaluate the common variance of these proxy records over time with a sliding window of 3000 years. The EPS is a measure of how well a finite number of proxy records represents an infinite population chronology53. The Rbar indicates the average correlation between all possible pairs of proxy records53. The persistent positive values of sliding EPS and Rbar quantitatively indicate that these proxy records respond to a common climate signal and the proxy-climate relationships remain stable over time (Fig. 2e). Our results indicate that the four proxy records from the Qingtongyang Maar Lake exhibit synchronous changes over the past 22,000 years, especially during the last deglaciation and mid-to-late Holocene when abrupt climate changes occurred, which is supported by higher EPS and Rbar values (Fig. 2). All proxy records show low values during the Bølling–Allerød (BA) interstadial and the early to middle Holocene, which qualitatively indicates anaerobic conditions in the lake and increased freshwater input (Supplementary Discussion).

Our rainfall reconstruction ensemble consists of 14,850 members that account for uncertainties in the reconstruction method, proxy measurement, proxy selection, calibration period, weighting scheme, and regression residual. We use nine statistical and machine learning methods to reconstruct rainfall: composite plus scaling (CPS) and principal component regression (PCR) and their variants, i.e., optimal information extraction (OIE) and point-by-point regression (PPR), partial least squares (PLS) regression and two regularization methods, i.e., ridge regression (RIG) and elastic net (ELN), as well as two machine learning methods, i.e., artificial neural network (ANN) and convolutional neural network (CNN) (Methods). All methods use the same proxy matrix as input and each generates 1650 members of reconstruction. Cross-validations indicate that these methods have excellent performance in reconstructing rainfall variability, with the regularization and machine learning methods showing superior performance and a narrower range of uncertainties (Supplementary Table 3 and Supplementary Fig. 7). The rainfalls reconstructed by these methods closely track station-based rainfall observations during the period 1953–2005 CE and the GPCC rainfall during the period 1891–2005 CE (Supplementary Fig. 8).

Abrupt shifts in the EASM rainfall since the LGM
The rainfall reconstructions obtained by the nine methods exhibit similar oscillations on the decadal to millennial timescales over the past 22,000 years (Fig. 3 and Supplementary Figs. 8, 9), indicating that our reconstruction is robust to the methodology used. Based on the framework of ensemble reconstruction, we consider as many uncertainties as possible when quantitatively reconstructing the EASM rainfall (Methods). The influence of age model uncertainty has also been considered. The age model uncertainty mainly affects the weights of proxy records and the timings of abrupt climate changes (Supplementary Fig. 4). Therefore, we addressed this uncertainty by randomly adjusting proxy weights during reconstruction and applying Bayesian age modeling to ascertain the timings of abrupt shifts (Methods). Even after accounting for these uncertainties, the amplitudes of rainfall changes in our reconstructions during the past 22,000 years are still larger than the total uncertainties (Fig. 3 and Supplementary Figs. 8, 9), demonstrating that the reconstructed rainfall changes are significantly better than the stochastic process.

Consistent with the physical and chemical interpretations of the proxy records (Fig. 2 and Supplementary Discussion), all the quantitative reconstructions indicate that the EASM rainfall was most abundant during the BA interstadial and the early-middle Holocene, with rainfall amounts being ~440 mm (28%) higher than the modern value of 1592 mm in the period 1961–1990 CE. The rainfall was substantially lower during the LGM, Heinrich Stadial 1 (HS1), Younger Dryas (YD), and late Holocene, with the lowest rainfall occurring during the HS1. The general trend of the EASM rainfall evolution revealed by our reconstructions is consistent with the pattern of summer insolation change in the Northern Hemisphere, which peaked during the early Holocene (Fig. 3). However, unlike the gradual change in summer insolation, our rainfall reconstructions strongly suggest that the EASM rainfall was characterized by contrasting stable states and abrupt transitions between these states. Based on the regime shift analysis54,55 (Supplementary Fig. 10), the evolution of the EASM rainfall over the past 22,000 years can be divided into six major climate states (or periods) that have different characteristics and internal dynamics, aligning with the LGM, HS1, BA, YD, early-middle Holocene, and late Holocene. The rainfall anomalies relative to the mean of 1592 mm in the period 1961–1990 CE during these periods are 141±97 mm, −1±98 mm, 447±86 mm, 53±89 mm, 440±75 mm, and 102±105 mm, respectively. Between two adjacent climate states are abrupt and irreversible shifts in rainfall amounts. Here, “abrupt” indicates that the change in a system is substantially faster than the rate of its past change, while “irreversible” indicates that the system does not recover to its previous state immediately after shift, but remains in a new state for a long time (>1000 years in this study) until its next shift10. Based on the Bayesian age modeling of our 14C dates from the Qingtongyang Maar Lake (Methods), the timings of these abrupt shifts are estimated to be centered at 18860±310, 14640±400, 12870±290, 11360±170, and 4370±60 yr BP (Supplementary Table 4). The transition from the HS1 to the BA is characterized by an abrupt increase in rainfall by 448 mm, while the transitions into and out of the YD are characterized by abrupt changes in rainfall by 394 mm and 387 mm, respectively. The transition from middle to late Holocene is characterized by a decrease in rainfall by 338 mm. The amplitude of rainfall decrease from the LGM to the HS1 is relatively small, with an amplitude of 142 mm, possibly because the rainfall amount during the LGM was low. Nonetheless, all rainfall changes are statistically significant at the 95% confidence level based on the two-tailed two-sample t-test (Supplementary Fig. 10). These abrupt shifts can be identified by another objective method, the Bayesian change-point analysis56, which further supports our finding based on the regime shift analysis (Supplementary Fig. 10). These two independent methods are consistent in revealing the timings and structures of abrupt shifts in the EASM rainfall in most cases, but the Bayesian method can provide more details about the abrupt shift during the mid-to-late Holocene transition. Based on the Bayesian change-point analysis, the mid-to-late Holocene transition is characterized by an abrupt drop of rainfall with an amplitude of ~140 mm at 4370±60 yr BP and a subsequent less abrupt but persistent decrease in rainfall with an amplitude of 416 mm until 3290±70 yr BP (Supplementary Fig. 10 and Supplementary Table 4). Therefore, our reconstruction reveals a two-step abrupt shift in the EASM rainfall during the mid-to-late Holocene transition, which has implications for a better understanding of the forcing mechanisms and cultural consequence of the well-known 4.2 ka event41,57.

To test the representativeness of our rainfall reconstruction, we compare our reconstruction with other proxy records from the same region, including the Huguangyan and Tianyang maar lakes (Supplementary Fig. 11). Although the proxy records from Huguangyan and Tianyang have much lower temporal resolution, millennial-scale events such as the HS1 and YD are clearly identified, supporting that the abrupt shifts in our rainfall reconstruction are regional phenomena rather than local noises. We also make a comprehensive compilation of paleoclimate records from the whole EASM domain (Methods), which consistently indicate that abrupt shifts in the rainfall regime are a continental-scale feature and intrinsic mode of the EASM rainfall variability since the LGM (Supplementary Tables 5–9 and Supplementary Figs. 12, 13). To help better understand the climatic significance of speleothem δ18O records in East Asia, we compare our rainfall reconstruction with the speleothem δ18O records from Dongge and Hulu caves in southern China (Supplementary Fig. 11d). The strong similarity between them during the last deglaciation indicates that the speleothem δ18O records in East Asia might be indeed related to hydroclimate changes during this period. After accounting for the effect of global ice volume change, the speleothem δ18O values during the LGM are slightly lighter than those during the late Holocene. This is consistent with our reconstruction indicating higher rainfall amounts during the LGM than in the late Holocene. However, the changes in speleothem δ18O records during the Holocene are much more gradual, indicating decoupled changes between speleothem δ18O records and the East Asian hydroclimate. Previously, the speleothem δ18O records in East Asia were interpreted as monsoon wind intensity, which is positively (negatively) correlated with monsoon rainfall in northern (southern) China58. However, the decoupled changes between speleothem δ18O records and our rainfall reconstruction, as well as other hydroclimate records in East Asia, indicate that speleothem δ18O records are not a direct signal of monsoon rainfall.

AMOC forcing of EASM rainfall during the last deglaciation
The EASM rainfall change can be driven by both dynamic and thermodynamic components59, with the former related to monsoon circulation and the latter related to atmospheric water vapor. The global mean temperature during the LGM was 6.1°C (95% confidence interval: 5.7–6.5°C) cooler than that during the late Holocene60, and the associated thermodynamic effect was expected to reduce atmospheric water vapor content and thus lead to less rainfall in East Asia. However, our reconstruction indicates that the rainfall amount in southern China during the LGM was similar to or even higher than that during the late Holocene (Fig. 3). This finding suggests that the global mean temperature change may not be the dominant control of the EASM rainfall anomaly during the LGM, at least in the core monsoon region. Moreover, the reconstructed rainfall changes during the last deglaciation did not follow the gradual increasing trend in the global mean temperature61,62 and atmospheric CO2 concentration63 (Fig. 4). Therefore, the observed variability in EASM rainfall is likely caused by other mechanisms, rather than a simple response to the thermodynamic effect of changes in the global mean temperature.

The only climatic forcing that shared similarity between the LGM and late Holocene is the local summer insolation at 20°N64 (Fig. 3). We, therefore, attribute the similar rainfall amounts to the dynamic effect of nearly identical summer insolation, which could compensate for the thermodynamic effect of global mean temperature through changes in the land-sea thermal contrast and latitudinal temperature gradient. However, insolation forcing alone cannot explain the EASM rainfall changes during the last deglaciation when the EASM rainfall was punctuated by two millennial-scale drought events (i.e., HS1 and YD) at times of high summer insolation (Fig. 3). In contrast, the oscillation of deglacial EASM rainfall was in parallel with the AMOC strength65 and the interhemispheric temperature gradient61 (Fig. 4), indicating a dominant role of the AMOC forcing. The AMOC has experienced tipping points for abrupt transition between strong states and weak states during the last deglaciation65, which could have strong impacts on the EASM rainfall15,66. For instance, remarkable weakening of the AMOC during the HS1 and YD65 reduced northward heat transport and generated a bipolar temperature seesaw pattern, with cooling in the Greenland67 and Northern Hemisphere and warming in the Antarctic68 and Southern Hemisphere, which diminished the interhemispheric temperature gradient61 and forced the ITCZ to move southward, reducing rainfall in East Asia (Fig. 4). Comparison of our rainfall reconstruction with hydroclimate records along a north-south transect across the East Asian-Australian monsoon system reveals an opposite rainfall pattern between the two hemispheres during the last deglaciation (Supplementary Fig. 14), supporting a southward migration of the ITCZ mean position during the HS1 and YD in response to the weakening of the AMOC during the HS1 and YD26,69. The cooling of the North Atlantic during the HS1 and YD can also influence the EASM rainfall via changes in the intensity and position of the westerlies15,66,70. However, the spatial pattern of EASM rainfall response to abrupt shifts in the AMOC during the last deglaciation remains ambiguous due to the lack of a systematic analysis of both proxy records and climate simulations. For instance, a recent study based on a trace-element record from Haozhu Cave and a hosing experiment under preindustrial conditions suggests that central-eastern China was wet (dry) during the times of a weak (strong) AMOC and the rainfall anomaly in East Asia exhibited a tripole mode15. The Haozhu record further indicates that the climate during the LGM was wetter than that of the early Holocene15. However, a wetter central-eastern China during North Atlantic cold events and LGM is challenged by other proxy records from the same region17–21, including two trace-element records from Heshang Cave17 and Hulu Cave21. Moreover, another hosing experiment under the LGM boundary conditions indicates rainfall decrease rather than increase in central-eastern China in response to the weakening of the AMOC66, which highlights the need for multi-model comparison. Previous studies mainly focused on the response of EASM rainfall to the weakening of the AMOC but overlooked the impact of AMOC strengthening. However, the responses of climate to the weakening and strengthening of the AMOC may not be a reversed process. In particular, abrupt shifts in the AMOC during the last deglaciation were accompanied by variations in other climatic forcings such as solar insolation, ice sheet, and greenhouse gas concentrations. Therefore, it is necessary to analyze the individual characteristics of EASM rainfall response during both the transitions into and out of HS1 and YD, which has not yet been done. Here we collect a large number of hydroclimate records in East Asia based on rigorous criteria (Methods, Supplementary Tables 5–8) and compare them with a state-of-the-science isotope-enabled transient climate experiment (iTRACE)16 and an earlier version without isotope module (TRACE-21ka)71 (Methods). Although continuous records with good age control and high temporal resolution are scarce for the last deglaciation, the general trend of hydroclimate change can be roughly inferred from these proxy fragments.

From LGM to HS1, the majority of proxy records, including the trace-element record from Haozhu Cave15, indicate a drying climate over East Asia, except for two proxy records from southeastern China72,73, which indicate a wetting climate (Fig. 5a–c and Supplementary Table 5). In contrast to the proxy records, the iTRACE simulation indicates that both annual and summer rainfalls increased during the HS1 relative to the LGM over Mongolia and eastern China (Fig. 5a, b). In the TRACE-21ka simulation, the annual and summer rainfalls increased in northern China and northeastern China but decreased in southern China (Supplementary Fig. 15a, b). Since it has been suggested that autumn rainfall could contribute to the changes in East Asian hydroclimate74, we also analyze the response to autumn rainfall. Indeed, the magnitude of autumn rainfall variability is comparable to that of summer rainfall and even larger (Fig. 5c). The decrease in autumn rainfall in northeastern China can partly explain the discrepancy between proxy records and model simulations in this region, but not for that in northern China and southern China. As revealed by sensitivity experiments, the summer rainfall variability in northern China is driven primarily by summer insolation and secondarily by meltwater, while summer rainfall variability in southern China is driven primarily by meltwater and secondarily by summer insolation16. In the iTRACE simulation, the increasing summer insolation during the HS1 tends to increase summer rainfall in northern China and northeastern China but decrease summer rainfall in southern China, while the weakened AMOC tends to decrease summer rainfall in northern China and northeastern China but increase summer rainfall in southern China16. As a consequence, summer rainfall in the whole of eastern China increased during the HS1 in the iTRACE simulation (Fig. 5b). However, such a response of EASM rainfall in the iTRACE simulation may not be realistic as it is inconsistent with the currently available proxy records (Supplementary Table 5) and there are large uncertainties in simulating rainfall response in East Asia to meltwater forcing75. We argue that the effect of AMOC slowdown on the global summer climate during HS1 might have been underestimated. For instance, summer cooling during the HS1 is limited to the North Atlantic in the iTRACE simulation, and the Eurasian continent experiences summer warming (Fig. 6b). The warming of the Eurasian continent may have been responsible for the simulated rainfall increase in East Asia during the HS1. However, a recent reconstruction of summer temperature in northeastern China based on the distribution of branched glycerol dialkyl glycerol tetraethers (brGDGTs) in the sediment core from the Arxan Lake indicates substantial summer cooling during the HS176.

From HS1 to BA, all proxy records, except for the trace-element record from Haozhu Cave15, indicate a wetting climate over the whole of East Asia (Fig. 5d–f and Supplementary Table 6). The iTRACE simulation shows that both annual and summer rainfalls, as well as autumn rainfall, increased in northeastern China and northern China but decreased in southern China during the BA period (Fig. 5d–f). Similar to the iTRACE simulation, annual and autumn rainfalls in southern China decreased during the BA in the TRACE-21ka simulation, but summer rainfall increased in both northern and southern China (Supplementary Fig. 15d–f). The global temperature during the transition from HS1 to BA is characterized by strong warming in the Northern Hemisphere with small cooling in the Southern Hemisphere and weak warming in Antarctica (Fig. 6c, d), while the transition from LGM to HS1 is characterized by strong cooling in the North Atlantic and strong warming in the Southern Hemisphere (Fig. 6a, b). These evidences suggest that the climatic transition from HS1 to BA is not a simple reverse of the transition from LGM to HS1.

From BA to YD, all proxy records, except for the trace-element record from Haozhu Cave15, indicate a drying climate over the whole of East Asia (Fig. 5g–i and Supplementary Table 7). The iTRACE simulation shows that annual, summer, and autumn rainfalls significantly increase in southern China, southwestern China, and North China Plain, but decrease in northeastern China, semi-arid northern China, and northeastern Tibetan Plateau (Fig. 5g–i). The simulated summer rainfall in the northern part of northeastern China also increased during the YD, but not significantly. The TRACE-21ka simulation shows rainfall increase in both northern and southern China (Supplementary Fig. 15g–i). The global temperature during the transition from BA to YD is characterized by strong cooling in the Northern Hemisphere and strong warming in the Southern Hemisphere (Fig. 6e, f). Unlike the temperature response during the HS1, the simulated summer cooling during the YD is not limited to North Atlantic, and the summer temperature in East Asia decreased during the YD in the iTRACE simulation (Fig. 6f). This difference may explain the better agreement between proxy records and model simulation during the YD than HS1.

From YD to early Holocene, all proxy records, except for a trace-element record from Haozhu Cave15 and a hopanoid record from Dajiuhu peatland77, indicate a wetting climate over the whole of East Asia (Fig. 5j–l and Supplementary Table 8). This pattern is generally consistent with the model results from the TRACE-21ka simulation (Supplementary Fig. 15j–l). The rainfall response at that time in the iTRACE simulation is in good agreement with proxy records from northeastern China, northern China, southwestern China, and the Leizhou Peninsula (Fig. 5j–l). However, the iTRACE simulation shows a rainfall decrease in southeastern China where various proxy records indicate an increasing rainfall. This divergence is robust and cannot be attributed to shortcomings in proxy records, because a well-dated speleothem carbon isotope (δ13C) record from Shennong Cave clearly indicates higher rainfall during the early Holocene than the YD78. The transition from YD to early Holocene in the iTRACE is characterized by warming in both the Northern and Southern Hemispheres, with the strongest warming occurring in the Northern Hemisphere (Fig. 6g, h). This pattern is generally consistent with global temperature reconstruction61.

Taken together, there are large discrepancies between proxy records and climate models simulations in reconstructing the hydroclimate of southern China during the last deglaciation, particularly for the iTRACE simulation running with the Community Earth System Model (CESM). We argue that these discrepancies may be attributable to the inability of climate models in correctly simulating the EASM rainfall pattern. It has been suggested that the coupled general circulation models (CGCMs) have large biases in simulating the mean state of the tropical Pacific with an excessive equatorial Pacific cold tongue when compared to observations79. As a result, current CGCMs are unable to accurately reproduce the spatial pattern of EASM rainfall during different seasonal stages, with the rain belt shifting northward too quickly80. For example, the position of the spring and pre-Mei-Yu rain belts are shifted northward compared to rainfall observations, resulting in deficient rainfall in southern China80.

Link to the end of Green Sahara and associated feedbacks
The climate boundary conditions during the mid-to late Holocene are largely constant81. Therefore, abrupt shifts in EASM rainfall under gradual insolation forcing must involve strongly nonlinear feedbacks. We attribute it to the end of the Green Sahara and associated vegetation-dust-cryosphere feedbacks (Fig. 7). However, the timing of the termination of the Green Sahara has not been well defined, with previous studies showing different results34,35. Therefore, our first goal is to objectively estimate this timing. Based on six high-resolution and well-dated leaf wax hydrogen isotopic composition (δDwax) records from the West African and East African monsoon regions (Supplementary Table 10), we generate a synthesis of African summer monsoon that considers uncertainties in both proxy measurement and age model (Methods, Fig. 7e and Supplementary Fig. 16). Based on regime shift analysis54,55 and Bayesian change-point analysis56, the timing of the termination of the Green Sahara is estimated to be 4620±180 yr BP, after which the African summer monsoon experienced a rapid decline with the rate of change exceeding two standard deviations during the Holocene (Supplementary Fig. 16). The timing of abrupt decline in the African summer monsoon is estimated to be centered at 4350±140 yr BP (Supplementary Fig. 16), statistically coinciding with abrupt rainfall decrease at 4370±60 yr BP in our reconstruction based on the reduced chi-square test82 (χ2=0.017; P=0.9). This timing also coincides with strong reductions in the tropical trees and Sahelian grassland cover in the Saharan region as evidenced by a continuous pollen record from Lake Yoa in northern Chad83 (Fig. 7g). Similarly, these events coincide with the onset of maximum dust flux recorded in the Kilimanjaro ice core84 and marine sediment cores from the Nile Delta85 and the Gulf of Oman86 (Fig. 7f), the timing of which was estimated to be centered at 4260±66 yr BP87 (χ2=0.78; P=0.46). These evidences suggest that the abrupt reduction in the EASM rainfall during the mid-to late Holocene is synchronous with the rapid degeneration of regional vegetation and abrupt increase in dust emissions across the Sahara in response to the termination of the Green Sahara. We further examine the climatic effect of changes in the Saharan vegetation and dust by performing two sets of mid-Holocene (MH) sensitivity experiments using the European Community Earth System Model (EC-Earth): one control experiment driven by orbital forcing only (MHORB), and an idealized experiment (MHGSRD) incorporating orbital forcing and a vegetated Sahara with reduced dust emissions (Methods).

The deterioration of vegetation at the end of the Green Sahara could influence the climate by increasing surface albedo88 and enhancing dust mobilization36, which jointly reduce the absorption of solar radiation and enhance radiative cooling of the troposphere. The climatic impacts of the end of the Green Sahara are not limited to northern Africa and can affect the global climate from the tropics to the poles89. Our sensitivity experiments suggest that the reduced vegetation cover and increased dust emissions in the Saharan region could cause global surface cooling during both summer and winter, with the most pronounced cooling occurring in northern Africa, West Asia, the Tibetan Plateau, and the northern high latitudes (Fig. 8a, c, e). The cooling of northern Africa and West Asia is directly related to the increase in local land-surface albedo and atmospheric dust concentrations36, while the cooling of the Tibetan Plateau and northern high latitudes is linked to the amplification of cryosphere feedbacks, such as changes in snow cover and sea ice extent39. The asymmetric cooling between the Northern and Southern Hemispheres due to different land-sea configurations could alter the interhemispheric temperature gradient (Fig. 7c)90 and shift the ITCZ southward, thereby reducing the rainfall in East Asia. Additionally, the cooling of the northern high latitudes can influence the EASM rainfall through changes in the North Atlantic Oscillation (NAO). The NAO is defined as the difference in surface sea level pressure between the subpolar Icelandic Low and the subtropical Azores High, which controls the strength and direction of the westerlies and storm tracks across the North Atlantic91. At 4384±23 yr BP (weighted mean), the reconstructed NAO index shifted from primarily positive phases to a weakly positive phase with intermittent strong negative phases92 (Fig. 7b), which closely tracks the change in the EASM rainfall (χ2=0.047; P=0.83), indicating an inherent dynamical connection between them. The shift of NAO into more negative phases in response to the termination of the Green Sahara can also be reproduced in our climate models simulations, which indicate increasing sea level pressure in the northern Europe and/or Nordic Seas together with decreasing sea level pressure in the Azores and/or southern Europe (Fig. 8b, d, f). The NAO persistently influences the surface temperature and snow cover over the Eurasian continent, sea ice and SST in the North Atlantic from winter to summer93. The NAO-induced SST anomalies in the North Atlantic can excite a wavetrain propagating toward East Asia, which affects the large-scale atmospheric circulation and, thus, regional rainfall pattern94. In addition, the negative NAO could weaken the wind-driven AMOC and reduce the poleward ocean heat transport in the North Atlantic39, further reducing the interhemispheric temperature gradient and the EASM rainfall. However, multiple lines of evidence suggest that the basin-scale AMOC remained fairly stable throughout the Holocene95. Therefore, the feedbacks associated with the Holocene AMOC likely play a minor role. In contrast, the Tibetan Plateau, known as the world’s third pole with an average elevation of 4000 m, may play a more significant role in shaping the EASM rainfall. The snow-albedo feedback may have amplified the temperature response on the Tibetan Plateau to the end of the Green Sahara. An alkenone-based summer temperature reconstruction from the Qinghai Lake on northeastern Tibetan Plateau indicates a substantial cooling (>5°C) between 5 and 3.5 ka96 (Fig. 7d). Since the Tibetan Plateau is the main heat source driving the Asian summer monsoon97, its abrupt cooling during the end of the Green Sahara could weaken land-sea thermal contrast and greatly reduce the EASM rainfall.

To reveal the spatial pattern of EASM rainfall response to the termination of the Green Sahara, we collect a wealth of hydroclimate records in East Asia based on rigorous criteria (Methods, Supplementary Table 9) and compare them with four climate model simulations (Fig. 9). In addition to our sensitivity experiments conducted with the EC-Earth model36 (Methods), similar experiments based on the Community Climate System Model version 4 (CCSM4)98, the water isotope-enabled Community Earth System Model (iCESM)99, and climate model from the Goddard Institute for Space Studies (GISS)100 were included for multi-model comparison. Consistent with the proxy records (Supplementary Table 9), nearly all climate models simulated decreasing annual and summer rainfalls in northeastern China, northern China, northeastern Tibetan Plateau, and southwestern China in response to the desertification of the Sahara (Fig. 9). However, the simulated increases in annual and summer rainfalls in southern China contrast with the available proxy records there (Supplementary Table 9 and Fig. 9). Similar to the last deglaciation, we argue that the biases in climate models are likely responsible for this model-proxy discrepancy. We also analyze the changes in autumn rainfall in climate simulations, which show a better agreement with proxy records (Fig. 9). In particular, the EC-Earth and CCSM4 models simulate a widespread autumn drying over the whole East Asia, closely resembling the pattern of proxy records (Fig. 9).

In summary, our quantitative reconstructions of the EASM rainfall amount changes over the past 22,000 years based on multi-methods and multi-proxies provide a benchmark for understanding the forcing and response of EASM rainfall in the context of long-term climate change. Our reconstructions suggest that the EASM rainfall since the LGM is characterized by alternative contrasting stable states and abrupt shifts between these states, indicating a strongly nonlinear nature of the EASM rainfall. In particular, we uncover and define five prominent tipping points in the EASM rainfall, which are linked to abrupt shifts in the AMOC and/or Saharan vegetation cover. In combination with comprehensive paleoclimate data compilations and multi-model simulations, we have drawn a picture of the EASM rainfall responses to abrupt shifts in the AMOC and Saharan vegetation. Our finding that the EASM rainfall was tightly coupled with known tipping elements in the Earth system during the past has practical implications for future change in the EASM rainfall under anthropogenic global warming. The potential “tipping” of tipping elements when global warming exceeds a specific threshold can trigger domino effects on the EASM rainfall and threaten the lives of billions of people.

Methods
Sitedescriptionandmodernclimate
TheQingtongyangMaarLake(110°10′E,20°34′N,elevation131m)islocated in the Leizhou Peninsula of Guangdong Province, southern China. It is a closed-basin crater lake that has a surface area of 8.81 km2 and a small catchment of 32.13 km2. Because of the small catchment and lack of river input, the lake level is mainly controlled by the rainfall amount falling in the catchment, which makes it an excellent measure of past rainfall changes. Also, this lake lies in the northern boundary of the present-day seasonal ITCZ migration, and thus, the lake status should have been sensitive to past changes in the position of the ITCZ. The surrounding bedrock of this lake is composed of volcanic basalt101. No limestone in the catchment minimizes the possibility of introducing old radiocarbon into the lake, which lends credit to 14C dating and ensures the reliability of age model.

Based on station data available from the China Meteorological Data Service Centre (https://data.cma.cn) during the period 1957–2018 CE (Supplementary Fig. 2), the long-term annual mean temperature (Tann) is 23.3°C with the mean coldest month (January) temperature of 15.6°C and the mean warmest month (July) temperature of 28.8°C. The annual total rainfall (Pann) is 1611 mm, with the driest month (December) rainfall of 24 mm and the wettest month (August) rainfall of 327 mm. The rainy season is from May to October under the influence of the EASM and northward migration of the ITCZ. Since more than 80% of the annual rainfall is contributed by summer rainfall, the variations of annual rainfall are almost identical to that of summer rainfall (Supplementary Fig. 2).

Drillingandlithology
The sediment cores were recovered from the center of Qingtongyang Maar Lake with piston corer in November 2016. Here we report the upper part of the core, which is coded as QTY-2-1 and spans the last ~22,000 years with firm age control. The lithology of the core is described as follows: 0–5 cm is a surface mixed layer, 5–8 cm is grayish black clay, 8–34.4 cm is grayish brown silty clay, 34.4–73 cm is grayish black clay, 73–82 cm is grayish yellow sandy clay, 82–89 cm is grayish black clay, 89–95.5 cm is grayish yellow sandy clay, 95.5–123 cm is gray-brown yellow sandy clay.

Datingandchronology
210Pb and 137Cs dating. Samples at 0.5-cm intervals from the top 16 cm of the sediment core were measured for the natural fallout radionuclide 210Pb and anthropogenic fallout radionuclides 137Cs. Their activities were measured by the gamma spectrometer using a low-background, hyper-pure germanium detector. Each sample was counted for at least 24 h to ensure that both detection limit and measurement error are acceptable. Prior to measurement, samples were freeze-dried, weighed, and stored in sealed containers for 3 weeks to allow for radioactive equilibrium between 214Pb and its parent radioisotope 226Ra102. The activities of 210Pb and 137Cs were determined by gamma emissions at 46.5 and 662 keV, respectively. 226Ra activity was measured by the gamma emissions of its daughter radioisotope 214Pb at 295 and 352 keV. The excess 210Pb activity (210Pbex) was calculated by subtracting the 226Ra activity from the total 210Pb activity. Because the sediment was mixed in the top 0–5 cm, we constructed the 210Pb age model based on the piecewise constant rate of supply (CRS) model49 constrained by the 137Cs peak at 1963 CE and two post-bomb 14C ages at 1955 CE and 2014 CE (Supplementary Fig. 6), which was implemented in the R package ‘serac’103.

AMS 14C dating. Sixteen charcoal samples from various depths of the sediment core were collected for AMS (accelerator mass spectrometry) 14C dating (Supplementary Table 1). Samples were chemically pretreated using standard acid-alkali-acid method104 and graphitized on a fully automated system105. 14C content was measured on the MICADAS system106 in the AMS laboratory of Nanjing University107. The positive 14C dates were calibrated to calendar years before present (BP), where present refers to 1950 CE, based on the IntCal20 curve108. The negative 14C dates were calibrated using monthly data of the Northern Hemisphere post-bomb calibration curve3109. Calibrations were implemented in the software OxCal110. For stability, the 14C age model is based on linear interpolation and extrapolation of the median values of the calibrated 14C dates, with two outliers excluded (Supplementary Fig. 5). According to the age model, the bottom of core QTY-2-1 is estimated to be 22,680 yr BP, but we limit our analysis to the last 22,000 years to account for potential disturbances at the core’s end during drilling.

Corescanningandcalibration
The color of sediment is an indicator of sediment composition changes and has been used successfully as a proxy for paleomonsoon precipitation changes111–115. Digital images of the sediment core surface were taken by a high-resolution line scan camera installed in the Geotech MSCL-XRF core scanner. Sediment color b*, representing the blue (negative) to yellow (positive) component, was extracted from digital images in a rectangle along the core direction with a width of ~15 mm and an unprecedented sampling interval of 100 μm, yielding an overall temporal resolution of ~1 years and 0.1-year resolution in the period overlapping with instrumental records.

Continuous and non-destructive X-ray fluorescence (XRF) element scanning was carried out at 1 mm resolution on the split sediment cores, corresponding to an overall sampling resolution of ~10 years and 1-year resolution in the period overlapping with instrumental records. The sediment core surfaces were covered with a 4 μm-thick thin film to avoid contamination of the scanning instruments. All XRF scanning measurements were conducted with a generator setting of 10 KV and a sampling time of 20 s. The irradiated dimensions of the core surface are 15 × 1 mm2. The split sediment cores were freeze-dried before scanning to avoid the influence of water content on the results of XRF116. The geochemical elements were calibrated using the immobile element titanium (Ti) to isolate the proportion of changes associated with climate:

Xcalib = Xraw × (Xraw + Ti) / Ti (1)

In Eq. (1), Xraw represents the raw value and Xcalib represents the calibrated value. From a physical and chemical perspective, Xcalib represents the proportion of element change attributable to redox conditions and/or salinity. A similar process was used to quantify the near-infrared reflectance of terrestrial vegetation based on satellite remote sensing117. This equation is an onlinear transform of the element ratio (Xraw/Ti), which is insensitive to dilution effects and can accommodate the inherent non-linearity of the relation between element intensity and concentration118. As a ratio, it has the advantage of minimizing some element-correlated noise and influences attributable to variations in direct/diffuse irradiance, surface morphology, and atmospheric attenuation116,118,119. To a lesser degree, the ratio can reduce instrument-related errors119. The distribution of iron after calibration was still highly skewed and it was further logarithmically transformed to make its response to climate more linearly.

Rainfallreconstruction
Instrumental target. Our reconstruction target is annual rainfall based on the calendar year that covers a full cycle of the rainy season. Monthly rainfall deviations from the climatological mean (1961–1990 CE) of three meteorological stations surrounding the Qingtongyang Maar Lake were averaged to represent regional rainfall variability (Supplementary Fig. 2). The stations from Hainan Island were not included because rainfall variability on Hainan Island is poorly correlated with that of the Leizhou Peninsula. Climate data were downloaded from the China Meteorological Data Service Centre (https://data.cma.cn). The period 1953–2005 CE were used for calibration because of the least missing values and overlapping with proxy data. The period prior to 1953 CE was not used for calibration because of large uncertainty in the climate data. Missing climate values were filled in using the regularized expectation maximization (RegEM) algorithm120 with truncated total least squares regression121. Grided rainfall data from the Global Precipitation Climatology Centre (GPCC) version 202244 were used as an additional independent validation of our reconstructions. To preserve the full amplitude of climate variability and avoid spurious signals, both climate data and proxy records were decadally smoothed with a cubic-smoothing spline45 to have roughly the same frequency (Supplementary Fig. 3), following the suggestion of previous studies46–48. Since smoothing can influence the degrees of freedom (DOF), we assumed DOF=51 for the raw data, which has an effective sample size (ESS) of 53 years, and DOF=3 for the smoothed data, which has an ESS of five decades47.

Proxy data. Proxy data measured on core QTY-2-1 from the Qingtongyang Maar Lake were used as predictors for rainfall reconstruction, including color b* and elemental concentrations of iron, manganese, and calcium. These proxy records are sensitive to redox conditions and/or salinity in the lake which are controlled by regional rainfall (Supplementary Discussion). We decided to choose these four proxy records based on three criteria. Firstly, from the statistical perspective, there should be a significant correlation (P < 0.02) between proxy and climate targets in their common time period. We utilized both correlation and partial correlation to analyze the relationship between our proxies and climate data, which shows that rainfall is the dominant controlling factor of changes in these proxy records (Supplementary Table 2). The impact of smoothing on the degrees of freedom was considered during statistical analysis. Secondly, from the physical and chemical perspective, the qualitative relationship between proxies and climate can be explained and the mechanisms have been relatively well recognized (Supplementary Discussion). Lastly, the temporal resolution of proxy records should be at least one year or finer during the calibration period to facilitate meaningful calibration with rainfall observations. To test the stability of the climate-proxy relationship over time, we calculated the running expressed population signal (EPS)53 and mean inter-series correlation (Rbar)53 using a sliding window of 3000 years (Fig. 2e).

Ensemble parameters. To address various sources of uncertainty, 1650 members of reconstructions were obtained for each reconstruction method by perturbing the following parameters:
(1) 0–25% of proxy records were randomly removed before reconstruction to test for the robustness of proxy selection122.
(2) To overcome bias in the choice of calibration period46 whilst retaining the autocorrelation structure123, all possible combinations of 5-year blocks and additional successive years with a length of 33/38 years were used for training, while the remaining 20/15 years were used for validation.
(3) Gaussian noises with standard deviations equal to measurement errors were added to the proxy data to simulate the influence of measurement uncertainties. This treatment can improve the generalization of the reconstruction model and avoid overfitting.
(4) The weight of each proxy is multiplied by a random scalar drawn from a uniform distribution in the interval between 1/1.5 and 1.5124. This treatment simulates the influence of proxy weights122 and age model uncertainties (Supplementary Fig. 4). Before weighting, all proxy records were normalized to have zero mean and standard deviation units over the calibration period.
(5) Red noise with its standard deviation equal to the root mean square error (RMSE) and its first-order autoregression (AR1) coefficient identical to that of regression residual was added to the result to represent the unexplained variance.

Composite plus scale (CPS). CPS is a classical method for quantitative climate index reconstruction and has been widely used to reconstruct mean climate of various spatial scale47,125,126. The proxy matrix is normalized over the calibration period and weighted averaged to form a single time series, which is then scaled to have the same mean and standard deviation as the reconstruction target. Here we perform CPS in a way similar to previous studies122,124: the normalized proxies were weighted by their non-detrended correlation with rainfall observation over the calibration period. Like many linear methods, CPS is prone to variance losses and the reconstructions tend to be biased toward zero, which are directly related to the correlation between proxies and the climate target127. CPS is computationally efficient, and its performance is often comparable to that of more complicated methods (Supplementary Fig. 7).

Principal component regression (PCR). PCR is a popular method commonly used in both climate index and climate field reconstructions47,122,124,125,128–130. PCR reduces the dimension of the proxy matrix by principal component analysis (PCA) and only retains the first few principal components (PC) as predictors in the regression. All PCs are linear combinations of the original proxies and orthogonal to each other, which copes well with the problem of multicollinearity. Here we perform PCA over the calibration period using the singular valued decomposition (SVD) algorithm131, retaining the first few PCs that explain up to 95% of the total variance to reduce noise. If PC1 explains more than 95% of the total variance, we retain only PC1. The number of retained PCs was selected in such a way because there exist many methods to truncate PCs without an objectively discernible best method, and the truncation is sensitive to the period over which the PCA is performed125.

Partial least squares (PLS). Though PLS regression is a basic tool widely used in chemometrics132, it is less seen in tree-ring-based climate reconstructions compared to CPS and PCR, except in some pollen- and diatom-based quantitative reconstructions133. PLS regression is of particular interest because it has the ability to analyze data with numerous noisy and collinear variables in two data matrices, X and Y132. Similar to the PCR method, PLS constructs components as linear combinations of the original predictor variables but has more advantages. In contrary to PCR, which constructs components to maximize the explained variance of predictors, PLS finds components that have a large covariance with the response variables. Therefore, PLS regression can achieve a parsimonious model with reliable predictive power. Here, we perform PLS regression using the SIMPLS algorithm134 and select the optimal number of PLS component based on two-sample t-test (P < 0.01)135 and require a reduction of RMSE in the validation period by at least 5%136.

Optimal information extraction (OIE). The OIE method is a variant of the CPS method, which integrates the features of local (LOC) method137, a Bayesian framework138, the generalized likelihood uncertainty estimation (GLUE)139,140, and ensemble reconstruction46,122. This method has the advantage of efficiently extracting low-frequency climate signals, since it accounts for the nonlinear response of proxy to climate based on Bayesian theory141. Instead of variance matching in the CPS, the OIE sets the prior distribution of regression coefficients as a uniform distribution in the ranges between direct and indirect regressions, and the posterior distribution is determined by the GLUE method139,140. Here we run the Markov Chain Monte Carlo (MCMC) simulation in the GLUE using the Metropolis–Hastings algorithm142 and calculate the acceptance probability of regression coefficients based on the RMSE of the validation period143. The MCMC sampling process generates 1000 realizations of reconstruction and the median is used as one of the ensemble members.

Point-by-point regression (PPR). PPR is a well-tested and easily interpreted PCR method that has been applied to create high-quality drought atlases worldwide129,144. PPR assumes that only those proxy records located relatively close to a given point are likely to be true predictors of climate at that location, where “true” indicates a causal relationship between proxy and climate that is stable through time. The PPR is initially designed for climate field reconstruction but can also be used for climate index reconstruction. Unlike PCR, which performs PCA over the calibration period, the PPR performs PCA over the whole period of proxy records and puts the orthogonal PCs into stepwise regression by the order of explanatory variance, which can avoid the insufficient representativeness of a short calibration period and extract the common variability linked to climate change. Here we determine the optimum number of PCs put into regression model based on two-sample t-test135 and a threshold of at least 5% reduction in the RMSE in the validation period136.

Ridge regression (RIG). Ridge regression is an approach for determining the coefficients of a linear model that contains linearly correlated predictors145. Estimations of coefficients for multiple linear regression models rely on the independence of model terms. However, when terms are correlated, and the columns of predictor matrix X have an approximate linear dependence, the matrix (XTX)-1 is close to singular. As a result, the least squares estimate is sensitive to random errors in the response variable Y, producing large variance, which is called multicollinearity. Ridge regression resolves the problem of multicollinearity by estimating regression coefficients using β = (XTX + kZ)-1XTY, where k is the ridge parameter, and Z is the identity matrix. Small, positive values of k improve the conditioning of the problem and reduce the variance of the estimates. While biased, the reduced variance of ridge estimates often results in a smaller mean squared error when compared to least squares estimates. Here, we set k as a random scalar evenly distributed in the range from 0 to 100, and its final distribution is determined using the GLUE method139,140 and other details on the MCMC sampling are identical to that of the OIE method.

Lasso and elastic net (ELN). Similar to the ridge regression, the lasso is a regularization method for performing linear regression146. However, lasso is a shrinkage estimator which includes a penalty term constraining the size of the estimated coefficients and thus forcing some coefficients to be zero. As the penalty term increases, lasso lets more coefficients to be zero. Elastic net is a hybrid of ridge regression and lasso regularization147. Like lasso, elastic net can generate reduced models by generating zero-valued coefficients, but it outperforms lasso on data with highly correlated predictors. Here, we perform elastic net by setting the weight of lasso versus ridge optimization as 0.5, and the regularization coefficient ‘lambda’ is determined by ten-fold cross-validation and equals the maximum “lambda” value such that the RMSE is within one standard error of the minimum cross-validation error.

Artificial neural network (ANN). ANN is a machine learning method that imitates biological neural network, a part of artificial intelligence technology148,149. It has the ability to learn and make inferences from data, mimicking the way that the brain processes information. ANN model consists of a series of layers, namely the input, hidden, and output layers. The first hidden layer has connections from input data, and each subsequent layer has connections from the previous layer. Each layer multiplies the input by a weight matrix and adds a bias vector. An activation function follows each hidden layer. The last layer produces the network’s output, namely predicted values. ANN with only one hidden layer can solve any problem of finite input-output mapping, given enough neurons in the hidden layer. The model used here is a fully connected feedforward backpropagation ANN with three hidden layers containing 100, 25, and 10 neurons, respectively. The linear transfer function is used. Data from the validation period are used to monitor training convergence every ten iterations.

Convolutional neural network (CNN). CNN is a deep learning method designed to process data in the form of multiple arrays150. CNN is commonly applied in tasks such as image detection, recognition, and classification, but can also be used for regression problems151. It is the convolutional layer and pooling layer that make CNN unique. Our CNN consists of seven layers in order: a sequence input layer, a one-dimensional (1D) convolutional layer with filter size of 3 and neurons of 32, a batch normalization layer to speed up training and reduce the sensitivity to model initialization152, a 1D average pooling layer with pool size equal to filter size, a Monte Carlo dropout layer with dropping probability of 50% to prevent the model from overfitting153, and a fully connected layer followed by a regression output layer. The weights of the model are initialized with the He initializer, which samples from a normal distribution with zero mean and variance of 2 divided by input size154. Padding was conducted so that the output has the same size as the input and reduces the influence of the edge effect. Data augmentation was performed by adding Gaussian noise to the input data to make the model more stable and avoid overfitting155. The model was trained with the Adam optimizer with a maximum epoch of 100, minimum batch size of 20%, and initial learning rate of 0.01. Gradient clipping was applied to prevent gradient exploding. Furthermore, we conducted independent validation every ten iterations and stop the training early with a patience of five epochs, which is another regularization method to prevent overfitting. The validation data was shuffled before each validation.

Verification statistics. We use four verification statistics to assess the reliability of the reconstruction models. These statistics evaluate the ability of the proxy to reconstruct climate over verification periods that are independent of calibration periods. (i) Pearson correlation coefficient (R), a powerful statistic for testing the relative association between two variables. (ii) Reduction of error (RE), which determines if a reconstruction is better than climatology, i.e., the mean of rainfall in the calibration period. There are no formal statistical significance tests for RE, but a RE > 0 indicates that the reconstruction is better than using the mean of the calibration period as forecast. RE is a rigorous verification statistic because it has no lower boundary. (iii) Coefficient of efficiency (CE). Like RE, there are no formal statistical significance test for CE, but a CE > 0 indicates that the reconstruction is better than using the mean of the verification period as forecast. CE is more rigorous than RE129. (iv) Root mean square error (RMSE), which measures the average difference between the reconstructed and observed rainfall.

Compilation of East Asian hydroclimate records
To test the robustness of abrupt shifts in East Asian hydroclimate, we compile proxy records from various climatic zones based on rigorous criteria: (i) Reflect local and/or regional summer hydroclimate changes. Because the climatic interpretation of speleothem δ18O records in East Asia remains controversial13, we only include the speleothem δ18O records that have been linked to hydroclimate change based on the correlation with instrumental records and/or independent rainfall-related proxy records such as δ13C records. (ii) The proxy records should at least cover the transition period from mid to late Holocene at 5–3 ka BP and/or one of the millennial-scale events during the last deglaciation (HS1, BA, YD) without hiatus longer than 500 years. (iii) The standard errors of age-control points should be less than 100 years in most cases and no larger than 800 years. (iv) At least three age-control points span or closely bracket the period of interest. The gap between successive ages should be smaller than 3000 years for the period 18–11 ka BP, and there should be at least one age-control point within the period 3.5–4.5 ka BP. (v) Has a mean resolution finer than 200 years. But this was relaxed to 300 years for regions where proxy records are scarce and 500 years for the lake level records. Details of the compiled proxy records are listed in Supplementary Tables 5–9 and some records are shown in Supplementary Figs. 12, 13. Most records are not continuous over the past 22,000 years and the source data are not publicly available. Therefore, it is impossible to make a continuous composite time series to describe the evolution of the East Asian hydroclimate since the LGM. However, these records can provide insights into hydroclimate changes during specific abrupt climate events, which were compared with model simulations.

Synthesis of African summer monsoon
At a total of six δDwax records (Supplementary Table 10) were combined into a composite time series to characterize the evolution of the African summer monsoon, encompassing both the West African and East African monsoons, during the Holocene. The proxy records were selected based on these criteria: (i) Reflect local and/or regional summer hydroclimate rather than winter, excluding proxy records from areas currently influenced by the Mediterranean climate. (ii) Span a minimum duration of 9000 years with no recorded hiatus during the period 11–0 ka BP. (iii) There should be at least one age-control point every 3000-year interval. (iv) Has a median resolution finer than 300 years in the period 6–2 ka BP; (v) Both the sampling depths and proxy data are publicly available for reanalysis of the age model.

There are two main sources of uncertainty in the paleoclimate data, i.e., proxy measurement uncertainty and age model uncertainty156. We generated 10,000 age models using the R package ‘rbacon’157 for each record and perturbed the proxy by adding random noise drawing from a normal distribution with a standard deviation equal to the measurement error. Then the perturbed proxy records were interpolated against the perturbed age models at 10-year resolution using the shape-preserving piecewise cubic interpolation method. Prior to calculating the composite mean, the interpolated δDwax records were corrected for ice volume change158 and subtracted by the mean and divided by the standard deviation during the period 12–2 ka BP that is common to all records to form a standardized time series. Furthermore, the standardized time series were multiplied by −1 so that greater proxy values indicate higher rainfall amounts.

Detection of abrupt shifts and their timings
Three different methods were applied to detect abrupt shifts associated with tipping points in the proxy records, including the regime shift analysis54,55, Bayesian change-point analysis56, and the rate of change. Regime shifts are defined as rapid reorganizations of a system from one relatively stable state to another54. The regime shift analysis54,55 splits a time series into multiple segments with significantly different mean and identifies the point at which abrupt changes in the mean occur. Regime shift analysis was performed using the function “ischange” in the software MATLAB with the maximum number of change points set to 5 for the past 22,000 years. To validate the results of regime shift analysis, a Bayesian ensemble algorithm56 was independently used to detect abrupt shifts in the proxy records. Unlike conventional criterion-based methods that select a single best model, the Bayesian paradigm embraces all candidate models, evaluates their probabilities, and synthesizes them into an average model56. These two independent methods perform well in detecting abrupt shifts in our rainfall reconstruction and show consistent results (Supplementary Fig. 10). Once the timings of abrupt shifts in our rainfall reconstruction were determined, the corresponding depths were calculated, and Bayesian age modeling157 was employed to assess age uncertainties. However, for low-resolution records, the timings of abrupt shifts may not align well with those detected by these two methods. In this case, the rate of change was used to determine the precise timings of abrupt shifts. The rate of change in the composite time series of the African summer monsoon was calculated, and the period with the rate of change exceeding two standard deviations of the period 11–0 ka BP was considered as the range of abrupt change in the time series (Supplementary Fig. 16). Since the composite time series of the African summer monsoon has been smoothed by averaging different realizations of uncertainties in both proxy measurement and age model, no additional smoothing was applied. The rate of change was not applied to detect abrupt changes in our rainfall reconstruction because the changes in ultra-high-resolution records are very abrupt, and the selection of smoothing parameters, such as window length, is quite subjective and cannot adapt universally to abrupt events with different amplitudes.

Climatemodelsimulations
TRACE-21ka and iTRACE. To reveal the spatial pattern of rainfall response in East Asia to abrupt shifts in the AMOC during the last deglaciation, we compare our rainfall reconstruction and the proxy records we compiled with model results from two transient climate experiments, i.e., the TRACE-21ka71 and iTRACE16. The TRACE-21ka simulated continuous climate evolution of the last 21,000 years using the Community Climate System Model version 3 (CCSM3) with a spatial resolution of T31 (3.75°)71. The transient climate forcings include orbitally-induced solar insolation variations159, atmospheric greenhouse gas concentrations160, meltwater fluxes71, and continental ice sheets161. The iTRACE was performed in the state-of-the-art water isotope-enabled version of the Community Earth System Model version 1.3 (iCESM1.3)16, which largely follows the strategy of the TRACE-21ka71. The iCESM1.3162 is composed of the Community Atmosphere Model (CAM5.3), Parallel Ocean Program (POP2), Los Alamos Sea Ice Model (CICE4), and Community Land Model (CLM4). The spatial resolution of atmosphere and land is 1.9° in latitude and 2.5° in longitude, with 30 vertical levels in the atmosphere. The horizontal resolution of ocean and sea ice is nominal 1° (gx1v6), with 60 vertical levels in the ocean.

All experiments in the iTRACE were integrated from the LGM (20 ka) to the early Holocene (11 ka). Three parallel transient sensitivity experiments (ICE, ICE+ORB, and ICE+ORB+GHG) were performed with three climatic forcing added one by one, i.e., the ice sheet (ICE-6G)163, the orbital forcing159, and the greenhouse gas forcing160. The baseline iTRACE simulation (i.e., ICE+ORB+GHG+MWF) was branched from the ICE+ORB+GHG run at 19 ka, with the imposed meltwater forcing (MWF) similar to that in the TRACE-21ka71. The ice sheet in the iTRACE was changed every thousand years from 19 ka to 11 ka. The ocean bathymetry was modified twice at 14 ka and 12 ka based on the ICE-6G reconstruction163.

GreenSaharaexperiments. The Green Sahara sensitivity experiments were performed using the European Community Earth System Model (EC-Earth) version 3.1164. The EC-Earth is a fully coupled Earth system model that includes several state-of-the-art components to describe the atmosphere, ocean, sea ice, land surface, dynamic vegetation, atmospheric composition, ocean biogeochemistry, and the Greenland Ice Sheet165. The atmospheric component is based on the Integrated Forecasting System (IFS cycle 36r4) of the European Centre for Medium-Range Weather Forecasts with a horizontal resolution of T159 (1.125°) and 62 vertical levels. The ocean component is based on the Nucleus for European Modelling of the Ocean (NEMO), which uses a tripolar grid with a horizontal resolution of 1° and 46 vertical levels. The EC-Earth has proven its success in simulating the Earth’s climate across past and future scenarios165,166. It shows good skills in capturing the nuanced spatial and temporal patterns of monsoon rainfall in present-day climate36,167.

To understand the response of rainfall in East Asia to the changes in Saharan vegetation during the end of the Green Sahara, we performed two sets of mid-Holocene (MH) experiments: one control experiment driven by orbital forcing only (MHORB), and an idealized experiment (MHGSRD) incorporating both orbital forcing and a vegetated Sahara with reduced dust emissions. The MHORB experiment was performed following the protocol of the Paleoclimate Modeling Intercomparison Project phase 3 (PMIP3)168 with orbital parameters fixed at 6 ka169. The concentrations of atmospheric CO2, CH4, and N2O were fixed at the preindustrial level, i.e., 280 ppm, 650 ppb, and 270 ppb, respectively. Other boundary conditions, such as land-sea distributions, ice sheets, topography, vegetation cover, and dust emission were the same as the preindustrial period. In the MHGSRD experiment, the vegetation type over the Sahara region (11°N–33°N and 15°W–35°E) was modified to shrub, and the dust emission from the Sahara Desert was reduced by 80%36 according to the proxy-based estimation of Saharan dust flux reduction during the mid-Holocene34,170. The change of vegetation cover from desert (MHORB) to shrub (MHGSRD) corresponds to a decrease in surface albedo from 0.3 to 0.15 and an increase in the leaf area index from 0.2 to 2.636. The dust reduction from the MHORB to the MHGSRD led to a decrease in the global dust aerosol optical depth (AOD) by about 60% and a total AOD of 0.0236. Such changes in vegetation cover and dust emissions were not intended to precisely replicate mid-Holocene Saharan conditions but were designed to provide insights into their potential climate feedback.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of the manuscript “Tipping point-induced abrupt shifts in East Asian hydroclimate since the Last Glacial Maximum”\nby Lu et al. [NCOMMS-24-28148]\nGeneral comments:\nThe East Asian summer monsoon (EASM), one of the key components of the Earth climate system, has been extensively\ninvestigated for both the processes and the forcing mechanisms on different timescales based on multi-proxy analyses of the\nChinese loess–paleosol sequences, stalagmites, and the other sedimentary archives. However, how the EASM behaves\nwhen global warming exceeds the specific threshold remains unknown due to the lack of quantitative records of the summer\nmonsoon rainfall in the geological past. In this study, the authors present a robust reconstruction of the EASM rainfall\nchanges in southern China over the past 22,000 years based on a multi-method, multi-proxy framework, and suggest that the\nEASM rainfall was characterized by abrupt shifts with a magnitude of ~500 mm (30%). Combined with climate simulations,\nmoreover, the authors propose that these abrupt shifts in the EASM rainfall would be linked to tipping points in the Atlantic\nMeridional Overturning Circulation (AMOC) and the Saharan vegetation cover. The data and inferences presented in this\nstudy are of great significance and would contribute to a deep understanding of the response of the EASM to future tipping\npoints. I recommend acceptance of this manuscript for publication in Nature Communications after minor revisions.\nSpecific comments:\n1. Lines 64–84: it is better to explain the physical link between the EASM and the AMOC / Saharan vegetation cover in this\npart rather than only to emphasize the importance of both in the Earth climate system.\n2. Lines 85–94: This part contains some of the results, discussion and conclusions, and should be moved separately to the\nfollowing sections ‘Results and discussion’ or ‘Conclusion’.\n3. Lines 201–207: Did the cooling in the North Atlantic caused by the end of the Green Sahara have any impact on the\nAMOC and then on the EASM rainfall?\n4. Lines 215–216: What is your hypothesis to explain such a conflict between model simulations and proxy records?\n5. Lines 237–239: I am confused about ‘a weakening AMOC’ and ‘a greening Sahara’ mentioned together with each other\nhere.\nMinor comments:\n1. Line 97: Change ‘sensitive’ to ‘closely related’.\n2. Line 110: Delete ‘consistently’.\n3. Line 111: Change ‘on’ to ‘of’.\n4. Line 132: Change ‘between states’ to ‘between the states’.\n5. Line 132: Change ‘According to’ to ‘Based on’.\n6. Line 141: Change ‘has’ to ‘had’.\n7. Line 142: Change ‘vast’ to ‘the whole’.\n8. Line 142: Move ‘the ’ before ‘observed rapid shifts’ to before ‘rainfall regime’.\n9. Lines 145–146: Change ‘the dynamic component’ to ‘the former’ and ‘the thermodynamic component’ to ‘the latter’.\n10. Line 151: Change ‘on’ to ‘of’.\n11. Figure 1: Show the latitude for the JJA insolation in Figure 1.\n12. Line 163: Change ‘for’ to ‘of’.\n13. Line 171: Insert ‘two’ before ‘hemispheres’.\n14. Lines 173–177: It may be better to remove these words about the relation between the ITCZ and the ENSO / the Walker\nCirculation in order to emphasize the relation between the AMOC and the ITCZ / the Westerlies.\n15. Line 224: Change ‘between states’ to ‘between the states’.", "reviewer_2": "The paper investigates precipitation change in the East Asian Monsoon domain during the last 22 kilo years. It does so by\nproviding new proxy records from a lake in South China. The records are post-processed and calibrated to estimate a\ntimeseries of precipitation, which the authors use as a proxy for East Asian Monsoon intensity. The authors particularly focus\non the abrupt changes that the timeseries exhibits. They link these changes to AMOC variability during the deglaciation, and\nto the end of the green Sahara during the mid-Holocene. They further compared their records to other proxy records and\nclimate simulations.\nI am not an expert on precipitation reconstruction from these particular proxy types. Yet, from a data analysis perspective, the\npost-processing and calibration methods used to retrieve the precipitation signal from the proxy records seem sound. In\nparticular, the large number of methods used to compute the reconstruction is impressive and valuable. This increases the\nconfidence in the reconstruction. This one also exhibits striking features that the authors discuss.\nThe authors also provided model data to compare their and other proxy records to. I find the analysis of these model outputs\nto be quite limited, and their added value insufficient. I provide suggestions regarding that matter below.\nFinally, the influence of AMOC and the end of the Sahara greening on the EASM is an already known feature, that has been\ndiscussed in previous studies both from a proxy and model perspective. The precipitation reconstruction from this study adds\nanother evidence to it.\nI have a few general suggestions that would help the reader better grasp what you are offering:\n- Based on your narrative \"We provide an ensemble reconstruction of the EASM based on a multi-proxy and multi-method\nframework\", I thought you were post-processing an ensemble of already published proxy records from the entire EASM\ndomain, and that your paper focused on the methods to derive the EASM reconstruction from these. I was therefore very\nconfused when I started to read the result section. I suggest that you make clear in the abstract and the introduction that you\nare publishing new proxy records from a single core/location and that you focus on retrieving the precipitation signal from\nthese records. I would also suggest further discussing how this location is representative of the entire EASM domain.\n- The evaluation of the uncertainty related to the specific method used is impressive. The discussion on the age model is\nalso quite clear. Is there however other known or potential limitation more specific to the proxy type and the calibration\nperiod used? In particular, how certain is it that the LGM precipitation was similar to the present day? Can the reconstruction\nof the long timescales be biased? Also, for a non-expert like me, I would welcome some references that explain how\nprecipitation influences the colour of the core and the concentration of iron, manganese, and calcium.\n- I find the comparison to other proxy records a bit lacking. For example, I did not find any information on the records used in\nFig. 4, and therefore cannot comment on it. The authors also show that the literature has been debating the heterogeneity of\nEASM response to remote forcing. However, the authors only mention once the \"tripole\" without any context to it, and mostly\ndismiss it as a possible feature. However, the tripole pattern is quite clear from Fig.4 b (although maybe a bit too\nSouth, possibly due to a biased location of the Meiyu front in the model), and is also the object of the paper that published the\nrecord you show in Fig. S4c (Shouldn't you flip the y-axis of this one, so that the top correspond to wetter and the bottom to\ndrier as for the other records?). I strongly suggest the author to further discuss this question of heterogeneity of the EASM\nresponse, and make use of the high resolution of their records to potentially discuss the timing and sequence of events\ncharacterising the abrupt changes. Another suggestion would also be to compare this reconstruction to other nearby proxy\nrecords, for example from Huguan Maar and Tianyang Maar, or even the Dongge cave. This could help to better understand\nregional features (or proxy limitations).\n- As it stands, I do not see the added value of the model data. I see three courses of action:\n- Entirely remove the model data, to better focus on the proxy records.\n- To compare the model data to previous studies, in order to bring a multi-model (or at least multi-simulation) approach to the\nanalysis. The authors could then discuss the reliability/robustness of their model data and potentially explain why we should\nbetter trust their higher resolution simulation.\n- The authors could then further investigate how AMOC / End of Sahara greening influences the EASM in their simulations.\nThis could also be done with sensitivity experiments to further understand which component of the forcing influence the\nEASM, and how. I am also not convinced about the use of PI boundary conditions in the hosing experiment. How different\nwould the AMOC impact on EASM be, if LGM boundary conditions were used?\n- I feel the focus on \"tipping points\" is either too extensive or too general. AMOC or End of Sahara greening are well known\ntipping points, and their impact on the EASM has already been discussed. In its current state, I do not see much added value\nfrom the paper that justifies this focus. However, I do see an opportunity, if the authors were to further discuss tipping points\nwithin the EASM domain, as part of a tipping point cascade. Neither Wand et al. 2023 nor Wunderling et al. 2024 recent\nreview on tipping points mention the EASM as a tipping element. The fact that the reconstructed time series exhibits a bi-\nstability could be a proof of the opposite. Simulations can be used to investigate this tipping point, and answer for example\nthe question: Is the EASM response linear to the freshwater hosing / AMOC slow down?\nNot all of these suggestions need to be implemented. It is up to the authors to shape their paper more as a publication of a\nnew record / reconstruction, as a comparison between proxy records and simulations, or as a discussion of tipping points in\nthe EASM domain.\nOther comments along the text:\n28: \"30%\". is this uncertainty or change compared to the mean value?\n30-31: \"the collapse of the Atlantic meridional overturning circulation [...] caused drying\": The re-establishment of the AMOC\nwas seemingly as abrupt with consequences on the EASM. Why not discuss that too?\n31-33: \"widespread drying\", \"drying in northern East Asia and southwestern China and possible\nwetting in southeastern China\". These are the raw results from the simulations. Either clarify that these are your model\nresults, or further discuss in the text why you think these results are better compared to other simulations.\n33-34: \"Our findings provide crucial evidence for projecting the response of East Asian hydroclimate to future tipping points.\"\nThe crucial evidence is not so evident to me. The current climate conditions are quite different to the AMOC collapse during\nthe deglaciation and the end of the Sahara greening.\n56-57: \"However, the interpretation of speleothem δ 18 O records in East Asia remains highly debated, excluding it as a\nlocal rainfall proxy\". You could also cite He et al. 2021 here. It quite clearly shows this. Does that mean you exclude all d18O\nrecords from your analysis? Could your new reconstruction help to understand other d18O records?\n63: \"particularly in southern China.\" Why there? How important is this area in the EASM?\n65: \"application of [...]\". Add here that you are using records from one lake.\n66: \"unprecedented accuracy\". How do you know the accuracy?\n65: \"to assess the impacts of climate tipping points on the EASM rainfall\". The phrasing of the goal is a bit strange. You\nfound abrupt changes in your timeseries, and you want to understand their cause. Comparison with other proxies helps you\nto formulate a hypothesis, that you then test with climate simulations.\n67: \"integrating\". You are just comparing the two. Integrating sounds like data assimilation, which you do not do.\n73: \"AMOC will collapse\". AMOC may collapse? That reflects better the uncertainty/controversy.\n74: \"climatic impact of AMOC collapse remains unclear.\" It is quite clear from a paleo perspective. What is unknown is what\nwill happen in the context of anthropogenic changes.\n82-84: These claims are far too generic. Maybe a better context would be the recent shifts / variability of the West African\nmonsoon, their relationship with vegetation in the Sahel, and potential teleconenction with the EASM?\n86: \"between states\". between these states\n85-96: Is this usual to have a paragraph summarising the results at the end of the introduction?\n101: \"annual rainfall and summer rainfall are equivalent\". They are only similar.\n102-105: First sentence on the fact that the core was recovered and the records measured for this study. Second sentence to\nsay that the records were used as predictors.\n105: \"The age model\". There are two outliers in the radiocarbon ages right around the End of Sahara greening, how does\nthat affect the reconstruction?\n111: \"rainfall is the dominant control on changes\". Yes, but the temperature is also significant, can it still be derived?\n124 \"our reconstruction is robust to methodology.\" Yes. Should you at this moment discuss other sources of uncertainties?\n126: \"the modern value.\" Please refer to it.\n127: \"substantially low\". substantially lower\n134: \"major states\". Major time periods? Where some time periods have the same state?\n135: \"Between two adjacent states are abrupt and irreversible shifts\". How do you define abruptness? Why not use an\nobjective criterion (specific rate of change after applying a smoothing)? In particular, I do not see a particular abrupt change\nbetween LGM and HS1: What do you base the division on, and the exact dating? Also, what does irreversible mean here?\n143: Fig S5f shows an interesting trend analysis describing a specific sequence of events. Why not mention it here or later?\n145: \"can be contributed\". can be driven\n147: \"The global mean temperature during the LGM was 6.1 °C cooler\". The authors only cite one study with no uncertainty\nrange. You could add Shakun et al. 2012, Annan et al. 2022 and Osman et al. 2021 (among others), for an estimate of\naround 4 to 6 degree cooler.\n150-151: \"This unexpected finding suggests that global mean temperature change may not be the dominant control on the\nEASM rainfall during the LGM,\". Well, none of the monsoon regimes follows the GMT, but rather the obliquity, at timescale >\n10 kyr. So, that is not at all unexpected.\n157-158: \"We therefore attribute the similar rainfall amounts to the dynamic effect of nearly identical summer insolation\ndespite [...]\". You could add that the dynamical effect needs to compensate for the thermodynamic one. And wonder if there\nis something that constrains the dynamical effect to match the thermodynamic one ... (e.g. latitudinal gradients?)\n163-165: \"The AMOC has experienced tipping points for abrupt transition between strong states and weak states during the\nlast deglaciation, which can have strong impacts on the EASM rainfall\". You can cite specific papers that looked at the\nrelationship between AMOC and EASM\n172-173: \"supporting a southward migration in the ITCZ position during the HS1 and YD in\nresponse to the weakening of the AMOC.\" Again, cite the paper on the topic.\n179-180: \"disagree on the direction of rainfall changes during North Atlantic cold events\". Can you develop? For example,\nyou can introduce the question of the tripole pattern here.\n181: \"we compile a large number of hydroclimate records\". References?\n185: \"consistent with the majority of the available hydroclimate records\". But you did say previously that the\nproxy records\ndid not agree with each other. Why does your selection give a more homogeneous pattern?\n186: \"but do not support the existence of a “tripole” rainfall pattern.\" Can you further discuss the differences with previous\nsimulations? Why is this simulation different to the previous one, and be considered closer to the truth? For example, I find\nthe simulated timeseries presented in He et al. 2021 to be very close to your reconstruction. However, this study does\nsupport the tripole pattern.\n199: \"with the strongest cooling occurring in the northern high latitudes owing to amplification of cryosphere feedbacks\". The\npaper cited here does not refer to the End of the Sahara greening. You look for other papers, such as Muschitiello et al.\n2015.\n200: \"altered interhemispheric temperature gradient\" I am a bit confused as to why the cooling should be confined to the\nNorthern Hemisphere? I understand the argument from Park et al. 2019 concerning the change in orbital forcing, but how\ndoes it relate to the WAM waning?\n203-209: \"which closely tracks rainfall change in East Asia, [...]\". Yes, because of how westerlies impact/drive EASM. I like\nthe discussion, you can maybe cite the EAWM records. Can the abruptness/the bistable state of the EASM change actually\nbe related to westerlies?\n208-209: \"the snow-albedo feedbacks could have magnified its climatic response to Saharan\nvegetation change\". This is also something you can check in your simulations.\n212-216: This discussion of the model results and comparison with proxy records is very short and inconclusive.\n226: \"comprehensive paleoclimate data integration\". Again, I do not see what the authors mean by \"comprehensive\nintegration\". You are rather comparing the results of an ensemble of proxy records to model data.\n227: \"we have drawn the hitherto most complete picture\". I would be happy to consider that your set of proxy recordings\ngives a complete picture if I knew what that set was and how you had analysed it. But as far as the model perspective is\nconcerned, I don't consider that a single simulation, whatever its resolution, provides a complete picture.\n230: \"[...] has explicit implications for future change in the EASM\". It is not wrong, but I don't think that your study really\nbrought these implications forward (as a new finding), hence the suggestion to discuss tipping points within the EASM\ndomain.\n232-233: \"For instance, collapse of the AMOC to a weak state will lead to long-lasting and widespread extreme droughts in\nEast Asia.\" But does this relationship hold for our current and future climate state? This state is quite different from glacial\nconditions.\n237-238: \"responses of the EASM to a [...] a greening Sahara\". First, you haven't really discussed the effect of Sahara\ngreening, which was less abrupt than its termination. Second, what makes you think that the Sahara will become greener\nunder anthropogenic changes? Maybe say that the dependency between WAM and EASM should be further investigated in\nlight of the recent WAM variability, particularly in terms of de/afforestation policies in the Sahel region and aerosol\nemissions.\n339-340 \"the qualitative relationship between proxy and climate can be explained and the\nmechanisms have been relatively well studied.\" Please give some references here.\n492-493: \"As the model resolution is very high, there are significant reductions in\nclimate biases and model drifts.\" high is enough (there are versions with higher resolutions). And the paper shows a large\nprecipitation bias over South Asia, with no improvement with increased resolution. I do not know how this bias can affect\nyour results (Would probably need to look into the seasonality)\n497: \"the control experiment as initial conditions\". It is not clearly indicated, but I surmise that the boundary conditions are\nthat of the PI. How would different GHG, ice sheet, and sea level boundaries affect your simulation?\n499: \"In our high-resolution model simulation\" This seems the same set-up as in Jackson and Wood 2018. Could you clarify\nif these are indeed the same simulations?\n515: \"we performed two sets of mid-Holocene (MH) experiments\". Similarly, can you clarify if this is the same simulation as\nin Pausata et al. 2016?\n582: \"Brown stars (this study) and circles indicate proxy records\". Please indicate which proxy records you used, and how\nyou determined whether they showed abrupt drying or wetting. I am sure they are also proxy records that exhibit no abrupt\nchanges (e.g. the d18O record from the Dongge cave).\n588: \"Source data are provided with this paper.\" Does that mean you will publish the reconstruction (and the proxy records)\nalong with the paper? How about the model data?\n589 \"Code Availability\". It would also be interesting to have the code for the various reconstruction methods. Or which\npackage you use to do it.\nThe literature I mention in the comments:\nC. He et al. ,Hydroclimate footprint of pan-Asian monsoon water isotope during the last deglaciation.Sci.\nAdv.7,eabe2611(2021).DOI:10.1126/sciadv.abe2611\nMuschitiello, F., Zhang, Q., Sundqvist, H. S., Davies, F. J. & Renssen, H. Arctic climate response to the termination of the\nAfrican Humid Period. Quat. Sci. Rev. 125, 91–97 (2015).", "reviewer_3": "The findings of the submitted manuscript by Lu et al. is mainly relies on climate proxy data from a sediment section from\nQingtongyang maar lake in South China, which seems to cover the last 22,000 years. However, no age-depth model is\nprovided for the studied sediment core, which makes it impossible to evaluate the quality of the chronology of the climate\nreconstructions and their discussion in context of climatic tipping points.\nFurthermore, information on the applied rainfall reconstruction method based mainly on the concentrations of three chemical\nelements is insufficient. There is only visual but no statistical information that shows the performance/reliability of the applied\nreconstruction method using the modern reference datasets. How do you explain the statistically detected relationships\nbetween element concentrations and rainfall variations? Are the chemical element–rainfall relationships applicable\nuniversally or do they only exist in the studied maar lake? It is also not clear which climate data was exactly used. It seems\nthat also climate data from the 19th century CE was used. Why was it used and how reliable are such climate observations?"}, "author_response": "1 Response to Reviewer # Manuscript NCOMMS-24-28148\n4\n7\n9 paper, we have carefully revised our manuscript according to the reviewer’s comments and suggestions,\n10 and the changes are marked in red in the main text. Followings are the point-by-point replies.\n11\n16\n19 been extensively investigated for both the processes and the forcing mechanisms on different timescales\n20 based on multi-proxy analyses of the Chinese loess–paleosol sequences, stalagmites, and the other\n21 sedimentary archives. However, how the EASM behaves when global warming exceeds the specific\n22 threshold remains unknown due to the lack of quantitative records of the summer monsoon rainfall in the\n23 geological past. In this study, the authors present a robust reconstruction of the EASM rainfall changes in\n24 southern China over the past 22,000 years based on a multi-method, multi-proxy framework, and suggest\n25 that the EASM rainfall was characterized by abrupt shifts with a magnitude of ~500 mm (30%). Combined\n26 with climate simulations, moreover, the authors propose that these abrupt shifts in the EASM rainfall would\n27 be linked to tipping points in the Atlantic Meridional Overturning Circulation (AMOC) and the Saharan\n28 vegetation cover. The data and inferences presented in this study are of great significance and would\n29 contribute to a deep understanding of the response of the EASM to future tipping points. I recommend\n30 acceptance of this manuscript for publication in Nature Communications after minor revisions.\n32\n34 1. Lines 64–84: it is better to explain the physical link between the EASM and the AMOC / Saharan\n35 vegetation cover in this part rather than only to emphasize the importance of both in the Earth climate\n36 system.\n38 “Through its impact on the interhemispheric thermal gradients, the AMOC can alter the position and\n1\n39 intensity of low latitude rainfall systems, including the intertropical convergence zone (ITCZ) and the\n41\n42 “Recent studies suggest that the changes in Saharan vegetation cover and dust emissions can not only\n43 influence local climate in northern Africa but also the distant climates ranging from the tropics to the Arctic.\n45 such as the upper-level Rossby wave train and the Pacific Walker circulation.”\n46\n47 2. Lines 85–94: This part contains some of the results, discussion, and conclusions, and should be moved\n48 separately to the following sections ‘Results and discussion’ or ‘Conclusion’.\n50\n51 3. Lines 201–207: Did the cooling in the North Atlantic caused by the end of the Green Sahara have any\n52 impact on the AMOC and then on the EASM rainfall?\n54 “In addition, the negative NAO could weaken the wind-driven AMOC and reduce the poleward ocean heat\n55 transport in the North Atlantic, further reducing the interhemispheric temperature gradient and the EASM\n56 rainfall. However, multiple lines of evidence suggest that the basin-scale AMOC remained fairly stable\n57 throughout the Holocene. Therefore, the feedbacks associated with the Holocene AMOC likely play a minor\n58 role.”\n59\n60 4. Lines 215–216: What is your hypothesis to explain such a conflict between model simulations and proxy\n61 records?\n63 “Taken together, there are large discrepancies between proxy records and climate model simulation in\n64 reconstructing the hydroclimate of southern China during the last deglaciation, particularly for the iTRACE\n65 simulation running with the Community Earth System Model (CESM). We argue that these discrepancies\n66 may be attributed to the inability of climate models in correctly simulating the EASM rainfall pattern. It\n67 has been suggested that the coupled general circulation models (CGCMs) have large biases in simulating\n68 the mean state of the tropical Pacific with an excessive equatorial Pacific cold tongue when compared to\n69 observations. As a result, current CGCMs are unable to accurately reproduce the spatial pattern of EASM\n70 rainfall during different seasonal stages, with the rain belt shifting northward too quickly. For example, the\n71 position of the spring and pre-mei-yu rain belts are shifted northward compared to rainfall observations,\n72 resulting in deficient rainfall in southern China.”\n73\n74 “However, the simulated increases in annual and summer rainfalls in southern China contrast with the\n75 available proxy records there (Supplementary Table 9, Fig. 9). Similar to the last deglaciation, we argue\n76 that the biases in climate models are likely responsible for this model-proxy discrepancy.”\n77\n78 5. Lines 237–239: I am confused about ‘a weakening AMOC’ and ‘a greening Sahara’ mentioned together\n79 with each other here.\n81\n2\n83 1. Line 97: Change ‘sensitive’ to ‘closely related’.\n85\n86 2. Line 110: Delete ‘consistently’.\n88\n89 3. Line 111: Change ‘on’ to ‘of’.\n91\n92 4. Line 132: Change ‘between states’ to ‘between the states’.\n94\n95 5. Line 132: Change ‘According to’ to ‘Based on’.\n97\n98 6. Line 141: Change ‘has’ to ‘had’.\n100\n101 7. Line 142: Change ‘vast’ to ‘the whole’.\n103\n104 8. Line 142: Move ‘the’ before ‘observed rapid shifts’ to before ‘rainfall regime’.\n106\n107 9. Lines 145–146: Change ‘the dynamic component’ to ‘the former’ and ‘the thermodynamic component’\n108 to ‘the latter’.\n110\n111 10. Line 151: Change ‘on’ to ‘of’.\n113\n114 11. Figure 1: Show the latitude for the JJA insolation in Figure 1.\n116\n117 12. Line 163: Change ‘for’ to ‘of’.\n119\n120 13. Line 171: Insert ‘two’ before ‘hemispheres’.\n122\n123 14. Lines 173–177: It may be better to remove these words about the relation between the ITCZ and the\n3\n127\n128 15. Line 224: Change ‘between states’ to ‘between the states’.\n130\n134 calibrated to estimate a timeseries of precipitation, which the authors use as a proxy for East Asian Monsoon\n135 intensity. The authors particularly focus on the abrupt changes that the time series exhibits. They link these\n136 changes to AMOC variability during the deglaciation, and to the end of the green Sahara during the mid-\n138\n140 analysis perspective, the post-processing and calibration methods used to retrieve the precipitation signal\n141 from the proxy records seem sound. In particular, the large number of methods used to compute the\n142 reconstruction is impressive and valuable. This increases the confidence in the reconstruction. This one also\n143 exhibits striking features that the authors discuss.\n145\n147 these model outputs to be quite limited, and their added value insufficient. I provide suggestions regarding\n148 that matter below.\n151 results.\n152\n154 feature, that has been discussed in previous studies both from a proxy and model perspective. The\n155 precipitation reconstruction from this study adds another evidence to it.\n157\n159\n160 - Based on your narrative \"We provide an ensemble reconstruction of the EASM based on a multi-proxy\n161 and multi-method framework\", I thought you were post-processing an ensemble of already published proxy\n162 records from the entire EASM domain, and that your paper focused on the methods to derive the EASM\n163 reconstruction from these. I was therefore very confused when I started to read the result section. I suggest\n164 that you make clear in the abstract and the introduction that you are publishing new proxy records from a\n165 single core/location and that you focus on retrieving the precipitation signal from these records.\n167\n4\n169 rainfall since the Last Glacial Maximum (LGM) using nine statistical and machine learning methods based\n170 on multi-proxy records from a maar lake in southern China.”\n171\n173 southern China and use it to detect tipping points in the EASM rainfall since the LGM and assess the\n174 cascading impacts of climate tipping points on the EASM rainfall. The first goal is achieved by applying\n175 nine statistical and machine learning methods to retrieve rainfall signal from multi-proxy records of a maar\n176 lake in southern China.”\n177\n180 based on empirical orthogonal function (EOF) decomposition of the GPCC rainfall dataset. Then we discuss\n181 in the main text:\n182 “Our analysis indicates that the rainfall variability in the Leizhou Peninsula is representative of rainfall\n183 variability in southern China (Supplementary Fig. 1).”\n184\n186 climatic zones in East Asia using rigorous criteria. We have added this part to Methods section and discussed\n187 it in the main text.\n188\n189 - The evaluation of the uncertainty related to the specific method used is impressive. The discussion on the\n190 age model is also quite clear.\n192\n194 period used?\n196 very small catchment and lacks river input, where regional rainfall influences the redox conditions and\n197 salinity of the lake, thereby affecting the concentration changes of elements. For other lakes where sediment\n198 transport and deposition are controlled by hydrodynamic conditions, the interpretation may differ. So, these\n199 proxy records should not be used blindly, and careful analyses are needed. We have added more discussion\n200 on the physical and chemical mechanisms of our proxy records to the Supplementary Discussion.\n201\n203 period as one of the ensemble parameters (see Methods):\n204 “To overcome bias in the choice of calibration period whilst retaining the autocorrelation structure, all\n205 possible combinations of 5-year blocks and additional successive years with a length of 33/38 years were\n206 used for training, while the remaining 20/15 years were used for validation.”\n207\n210 with other proxy records from the same region, including the Huguangyan and Tianyang maar lakes\n211 (Supplementary Fig. 11). They are similar. We also compare our rainfall reconstruction with the speleothem\n212 δ18O records from Dongge and Hulu caves in southern China (Supplementary Fig. 11). After accounting\n5\n213 for the effect of global ice volume change, the speleothem δ18O values during the LGM are slightly lighter\n214 than those during the late Holocene. This is consistent with our reconstruction indicating similar or higher\n215 rainfall amounts during the LGM than the late Holocene. Now we have discussed these in the main text.\n216\n219 biased toward zero, which are directly related to the correlation between proxies and the climate target. The\n220 biases are inevitable, but our use of multi-proxy and multi-method can reduce such biases.\n221\n224 reconstructions of the past two millennia. Reviews of Geophysics 55, 40-96 (2017).\n227\n229 precipitation influences the colour of the core and the concentration of iron, manganese, and calcium.\n231 the physical and chemical mechanisms of our proxy records, explaining how precipitation influences the\n232 colour of the core and the concentration of iron, manganese, and calcium.\n233\n234 - I find the comparison to other proxy records a bit lacking. For example, I did not find any information on\n235 the records used in Fig. 4, and therefore cannot comment on it.\n238 information of the collected hydroclimate records in the Supplementary Tables 5 to 9.\n239\n241 forcing. However, the authors only mention once the \"tripole\" without any context to it, and mostly dismiss\n242 it as a possible feature. However, the tripole pattern is quite clear from Fig.4 b (although maybe a bit too\n244 that published the record you show in Fig. S4c (Shouldn't you flip the y-axis of this one, so that the top\n245 correspond to wetter and the bottom to drier as for the other records?). I strongly suggest the author to\n246 further discuss this question of heterogeneity of the EASM response.\n249 deglaciation remains ambiguous due to the lack of a systematic analysis of both proxy records and climate\n250 simulations. For instance, a recent study based on a trace-element record from Haozhu Cave and a hosing\n251 experiment under preindustrial conditions suggests that the central eastern China was wet (dry) during the\n252 times of a weak (strong) AMOC and the rainfall anomaly in East Asia exhibited a “tripole” mode 1. The\n255 proxy records from the same region 2-6, including two trace-element records from Heshang Cave 2 and Hulu\n6\n257 decrease rather than increase in central eastern China in response to the weakening of the AMOC 7, which\n258 highlights the need for multi-model comparison. Previous studies mainly focused on the response of EASM\n259 rainfall to the weakening of the AMOC but overlooked the impact of AMOC strengthening. However, the\n260 responses of climate to the weakening and strengthening of the AMOC may not be a reversed process. In\n261 particular, abrupt shifts in the AMOC during the last deglaciation were accompanied by variations in other\n262 climatic forcings such as solar insolation, ice sheet, and greenhouse gas concentrations. Therefore, it is\n263 necessary to analyze the individual characteristics of EASM rainfall response during both the transitions\n264 into and out of HS1 and YD, which has not yet been done. Here we collect a large number of hydroclimate\n265 records in East Asia based on rigorous criteria (Methods, Supplementary Tables 5 to 8) and compare them\n266 with a state-of-the-science isotope-enabled transient climate experiment (iTRACE) 8 and an earlier version\n267 without isotope module (TRACE-21ka) 9 (Methods). Although continuous records with good age control\n268 and high temporal resolution are scarce for the last deglaciation, the general trend of hydroclimate change\n269 can be roughly inferred from these proxy fragments.\n271 1 Zhang, H. et al. East Asian hydroclimate modulated by the position of the westerlies during Termination\n273 2 Hu, C. et al. Adsorbed silica in stalagmite carbonate and its relationship to past rainfall. Geochimica et\n275 3 Wang, Q. et al. Orbital- and millennial-scale hydroclimate changes in central China during the last\n276 glacial period. Quaternary Science Reviews 337, 108802 (2024).\n277 4 Han, Z. et al. Extreme monsoon aridity episodes recorded in South China during Heinrich events.\n279 5 Zhou, X. et al. Catastrophic drought in East Asian monsoon region during Heinrich event 1. Quaternary\n281 6 Zhang, W., Wang, Y., Wu, J. & Duan, F. Last deglacial climate variations inferred from trace elements\n282 in a stalagmite from Hulu Cave, Nanjing. Quaternary Sciences 34, 1227-1237 (2014).\n283 7 Sun, Y. et al. Influence of Atlantic meridional overturning circulation on the East Asian winter monsoon.\n285 8 He, C. et al. Hydroclimate footprint of pan-Asian monsoon water isotope during the last deglaciation.\n287 9 Liu, Z. et al. Transient simulation of last deglaciation with a new mechanism for Bølling-Allerød\n288 warming. Science 325, 310-314 (2009).\n289\n291 characterizing the abrupt changes.\n293\n295 “Based on the Bayesian age modeling of our 14C dates from the Qingtongyang Maar Lake (Methods), the\n296 timings of these abrupt shifts are estimated to be centered at 18860 ± 310, 14640 ± 400, 12870 ± 290, 11360\n297 ± 170, and 4370 ± 60 yr BP (Supplementary Table 4).”\n298\n300 “These two independent methods are consistent in revealing the timings and structures of abrupt shifts in\n7\n301 the EASM rainfall in most case, but the Bayesian method can provide more details about the abrupt shift\n302 during the mid- to late Holocene transition. Based on the Bayesian change-point analysis, the mid- to late\n304 60 yr BP and a subsequent less abrupt but persistent decrease in rainfall with an amplitude of 416 mm until\n305 3290 ± 70 yr BP (Supplementary Fig. 10 and Supplementary Table 4). Therefore, our reconstruction reveals\n306 a two-step abrupt shift in the EASM rainfall during the mid- to late Holocene transition, which has\n307 implications for better understanding of the forcing mechanisms and cultural consequence of the well-\n308 known 4.2-ka event.”\n309\n311 from Huguangyan Maar and Tianyang Maar, or even the Dongge cave. This could help to better understand\n312 regional features (or proxy limitations).\n315 and discussion in the main text).\n316\n317 - As it stands, I do not see the added value of the model data. I see three courses of action:\n318 - Entirely remove the model data, to better focus on the proxy records.\n319 - To compare the model data to previous studies, in order to bring a multi-model (or at least multi-simulation)\n320 approach to the analysis. The authors could then discuss the reliability/robustness of their model data and\n321 potentially explain why we should better trust their higher resolution simulation.\n322 - The authors could then further investigate how AMOC / End of Sahara greening influences the EASM in\n323 their simulations. This could also be done with sensitivity experiments to further understand which\n324 component of the forcing influence the EASM, and how.\n326 in the AMOC, we compare proxy records with the iTRACE and TRACE-21ka simulations. For the end of\n327 the Green Sahara, we compare with climate simulations using the EC-Earth, CCSM4, iCESM, and GISS\n328 models. We also analyze temperature and sea level pressure in simulations to see how AMOC / End of\n330\n332 would the AMOC impact on EASM be, if LGM boundary conditions were used?\n334 experiment based on PI boundary conditions. The different rainfall response in East Asia under different\n335 boundary conditions has also been mentioned in the main text.\n336\n337 - I feel the focus on \"tipping points\" is either too extensive or too general. AMOC or End of Sahara greening\n338 are well known tipping points, and their impact on the EASM has already been discussed. In its current\n339 state, I do not see much added value from the paper that justifies this focus. However, I do see an opportunity,\n340 if the authors were to further discuss tipping points within the EASM domain, as part of a tipping point\n341 cascade. Neither Wang et al. 2023 nor Wunderling et al. 2024 recent review on tipping points mention the\n343 proof of the opposite. Simulations can be used to investigate this tipping point, and answer for example the\n344 question: Is the EASM response linear to the freshwater hosing / AMOC slow down?\n8\n347\n349 publication of a new record / reconstruction, as a comparison between proxy records and simulations, or as\n350 a discussion of tipping points in the EASM domain.\n352\n354\n355 28: \"30%\". is this uncertainty or change compared to the mean value?\n357 period 1961–1990 CE. In this version of our paper, this value has been revised using more objective method,\n358 i.e., regime shift analysis, and the median value was estimated to be 387 ± 73 mm (24 ± 5 %).\n359\n360 30-31: \"the collapse of the Atlantic meridional overturning circulation [...] caused drying\": The re-\n361 establishment of the AMOC was seemingly as abrupt with consequences on the EASM. Why not discuss\n362 that too?\n364\n365 31-33: \"widespread drying\", \"drying in northern East Asia and southwestern China and possible wetting in\n366 southeastern China\". These are the raw results from the simulations. Either clarify that these are your model\n367 results, or further discuss in the text why you think these results are better compared to other simulations.\n369 between proxy records and Green Sahara simulations in the main text.\n370\n371 33-34: \"Our findings provide crucial evidence for projecting the response of East Asian hydroclimate to\n372 future tipping points.\" The crucial evidence is not so evident to me. The current climate conditions are quite\n373 different to the AMOC collapse during the deglaciation and the end of the Sahara greening.\n375 and its coupling with other tipping elements.”\n376\n377 56-57: \"However, the interpretation of speleothem δ18O records in East Asia remains highly debated,\n378 excluding it as a local rainfall proxy\". You could also cite He et al. 2021 here.It quite clearly shows this.\n380\n383 “Because the climatic interpretation of speleothem δ18O records in East Asia remains controversial, we only\n384 include the speleothem δ18O records that have been linked to hydroclimate change based on correlation\n385 with instrumental records and/or independent rainfall-related proxy records such as δ13C records.”\n386\n9\n389 “To help for better understanding the climatic significance of speleothem δ18O records in East Asia, we\n390 compare our rainfall reconstruction with the speleothem δ18O records from Dongge and Hulu caves in\n391 southern China (Supplementary Fig. 11). The strong similarity between them during the last deglaciation\n392 indicates that the speleothem δ18O records in East Asia might be indeed related to hydroclimate changes\n393 during this period. After accounting for the effect of global ice volume change, the speleothem δ18O values\n394 during the LGM are slightly lighter than those during the late Holocene. This is consistent with our\n395 reconstruction indicating higher rainfall amounts during the LGM than the late Holocene. However, the\n396 changes in speleothem δ18O records during the Holocene are much more gradual, indicating decoupled\n397 changes between speleothem δ18O records and East Asian hydroclimate. Previously, the speleothem δ18O\n398 records in East Asia were interpreted as monsoon wind intensity, which is positively (negatively) correlated\n399 with monsoon rainfall in northern (southern) China. However, the decoupled changes between speleothem\n400 δ18O records and our rainfall reconstruction as well as other hydroclimate records in East Asia indicates\n401 that speleothem δ18O records are not a pure signal of monsoon rainfall.”\n402\n403 63: \"particularly in southern China.\" Why there? How important is this area in the EASM?\n405 “Absent from these debates is an independent quantitative reconstruction of past rainfall changes in the\n407 annual rainfall exceeds 800 mm and rainfall variability is the largest (Fig. 1).”\n408\n409 65: \"application of [...]\". Add here that you are using records from one lake.\n411\n412 66: \"unprecedented accuracy\". How do you know the accuracy?\n414\n415 65: \"to assess the impacts of climate tipping points on the EASM rainfall\". The phrasing of the goal is a bit\n416 strange. You found abrupt changes in your timeseries, and you want to understand their cause. Comparison\n417 with other proxies helps you to formulate a hypothesis, that you then test with climate simulations.\n419 “In this study, we aim to provide a robust reconstruction of the EASM rainfall in southern China and use it\n420 to detect tipping points in the EASM rainfall since the LGM and assess the cascading impacts of climate\n421 tipping points on the EASM rainfall.”\n422\n423 67: \"integrating\". You are just comparing the two. Integrating sounds like data assimilation, which you do\n424 not do.\n426\n427 73: \"AMOC will collapse\". AMOC may collapse? That reflects better the uncertainty/controversy.\n429\n430 74: \"climatic impact of AMOC collapse remains unclear.\" It is quite clear from a paleo perspective. What\n431 is unknown is what will happen in the context of anthropogenic changes.\n10\n433 “However, the climatic impact of AMOC collapse in the context of anthropogenic warming remains\n434 unclear.”\n435\n436 82-84: These claims are far too generic. Maybe a better context would be the recent shifts / variability of\n437 the West African monsoon, their relationship with vegetation in the Sahel, and potential teleconnection with\n438 the EASM?\n440 “Given that the Sahara/Sahel vegetation and the West African monsoon have been identified as a tipping\n441 element that may experience abrupt shifts under future anthropogenic warming, investigating the effect of\n443 vegetation feedback on future climate.”\n444\n445 86: \"between states\". between these states\n447\n448 85-96: Is this usual to have a paragraph summarizing the results at the end of the introduction?\n450\n451 101: \"annual rainfall and summer rainfall are equivalent\". They are only similar.\n453 “Because annual rainfall in East Asia is dominated by summer rainfall and their changes are highly\n454 correlated with similar amplitudes (Supplementary Fig. 2), we only reconstruct the annual rainfall.”\n455\n456 102-105: First sentence on the fact that the core was recovered and the records measured for this study.\n459\n460 105: \"The age model\". There are two outliers in the radiocarbon ages right around the End of Sahara\n461 greening, how does that affect the reconstruction?\n464 and extrapolation of the median values of sixteen AMS (accelerator mass spectrometry) 14C dates with two\n465 outliers excluded.”\n466\n467 111: \"rainfall is the dominant control on changes\". Yes, but the temperature is also significant, can it still\n468 be derived?\n470 the correlation with temperature are not significant (Supplementary Table 2), so it cannot be derived from\n471 these proxies faithfully.\n472\n473 124 \"our reconstruction is robust to methodology.\" Yes. Should you at this moment discuss other sources\n474 of uncertainties?\n476 “The rainfall reconstructions obtained by the nine methods exhibit similar oscillations on the decadal to\n11\n477 millennial timescales over the past 22,000 years (Fig. 3, Supplementary Fig. 8 and 9), indicating that our\n478 reconstruction is robust to methodology. This is the first quantitative reconstruction of the EASM rainfall\n479 that considers as many uncertainties as possible (Methods). The influence of age model uncertainty has also\n480 been considered. The age model uncertainty mainly affects the weights of proxy records and the timings of\n481 abrupt climate changes (Supplementary Fig. 4). Therefore, we addressed this uncertainty by randomly\n482 adjusting proxy weights during reconstruction and applying Bayesian age modeling to ascertain the timings\n483 of abrupt shifts (Methods). Even after accounting for these uncertainties, the amplitudes of rainfall changes\n484 in our reconstructions during the past 22,000 years are still larger than the total uncertainties (Fig. 3,\n486 than stochastic process.”\n487\n488 126: \"the modern value.\" Please refer to it.\n490\n491 127: \"substantially low\". substantially lower\n493\n494 134: \"major states\". Major time periods? Where some time periods have the same state?\n496 same climate state. Now we use “climate states” to make it clear:\n497 “The evolution of the EASM rainfall over the past 22,000 years can be divided into six major climate states\n498 (or periods) that have different characteristics and internal dynamics, aligning with the LGM, HS1, BA,\n500\n501 135: \"Between two adjacent states are abrupt and irreversible shifts\". How do you define abruptness?\n503 change.\n504\n507 proxy records, including the regime shift analysis, Bayesian change-point analysis, and the rate of change\n508 (see Methods). However, the rate of change was not applied to detect abrupt changes in our rainfall\n509 reconstruction because the changes in ultra-high-resolution records are very abrupt and the selection of\n510 smoothing parameters such as window length is quite subjective and cannot adapt universally to all\n511 situations of different abrupt events.\n512\n514 on, and the exact dating?\n516 of 142 mm, possibly because the rainfall amount during the LGM had been low. Nonetheless, this rainfall\n517 change is statistically significant at the 95% confidence level based on the two-tailed two-sample t-test. The\n518 regime shift analysis and Bayesian change-point analysis were used to detect this abrupt change from LGM\n519 to HS1 (Supplementary Fig. 10). Bayesian age modeling was employed to assess age uncertainties\n520 (Supplementary Table 4).\n12\n521\n524 shift, but remains in a new state for a long time (>1000 years in this study) until its next shift.\n525\n526 143: Fig S5f shows an interesting trend analysis describing a specific sequence of events. Why not mention\n527 it here or later?\n529 “Based on the Bayesian change-point analysis, the mid- to late Holocene transition is characterized by an\n530 abrupt drop of rainfall with an amplitude of ~140 mm at 4370 ± 60 yr BP and a subsequent less abrupt but\n531 persistent decrease in rainfall with an amplitude of 416 mm until 3290 ± 70 yr BP (Supplementary Fig. 10\n532 and Supplementary Table 4). Therefore, our reconstruction reveals a two-step abrupt shift in the EASM\n533 rainfall during the mid- to late Holocene transition, which has implications for better understanding of the\n534 forcing mechanisms and cultural consequence of the well-known 4.2-ka event.”\n535\n536 145: \"can be contributed\". can be driven\n538\n539 147: \"The global mean temperature during the LGM was 6.1 °C cooler\". The authors only cite one study\n540 with no uncertainty range. You could add Shakun et al. 2012, Annan et al. 2022 and Osman et al. 2021\n541 (among others), for an estimate of around 4 to 6 degree cooler.\n543 “The global mean temperature during the LGM was 6.1°C (95% confidence interval: 5.7 to 6.5 °C) cooler\n544 than that during the late Holocene.”\n545\n547 studies indicate the difference between LGM and preindustrial.\n548\n549 150-151:"}, {"round": 2, "reviewer_comments": {"reviewer_1": "All my concerns have been addressed in the revisions, and I have no further comments on the revised manuscript. I also looked into the reviewer reports provided by the other two reviewers, and found that the authors have made careful and acceptable replies to all the comments. I think that the revised manuscript has been greatly improved. Therefore I recommend acceptance of the manuscript in the present form for publication in Nature Communications.", "reviewer_2": "I am positively impressed by this version of the manuscript. The authors have performed a significant amount of work, and the improvements of the manuscript have overcome my expectations. All the reviewer's comments have been successfully addressed. Furthermore, the scope is now clear, and the readability satisfactory. I have enjoyed reading the manuscript. The current state of the manuscript is suitable for publication.\nI just have a few very minor comments:\nL.74 \"in the Leizhou Peninsula\": it is here only implicit that the lake is located there. I would suggest \"at the lake location\" or \" in the Leizhou Peninsula, where the lake is located\".\nL. 174 \"is robust to methodology.\": \"is robust to the methodology used\"\nL. 208 \"the rainfall amount during the LGM had been low\": \"the rainfall amount during the LGM was low\".\nL. 209 \"with a median of 387\": I don't see the added value of this information.\nL.229 \"To help for better understanding\": \"To help better understand\"\nL.258-259 \"The dynamic effect of summer insolation could compensate for the thermodynamic effect of global mean temperature through changes in the land-sea thermal contrast and latitudinal temperature gradient.\": I find the sentence a bit disconnected from the previous one. I suggest a slight reformulation.\n\"Yet, the thermodynamic effect of global mean temperature needs to be compensated by dynamic effects, maybe related to changes in the land-sea thermal contrast and latitudinal temperature gradient.\"\nL. 271 \"between two hemispheres\": \"between the two hemispheres\"\nL. 391 \"Therefore, it is convinced\": \"Therefore, it appears\" or \"These evidences suggest\"\nL. 441 \"In consistent with\": \"Consistent with\"\nL. 861 \"in the latitude of 20°N\": \"at 20°N\""}, "author_response": "Response to Reviewer # Manuscript NCOMMS-24-28148A\nTipping point-induced abrupt shifts in East Asian hydroclimate since\nthe Last Glacial Maximum\nFuzhi Lu, Huayu Lu*, Yao Gu, Pengyu Lin, Zhengyao Lu, Qiong Zhang, Hongyan\nZhang, Fan Yang, Xiaoyi Dong, Shuangwen Yi, Deliang Chen, Francesco S. R. Pausata,\nMaya Ben-Yami, Jennifer V. Mecking\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nAll my concerns have been addressed in the revisions, and I have no further comments\non the revised manuscript. I also looked into the reviewer reports provided by the other\ntwo reviewers, and found that the authors have made careful and acceptable replies to\nall the comments. I think that the revised manuscript has been greatly improved.\nTherefore, I recommend acceptance of the manuscript in the present form for\npublication in Nature Communications.\nReply: Thank you for your careful review and supportive comments.\nReviewer #2 (Remarks to the Author):\nI am positively impressed by this version of the manuscript. The authors have\nperformed a significant amount of work, and the improvements of the manuscript have\novercome my expectations. All the reviewer's comments have been successfully\naddressed. Furthermore, the scope is now clear, and the readability satisfactory. I have\nenjoyed reading the manuscript. The current state of the manuscript is suitable for\npublication.\nReply: Thank you for your careful review and supportive comments.\nI just have a few very minor comments:\nL.74 \"in the Leizhou Peninsula\": it is here only implicit that the lake is located there. I\nwould suggest \"at the lake location\" or \" in the Leizhou Peninsula, where the lake is\nlocated\".\nReply: Done.\nL. 174 \"is robust to methodology.\": \"is robust to the methodology used\"\nReply: Done.\nL. 208 \"the rainfall amount during the LGM had been low\": \"the rainfall amount during\nthe LGM was low\".\nReply: Done.\nL. 209 \"with a median of 387\": I don't see the added value of this information.\nReply: Deleted.\nL.229 \"To help for better understanding\": \"To help better understand\"\nReply: Done.\nL.258-259 \"The dynamic effect of summer insolation could compensate for the\nthermodynamic effect of global mean temperature through changes in the land-sea\nthermal contrast and latitudinal temperature gradient.\": I find the sentence a bit\ndisconnected from the previous one. I suggest a slight reformulation.\n\"Yet, the thermodynamic effect of global mean temperature needs to be compensated\nby dynamic effects, maybe related to changes in the land-sea thermal contrast and\nlatitudinal temperature gradient.\"\nReply: We merge these two sentences to one sentence:\n“We therefore attribute the similar rainfall amounts to the dynamic effect of nearly\nidentical summer insolation, which could compensate for the thermodynamic effect of\nglobal mean temperature through changes in the land-sea thermal contrast and\nlatitudinal temperature gradient.”\nL. 271 \"between two hemispheres\": \"between the two hemispheres\"\nReply: Done.\nL. 391 \"Therefore, it is convinced\": \"Therefore, it appears\" or \"These evidences\nsuggest\"\nReply: Done.\nL. 441 \"In consistent with\": \"Consistent with\"\nReply: Done.\nL. 861 \"in the latitude of 20°N\": \"at 20°N\"\nReply: Done.\nAgain, thank you very much. We hope that we have addressed all you concerns."}]

Decision letter