100-kyr climate cycles caused by 2.4-Myr eccentricity-modulated carbon cycles

Domain: earth_science
Published: 2025-08-28
Editorial outcome: major_revision

Abstract
Earth’s climate has been dominated by ~100-kyr glacial cycles over the past ~800 ka, yet the mechanism remains debated. Here, we present correlation analyses of spectral power ratios of global records spanning the past 2.7 Ma, revealing a persistent anticorrelation between ~21-kyr and ~100-kyr power ratios, but no significant relationship between ~41-kyr and ~100-kyr power ratios. This suggests that ~100-kyr climate cycles are more related to eccentricity-modulated precession than to obliquity. Phase analyses of benthic δ18O/ice volume and δ13C (carbon cycle) since Antarctic glaciation onset (~34 Ma) show that strong ~100-kyr cycles emerged only when these proxies were phase-coupled. Such coupling recurred at ~2.4-Myr eccentricity maxima during the unipolar regime (before 7.5 Ma) and minima during the bipolar regime (after 4 Ma), explaining the persistent ~21-kyr/~100-kyr anticorrelation because eccentricity modulates precession amplitude. We propose that internal carbon cycle dynamics and ~2.4-Myr eccentricity-modulated δ¹⁸O/ice volume–δ¹³C coupling amplified ~100-kyr climate cycles not only over the past ~800 ka but since 34 Ma. Given that eccentricity will remain low for the next 400 kyr, ~100-kyr periodicities may continue to dominate future climate variability, assuming Earth remains in a bipolar regime.

Introduction
As global cooling progressed through the Cenozoic, Arctic ice sheets began developing in the late Miocene, and large Northern Hemisphere ice sheets first emerged at the onset of the Pleistocene (~2.6-million-year age; Ma), marking Earth’s shift from a unipolar to a bipolar ice-sheet regime. Since then, Earth’s climate has been characterized by periodical glacial-interglacial cycles. The alignment of glacial cycles with astronomical periodicities leads to linking these climate cycles to variations in Earth’s orbital geometry and axial inclination that cause the insolation received on Earth to fluctuate at the frequency of obliquity (~1/41 kyr−1), precession (~1/19 and 1/23 kyr−1), and eccentricity (~1/100 and ~1/400 kyr−1). Obliquity alters the latitudinal distribution of insolation and seasonal contrast; precession affects the timing and intensity of seasons relative to Earth’s position in its elliptical orbit; and eccentricity influences the Earth–Sun distance, modulating the total insolation received by Earth. However, during the period after the Mid-Pleistocene Transition (MPT) (~800–1200 ka), the glacial-interglacial cycles display a dominant ~100-kyr cyclicities, no astronomical parameter could directly produce such strong ~100-kyr oscillations. This is because precession and obliquity cannot directly generate periodicities at ~100-kyr band, and eccentricity’s effect on total insolation is too weak to explain the strong ~100-kyr glacial power, despite eccentricity can produce ~100-kyr cycles. Amplitude mismatch between ~100 kyr cycles in eccentricity and climate, and the absence of the 405-kyr power in benthic foraminiferal oxygen isotope (δ18O), complicate the understanding of the ~100-kyr glacial cycles.

Before the MPT, glacial cycles were directly attributed to obliquity due to dominant ~40-kyr cyclicities, though interpretations for the dominance remain diverse. This idea was extended to explain the post-MPT ~100-kyr glacial cycles as a quantum of 2–3 obliquity cycles. Indeed, convincing evidence confirms obliquity’s role in late Pleistocene glacial terminations. However, the precession and obliquity phases of glacial terminations based on 230Th-dated speleothem chronology reveal that late Pleistocene glacial terminations have a closer relationship with precession than obliquity. This finding supports interpreting ~100-kyr cycles as a quantum of 4–6 precession cycles or a result of a beat between the main precessional frequencies. The research on the morphology of deglaciation/termination of glacial cycles further unveils that precession has more influence on deglacial onset, while obliquity is more important for attainment of peak interglacial conditions and glacial inception.

Despite these phase analyses confirm the role of astronomical forcing in driving glacial cycles, the mechanism behind dominant ~100-kyr power remains unclear. A comparison based on spectral analysis suggests that variations in inclination of Earth’s orbital plane could explain late Pleistocene ~100-kyr cycles. However, this is contested as an artifact of spectral methods, because unlike inclination’s singular ~100-kyr peak, glacial cycles exhibit variable ~80–120 kyr coeval periodicities.

In contrast to the external astronomical origin, late Pleistocene ~100-kyr glacial cycles were thought as the self-sustained oscillation driven by internal climate feedbacks. However, the solely internal mechanisms produce arbitrary phase, conflicting with the observed phase-locking between ~100-kyr cycles of eccentricity and climate. Thus, the plausible explanation for the late-Pleistocene ~100-kyr glacial cycle lies on the models that express these cycles as a nonlinear interaction between orbitally forced response (precession and obliquity) and the internal dynamics of atmosphere, oceans, ice sheets and lithosphere. This nonlinearity invokes a set of internal feedbacks to channel some of the system’s potential power into the ~100-kyr band. It is accepted that the ~100-kyr glacial cycles were phase-locked to eccentricity over the last 1.2 Ma, thus, the key to decoding the strong ~100-kyr glacial cycles depends on finding out the power source that amplifies the ~100-kyr spectral power. Both nonlinear response to precession and obliquity can produce ~100-kyr cycles, potentially amplifying ~100-kyr glacial power. Thus, confirming a potential link between the strong ~100-kyr glacial power and eccentricity-modulated precession or obliquity is pivotal to understanding glacial cycles. In this regard, analyzing the spectral power ratio relationship among the ~100-kyr cycles, precession, and obliquity components could provide clues.

The spectral power ratio relationship solely based on benthic δ18O lack statistical significance because the reverse relationship between the ~21-kyr and ~100-kyr power ratios in δ18O is masked by a nearly antiphase relationship between the ~41-kyr and ~100-kyr power ratio over the last 800 ka, which obscures the roles of precession and obliquity in amplifying ~100-kyr power. Nevertheless, as mentioned above, given that the ~100-kyr glacial cycles were closely related to internal processes, the 100-kyr cycles across climate components should be highly correlated with each other through a series of positive and negative feedbacks, despite potential time lags or leads between the forcing and the climate feedback. For example, the expansion and retreat of Pleistocene ice sheets could perturbate sea level and atmospheric-oceanic circulations via meltwater pulses and reshape Earth’s energy budget through albedo feedback. Consequently, the spectral power ratio relationship between precession/obliquity and the ~100-kyr components, preserved in ice volume or δ18O, should be projected onto other climate archives that are significantly affected by ice sheets. We therefore examine the spectral power ratio relationship among various paleoclimate proxies, including benthic δ18O, to explore the ~100-kyr glacial cycles. The spectral power ratio serves as a proxy for estimating spectral strength or variance within a specific frequency band (see Methods).

Given the shift from dominant ~40-kyr to ~100-kyr glacial cycles after the MPT and the current availability of high-resolution paleoclimate data, here, we statistically analyze spectral power ratio relationship among the ~21-kyr, ~41-kyr, and ~100-kyr components within datasets across 57 globally distributed paleoclimate archives (including ocean, loess, ice core, and other archives) from 33 sites. These continuous records span pre-MPT (1200–2700 ka), MPT (800–1200 ka), and post-MPT (0–800 ka) intervals, all with robust age models and a high temporal resolution of at least a 2-kyr sampling rate. The ~100-kyr power ratios display persistent negative correlation with ~21-kyr power ratios across these three intervals, while no relationship with ~41-kyr power ratios. Notably, strong/amplified ~100-kyr climate cycles also occurred before the Pleistocene, such as during the Pliocene, Miocene, and Oligocene. To test whether these earlier ~100-kyr climate cycles share a similar mechanism with those of the past 800 ka, we extend the analysis to the Eocene–Oligocene Transition (EOT; ~34 Ma), when the Antarctic ice sheet began to develop, by exploring the phase relationship between ice volume/benthic δ18O and benthic δ13C. It turns out that ~100-kyr climate cycles were periodically amplified, coinciding with δ18O/ice volume–δ13C phase coupling, which recurred at ~2.4-Myr eccentricity maxima (before 7.5 Ma) and minima (after 4 Ma).

Results
The power decomposition analysis (PDA) (see Methods) was used to quantitatively estimate the spectral strength at frequencies corresponding to precession, obliquity, and ~100-kyr bands in paleoclimate records. Figure 2 presents the summarized PDA results for the collected 57 datasets, each based on their original age models, while detailed PDA results are shown in Figs. S2–S5. For comparisons, the overall power ratios of the ~100-kyr, ~41-kyr, and ~21-kyr cycles are shown using box whisker plots (Fig. 2A–C). To ensure significance, only PDA results with spectral power exceeding the 90% significance level at ~1/21 kyr−1, ~1/41 kyr−1, and ~1/100 kyr−1 were included. Detailed spectral estimations of these datasets are provided in Figs. S6–S9.

Post-MPT ~100-kyr power ratios exceed those of the pre-MPT interval in all datasets, while the opposite trend occurs for ~41-kyr power ratios. These shifts align with the transition from dominant ~41-kyr to ~100-kyr glacial cycles, confirming glacial influence on these records (Fig. 2A, C). Most paleoclimate records over the past ~800 ka show an inverse correlation between the ~100-kyr and ~21-kyr power ratios (Fig. 2A, B), a relationship that persists during and before the MPT (Fig. 2A, B). In contrast, the ~100-kyr power ratio shows no significant correlation with the ~41-kyr power ratio (Fig. 2A, C). Further analysis using the median power ratio values for the ~21-kyr, ~41-kyr, and ~100-kyr bands, represented by the horizontal lines in the box whisker plots, confirms a consistent negative correlation between the ~100-kyr and ~21-kyr power ratios (Fig. 2D, F, H), and no correlation between the ~100-kyr and ~41-kyr power ratios (Fig. 2E, G, I). Heterogeneity in dataset length, age models, and orbital tuning do not influence these relationship (Fig. S10C–F). Additionally, adjusting the frequency bandwidth setting in the PDA does not alter these relationships among the ~21-kyr, ~41-kyr, and ~100-kyr power ratios (Fig. S11).

Methods
Extracting precession, obliquity and ~100-kyr spectral power ratios
To detect long-term evolution of changes in the spectral strength (or variance) of precession, obliquity, and ~100-kyr cycles in the collected climate archives, the evolutionary spectral power ratios at the precession, obliquity, and ~100-kyr bands to the total spectral power of the climate datasets were estimated. The estimation includes three steps. (1) Initially, the evolutionary power spectra of the datasets were calculated using the 2π multi-taper method (MTM) with a sliding time window (such as Fig. S1B). For the paleoclimate datasets spanning 0–800 ka, 800–1200 ka, and 1200–2700 ka, the sliding windows were set at 500 kyr, 200 kyr, and 500 kyr, respectively. (2) Then, the time-dependent amplitude modulations (power) at these three frequency bands were quantitatively calculated by integrating the evolutionary power at corresponding frequency bands. (3) Thereafter, to facilitate comparisons among different paleoclimate datasets, the integrated evolutionary spectral power at each specific band was normalized by the total evolutionary power (integrated from 0 to the Nyquist frequency) to get the evolution spectral power ratio. This is because the unit of different kind of paleoclimate proxies is different.

Over the last 2700 ka, the fundamental frequencies for precession and obliquity are 1/21 kyr-1 and 1/41 kyr-1. Thus, during the power integrating process, the frequency bands for the precession, obliquity, and ~100-kyr components were set to 1/18–1/26 kyr−1, 1/33.4–1/50.3 kyr-1, and 1/77–1/175 kyr-1, respectively.

Two software tools can be used to calculate the spectral power ratios: (1) the integratePower function in R package Astrochron, and (2) the Power Decomposition Analysis (PDA) function in Acycle software. The integratePower function processes both spectrograms (time-frequency or evolutionary spectra) and single spectra. Thus, for estimating the overall spectral power ratios of spectific frequenct bands in the clipped 0E0T-P series, we used integratePower function (e.g., Fig. 3J). Nevertheless, integratePower causes windows loss when estimating the evolutionary spectral power ratios, an issue resolved by the PDA function. Therefore, we used the PDA to calculate the evolutionary spectral power ratios for the 57 collected paleoclimate datasets.

Robustness estimation for the spectral power ratio results
Extracting evolutionary spectral power ratios is a non-stationary method and is relatively insensitive to age models because slight variations in the age model do not influence the long-term evolution trend of the spectral power ratio at a specific frequency. As a test, the evolution spectral power ratios at the precession and obliquity frequency bands for three different benthic δ18O stacks were examined (Fig. S1D, F, G). These stacks include the LR04 stack with LR04 age model, the LR04 stack with LR10 age model, and the Prob_stack, (Fig. S1D, F). The LR04 age model was established by tuning the δ18O stack to a target derived from an ice sheet model driven by 65°N summer insolation. The LR10 age model was based on the average sedimentation rate between tie points, independent of orbital forcing assumptions. The Prob_stack, compiled from 180 globally distributed benthic δ18O records and developed using a profile hidden Markov model, is more representative than the LR04 stack. Nevertheless, the Prob_stack adopted the LR04 age model as an initial condition during the stacking process, thereby inheriting the LR04 age model.

The calculated spectral power ratios of these stacks all show similar evolution trends, which also align with those acquired through the ENOF method (Fig. S1C, E). The ENOF is a non stationary method and employed to estimate the temporal amplitude of precession and obliquity in the benthic δ18O over the last 2700 ka. Therefore, the PDA method is also suitable for non-stationary analysis and can be regarded as an alternative to the ENOF method.

Heterogeneity in dataset length, age models, and orbital tuning may distort PDA results, potentially altering the spectral power ratio correlations. To mitigate these impacts, only PDA results from paleoclimate archives longer than 600 kyr over the past 800 ka were considered (Fig. S10C). Additionally, archives subjected to orbital tuning were separately assessed (Fig. S10D). To further reduce the effects of orbital tuning, datasets primarily constrained by the LR04 age model (Fig. S10E, F; Table S2) were tuned to the LR10 age model. The LR10 age model is independent of orbital tuning. Across all scenarios, the ~21-kyr power ratios negatively correlate with ~100-kyr power ratios, with no observed relationship between the ~41-kyr and ~100-kyr power ratios (Fig. S10). Additionally, the frequency bandwidth setting for the PDA does not influence the relationship among the ~21-kyr, ~41-kyr, and ~100-kyr power ratios (Fig. S11).

Phase coupling estimation between benthic δ13C and δ18O/ice volume
To obtain the phase relationship between benthic δ13C and δ18O/ice volume at the ~100-kyr timescale, the ~100-kyr components were initially extracted from benthic δ13C and δ18O/ice volume using a Gaussian bandpass filter with frequency bands setting of ~1/77–1/175 kyr-1. The Cross-Wavelet Transform (XWT), based on these filters was then estimated using the MATLAB Cross Wavelet and Wavelet Coherence Toolbox. To further quantitatively estimate their coupling relationship, The evolutionary correlation coefficient between the ~100-kyr filters of benthic δ13C and δ18O/ice volume was calculated with a 1000-kyr sliding windows. A correlation coefficient above 0.5 or below -0.5 indicates in-phase or antiphase coupling, respectively. To keep the consistency between δ18O and ice volume, the sea-level equivalents were reversed before these calculations.

Lead-lag relationship analysis
Lead-lag analysis was used to estimate the lead/lag relationship between two datasets, reference x and target y series. When target y lags behind reference x, then a positive lag value will be return, which means that reference x occurred earlier than y. The lead-lag relationship analysis was realized using Acycle software.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper tackles a fundamental problem in understanding ice age behavior in the Pleistocene, which records several mysterious changes in behavior. No one doubts the critical influence of orbital cycles in controlling the timing of the ice age cycles, but no linear model comes close to capturing the geological observations. This paper attempts a statistical approach to the problem, using a variety of previously published archives of climate change.\nUnfortunately, the paper is fundamentally flawed in its present form. I list below briefly major issues which I urge the authors to evaluate.\n1. Although the role of CO2 is invoked frequently, high resolution and replicated estimates are only available for the past ~800,000 years. The record is being extended, most successfully in my opinion using the d11B proxy, but is incomplete. Critically, the proxy data do not support the idea of a gradual drawdown of CO2 at around the tijme of the MPT. Neither do temperature records show such a behavior. I do not see citations to these CO2 or SST reconstructions in the manuscript. Yet this argument of interglacial CO2 drawdown and cooling “pre-conditioning” ice sheet behavior during the subsequent glacial cycle is central to the interpretation here. The CO2 estimates used here consist of model results, using the d18O proxy and cannot be considered as independent of the isotopic data. They have not been validated with actual palo-CO2 estimates.\n2. The authors assume that the isotope curve linearly reflects ice volume. Work over the last decade emphasizes the fact that there may be a significant, and evolving overprint of deep ocean temperature. Without separating temperature and ice volume, the interpretations presented here, which focus on the behavior of ice sheets, is incomplete and possibly incorrect. I’d urge them to look at the recent paper by Clark et al. in Science (2024) to see how much room for debate there is in interpreting what d18O is actually telling us.\n3. The paper does not consider issues of age models and sampling resolution, both of which can be critical to isolating the precessional component in geological time series. Given “real-world” issues of sampling resolution and changes in sedimentation rate, accurately capturing precessional signals is quite difficult. A particular problem has to do with the potential circularity in assessing the phasing of precession relative to obliquity in records tuned to obliquity. As an example, the IODP 1385 age model referred to was “retuned” from the standard isotopic age model using a precessional signal by Liautaud et al. This retuning significantly increase the amplitude of the precessional signal, although it comes with its own assumptions.\n4. The time series methods need to be better described. How do the authors treat non-stationarity? With moving windows? Given that precession is rarely resolved well in spectra of “real” data, how do they estimate uncertainties? And finally, how the “phase” is derived robustly is very difficult to determine, yet central to the conclusions.\nDetection of significant climatic precession variability in early Pleistocene glacial cycles\nLiautaud, P. R. Liautaud, D. A. Hodell and P. J. Huybers, Earth and Planetary Science Letters 2020 Vol. 536 Pages 116137\nDOI: https://doi.org/10.1016/j.epsl.2020.116137", "reviewer_2": "1 Response Letter\n2\n3 Response to Reviews\n4 Reviewer #1:\n5 (General comments) This paper tackles a fundamental problem in understanding ice\n6 age behavior in the Pleistocene, which records several mysterious changes in behavior.\n7 No one doubts the critical influence of orbital cycles in controlling the timing of the ice\n8 age cycles, but no linear model comes close to capturing the geological observations.\n9 This paper attempts a statistical approach to the problem, using a variety of previously\n10 published archives of climate change.\n11 Unfortunately, the paper is fundamentally flawed in its present form. I list below\n12 briefly major issues which I urge the authors to evaluate.\n13\n14 Major points:\n15 1.1. Although the role of CO is invoked frequently, high resolution and replicated\n2\n16 estimates are only available for the past ~800,000 years. The record is being extended,\n17 most successfully in my opinion using the δ11B proxy, but is incomplete. Critically, the\n18 proxy data do not support the idea of a gradual drawdown of CO at around the time of\n2\n19 the MPT. Neither do temperature records show such a behavior. I do not see citations\n20 to these CO or SST reconstructions in the manuscript. Yet this argument of interglacial\n2\n21 CO drawdown and cooling “pre-conditioning” ice sheet behavior during the\n2\n22 subsequent glacial cycle is central to the interpretation here. The CO estimates used\n2\n23 here consist of model results, using the δ18O proxy and cannot be considered as\n24 independent of the isotopic data. They have not been validated with actual palo-CO\n2\n25 estimates\n26 Our reply: Thank you for this comment. In the revised manuscript, we have\n27 incorporated the reconstructed CO records, including those based on δ11B (Foster,\n2\n28 2008; Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011; Henehan et al., 2013;\n29 Chalk et al., 2017; Dyez et al., 2018), reconstructions derived from Allan Hills Blue Ice\n30 (Yan et al., 2019) and EPICA ice core (0-800ka; Bereiter et al., 2015), and a continuous\n31 CO reconstruction based on the leaf-wax δ13C proxy spanning the last 1500 ka\n2\n32 (Yamamoto et al., 2022) (Fig. 5A). Furthermore, we included the new stack of global\n33 mean surface temperature (GMST) (Clark et al., 2024) (Fig. 5B).\n34 The proxy-based CO reconstructions do show a glacial-stage decline across the\n2\n35 MPT as evidenced by CO reconstructions based on δ11B (Hönisch et al., 2009; Chalk\n2\n36 et al., 2017; Dyez et al., 2018), ice core (Bereiter et al., 2015; Yan et al., 2019), and\n37 leaf-wax δ13C proxies (Yamamoto et al., 2022) (Fig. 5A). As glacial-stage temperature\n38 determines the maximal extension of ice sheets, the glacial-stage CO decline could be\n2\n39 critical for understanding the changes in ice sheet dynamics across the MPT.\n40 Moreover, the stacked global mean surface temperature (GMST) data (Clark et al.,\n41 2024) clearly demonstrate a long-term cooling trend across the MPT (Fig. 3D; Fig. 5B).\n42 All these lines of evidence support a long-term cooling trend across the MPT. These\n43 findings strengthen our interpretation of the significant role that background\n44 temperature changes play in modulating the glacial response to precession.\n45\n46 1.2. The authors assume that the isotope curve linearly reflects ice volume. Work over\n47 the last decade emphasizes the fact that there may be a significant, and evolving\n48 overprint of deep ocean temperature. Without separating temperature and ice volume,\n49 the interpretations presented here, which focus on the behavior of ice sheets, is\n50 incomplete and possibly incorrect. I’d urge them to look at the recent paper by Clark et\n51 al. in Science (2024) to see how much room for debate there is in interpreting what\n52 d18O is actually telling us.\n53 Our reply: Thank you for your comment. We agree that benthic δ18O could be affected\n54 by local deep-sea temperature, but it remains to be a valuable proxy for discussing\n55 global ice volume changes as has been done and continues to be done in numerous\n56 studies (e.g. Huybers and Wunsch, 2005; Huybers, 2006; Drysdale et al., 2009; Lisiecki,\n57 2010; Huybers, 2011; Tzedakis et al., 2017; Liautaud et al., 2020; Hobart et al., 2023),\n58 especially when many δ18O records from different oceans are used like in our study,\n59 which is expected to mitigate the effect of local deep-sea temperature and provide a\n60 more accurate reflection of global ice volume changes.\n61 In our manuscript, the definition of the glacial cycles follows the scheme proposed\n62 by Tzedakis et al (2017), which is based on the detrended LR04 stack. The purpose of\n63 detrending is to remove the temperature contribution to benthic δ18O (Tzedakis et al.,\n64 2017). As depicted in Fig. A1D, the detrended LR04 stack closely aligns with the\n65 separate ice volume proxy (Elderfield et al., 2012). Thus, the estimated precession and\n66 obliquity phases are minimally affected by temperature influences on benthic δ18O.\n67 Nevertheless, we also considered the δ18Ow record from Site ODP1123 for further\n68 validation of our results because this record represents the contribution of ice volume\n69 to benthic δ18O (Elderfield et al., 2012). As shown in Fig. A1, the estimated precession\n70 and obliquity phases display similar trends to those derived from the LR04 stack (Fig.\n71 A1B, C, H&I). Moreover, the extracted ~21-kyr spectral power ratio of δ18Ow shows a\n72 negative correlation with the ~100-kyr spectral power ratio after the MPT (r=-0.89),\n73 whereas their relationship is positive before the MPT (r=0.66) (Fig. A2E). This suggests\n74 a shift in the response of ice sheets to precession from early linearity to later nonlinearity,\n75 which is consistent with the results acquired from the H07&NA stack (Fig. 5E, G). Thus,\n76 our estimated precession and obliquity phases are minimally affected by temperature\n77 influences on benthic δ18O.\n78\n79 Figure A1. The precession and obliquity phase of glacial terminations, based on ice volume proxy\n80 of δ18Ow in Site ODP1123 (Elderfield et al., 2012). A, the detrended LR04 benthic δ18O stack\n81 (Lisiecki and Raymo, 2005; Tzedakis et al., 2017), the red and green dots denote the glacial minima\n82 (or interglacial maxima) and maxima within glacial cycles, the orange dotes refer to the timing of\n83 the local maximal changing rate of benthic oxygen isotope within a given glacial cycle. The black\n84 numbers above the δ18O stack mark the marine isotope stages (MIS). The orange tags and numbers\n85 bellow the δ18O stack label the Termination numbers (T) and corresponding ages. The horizontal\n86 dotted lines are the threshold used to define the glacial cycles (Tzedakis et al., 2017). B&C, the\n87 precession and obliquity phase during the maximal changing rate of LR04 benthic δ18O stack,\n88 as well as their corresponding 5-points moving average (5pma; in blue). D, the δ18Ow in Site ODP1123,\n89 which represents the ice volume component extracted from the benthic δ18O (Elderfield et al., 2012).\n90 Same as in A, the red and green dots denote the glacial minima and maxima within glacial cycles,\n91 the orange dotes refer to the timing of the local maximal changing rate of δ18Ow. The rosters of\n92 glacial cycles refer to scheme provided by Tzedakis et al. (2017). E, the derivative of δ18Ow, the\n93 orange dots label the maximal changing rate of δ18Ow within glacial cycles. F&G, the precession\n94 and obliquity parameters (Berger and Loutre, 1991). H&I, the precession and obliquity phase during\n95 the maximal changing rate of δ18Ow (in thin black curve), as well as their corresponding 5-points\n96 moving average (5pma; in thick black curve).\n97\n98\n99 Figure A2. The calculated ~100-kyr (A) and ~21-kyr (B) spectral power ratios, as well as the 5-\n100 points moving average of the estimated precession phase (C) based on δ18Ow in Site ODP1123\n101 (Elderfield et al., 2012). The correlation plots among these variabilities were also created for each\n102 pair (D, E, and F).\n103\n104 1.3. The paper does not consider issues of age models and sampling resolution, both of\n105 which can be critical to isolating the precessional component in geological time series.\n106 Given “real-world” issues of sampling resolution and changes in sedimentation rate,\n107 accurately capturing precessional signals is quite difficult. A particular problem has to\n108 do with the potential circularity in assessing the phasing of precession relative to\n109 obliquity in records tuned to obliquity. As an example, the IODP 1385 age model\n110 referred to was “retuned” from the standard isotopic age model using a precessional\n111 signal by Liautaud et al. This retuning significantly increases the amplitude of the\n112 precessional signal, although it comes with its own assumptions.\n113 Our reply: Thank you for this comment. For all of our complied datasets, the sampling\n114 resolution is higher than 2 samples per kyr. According to the Nyquist frequency, this\n115 resolution is sufficient to detect the precession signal. Moreover, in the revised\n116 manuscript, we tested the confidence level (CL) of the power peaks at frequencies of\n117 ~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1 within the datasets (Figs. S6-S9). Only the\n118 power ratios with corresponding power peaks above 90% CL were considered in\n119 assessing relationship among spectral power ratios (Fig. S10A, B).\n120 Additionally, to establish a robust relationship in spectral power ratios, we utilized\n121 datasets spanning the past 800 ka and discussed six different scenarios regarding the\n122 influence of significant level of power peaks (Fig. S10A, B), datasets length (Fig. S10C),\n123 age models used (Fig. S10D), as well as the effects of orbital tuning (Fig. S10E, F) on\n124 the correlation of power ratios at ~21-kyr, ~41-kyr. and ~100-kyr frequency bands.\n125 Specifically, to eliminate the influence of orbital tuning, we tuned datasets mainly\n126 constrained by the LR04 age model to the H07&NA age model (Fig. 10E, F; Figs. S12).\n127 The H07&NA stack combines the depth-derived H07 stack (Huybers, 2007) with its\n128 last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack) (Hobart et al.,\n129 2023) (Fig. S11; Table S2). The age model for the NA stack was based on the\n130 speleothem chronology (Cheng et al., 2016), making the H07&NA stack independent\n131 of orbital tuning. Across all scenarios, the ~21-kyr power ratios negatively correlate\n132 with the ~100-kyr power ratios, while there is no significant correlation between the\n133 ~41-kyr and ~100-kyr power ratios (Fig. S10), consistent with our previous findings.\n134 Regarding the potential issue of circular reasoning in assessing the phasing of\n135 precession and obliquity. To avoid this, in the revised manuscript, we recalculated the\n136 precession and obliquity phases of glacial terminations based on the H07&NA stack,\n137 featuring an age model independent of orbital tuning. The orbital phases calculated\n138 from this stack align with the evolutionary trends obtained from the LR04 stack and\n139 Prob_stack (Fig. 3; Fig. S14–16), thereby confirming our previous results.\n140\n141 1.4. The time series methods need to be better described. How do the authors treat non-\n142 stationarity? With moving windows? Given that precession is rarely resolved well in\n143 spectra of “real” data, how do they estimate uncertainties?\n144 Our reply: In terms of the time series methods, we used the power decomposition\n145 analysis (PDA) (Li et al., 2019). This method allows for the quantitative estimation of\n146 the influence of precession, obliquity, and ~100-kyr cycles on paleoclimate archives.\n147 Analogous to the Empirical Nonlinear Orbital Fitting (ENOF) method (Liautaud et al.,\n148 2020), this non-stationary technique, PDA, employs a sliding time window to calculate\n149 evolutionary power ratios within specific frequency bands (Fig. S1). For the\n150 paleoclimate datasets spanning 0–800 ka, 800–1200 ka, and 1200–2700 ka, the sliding\n151 windows were set at 500 kyr, 200 kyr, and 500 kyr, respectively.\n152 In the revised manuscript, we have added more details to introduce the PDA\n153 method (see Method section) and rewrote the first section of the Method to improve\n154 clarity on the principles and processing steps of PDA. The revised description is as\n155 follows: (1) Initially, the evolutionary power spectra of the datasets were calculated\n156 using the 2π multi-taper method (MTM) (Thomson, 1982) with a sliding time window\n157 (Meyers et al., 2012) (Fig. S1B). (2) Next, the time-dependent amplitude modulations\n158 (power) at these three frequency bands (~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1) were\n159 quantitatively calculated by integrating the evolutionary power at the corresponding\n160 frequency bands (Meyers et al., 2012; Li et al., 2019). (3) Thereafter, to facilitate\n161 comparisons among different paleoclimate datasets, the integrated evolutionary spectral\n162 power at each specific band was normalized by the total evolutionary power (integrated\n163 from 0 to the Nyquist frequency) to obtain the evolution spectral power ratio. This\n164 process is necessary because different paleoclimate proxies have varying units.\n165 Moreover, to better illustrate this method, the evolution spectral power ratios at the\n166 precession and obliquity frequency bands for three different benthic δ18O stacks were\n167 examined. These include the LR04 stack (Lisiecki and Raymo, 2005), the Prob_stack\n168 (Ahn et al., 2017), and the H07&NA stack (Huybers, 2007) (Fig. S1D, F). The\n169 calculated spectral power ratios from these stack exhibit similar evolution trends,\n170 consistent with the results obtained through the Empirical Nonlinear Orbital Fitting\n171 (ENOF) (Fig. S1C, E) (Liautaud et al., 2020).\n172\n173 1.5. And finally, how the “phase” is derived robustly is very difficult to determine, yet\n174 central to the conclusions.\n175 Our reply: We took several steps to ensure the robustness of the phase estimation.\n176 Initially, several different benthic δ18O stack and age models were used to estimate the\n177 precession and obliquity phases (Fig. 3B, C; Figs. S14-16); Secondly, to minimize the\n178 influence of errors from individual precession or obliquity phases at specific glacial\n179 terminations on the overall evolution trend of phases, we applied a 5-point moving\n180 average (5 pma). Thirdly, the H07&NA stack was employed to estimate the\n181 astronomical phases of glacial terminations, therefore mitigating the circular reasoning.\n182 The 5 pma values of the estimated precession and obliquity phases, derived from three\n183 different benthic δ18O stacks, display consistent evolution trends, demonstrating the\n184 robustness of our phase results.\n185\n186 Minor points:\n187 1.6. line 23: repetition of \"transition\"\n188 Our reply: Thank you. In the revised manuscript, we have reorganized the abstract,\n189 and eliminated redundancy.\n190\n191 1.7. line 24-30: since there are no continuous, orbitally resolved pCO measurements >\n2\n192 800 ka, this can't be supported as written. Note also that the δ18O is equated with ice\n193 volume- this is increasingly questioned with recent work, e.g. Clark et al. 2024. Also it\n194 isn't clear what \"phase covaries with pCO … dynamics...\"- means concretely in terms\n2\n195 of the data and/or statistical analysis\n196 Our reply: Thank you. In terms of the usage of pCO as a reference for comparison\n2\n197 with precession phase. We have reorganized the Abstract and emphasized the role of\n198 gradually-cooling temperature in modulating glacial response to precession (Fig. 5B).\n199 This emphasis is supported by the recent publication of a continuous stack of global\n200 mean surface temperatures, facilitating a direct correlation with reconstructed\n201 precession phase.\n202 A continuous CO reconstruction throughout the Pleistocene is still lacking.\n2\n203 However, the proxy-based CO reconstructions do indicate a glacial-stage decline trend\n2\n204 across the MPT (Fig. 5A). Therefore, the decline in glacial-stage pCO (Yan et al., 2019)\n2\n205 suggests a potential relationship with the temperature decline observed across the MPT.\n206 For further details, please refer to our response to the first major comment (1.1).\n207 Regarding the concern about using δ18O as a proxy for ice volume, please see our\n208 response to the second major comment (1.2).\n209\n210 1.8. Line 52: signal usually refers to a representation in the earth system (e.g. proxy\n211 evidence) not to a pacing/forcing such as orbital cycles.\n212 Our reply: Thank you. Now we used the astronomical ‘parameters’ to substitute the\n213 ‘signals.’\n214\n215 1.9. Line 59: without clarification, the reader can't follow what is meant by \"only\n216 specific instances\"\n217 Our reply: Thank you. In the revised manuscript, we have reorganized the\n218 ‘Introduction’, and this sentence has been deleted.\n219\n220 1.10. Line 70-73: since any CO reconstruction at this point is indirect, I wouldn't use\n2\n221 this phrase.\n222 Our reply: Thank you. In the revised manuscript, we have changed this narrative as\n223 follows: These evolutionary phases enable a comparison with previously simulated and\n224 reconstructed pCO (Foster, 2008; Hönisch et al., 2009; Seki et al., 2010; Bartoli et al.,\n2\n225 2011; Henehan et al., 2013; Bereiter et al., 2015; Chalk et al., 2017; Dyez et al., 2018;\n226 Yan et al., 2019; Willeit et al., 2019; Yamamoto et al., 2022), global mean surface\n227 temperature (ΔGMST) (Clark et al., 2024), integrated glacial latitudinal extension\n228 (Köhler and van de Wal, 2020), and ice volume (de Boer et al., 2014; Willeit et al., 2019;\n229 Berends et al., 2021-CP), thereby allowing for the establishment of a direct and\n230 successive correlation between internal and external factors.\n231\n232 1.11. Line 73-74: this sentence (These analyses contribute to piecing together a jigsaw\n233 for understanding the MPT and the ~100-kyr glacial cycles.) does not contain any useful\n234 information.\n235 Our reply: Thank you. We have deleted this sentence in the revised manuscript.\n236\n237 1.12. Line 78-83: some minimal discussion of the spectral technique is needed- is some\n238 kind of evolutive approach used to capture the time-dependent spectral characteristics\n239 (e.g. non-stationarity)? Also, since the records in many cases are not of the same length,\n240 how is this heterogeneity handled?\n241 Our reply: For the first suggestion, we have enhanced the introduce of the PDA method\n242 with additional details. In the revised manuscript, the Result section begins with ‘The\n243 power decomposition analysis (PDA) (Li et al., 2019-ESR) was employed to\n244 quantitatively estimate the influence of precession, obliquity, and ~100-kyr cycles on\n245 paleoclimate archives. Analogous to the Empirical Nonlinear Orbital Fitting (ENOF)\n246 method (Liautaud et al., 2020) …’ Moreover, we have included a supplementary figure\n247 (Fig. S1) to visually explain the principles and processing steps of the PDA method.\n248 The variability in dataset length may affect the correlation of spectral power ratios.\n249 To mitigate this influence, only PDA results from paleoclimate archives extending\n250 longer than 600 kyr within the last 800 ka were considered (Fig. S10C). The analysis\n251 shows a negative correlation between the ~21-kyr and ~100-kyr power ratios, while no\n252 significant relationship is observed between the ~41-kyr and ~100-kyr power ratios.\n253 Notably, for datasets spanning the 800-1200 ka, their datasets lengths are same and as\n254 long as 400 kyr, their medians of ~21-kyr power ratios still show a strong positive\n255 correlation with those at ~100-kyr bands (Fig. 2F, G). Thus, we conclude that, despite\n256 the heterogeneity in dataset length, our findings remain robust and reliable.\n257\n258 1.13. Line 95-96: Prob stack is still not independent of an underlying orbital tuning for\n259 its timing, as shown by its similarity to LR04. It isn't clear how the authors avoid either\n260 the circularity in tuning of the isotope record or the inherent uncertainties in any age\n261 model. The 1385 hybrid model represents a different approach from LR04/prob stack\n262 in that it is tuned to a lithological precession component, so although tuned, it does use\n263 the isotope record in the tuning. Ref. 23 uses a similar approach. Note that neither the\n264 1385 or Ref. 23 approach extend significantly before the MPT.\n265 Our reply: Thank you for your comment. We acknowledge that Prob_stack has adopted\n266 the LR04 age model as an initial condition during the stacking process, which resulted\n267 in our previous phase estimations being subject to circular reasoning. To address this\n268 issue in the revised manuscript, we used the depth-derived H07 stack (Huybers, 2007)\n269 with its last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack) (Hobart\n270 et al., 2023) (Fig. S11; Table S2). The age model for the NA stack was established based\n271 on the speleothem age model (Cheng et al., 2016). The tuned H07 stack, termed\n272 H07&NA stack, is thus independent of orbital tuning. Subsequently, the H07&NA stack\n273 was used to recalculate the precession and obliquity phases of glacial terminations (Fig.\n274 S16). The orbital phases calculated from this stack align with the evolutionary trends\n275 obtained using the LR04 stack and Prob_stack (Fig. 3; Fig. S14–16), confirming our\n276 previous results.\n277 Moreover, in the revised manuscript, the LR04 stack with U1385 hybrid age model\n278 (Hodell et al., 2023) and the LR04 stack with Konijnendijk age model (Konijnendijk et\n279 al., 2015) have not been considered for calculating precession and obliquity phases\n280 because of their shorter length and their age models predominantly inheriting the LR04\n281 age model.\n282\n283 1.14. Line 111-118: given that the data sets represent different aspects of climate, should\n284 we expect similar power ratios and/or \"phases\". How do the authors handle the\n285 dispersion of these quantities shown in Figure 2- in other words, what is the confidence\n286 level of the estimates? I also find the figure visually confusing.\n287 Our reply: Thank you for the comment. Although the climate implications of the\n288 collected datasets are different, it is important to note that the climate subsystems\n289 interact extensively. Thus, the relationship in spectral power ratios within these\n290 subsystems can propagate through the entire climate system via a series of positive and\n291 negative feedbacks (Clark et al., 1999; Denton et al., 2010). The expansion and retreat\n292 of Pleistocene ice sheets respond not only to external orbital forcing but also act as\n293 significant internal forcings, driving and influencing coeval regional and global\n294 climates. For example, they affect sea-level variations and atmospheric-oceanic\n295 circulations through the release of meltwater (Imbrie et al., 1993; Clark et al., 1999;\n296 Toggweiler et al., 2006; Denton et al., 2010; Lambert et al., 2012; Sun et al., 2012;\n297 Cheng et al., 2016) and alter Earth’s energy budget through albedo feedback.\n298 Consequently, the spectral power ratio relationship between precession/obliquity and\n299 the ~100-kyr components, retained in ice volume, can be projected onto other climate\n300 archives. Thus, regarding the internal nonlinear feedback mechanism that causes the\n301 transfer of power and subsequently amplifies the ~100-kyr cycles, regardless of\n302 whether it is a feature exclusive to the glacial subsystem or intrinsic to the entire climate\n303 system, an analysis of the spectral power ratio relationship between the ~100-kyr cycles,\n304 precession, and obliquity components within ice volume and other contemporaneous\n305 paleoclimate archives can determine whether the amplified power of ~100-kyr glacial\n306 cycles is derived from precession or obliquity.\n307 To ensure the robustness of our findings, we applied a statistical method to\n308 demonstrate a reversal correlation between power ratios at ~21-kyr and ~100-kyr bands\n309 in the climate system, rather than relying solely on benthic δ18O data. Moreover, we\n310 restricted our analysis to PDA results from the paleoclimate archives where spectral\n311 power at frequencies of ~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1 exceeding the 90%\n312 significance level. In terms of Fig. 2, in the revised manuscript, we have reorganized\n313 the presentation for clarity. For a better comparability of the power ratio changes across\n314 three intervals (0-800 ka; 800-1200 ka; 1200-2700 ka), we merged three panels into one\n315 and standardized their axis ranges. Furthermore, the scatter plots have been magnified\n316 for enhancing detail.\n317\n318 1.15. Line 121-125: the modelled CO in fact comes from an interpretation of the stable\n2\n319 isotope record as a forcing function, so it's hard to see how these are independent\n320 estimates of different aspects of earth climate. Note that there is no evidence of\n321 \"progressive decrease\" of CO from proxy evidence such as δ11B, if \"progressive\n2\n322 decrease\" means averaged over the period of several orbital cycles. The lack of trend in\n323 obliquity phase is likely to come from the assumptions of the age models that align the\n324 signals to obliquity.\n325 Our reply: In the revised manuscript, we have supplemented the proxy-based CO\n2\n326 reconstructions, including δ11B-based reconstruction (Chalk et al., 2017; Seki et al.,\n327 2010; Foster, 2008; Henehan et al., 2013; Hönisch et al., 2009; Dyez et al., 2018; Bartoli\n328 et al., 2011), reconstructions derived from Allan Hills Blue Ice (Yan et al., 2019) and\n329 EPICA ice core (0-800ka; Bereiter et al., 2015), and a continuous CO reconstruction\n2\n330 based on the leaf-wax δ13C proxy spanning the last 1500 ka (Yamamoto et al., 2022)\n331 (Fig. 5A).\n332 Despite the absence of continuous CO reconstruction throughout the past 2700 ka,\n2\n333 the reconstructions based on δ11B and ice core archives consistently display a glacial-\n334 stage decline in CO across the MPT. Moreover, the continuous CO reconstruction\n2 2\n335 based on the leaf-wax δ13C proxy over the past 1500 ka also indicates a gradual decline\n336 in glacial-stage CO until the MPT (Fig. 5A). We also included the new stack of global\n2\n337 mean surface temperature (GMST) (Clark et al., 2024) (Fig. 5B), which confirms a\n338 long-term cooling trend across the MPT (Fig. 3D; Fig. 5B). Therefore, various lines of\n339 evidence show that there is a consistent decline in the global mean temperature and\n340 glacial-stage pCO across the MPT.\n2\n341 Regarding the comment ‘The lack of trend in obliquity phase is likely to come from\n342 the assumptions of the age models that align the signals to obliquity.’\n343 In our revised manuscript, we present three obliquity phase reconstrutions based\n344 on the LR04 stack, Prob_Stack, and H07&NA stack, revealing a similar evolutionary\n345 trend. The H07&NA stack is indenpend of orbtial tuning, confirming that the evolution\n346 trends of constructed precession and obliquity phases are not induced by the orbital\n347 tuning. The estimated obliquity phases show an increasing trend from 2700 ka to around\n348 800 ka, stabilizing thereafter (Fig. 3B). Nevertheless, this evolutionary trend is less\n349 prominent. In contrast, the precession phase follows an opposite trend, suggesting an\n350 increasing glacial sensitivity to precession (Feng and Bailer-Jones, 2015; Liautaud et\n351 al., 2020; Barker et al., 2022). Menwhile, obliquity’s role becomes comparably less\n352 significant, even though it remains indispensable for deglaciations (Feng and Bailer-\n353 Jones, 2015).\n354\n355 1.16. Line 135-137: It's not clear why \"phase\" is connected with sensitivity. The phase\n356 offset is usually a measure of the response time of the sytem relative to the forcing, and\n357 response time != sensitivity (or sensitivity needs to be defined so that there is a clear\n358 connection to the \"phase\")\n359 Our reply: In our manuscript, the obliquity and precession phases are defined as\n360 the absolute age offsets between the timing of glacial terminations and the timing at the\n361 nearest precession minima or obliquity maxima. Following Huybers (2011), a glacial\n362 termination is pinpointed at the maximal changing rate of benthic δ18O within a glacial\n363 cycle. Thus, the phase here estimates how quickly the ice sheets respond to precession-\n364 or obliquity-induced insolation changes. A smaller phase implies a shorter response\n365 time of ice sheets to precession or obliquity, indicating greater susceptibility of the ice\n366 sheets to these forcing mechanisms. Building on prior studies that discussed whether\n367 precession or obliquity has a greater influence on ~100-kyr glacial cycles (e.g., Hobart\n368 et al., 2023), we use ‘sensitivity’ to express phase, with a smaller phase indicating\n369 higher sensitivity. In the revised manuscript, we have added this explanation in the\n370 Method section.\n371\n372 1.17. Line 150-154: this is asserted, but isn't born out so far by proxy CO estimates.\n2\n373 The \"close correlations\" mentiuoned here largely come from using δ18O as a surrogate\n374 for estimating CO2 so there is some circularity involved.\n375 Our reply: We agree with your comment and therefore have deleted this sentence.\n376 However, as clarified in our replies to comments 1, 7, 10, and 15, the evolution of CO\n2\n377 does not affect our major conclusions. This is supported by the stacked global mean\n378 surface temperature (Clark et al., 2024), which clearly display a cooling trend across\n379 the MPT, consistent with our interpretation. Moreover, our estimated precession phase\n380 displays a strong correlation with global mean surface temperature (r=0.88) (Fig. 5B,\n381 F; Fig. S17F). Thus, in the revised manuscript, we emphasized the role of background\n382 temperature in modulating glacial sensitivity to astronomical forcing, rather than\n383 establishing a direct connection between phase and CO .\n2\n384\n385 1.18. Line 183-186: Clarify whzt you mean specifically by \"evolution'. A number of\n386 points were made above and it's hard for the reader to follow what you are comparing\n387 the obliquity evolution to.\n388 Our reply: In the revised manuscript, we have revised this statement to: ‘The obliquity\n389 phase demonstrates an increasing trend from 2700 ka to around 800 ka, then stabilizes\n390 (Fig. 3B). Conversely, the precession phase follows an opposite trend (Fig. 3C),\n391 indicating a progressive increase in ice sheet sensitivity to precession until around 800\n392 ka, after which it remained stable.’\n393\n394 1.19. Line 233-235: since the power ratio is relative, isn't it possibl that the actual\n395 importance of precession increased toward the present, but that the 100 kyr grew more?\n396 This wouldn't mean that the 100 kyr increeases at the \"expensee\" of precession, but\n397 rather that it increased by a non-linear process as precession became more important.\n398 Our reply:\n399\n400 Figure A3. The ~100-kyr power ratio and ~21-kyr power ratio shows a positive correlation when\n401 they grow over time in the same direction.\n402\n403 As you supposed, the importance of ~100-kyr component may grow more rapidly\n404 compared to the ~21-kyr component. We demonstrated this situation using two simple\n405 linear functions, as shown in Fig. A3A. This situation will still result in a positive\n406 correlation between ~21-kyr and ~100-kyr power ratios, rather a negative one (Fig.\n407 A3B). In fact, as long as the ~100-kyr power ratio shows a similar increasing trend as\n408 the ~21-kyr power ratio, they will display a positive correlation. However, our statistical\n409 analysis indicates a negative correlation between the ~21-kyr and ~100-kyr power ratio.\n410 Thus, the only scenario aligning with the observed correlation is that the ~100-kyr\n411 power ratio decreases over time.\n412\n413 1.20. Line 273-276: However, there is not indication that interglacial tmperatures or\n414 CO were in fact decreasing. In fact, there is basically no trend.\n2\n415 Our reply: Although there is no significant trend in the interglacial CO levels across\n2\n416 the MPT; the glacial-stage CO displays a decreasing trend across the MPT (Fig. 5A).\n2\n417 Besides, the stacked global mean surface temperature (ΔGMST) shows cooling trend\n418 across the MPT (Fig. 3D) (Clark et al., 2024). Notably, during the MPT period, the\n419 cooling of ΔGMST displays an accelerating trend (Fig. 1D&E in Clark et al., 2024).\n420 These conditions favor the development of larger continental ice sheets during glacial\n421 periods after the MPT and the merging of several insolation peaks into ~100-kyr\n422 periodicities. Therefore, these observations explain the growing influence of precession\n423 on ice sheets across the MPT, as the ‘over’ developed ice sheets could extend southward\n424 gradually where precessional influence on daily insolation increases (Fig. 3C; Fig. 5F).\n425 Moreover, the estimated precession phase displays a strong correlation with\n426 temperature (r=0.88) (Fig. 5B, F; Fig. S17F). Thus, in the revised manuscript, we\n427 emphasized the role of background temperature in modulating glacial sensitivity to\n428 astronomical forcing, rather than establishing a direct connection between phase and\n429 CO . Please also see our replies to comments 1, 7, 10, 15, and 17.\n2\n430\n431 1.21. Line 314-317: is the \"phase\" analysis here strong enough to support the claim of\n432 low latitude processes?\n433 Our reply: In our revised manuscript, we have removed this section. To emphasize our\n434 new findings that the amplified ~100-kyr glacial power is most likely derived from\n435 precession rather than obliquity, we delve deeper into the reasons behind this conclusion\n436 (please refer to the first part of the revised Discussion) and discuss its implications for\n437 interpreting glacial cycles before and after the MPT (please refer to the second part of\n438 the revised Discussion).\n439\n440 Response to Reviews\n441 Reviewer #2:\n442 (General comments) The cause of the Mid-Pleistocene Transition (MPT) is one of the\n443 most studied and enduring topics in the paleoclimate community. A better\n444 understanding of the MPT will lead to advances in untangling the internal feedback\n445 mechanisms that govern the climate pacing, improving the long-term predictability of\n446 future change. Drawing from analyses done on 57 paleoclimate records as well as the\n447 insolation, Zhang et al. propose that the gradual shift in atmospheric CO and ice sheet\n2\n448 extent caused a sudden shift in climate from a linear to a nonlinear system in response\n449 to precession around 1.5 million years ago. The meticulous investigation of a\n450 comprehensive compilation of datasets produced a novel and mostly convincing theory\n451 on the unfolding of the MPT. In short, the manuscript should be accepted for publication\n452 pending minor revisions to enhance its conclusions.\n453 The strength of Zhang et al.’s argument lies in their use of a vast array of climate\n454 records in testing their hypothesis, as exemplified in their Figure 2. To my knowledge,\n455 an investigation on this scale has not been attempted before, and I must commend the\n456 authors for the tremendous effort. The coherent power decomposition results for a\n457 diverse collection of records enhance the conclusion that precession and eccentricity\n458 powers show a complementary relationship pre- and post-MPT. This relationship\n459 between precession and eccentricity was crucial for the proposed theory that the system\n460 transitioned from a linear to nonlinear response to insolation during the MPT.\n461\n462 Major points:\n463 2.1 With Figure 5, the manuscript makes an argument that before the MPT, the\n464 precession power in LR04 is positively related to eccentricity power and inversely\n465 related to the precession phase. After the MPT, the two relationships flip, signifying a\n466 nonlinear glacial response to insolation. The authors acknowledge that LR04 is orbitally\n467 tuned. To avoid circular reasoning, they did the same precession phase calculation on\n468 Prob_stack (Fig. 3D), stating, “compiled from 180 globally distributed benthic δ18O\n469 records using a profile hidden Markov model, the age model of Prob_stack can shift\n470 freely according to the alignments of benthic δ18O records, thus remaining independent\n471 of orbital tuning.” While theoretically, Prob_stack ages are allowed to shift freely, it\n472 largely inherits the ages of LR04, its alignment target. This can be clearly seen in Figure\n473 4 of Ahn et al. (2017), where the timing of glacial terminations of the two stacks mostly\n474 tracks each other. Compared to LR04, Prob_stack includes more records and quantifies\n475 uncertainty, but it is not immune from the “leak” of orbital tuning in LR04 as its\n476 alignment target. I recommend the authors to try the same tests on untuned benthic\n477 δ18O stacks such as those from Lisiecki (2010) or Huybers (2007). These untuned,\n478 depth-derived stacks are constructed using geomagnetic reversals only, assuming a\n479 constant sedimentation rates between the geomagnetic age control points. Because\n480 geomagnetic reversals are widely dated using a variety of techniques, orbital\n481 assumptions are minimized.\n482 Our reply: Thank you for your comments and suggestions. In the revised manuscript,\n483 we eatablished the H07&NA stack, derived from the depth-derived H07 stack (Huybers,\n484 2007) with its last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack)\n485 (Hobart et al., 2023) (Fig. S11; Table S3). The age model for the NA stack was\n486 established based on the speleothem age model (Cheng et al., 2016), making the\n487 H07&NA stack independent of orbital tuning. The estimated precession and obliquity\n488 phases based on the H07&NA stack align with the evolutionary trends of phases\n489 obtained using the LR04 stack and Prob_stack (Fig. 3; Fig. S14–S16), confirming our\n490 previous results. In the discussion on the reasons causing different glacial response to\n491 precession before and after the MPT, we only considered the astronomical phases of\n492 glacial terminations obtained from the H07&NA stack (Fig. 5).\n493\n494 2.2. Regarding the data presented in Figure 2, tests with different age models are done\n495 for two records in the Supplementary Materials. That is fine. However, it is not clear\n496 what age models are used for the rest of the records in the compilation. Are they using\n497 the age model provided by the respective study? Are they aligned to LR04 and, if so,\n498 do they potentially inherit the orbital tuning of LR04? The authors should discuss this,\n499 as it is likely important for the orbital power decomposition work done here.\n500 Our reply: During the Power Decomposition Analysis (PDA) for calculating the\n501 spectral power ratios, we directly used the original age models of the collected\n502 paleoclimate datasets. In the revised manuscript, we have added supplementary\n503 material introducing the chronologies of these datasets. Please refer to Table S1 and the\n504 section ‘Introduction to chronologies of the collected paleoclimate datasets’ in the\n505 supplementary materials. In summary, only several datasets have their own\n506 chronologies, while the majority have undergone orbital tuning.\n507 To eliminate the influence of orbital tuning on the PDA results, the datasets\n508 subjected to orbital tuning were separately assessed (Fig. S10D). Moreover, datasets\n509 primarily constrained by the LR04 age model (Lisiecki and Raymo, 2005) (Fig.\n510 S10E&F; Table S2) were tuned to the H07&NA age model, based on tie points of glacial\n511 maxima and minima (Figs. S12). As elaborated above, the age model of the H07&NA\n512 is independent of orbital tuning. In these occasions, the ~21-kyr power ratios negatively\n513 correlate with ~100-kyr power ratios, with no observed relationship between ~41-kyr\n514 and ~100-kyr power ratios (Fig. S10). Thus, we confirm that the orbital tuning has no\n515 influence on our results.\n516\n517 2.3. The visualization choices of Figure 2 can be improved. As it stands, it is difficult\n518 to compare the power ratio changes between the three panels. It does not help that the\n519 y-axis range is different among the panels. I encourage the authors to explore ways to\n520 improve Figure 2 so that important features stand out and are easier for readers to\n521 recognize. One possibility is to combine the three panels into one so that the power ratio\n522 boxes of different periods are shown side by side. See what I mean below. Then the\n523 correlation scatter plots can be enlarged, with colors of the scatter plots corresponding\n524 to the same color scheme.\n525 Our reply: It’s a very constructive suggestion., Fig. 2 in the revised manuscript has\n526 undergone major revisions following this suggestion.\n527\n528 Minor points:\n529 2.4. line 20: the full terms should be shown in the first appearances of “kyr” and “ka.”\n530 Also, if the manuscript argues that the MPT occurred ~1500 ka, perhaps avoid stating\n531 that it occurred “around ~800-1200 ka” in the first sentence of Abstract. Consider\n532 changing “around ~800-1200 ka” to “around 1000 ka.”\n533 Our reply: Thank you. We have completely reorganized the abstract. The previous\n534 expression, ‘The transition of glacial cycles from ~40 kyr to ~100 kyr around ~800-\n535 1200 ka, without evident astronomical changes, underscores the interplay of Earth's\n536 internal and external factors in this transition,’ has been deleted. In the revised abstract,\n537 we begin with, ‘The Middle-Pleistocene Transition (MPT) represents a significant\n538 climate change over the Quaternary, but its cause remains not fully understood.’.\n539\n540 2.5. line 47: this sentence is too long and awkward. Consider breaking it apart.\n541 Our reply: Thank you. We have broken it into two separate sentences as follows:\n542 Internal feedback mechanisms have been suggested as essential for transitioning to\n543 longer glacial cycles (Berger et al., 1999; Clark et al., 2006; Bintanja and van de Wal,\n544 2008; Willeit et al., 2019). These mechanisms include the long-term decline in\n545 temperature due to pCO decrease and the regolith removal beneath ice sheets (Clark\n2\n546 et al., 2006; Willeit et al., 2019), which enable the excessive development of ice sheets,\n547 thus gradually extending the glacial cycles into ~100-kyr periodicities (Legrain et al.,\n548 2023).\n549\n550 2.6. line 52: “phase-locking.”\n551 Our reply: Thanks. We have changed the ‘phased-locking’ to ‘phase-locking’.\n552\n553 2.7. line 68: remove “pioneeringly.”\n554 Our reply: Thanks. We have removed this adverb, and it has been deleted throughout\n555 the revised manuscript.\n556\n557 2.8. line 263: remove “both.”\n558 Our reply: We have removed this preposition.\n559\n560 2.9. line 295: linear/nonlinear what? Maybe add “system” after “nonlinear” or another\n561 word as the authors see fit.\n562 Our reply: Thanks for this reminder. We have added the word ‘mode’ after the\n563 ‘nonlinear’.\n564\n565 Response to Reviews\n566 Reviewer #3:\n567 (General comments) Zhang et al. provide new and compelling data analyses which\n568 show that the 100-kyr periodicity in Late Pleistocene climate archives likely derived\n569 from a change in the climate sensitivity to precession forcing prior to the Mid-\n570 Pleistocene transition (MPT). According to the authors, the link between precession and\n571 ice sheets is linked to the growth of ice sheets consecutively reaching the lower latitudes,\n572 where precession forcing becomes more and more influential on the glacial mass\n573 balance. This is a new and exciting discovery which will advance the field of both data\n574 and model paleo-science. I am convinced by the author’s data showing the change in\n575 sensitivity of climate to the precession forcing, and the subsequent hypothesis about\n576 gradual southward expansion of ice sheets against the backdrop of long-term CO\n2\n577 decline. (1) However, I currently disagree with the arguments in the last part of the\n578 manuscript which identify the onset of the MPT, and try to explain this mechanism to\n579 not just be a driver, but an instigator, of the MPT. In my opinion, shortening this section,\n580 and focusing on the manuscript’s main data analyses will elevate the manuscript and\n581 highlight the new findings that the analyses identify. (2) Additionally, the manuscript\n582 would benefit from including comparison of their findings and methods to other,\n583 already published, studies, such as work created by Huybers et al., or Berger et al. who\n584 either promote obliquity to control glacial-interglacial cycle development across the\n585 Pleistocene (both 41 and 100kyr worlds), or find precession to influence either the\n586 41kyr or 100kyr-driven climate cycles. It would be very interesting to read the opinion\n587 of the authors on those findings, and how their own results fit into the already existing\n588 narratives. (3) A final minor point, as I am not a specialist on the field of spectral\n589 analyses, but on data collection and MPT mechanisms, I sometimes found it difficult to\n590 follow the terminology in the current manuscript. I therefore included a list of\n591 suggestions below that hopefully help to clarify some questions and make the text more\n592 accessible to a wider audience. This includes minor changes, such as including\n593 definitions for certain spectral analysis vocabulary, and adding mathematical\n594 expressions of these definitions in the supplement, as well as adapting the figures for\n595 easier access.\n596\n597 Major points:\n598 3.1. However, I currently disagree with the arguments in the last part of the manuscript\n599 which identify the onset of the MPT, and try to explain this mechanism to not just be a\n600 driver, but an instigator, of the MPT. In my opinion, shortening this section, and\n601 focusing on the manuscript’s main data analyses will elevate the manuscript and\n602 highlight the new findings that the analyses identify.\n603 Our reply: Thanks for your constructive comment. In the revised manuscript, we have\n604 reorganized and merged the second and third parts of the Discussion into a single\n605 section (the second section of the Discussion in the revised manuscript). We no longer\n606 discuss the timing of the onset of the MPT; instead, we propose that the difference in\n607 periodicity before and after the MPT indicates a shift from a linear to a nonlinear glacial\n608 response to precession. This shift in response mode to precession is closely linked to\n609 the cooling-induced increase in deglacial threshold. Additionally, we discuss the links\n610 between cooling, the regolith removal beneath ice sheets, as well as the decline in CO\n2\n611 level.\n612\n613 3.2. Additionally, the manuscript would benefit from including comparison of their\n614 findings and methods to other, already published, studies, such as work created by\n615 Huybers et al., or Berger et al. who either promote obliquity to control glacial-\n616 interglacial cycle development across the Pleistocene (both 41 and 100kyr worlds), or\n617 find precession to influence either the 41kyr or 100kyr-driven climate cycles. It would\n618 be very interesting to read the opinion of the authors on those findings, and how their\n619 own results fit into the already existing narratives.\n620 Our reply: Thanks for your suggestion. We have compared our results and mains\n621 findings to previously published work, such as Berger (1978), Hagelberg et al (1991),\n622 Huybers& Wunsch (2005), Berger and Loutre (2010), Feng and Bailer-Jones (2015). In\n623 the revised manuscript, we wrote as follows:\n624 Our statistical analysis of PDA results indicates that the amplified ~100-kyr glacial\n625 power is most likely derived from precession, with a less likely an obliquity\n626 contribution. Nevertheless, given that the main mass of ice sheets remains in polar\n627 regions where daily insolation variation is dominated by obliquity variance, obliquity\n628 should persistently exert a strong influence on ice sheets (Huybers& Wunsch, 2005;\n629 Drysdale et al., 2009; Huybers, 2011; Feng & Bailer-Jones, 2015; Bajo et al., 2020).\n630 The spectral power transfer between precession and ~100 kyr cycles may stem\n631 from the intrinsic relationship between precession and eccentricity, as all the periods of\n632 eccentricity are related to precessional ones via the heterodyne mechanism (Berger,\n633 1978). The nonlinear resonance between the frequencies of two main precession cycles,\n634 producing a beat at ~100-kyr frequency, can amplify the ~100-kyr glacial power\n635 (Berger, 1977; Hagelberg et al., 1991; Berger and Loutre, 2010). However, it is difficult\n636 to see the physical significance of such resonance (Imbrie et al., 1993).\n637 Moreover, we also compared our main findings to the work by Huybers (2006) and\n638 Raymo et al (2006). Our results are consistent with their conclusions (please refer to\n639 the last part of the discussion).\n640\n641 3.3. A final minor point, as I am not a specialist on the field of spectral analyses, but on\n642 data collection and MPT mechanisms, I sometimes found it difficult to follow the\n643 terminology in the current manuscript. I therefore included a list of suggestions below\n644 that hopefully help to clarify some questions and make the text more accessible to a\n645 wider audience. This includes minor changes, such as including definitions for certain\n646 spectral analysis vocabulary, and adding mathematical expressions of these definitions\n647 in the supplement, as well as adapting the figures for easier access.\n648 Our reply: Thank you. In terms of spectral power ratio, it serves as a proxy for\n649 estimating spectral strength within a specific frequency band. We have included this\n650 clarification in the last paragraph of the Introduction of the revised manuscript.\n651 Regarding the obliquity and precession phases, they refer to the absolute age offsets\n652 between the timing of the local maximal change rate of the benthic δ18O stack within a\n653 given glacial cycle and the timing at the nearest precession minimum or obliquity\n654 maximum. This definition has now been supplemented in the last paragraph of the\n655 Introduction and in the Method section. Moreover, the definitions for the precession\n656 and obliquity phases are now expressed with two mathematical equations (see the\n657 Method section).\n658\n659 Minor points:\n660 3.4. line 1: As I currently cannot support the final chapter of the manuscript (see\n661 comments below), I suggest to change the title to “CO2-forced change in glacial\n662 response to precession drives the Middle-Pleistocene Transition”\n663 Our reply: Thank you. Given the lack of a continuous CO reconstruction over the last\n2\n664 2700 ka, the title of the manuscript has been revised to ‘Cooling-induced shift in glacial\n665 response to precession drives the Middle-Pleistocene Transition’. This revision aims to\n666 mitigate unnecessary debate and better align with our results because our estimated\n667 precession phase displays a generally consistent trend with the stacked global mean\n668 surface temperature (Fig. 3C&D; Fig. 5B&F) (Clark et al., 2024).\n669\n670 3.5. line 38: are there any references that confirm the frequency of obliquity, precession,\n671 and eccentricity for a wider audience?\n672 Our reply: The original paper (Berger, 1978) for these frequencies is now cited.\n673\n674 3.6. line 52: Other references state that the 100 kyr cycle is possibly not phase locked\n675 to eccentricity, but expresses multitides of precession and obliquity with internal\n676 feedback phase locking (e.g. Berger & Loutre, 2010, Global & Planetary change;\n677 Huybers & Wunsch, 2005, nature).\n678 Our reply: Thank you. These references are included in our revised introduction.\n679 The major contribution of paper by Huybers and Wunsch (2005, nature) is their\n680 validation of the obliquity influence on late Pleistocene glacial cycles through a\n681 statistical test analysis. They proposed that the climate state of skipping one or two\n682 obliquity beats before deglaciation can yield quantized glacial-cycle durations of either\n683 80 or 120 kyr through either a stochastic or deterministic process. Berger and Loutre\n684 (2010) proposed that the ~100-kyr glacial cycles are related to a non-linear\n685 amplification of the astronomical forcing involving a beat effect from the two main\n686 precession periodicities.\n687 Now we state that “The persistent influence of precession and obliquity on glacial\n688 cycles (Berger et al., 1999; Huybers and Wunsch, 2005; Drysdale et al., 2009; Cheng\n689 et al., 2009; Huybers, 2011; Cheng et al., 2016; Bajo et al., 2020; Liautaud et al., 2020;\n690 Barker et al., 2022; Hobart et al., 2023), along with the phase-locking between ~100-\n691 kyr glacial cycles and eccentricity (Imbrie et al., 1993; Tziperman et al., 2006; Lisiecki,\n692 2010), and the non-linear amplification of the astronomical forcing related to\n693 eccentricity and precession (Berger and Loutre, 2010) suggest that astronomically-\n694 forced variations in insolation must be indispensable for the glacial cycles\n695 (Milankovitch, 1941; Hays et al., 1976)”.\n696\n697 3.7. line 63: It would be helpful to include definitions of the main concepts that will be\n698 used throughout the manuscript. These include spectral power, spectral phase, spectral\n699 power ratio, power shift. It would be supportive to have mathematical expressions of\n700 those definitions in the supplement as well.\n701 Our reply: We have included additional definitions for spectral power and\n702 obliquity/precession phase in the revised manuscript. For more details, please refer to\n703 our reply to comment 3.3 above.\n704\n705 3.8. line 66: Briefly describe what the 57 paleoclimate archives are for a wider audience.\n706 I saw this information in the supplement, but briefly include why these 57 paleoclimate\n707 archives were chosen.\n708 Our reply: In the revised introduction of the manuscript, we have included the\n709 following sentence ‘These datasets, globally derived from ocean sediments, loess\n710 deposits, ice cores, and other archives (Fig. 1), are all continuous, with well-constrained\n711 age models (Table S1) and possess a high temporal resolution with at least a 2-kyr\n712 sampling rate.’\n713\n714 3.9. line 82: The authors state in line 82 that the relationship between the 21kyr power\n715 ratio and the 100kyr power ratio holds during and before the MPT. However, the cross-\n716 plot in Figure 2 panel A has the highest r of 0.7. That surely shows that the relationship\n717 holds throughout? Line 83: Can the authors comment on why their results show no link\n718 with obliquity, against what published studies have found (e.g. Feng & Bailer-Jones,\n719 2015, QSR).\n720 Our reply: In the revised manuscript, we discussed several scenarios to assess the\n721 robustness of the relationship between spectral power ratios of ~21-kyr, ~41-kyr, and\n722 ~100-kyr components (Fig. S10). These scenarios considered the potential influence of\n723 significance level of spectral power peaks, dataset length, age models used, and orbital\n724 tuning on the PDA results and their relationship. The results indicate that across all\n7", "reviewer_3": "Overall remarks\nZhang et al. provide new and compelling data analyses which show that the 100-kyr periodicity in Late Pleistocene climate archives likely derived from a change in the climate sensitivity to precession forcing prior to the Mid-Pleistocene transition (MPT). According to the authors, the link between precession and ice sheets is linked to the growth of ice sheets consecutively reaching the lower latitudes, where precession forcing becomes more and more influential on the glacial mass balance. This is a new and exciting discovery which will advance the field of both data and model paleo-science. I am convinced by the author’s data showing the change in sensitivity of climate to the precession forcing, and the subsequent hypothesis about gradual southward expansion of ice sheets against the backdrop of long-term CO2 decline. However, I currently disagree with the arguments in the last part of the manuscript which identify the onset of the MPT, and try to explain this mechanism to not just be a driver, but an instigator, of the MPT. In my opinion, shortening this section, and focusing on the manuscript’s main data analyses will elevate the manuscript and highlight the new findings that the analyses identify.\nAdditionally, the manuscript would benefit from including comparison of their findings and methods to other, already published, studies, such as work created by Huybers et al., or Berger et al. who either promote obliquity to control glacial-interglacial cycle development across the Pleistocene (both 41 and 100kyr worlds), or find precession to influence either the 41kyr or 100kyr-driven climate cycles. It would be very interesting to read the opinion of the authors on those findings, and how their own results fit into the already existing narratives. A final minor point, as I am not a specialist on the field of spectral analyses, but on data collection and MPT mechanisms, I sometimes found it difficult to follow the terminology in the current manuscript. I therefore included a list of suggestions below that hopefully help to clarify some questions and make the text more accessible to a wider audience. This includes minor changes, such as including definitions for certain spectral analysis vocabulary, and adding mathematical expressions of these definitions in the supplement, as well as adapting the figures for easier access.\nLine by line comments\nMain text:\nLine 1: As I currently cannot support the final chapter of the manuscript (see comments below), I suggest to change the title to “CO2-forced change in glacial response to precession drives the Middle-Pleistocene Transition”\nLine 38: are there any references that confirm the frequency of obliquity, precession, and eccentricity for a wider audience?\nLine 52: Other references state that the 100kyr cycle is possibly not phase locked to eccentricity, but expresses multitides of precession and obliquity with internal feedback phase locking (e.g. Berger & Loutre, 2010, Global & Planetary change; Huybers & Wunsch, 2015, nature).\nLine 63: It would be helpful to include definitions of the main concepts that will be used throughout the manuscript. These include spectral power, spectral phase, spectral power ratio, power shift. It would be supportive to have mathematical expressions of those definitions in the supplement as well.\nLine 66: Briefly describe what the 57 paleoclimate archives are for a wider audience. I saw this information in the supplement, but briefly include why these 57 paleoclimate archives were chosen.\nLine 82: The authors state in line 82 that the relationship between the 21kyr power ratio and the 100kyr power ratio holds during and before the MPT. However, the cross-plot in Figure 2 panel A has the highest r of 0.7. That surely shows that the relationship holds throughout?\nLine 83: Can the authors comment on why their results show no link with obliquity, against what published studies have found (e.g. Feng & Bailer-Jones, 2015, QSR)\nLine 94: can the authors confirm why the analysis was conducted on LR04, and not the 57 climate records which are plotted in Fig. 1?\nFigure 3:\na) Can the authors include the mathematical expression in the supplement or reference the software, of how the derivative was calculated?\nb) Also I suggest to change the y-axes labels to something more intuitive for wider audience accessibility, especially for panels C, D, E, F, K. I suggest the authors to consider a re-organisation of the panels. The small cross-plots (H, I, J) on the right invoke connection to panels E, F, G, which I believe is not the intention of the authors?\nc) Considering B: can the authors comment on why this sea level record was chosen? There are other sea level records that could be used (e.g. Wang et al., 2010; de Boer et al., 2014, Sosdian et al. 2009, and others).\nd) Have the authors considered looking at other CO2 records, such as Yamamoto et al., 2022? If so, why was the record in G chosen?\ne) drawing the vertical dotted lines through the entire graph might be more intuitive for the reader\nLine 186: I would be interested to see the author’s comment on the following publication here: Huybers et al., 2006, Science.\nLine 200: It might be worth mentioning that ice sheets do not only drive internal variability, but can be affected by other internal feedback mechanisms themselves.\nFigure 4. In general, I did not find this figure very intuitive and it took me some time to understand what the authors are trying to express. I kindly listed my issues and questions below.\na) Panel A: a wider audience will not be familiar with the term 1E1T-P. I assume this refers to the sum of all power spectra defined by the authors as Precession, Obliquity, Ecentricity. If so, I suggest to clarify this in the figure caption. I also suggest a small comment in the supplement, why this definition is chosen, and where it comes from. Maybe it is a common reference, but clarifying that would be useful.\nb) Since the term power ratio is not necessarily well known to a wider audience, it would be useful to either have a clear definition in the text, or re-word the definition very brief in the figure or y-axis figure caption. Sometimes annotations may also help.\nc) Is there a plot that shows the behaviour of obliquity in comparison to the plot in panel C? Also, for clarification, I assume that panel C is a cross-plot of panel B, is that correct?\nd) Is the threshold in panel D just an example of a threshold, given that panel H and I show multiple thresholds on the x-axis?\ne) I wonder whether it would be helpful to include a plot comparable to panel G, where each threshold’s power spectrum is plotted in the same plot, maybe using grey scales. Then the reader could see that changing the threshold will cause the appearance of the 95, 125, and 405 spectra.\nf) The jump from panel D to panel F is counter-intuitive, since the y-axis are flipped on panel F. This needs to be clearly stated, or highlighted on the y-axis.\ng) Using more descriptive titles on panels H, I, J would be useful.\nh) The figure caption must include a description of what the whisker plots represent. I assume they represent the amplitude of the spectrum at each threshold? What are the error bars on the whisker plots?\ni) In general, figure 4 is central to the hypothesis of the authors and shows the link between changes in climate sensitivity. Next to figure 2, it seems to be the most important figure. As such, I recommend the authors to consider re-organising the panels, or re-writing the figure caption, to help the general reader follow this information-heavy figure more easily.\nLine 233: Can the authors comment on why they do not find any influence of obliquity on climate in comparison to other studies? (e.g. Huybers & Wunsch, 2015; Feng & Bailer-Jones, 2015; Berger et al., 2005, Paleoceanography). How do the author’s results compare to other studies finding precession to be a dominant driver (e.g. Liautaud et al., 2020, EPSL)?\nLine 234: I was unable to follow the statement made in this sentence. Can the authors clarify why they claim that “the nonlinear amplification of 100kyr power is an intrinsic characteristic of the entire climate system, not just of ice volume?\nLine 237: can the authors elaborate why multiple thresholds are needed to induce sensitivity change, and not just one threshold?\nLine 240: I believe this analysis was conducted on data from LR04 correct? As such, why do the authors use LR04 as a proxy for global climate instead of ice volume, even though they have a perfect collection of 57 cores that are better at representing global climate than LR04?\nLine 244: I cannot quite follow the argument in the following lines. Are the authors suggesting that the 100kyr cycle in LR04 occurs due to the 400kyr power in eccentricity, and therefore disagree with a collection of publications claiming that the 100kyr is a “quasi 100kyr” cycle modulated by precession and obliquity due to internal feedbacks? If so, this would be a new suggestion, and requires more in-depth discussion.\nLine 262: At the moment, I cannot agree with this statement, since the authors receive different timings depending on whether their analysis uses LR04, or the probabilistic stack. Previously in the manuscript the authors argue that the use of the probabilistic stack is better, due to a more robust age model. I suggest that the authors stick with this logical approach, and present data from the probabilistic stack. The fact that all the other analyses show no sensitivity to different age models allows me to believe their results. In this case, I believe that the analysis shows sensitivity to age modelling, which in itself is an interesting result, but should not be over interpreted in terms of climate dynamics.\nLine 259 and following section: In the following lines the authors present a hypothesis for how precession may drive the MPT. At this point, I am unable to agree with the hypothesis, due to (i) a lack of robust data which shows that the shift in sensitivity to the precession cycle, and re-location of power to 100kyr cycle is due to long-term reduction in pCO2 (see line 274). (ii) a mis-representation of reference 7 (see line 271). Reference 7 reviews processes that influence the MPT, but does not state significant hypotheses on how the MPT was instigated. In contrast, ref 7 proposes a collection of research questions which would be beneficial for the field to solve. (iii) a misrepresentation of refs 16, 45 (see line 278) which consider the late Pleistocene, not the early Pleistocene prior to the MPT. As such, their evidence needs to be used after large-scale ice sheets are already occurring (see line 281). (iv) the statement that “bedrock-mass-balance feedback and increased basal-sliding influences ice sheet dynamics” would benefit from more than ref 9, as this process has been suggested to influence ice sheets on multiple timescales, ie during deglaciations, or as an instigator for the MPT (see Regolith hypothesis, Clark et al., 2006). (v) a misrepresentation of ref 11 which only tested the change in long-term CO2 coupled to regolith removal, vs. only regolith removal. It did not test whether abrupt or gradual changes in CO2 cause the MPT.\nLine 294: This statement is currently unsupported due to the above comments.\nLine 297: Given that the probabilistic stack analyses of the authors show an onset at 1.2Ma, the literature seems to support the findings of the authors with respect to the probabilistic stack. Although I would suggest not putting too much influence on this analysis, as it seems highly variable to changes in age models, as the authors admit themselves in line 300.\nLine 301: This sentence might need reworking.\nLine 303: “High-amplitude” glacial cycles are generally believed to occur after 800kyr as they describe the 100kyr world glacial cycles, which occur after the MPT. Are the authors stating that they suggest an alternative idea?\nLine 305: ref 47 needs to be discussed in line with more recent publications on this topic (some already listed throughout the review).\nLine 308: At this point, I cannot agree with the conclusion that the MPT is a gradual process. The authors do not provide robust evidence for that, and the list of arguments is lacking. In fact, I am unsure whether the authors are able to supply this information with the data they have at this stage. I therefore suggest to shorten the last section, which should not take away from the novelty of the work presented before.\nFigure 5:\na) The timing of 1452kyr seems dependent on the age model used which makes this graph hard to accept as robust\nb) Improvement of the y-axis label would help with the understanding of the graph.\nc) To help a wider audience, it might be useful to guide the reader through the cross plots and explain what the change in these patterns means for global climate. The correlations are compelling and I believe the analysis should focus more on the change in the correlation direction, which stays robust throughout different age models, rather than focusing on the timing of the MPT.\nd) In the cross plot in panel D, E, the yellow line is annotated with p=0.86. Should this be r=0.86 or does it refer to a p-value? If so, maybe also report r.\nMethods:\nLine 439: Can this be expressed as a mathematical equation?\nLine 442: It would be useful to include a detailed and clear definition for the different terms used throughout the manuscript, such as spectral power, spectral strength, spectral power ratio, spectral phase etc\nLine 444: I wonder how different the age models are between these different approaches. Given the recent improvement in age modelling, the age models tested here might not vary greatly which could lead to internal age model errors per age model being larger than between- age model errors. Could the authors include a plot with one record on all age models so the variability can be compared? Alternatively, the authors could consider making a plot like S1 or S2, but using the maximum error values for one analysis and comparing to the minimum error values from the same age model from a second analysis. This should also show the influence of internal age model errors on the analyses.\nLine 458: What do the authors mean by “robust” age models?\nLine 476: Can this be expressed as a mathematical equation?\nLine 480 and following: Some of these definitions should be included in the main text to help a wider audience understand the analyses more easily."}, "author_response": "1 Response Letter\n2\n5 (General comments) This paper tackles a fundamental problem in understanding ice\n6 age behavior in the Pleistocene, which records several mysterious changes in behavior.\n8 age cycles, but no linear model comes close to capturing the geological observations.\n10 published archives of climate change.\n12 briefly major issues which I urge the authors to evaluate.\n13\n15 1.1. Although the role of CO is invoked frequently, high resolution and replicated\n2\n16 estimates are only available for the past ~800,000 years. The record is being extended,\n17 most successfully in my opinion using the δ11B proxy, but is incomplete. Critically, the\n18 proxy data do not support the idea of a gradual drawdown of CO at around the time of\n2\n19 the MPT. Neither do temperature records show such a behavior. I do not see citations\n20 to these CO or SST reconstructions in the manuscript. Yet this argument of interglacial\n2\n2\n22 subsequent glacial cycle is central to the interpretation here. The CO estimates used\n2\n23 here consist of model results, using the δ18O proxy and cannot be considered as\n24 independent of the isotopic data. They have not been validated with actual palo-CO\n2\n25 estimates\n27 incorporated the reconstructed CO records, including those based on δ11B (Foster,\n2\n28 2008; Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011; Henehan et al., 2013;\n30 (Yan et al., 2019) and EPICA ice core (0-800ka; Bereiter et al., 2015), and a continuous\n2\n32 (Yamamoto et al., 2022) (Fig. 5A). Furthermore, we included the new stack of global\n33 mean surface temperature (GMST) (Clark et al., 2024) (Fig. 5B).\n2\n2\n36 et al., 2017; Dyez et al., 2018), ice core (Bereiter et al., 2015; Yan et al., 2019), and\n37 leaf-wax δ13C proxies (Yamamoto et al., 2022) (Fig. 5A). As glacial-stage temperature\n38 determines the maximal extension of ice sheets, the glacial-stage CO decline could be\n2\n39 critical for understanding the changes in ice sheet dynamics across the MPT.\n41 2024) clearly demonstrate a long-term cooling trend across the MPT (Fig. 3D; Fig. 5B).\n43 findings strengthen our interpretation of the significant role that background\n44 temperature changes play in modulating the glacial response to precession.\n45\n46 1.2. The authors assume that the isotope curve linearly reflects ice volume. Work over\n47 the last decade emphasizes the fact that there may be a significant, and evolving\n48 overprint of deep ocean temperature. Without separating temperature and ice volume,\n49 the interpretations presented here, which focus on the behavior of ice sheets, is\n50 incomplete and possibly incorrect. I’d urge them to look at the recent paper by Clark et\n51 al. in Science (2024) to see how much room for debate there is in interpreting what\n52 d18O is actually telling us.\n54 by local deep-sea temperature, but it remains to be a valuable proxy for discussing\n55 global ice volume changes as has been done and continues to be done in numerous\n56 studies (e.g. Huybers and Wunsch, 2005; Huybers, 2006; Drysdale et al., 2009; Lisiecki,\n57 2010; Huybers, 2011; Tzedakis et al., 2017; Liautaud et al., 2020; Hobart et al., 2023),\n58 especially when many δ18O records from different oceans are used like in our study,\n59 which is expected to mitigate the effect of local deep-sea temperature and provide a\n60 more accurate reflection of global ice volume changes.\n62 by Tzedakis et al (2017), which is based on the detrended LR04 stack. The purpose of\n63 detrending is to remove the temperature contribution to benthic δ18O (Tzedakis et al.,\n64 2017). As depicted in Fig. A1D, the detrended LR04 stack closely aligns with the\n65 separate ice volume proxy (Elderfield et al., 2012). Thus, the estimated precession and\n66 obliquity phases are minimally affected by temperature influences on benthic δ18O.\n68 validation of our results because this record represents the contribution of ice volume\n69 to benthic δ18O (Elderfield et al., 2012). As shown in Fig. A1, the estimated precession\n70 and obliquity phases display similar trends to those derived from the LR04 stack (Fig.\n72 negative correlation with the ~100-kyr spectral power ratio after the MPT (r=-0.89),\n73 whereas their relationship is positive before the MPT (r=0.66) (Fig. A2E). This suggests\n74 a shift in the response of ice sheets to precession from early linearity to later nonlinearity,\n75 which is consistent with the results acquired from the H07&NA stack (Fig. 5E, G). Thus,\n76 our estimated precession and obliquity phases are minimally affected by temperature\n77 influences on benthic δ18O.\n78\n80 of δ18Ow in Site ODP1123 (Elderfield et al., 2012). A, the detrended LR04 benthic δ18O stack\n81 (Lisiecki and Raymo, 2005; Tzedakis et al., 2017), the red and green dots denote the glacial minima\n82 (or interglacial maxima) and maxima within glacial cycles, the orange dotes refer to the timing of\n83 the local maximal changing rate of benthic oxygen isotope within a given glacial cycle. The black\n84 numbers above the δ18O stack mark the marine isotope stages (MIS). The orange tags and numbers\n85 bellow the δ18O stack label the Termination numbers (T) and corresponding ages. The horizontal\n86 dotted lines are the threshold used to define the glacial cycles (Tzedakis et al., 2017). B&C, the\n87 precession and obliquity phase during the maximal changing rate of LR04 benthic δ18O stack, as\n88 well as their corresponding 5-points moving average (5pma; in blue). D, the δ18Ow in Site ODP1123,\n89 which represents the ice volume component extracted from the benthic δ18O (Elderfield et al., 2012).\n91 the orange dotes refer to the timing of the local maximal changing rate of δ18Ow. The rosters of\n92 glacial cycles refer to scheme provided by Tzedakis et al. (2017). E, the derivative of δ18Ow, the\n93 orange dots label the maximal changing rate of δ18Ow within glacial cycles. F&G, the precession\n94 and obliquity parameters (Berger and Loutre, 1991). H&I, the precession and obliquity phase during\n95 the maximal changing rate of δ18Ow (in thin black curve), as well as their corresponding 5-points\n96 moving average (5pma; in thick black curve).\n97\n98\n100 points moving average of the estimated precession phase (C) based on δ18Ow in Site ODP1123\n101 (Elderfield et al., 2012). The correlation plots among these variabilities were also created for each\n102 pair (D, E, and F).\n103\n104 1.3. The paper does not consider issues of age models and sampling resolution, both of\n105 which can be critical to isolating the precessional component in geological time series.\n107 accurately capturing precessional signals is quite difficult. A particular problem has to\n108 do with the potential circularity in assessing the phasing of precession relative to\n109 obliquity in records tuned to obliquity. As an example, the IODP 1385 age model\n110 referred to was “retuned” from the standard isotopic age model using a precessional\n111 signal by Liautaud et al. This retuning significantly increases the amplitude of the\n112 precessional signal, although it comes with its own assumptions.\n114 resolution is higher than 2 samples per kyr. According to the Nyquist frequency, this\n115 resolution is sufficient to detect the precession signal. Moreover, in the revised\n116 manuscript, we tested the confidence level (CL) of the power peaks at frequencies of\n117 ~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1 within the datasets (Figs. S6-S9). Only the\n118 power ratios with corresponding power peaks above 90% CL were considered in\n119 assessing relationship among spectral power ratios (Fig. S10A, B).\n121 datasets spanning the past 800 ka and discussed six different scenarios regarding the\n122 influence of significant level of power peaks (Fig. S10A, B), datasets length (Fig. S10C),\n123 age models used (Fig. S10D), as well as the effects of orbital tuning (Fig. S10E, F) on\n124 the correlation of power ratios at ~21-kyr, ~41-kyr. and ~100-kyr frequency bands.\n126 constrained by the LR04 age model to the H07&NA age model (Fig. 10E, F; Figs. S12).\n128 last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack) (Hobart et al.,\n129 2023) (Fig. S11; Table S2). The age model for the NA stack was based on the\n130 speleothem chronology (Cheng et al., 2016), making the H07&NA stack independent\n131 of orbital tuning. Across all scenarios, the ~21-kyr power ratios negatively correlate\n132 with the ~100-kyr power ratios, while there is no significant correlation between the\n133 ~41-kyr and ~100-kyr power ratios (Fig. S10), consistent with our previous findings.\n135 precession and obliquity. To avoid this, in the revised manuscript, we recalculated the\n136 precession and obliquity phases of glacial terminations based on the H07&NA stack,\n137 featuring an age model independent of orbital tuning. The orbital phases calculated\n138 from this stack align with the evolutionary trends obtained from the LR04 stack and\n140\n141 1.4. The time series methods need to be better described. How do the authors treat non-\n142 stationarity? With moving windows? Given that precession is rarely resolved well in\n143 spectra of “real” data, how do they estimate uncertainties?\n145 analysis (PDA) (Li et al., 2019). This method allows for the quantitative estimation of\n146 the influence of precession, obliquity, and ~100-kyr cycles on paleoclimate archives.\n148 2020), this non-stationary technique, PDA, employs a sliding time window to calculate\n149 evolutionary power ratios within specific frequency bands (Fig. S1). For the\n150 paleoclimate datasets spanning 0–800 ka, 800–1200 ka, and 1200–2700 ka, the sliding\n151 windows were set at 500 kyr, 200 kyr, and 500 kyr, respectively.\n153 method (see Method section) and rewrote the first section of the Method to improve\n154 clarity on the principles and processing steps of PDA. The revised description is as\n155 follows: (1) Initially, the evolutionary power spectra of the datasets were calculated\n156 using the 2π multi-taper method (MTM) (Thomson, 1982) with a sliding time window\n157 (Meyers et al., 2012) (Fig. S1B). (2) Next, the time-dependent amplitude modulations\n158 (power) at these three frequency bands (~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1) were\n159 quantitatively calculated by integrating the evolutionary power at the corresponding\n160 frequency bands (Meyers et al., 2012; Li et al., 2019). (3) Thereafter, to facilitate\n161 comparisons among different paleoclimate datasets, the integrated evolutionary spectral\n162 power at each specific band was normalized by the total evolutionary power (integrated\n163 from 0 to the Nyquist frequency) to obtain the evolution spectral power ratio. This\n164 process is necessary because different paleoclimate proxies have varying units.\n166 precession and obliquity frequency bands for three different benthic δ18O stacks were\n167 examined. These include the LR04 stack (Lisiecki and Raymo, 2005), the Prob_stack\n168 (Ahn et al., 2017), and the H07&NA stack (Huybers, 2007) (Fig. S1D, F). The\n169 calculated spectral power ratios from these stack exhibit similar evolution trends,\n170 consistent with the results obtained through the Empirical Nonlinear Orbital Fitting\n171 (ENOF) (Fig. S1C, E) (Liautaud et al., 2020).\n172\n173 1.5. And finally, how the “phase” is derived robustly is very difficult to determine, yet\n174 central to the conclusions.\n177 precession and obliquity phases (Fig. 3B, C; Figs. S14-16); Secondly, to minimize the\n178 influence of errors from individual precession or obliquity phases at specific glacial\n179 terminations on the overall evolution trend of phases, we applied a 5-point moving\n180 average (5 pma). Thirdly, the H07&NA stack was employed to estimate the\n181 astronomical phases of glacial terminations, therefore mitigating the circular reasoning.\n183 different benthic δ18O stacks, display consistent evolution trends, demonstrating the\n184 robustness of our phase results.\n185\n187 1.6. line 23: repetition of \"transition\"\n189 and eliminated redundancy.\n190\n191 1.7. line 24-30: since there are no continuous, orbitally resolved pCO measurements >\n2\n192 800 ka, this can't be supported as written. Note also that the δ18O is equated with ice\n193 volume- this is increasingly questioned with recent work, e.g. Clark et al. 2024. Also it\n194 isn't clear what \"phase covaries with pCO … dynamics...\"- means concretely in terms\n2\n195 of the data and/or statistical analysis\n2\n197 with precession phase. We have reorganized the Abstract and emphasized the role of\n198 gradually-cooling temperature in modulating glacial response to precession (Fig. 5B).\n200 mean surface temperatures, facilitating a direct correlation with reconstructed\n201 precession phase.\n2\n2\n204 across the MPT (Fig. 5A). Therefore, the decline in glacial-stage pCO (Yan et al., 2019)\n2\n205 suggests a potential relationship with the temperature decline observed across the MPT.\n208 response to the second major comment (1.2).\n209\n210 1.8. Line 52: signal usually refers to a representation in the earth system (e.g. proxy\n211 evidence) not to a pacing/forcing such as orbital cycles.\n213 ‘signals.’\n214\n215 1.9. Line 59: without clarification, the reader can't follow what is meant by \"only\n216 specific instances\"\n218 ‘Introduction’, and this sentence has been deleted.\n219\n220 1.10. Line 70-73: since any CO reconstruction at this point is indirect, I wouldn't use\n2\n221 this phrase.\n223 follows: These evolutionary phases enable a comparison with previously simulated and\n224 reconstructed pCO (Foster, 2008; Hönisch et al., 2009; Seki et al., 2010; Bartoli et al.,\n2\n225 2011; Henehan et al., 2013; Bereiter et al., 2015; Chalk et al., 2017; Dyez et al., 2018;\n227 temperature (ΔGMST) (Clark et al., 2024), integrated glacial latitudinal extension\n228 (Köhler and van de Wal, 2020), and ice volume (de Boer et al., 2014; Willeit et al., 2019;\n230 successive correlation between internal and external factors.\n231\n232 1.11. Line 73-74: this sentence (These analyses contribute to piecing together a jigsaw\n233 for understanding the MPT and the ~100-kyr glacial cycles.) does not contain any useful\n234 information.\n236\n237 1.12. Line 78-83: some minimal discussion of the spectral technique is needed- is some\n238 kind of evolutive approach used to capture the time-dependent spectral characteristics\n239 (e.g. non-stationarity)? Also, since the records in many cases are not of the same length,\n240 how is this heterogeneity handled?\n242 with additional details. In the revised manuscript, the Result section begins with ‘The\n243 power decomposition analysis (PDA) (Li et al., 2019-ESR) was employed to\n244 quantitatively estimate the influence of precession, obliquity, and ~100-kyr cycles on\n245 paleoclimate archives. Analogous to the Empirical Nonlinear Orbital Fitting (ENOF)\n246 method (Liautaud et al., 2020) …’ Moreover, we have included a supplementary\n247 figure (Fig. S1) to visually explain the principles and processing steps of the PDA method.\n250 longer than 600 kyr within the last 800 ka were considered (Fig. S10C). The analysis\n251 shows a negative correlation between the ~21-kyr and ~100-kyr power ratios, while no\n252 significant relationship is observed between the ~41-kyr and ~100-kyr power ratios.\n254 long as 400 kyr, their medians of ~21-kyr power ratios still show a strong positive\n255 correlation with those at ~100-kyr bands (Fig. 2F, G). Thus, we conclude that, despite\n256 the heterogeneity in dataset length, our findings remain robust and reliable.\n257\n258 1.13. Line 95-96: Prob stack is still not independent of an underlying orbital tuning for\n259 its timing, as shown by its similarity to LR04. It isn't clear how the authors avoid either\n260 the circularity in tuning of the isotope record or the inherent uncertainties in any age\n261 model. The 1385 hybrid model represents a different approach from LR04/prob stack\n262 in that it is tuned to a lithological precession component, so although tuned, it does use\n263 the isotope record in the tuning. Ref. 23 uses a similar approach. Note that neither the\n264 1385 or Ref. 23 approach extend significantly before the MPT.\n266 the LR04 age model as an initial condition during the stacking process, which resulted\n267 in our previous phase estimations being subject to circular reasoning. To address this\n268 issue in the revised manuscript, we used the depth-derived H07 stack (Huybers, 2007)\n269 with its last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack) (Hobart\n270 et al., 2023) (Fig. S11; Table S2). The age model for the NA stack was established based\n271 on the speleothem age model (Cheng et al., 2016). The tuned H07 stack, termed\n273 was used to recalculate the precession and obliquity phases of glacial terminations (Fig.\n275 obtained using the LR04 stack and Prob_stack (Fig. 3; Fig. S14–16), confirming our\n276 previous results.\n278 (Hodell et al., 2023) and the LR04 stack with Konijnendijk age model (Konijnendijk et\n279 al., 2015) have not been considered for calculating precession and obliquity phases\n280 because of their shorter length and their age models predominantly inheriting the LR04\n281 age model.\n282\n283 1.14. Line 111-118: given that the data sets represent different aspects of climate, should\n284 we expect similar power ratios and/or \"phases\". How do the authors handle the\n285 dispersion of these quantities shown in Figure 2- in other words, what is the confidence\n286 level of the estimates? I also find the figure visually confusing.\n288 collected datasets are different, it is important to note that the climate subsystems\n289 interact extensively. Thus, the relationship in spectral power ratios within these\n290 subsystems can propagate through the entire climate system via a series of positive and\n291 negative feedbacks (Clark et al., 1999; Denton et al., 2010). The expansion and retreat\n292 of Pleistocene ice sheets respond not only to external orbital forcing but also act as\n293 significant internal forcings, driving and influencing coeval regional and global\n294 climates. For example, they affect sea-level variations and atmospheric-oceanic\n295 circulations through the release of meltwater (Imbrie et al., 1993; Clark et al., 1999;\n299 the ~100-kyr components, retained in ice volume, can be projected onto other climate\n300 archives. Thus, regarding the internal nonlinear feedback mechanism that causes the\n301 transfer of power and subsequently amplifies the ~100-kyr cycles, regardless of\n302 whether it is a feature exclusive to the glacial subsystem or intrinsic to the entire climate\n303 system, an analysis of the spectral power ratio relationship between the ~100-kyr cycles,\n304 precession, and obliquity components within ice volume and other contemporaneous\n305 paleoclimate archives can determine whether the amplified power of ~100-kyr glacial\n306 cycles is derived from precession or obliquity.\n308 demonstrate a reversal correlation between power ratios at ~21-kyr and ~100-kyr bands\n309 in the climate system, rather than relying solely on benthic δ18O data. Moreover, we\n310 restricted our analysis to PDA results from the paleoclimate archives where spectral\n311 power at frequencies of ~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1 exceeding the 90%\n312 significance level. In terms of Fig. 2, in the revised manuscript, we have reorganized\n313 the presentation for clarity. For a better comparability of the power ratio changes across\n314 three intervals (0-800 ka; 800-1200 ka; 1200-2700 ka), we merged three panels into one\n315 and standardized their axis ranges. Furthermore, the scatter plots have been magnified\n316 for enhancing detail.\n317\n318 1.15. Line 121-125: the modelled CO in fact comes from an interpretation of the stable\n2\n319 isotope record as a forcing function, so it's hard to see how these are independent\n320 estimates of different aspects of earth climate. Note that there is no evidence of\n321 \"progressive decrease\" of CO from proxy evidence such as δ11B, if \"progressive\n2\n322 decrease\" means averaged over the period of several orbital cycles. The lack of trend in\n323 obliquity phase is likely to come from the assumptions of the age models that align the\n324 signals to obliquity.\n2\n326 reconstructions, including δ11B-based reconstruction (Chalk et al., 2017; Seki et al.,\n327 2010; Foster, 2008; Henehan et al., 2013; Hönisch et al., 2009; Dyez et al., 2018; Bartoli\n328 et al., 2011), reconstructions derived from Allan Hills Blue Ice (Yan et al., 2019) and\n2\n330 based on the leaf-wax δ13C proxy spanning the last 1500 ka (Yamamoto et al., 2022)\n331 (Fig. 5A).\n2\n333 the reconstructions based on δ11B and ice core archives consistently display a glacial-\n334 stage decline in CO across the MPT. Moreover, the continuous CO reconstruction\n2 2\n335 based on the leaf-wax δ13C proxy over the past 1500 ka also indicates a gradual decline\n336 in glacial-stage CO until the MPT (Fig. 5A). We also included the new stack of global\n2\n337 mean surface temperature (GMST) (Clark et al., 2024) (Fig. 5B), which confirms a\n338 long-term cooling trend across the MPT (Fig. 3D; Fig. 5B). Therefore, various lines of\n339 evidence show that there is a consistent decline in the global mean temperature and\n340 glacial-stage pCO across the MPT.\n2\n342 the assumptions of the age models that align the signals to obliquity.’\n344 on the LR04 stack, Prob_Stack, and H07&NA stack, revealing a similar evolutionary\n345 trend. The H07&NA stack is indenpend of orbtial tuning, confirming that the evolution\n346 trends of constructed precession and obliquity phases are not induced by the orbital\n347 tuning. The estimated obliquity phases show an increasing trend from 2700 ka to around\n348 800 ka, stabilizing thereafter (Fig. 3B). Nevertheless, this evolutionary trend is less\n349 prominent. In contrast, the precession phase follows an opposite trend, suggesting an\n350 increasing glacial sensitivity to precession (Feng and Bailer-Jones, 2015; Liautaud et\n351 al., 2020; Barker et al., 2022). Menwhile, obliquity’s role becomes comparably less\n352 significant, even though it remains indispensable for deglaciations (Feng and Bailer-\n354\n355 1.16. Line 135-137: It's not clear why \"phase\" is connected with sensitivity. The phase\n356 offset is usually a measure of the response time of the sytem relative to the forcing, and\n357 response time != sensitivity (or sensitivity needs to be defined so that there is a clear\n358 connection to the \"phase\")\n360 the absolute age offsets between the timing of glacial terminations and the timing at the\n361 nearest precession minima or obliquity maxima. Following Huybers (2011), a glacial\n362 termination is pinpointed at the maximal changing rate of benthic δ18O within a glacial\n363 cycle. Thus, the phase here estimates how quickly the ice sheets respond to precession-\n364 or obliquity-induced insolation changes. A smaller phase implies a shorter response\n365 time of ice sheets to precession or obliquity, indicating greater susceptibility of the ice\n366 sheets to these forcing mechanisms. Building on prior studies that discussed whether\n367 precession or obliquity has a greater influence on ~100-kyr glacial cycles (e.g., Hobart\n368 et al., 2023), we use ‘sensitivity’ to express phase, with a smaller phase indicating\n369 higher sensitivity. In the revised manuscript, we have added this explanation in the\n371\n372 1.17. Line 150-154: this is asserted, but isn't born out so far by proxy CO estimates.\n2\n374 for estimating CO2 so there is some circularity involved.\n2\n377 does not affect our major conclusions. This is supported by the stacked global mean\n378 surface temperature (Clark et al., 2024), which clearly display a cooling trend across\n379 the MPT, consistent with our interpretation. Moreover, our estimated precession phase\n380 displays a strong correlation with global mean surface temperature (r=0.88) (Fig. 5B,\n382 temperature in modulating glacial sensitivity to astronomical forcing, rather than\n383 establishing a direct connection between phase and CO .\n2\n384\n385 1.18. Line 183-186: Clarify whzt you mean specifically by \"evolution'. A number of\n386 points were made above and it's hard for the reader to follow what you are comparing\n387 the obliquity evolution to.\n389 phase demonstrates an increasing trend from 2700 ka to around 800 ka, then stabilizes\n390 (Fig. 3B). Conversely, the precession phase follows an opposite trend (Fig. 3C),\n391 indicating a progressive increase in ice sheet sensitivity to precession until around 800\n392 ka, after which it remained stable.’\n393\n394 1.19. Line 233-235: since the power ratio is relative, isn't it possibl that the actual\n395 importance of precession increased toward the present, but that the 100 kyr grew more?\n397 rather that it increased by a non-linear process as precession became more important.\n399\n401 they grow over time in the same direction.\n402\n404 compared to the ~21-kyr component. We demonstrated this situation using two simple\n405 linear functions, as shown in Fig. A3A. This situation will still result in a positive\n406 correlation between ~21-kyr and ~100-kyr power ratios, rather a negative one (Fig.\n408 the ~21-kyr power ratio, they will display a positive correlation. However, our statistical\n409 analysis indicates a negative correlation between the ~21-kyr and ~100-kyr power ratio.\n411 power ratio decreases over time.\n412\n413 1.20. Line 273-276: However, there is not indication that interglacial tmperatures or\n2\n2\n416 the MPT; the glacial-stage CO displays a decreasing trend across the MPT (Fig. 5A).\n2\n418 across the MPT (Fig. 3D) (Clark et al., 2024). Notably, during the MPT period, the\n419 cooling of ΔGMST displays an accelerating trend (Fig. 1D&E in Clark et al., 2024).\n421 periods after the MPT and the merging of several insolation peaks into ~100-kyr\n422 periodicities. Therefore, these observations explain the growing influence of precession\n423 on ice sheets across the MPT, as the ‘over’ developed ice sheets could extend southward\n424 gradually where precessional influence on daily insolation increases (Fig. 3C; Fig. 5F).\n426 temperature (r=0.88) (Fig. 5B, F; Fig. S17F). Thus, in the revised manuscript, we\n427 emphasized the role of background temperature in modulating glacial sensitivity to\n428 astronomical forcing, rather than establishing a direct connection between phase and\n2\n430\n431 1.21. Line 314-317: is the \"phase\" analysis here strong enough to support the claim of\n432 low latitude processes?\n434 new findings that the amplified ~100-kyr glacial power is most likely derived from\n435 precession rather than obliquity, we delve deeper into the reasons behind this conclusion\n436 (please refer to the first part of the revised Discussion) and discuss its implications for\n437 interpreting glacial cycles before and after the MPT (please refer to the second part of\n438 the revised Discussion).\n439\n442 (General comments) The cause of the Mid-Pleistocene Transition (MPT) is one of the\n443 most studied and enduring topics in the paleoclimate community. A better\n444 understanding of the MPT will lead to advances in untangling the internal feedback\n445 mechanisms that govern the climate pacing, improving the long-term predictability of\n446 future change. Drawing from analyses done on 57 paleoclimate records as well as the\n447 insolation, Zhang et al. propose that the gradual shift in atmospheric CO and ice sheet\n2\n448 extent caused a sudden shift in climate from a linear to a nonlinear system in response\n449 to precession around 1.5 million years ago. The meticulous investigation of a\n450 comprehensive compilation of datasets produced a novel and mostly convincing theory\n451 on the unfolding of the MPT. In short, the manuscript should be accepted for publication\n452 pending minor revisions to enhance its conclusions.\n454 records in testing their hypothesis, as exemplified in their Figure 2. To my knowledge,\n455 an investigation on this scale has not been attempted before, and I must commend the\n456 authors for the tremendous effort. The coherent power decomposition results for a\n457 diverse collection of records enhance the conclusion that precession and eccentricity\n458 powers show a complementary relationship pre- and post-MPT. This relationship\n459 between precession and eccentricity was crucial for the proposed theory that the system\n460 transitioned from a linear to nonlinear response to insolation during the MPT.\n461\n463 2.1 With Figure 5, the manuscript makes an argument that before the MPT, the\n464 precession power in LR04 is positively related to eccentricity power and inversely\n465 related to the precession phase. After the MPT, the two relationships flip, signifying a\n466 nonlinear glacial response to insolation. The authors acknowledge that LR04 is orbitally\n467 tuned. To avoid circular reasoning, they did the same precession phase calculation on\n469 records using a profile hidden Markov model, the age model of Prob_stack can shift\n470 freely according to the alignments of benthic δ18O records, thus remaining independent\n471 of orbital tuning.” While theoretically, Prob_stack ages are allowed to shift freely, it\n472 largely inherits the ages of LR04, its alignment target. This can be clearly seen in Figure\n473 4 of Ahn et al. (2017), where the timing of glacial terminations of the two stacks mostly\n474 tracks each other. Compared to LR04, Prob_stack includes more records and quantifies\n475 uncertainty, but it is not immune from the “leak” of orbital tuning in LR04 as its\n476 alignment target. I recommend the authors to try the same tests on untuned benthic\n477 δ18O stacks such as those from Lisiecki (2010) or Huybers (2007). These untuned,\n478 depth-derived stacks are constructed using geomagnetic reversals only, assuming a\n479 constant sedimentation rates between the geomagnetic age control points. Because\n480 geomagnetic reversals are widely dated using a variety of techniques, orbital\n481 assumptions are minimized.\n483 we eatablished the H07&NA stack, derived from the depth-derived H07 stack (Huybers,\n484 2007) with its last 640 kyr fitted to the North Atlantic benthic δ18O stack (NA stack)\n485 (Hobart et al., 2023) (Fig. S11; Table S3). The age model for the NA stack was\n486 established based on the speleothem age model (Cheng et al., 2016), making the\n488 phases based on the H07&NA stack align with the evolutionary trends of phases\n489 obtained using the LR04 stack and Prob_stack (Fig. 3; Fig. S14–S16), confirming our\n490 previous results. In the discussion on the reasons causing different glacial response to\n491 precession before and after the MPT, we only considered the astronomical phases of\n492 glacial terminations obtained from the H07&NA stack (Fig. 5).\n493\n494 2.2. Regarding the data presented in Figure 2, tests with different age models are done\n495 for two records in the Supplementary Materials. That is fine. However, it is not clear\n496 what age models are used for the rest of the records in the compilation. Are they using\n497 the age model provided by the respective study? Are they aligned to LR04 and, if so,\n498 do they potentially inherit the orbital tuning of LR04? The authors should discuss this,\n499 as it is likely important for the orbital power decomposition work done here.\n501 spectral power ratios, we directly used the original age models of the collected\n502 paleoclimate datasets. In the revised manuscript, we have added supplementary\n503 material introducing the chronologies of these datasets. Please refer to Table S1 and the\n504 section ‘Introduction to chronologies of the collected paleoclimate datasets’ in the\n505 supplementary materials. In summary, only several datasets have their own\n506 chronologies, while the majority have undergone orbital tuning.\n508 subjected to orbital tuning were separately assessed (Fig. S10D). Moreover, datasets\n509 primarily constrained by the LR04 age model (Lisiecki and Raymo, 2005) (Fig.\n511 maxima and minima (Figs. S12). As elaborated above, the age model of the H07&NA\n512 is independent of orbital tuning. In these occasions, the ~21-kyr power ratios negatively\n513 correlate with ~100-kyr power ratios, with no observed relationship between ~41-kyr\n514 and ~100-kyr power ratios (Fig. S10). Thus, we confirm that the orbital tuning has no\n515 influence on our results.\n516\n517 2.3. The visualization choices of Figure 2 can be improved. As it stands, it is difficult\n518 to compare the power ratio changes between the three panels. It does not help that the\n519 y-axis range is different among the panels. I encourage the authors to explore ways to\n520 improve Figure 2 so that important features stand out and are easier for readers to\n521 recognize. One possibility is to combine the three panels into one so that the power ratio\n522 boxes of different periods are shown side by side. See what I mean below. Then the\n523 correlation scatter plots can be enlarged, with colors of the scatter plots corresponding\n524 to the same color scheme.\n526 undergone major revisions following this suggestion.\n527\n529 2.4. line 20: the full terms should be shown in the first appearances of “kyr” and “ka.”\n531 that it occurred “around ~800-1200 ka” in the first sentence of Abstract. Consider\n532 changing “around ~800-1200 ka” to “around 1000 ka.”\n534 expression, ‘The transition of glacial cycles from ~40 kyr to ~100 kyr around ~800-\n535 1200 ka, without evident astronomical changes, underscores the interplay of Earth's\n536 internal and external factors in this transition,’ has been deleted. In the revised abstract,\n537 we begin with, ‘The Middle-Pleistocene Transition (MPT) represents a significant\n538 climate change over the Quaternary, but its cause remains not fully understood.’.\n539\n540 2.5. line 47: this sentence is too long and awkward. Consider breaking it apart.\n543 longer glacial cycles (Berger et al., 1999; Clark et al., 2006; Bintanja and van de Wal,\n544 2008; Willeit et al., 2019). These mechanisms include the long-term decline in\n545 temperature due to pCO decrease and the regolith removal beneath ice sheets (Clark\n2\n546 et al., 2006; Willeit et al., 2019), which enable the excessive development of ice sheets,\n547 thus gradually extending the glacial cycles into ~100-kyr periodicities (Legrain et al.,\n548 2023).\n549\n550 2.6. line 52: “phase-locking.”\n552\n553 2.7. line 68: remove “pioneeringly.”\n555 the revised manuscript.\n556\n557 2.8. line 263: remove “both.”\n559\n560 2.9. line 295: linear/nonlinear what? Maybe add “system” after “nonlinear” or another\n561 word as the authors see fit.\n563 ‘nonlinear’.\n564\n567 (General comments) Zhang et al. provide new and compelling data analyses which\n568 show that the 100-kyr periodicity in Late Pleistocene climate archives likely derived\n569 from a change in the climate sensitivity to precession forcing prior to the Mid-\n571 ice sheets is linked to the growth of ice sheets consecutively reaching the lower latitudes,\n572 where precession forcing becomes more and more influential on the glacial mass\n573 balance. This is a new and exciting discovery which will advance the field of both data\n574 and model paleo-science. I am convinced by the author’s data showing the change in\n575 sensitivity of climate to the precession forcing, and the subsequent hypothesis about\n576 gradual southward expansion of ice sheets against the backdrop of long-term CO\n2\n577 decline. (1) However, I currently disagree with the arguments in the last part of the\n578 manuscript which identify the onset of the MPT, and try to explain this mechanism to\n579 not just be a driver, but an instigator, of the MPT. In my opinion, shortening this section,\n580 and focusing on the manuscript’s main data analyses will elevate the manuscript and\n581 highlight the new findings that the analyses identify. (2) Additionally, the manuscript\n582 would benefit from including comparison of their findings and methods to other,\n583 already published, studies, such as work created by Huybers et al., or Berger et al. who\n584 either promote obliquity to control glacial-interglacial cycle development across the\n586 41kyr or 100kyr-driven climate cycles. It would be very interesting to read the opinion\n587 of the authors on those findings, and how their own results fit into the already existing\n588 narratives. (3) A final minor point, as I am not a specialist on the field of spectral\n589 analyses, but on data collection and MPT mechanisms, I sometimes found it difficult to\n590 follow the terminology in the current manuscript. I therefore included a list of\n591 suggestions below that hopefully help to clarify some questions and make the text more\n592 accessible to a wider audience. This includes minor changes, such as including\n593 definitions for certain spectral analysis vocabulary, and adding mathematical\n594 expressions of these definitions in the supplement, as well as adapting the figures for\n595 easier access.\n596\n598 3.1. However, I currently disagree with the arguments in the last part of the manuscript\n599 which identify the onset of the MPT, and try to explain this mechanism to not just be a\n600 driver, but an instigator, of the MPT. In my opinion, shortening this section, and\n601 focusing on the manuscript’s main data analyses will elevate the manuscript and\n602 highlight the new findings that the analyses identify.\n604 reorganized and merged the second and third parts of the Discussion into a single\n605 section (the second section of the Discussion in the revised manuscript). We no longer\n606 discuss the timing of the onset of the MPT; instead, we propose that the difference in\n607 periodicity before and after the MPT indicates a shift from a linear to a nonlinear glacial\n608 response to precession. This shift in response mode to precession is closely linked to\n609 the cooling-induced increase in deglacial threshold. Additionally, we discuss the links\n610 between cooling, the regolith removal beneath ice sheets, as well as the decline in CO\n2\n611 level.\n612\n613 3.2. Additionally, the manuscript would benefit from including comparison of their\n614 findings and methods to other, already published, studies, such as work created by\n616 interglacial cycle development across the Pleistocene (both 41 and 100kyr worlds), or\n617 find precession to influence either the 41kyr or 100kyr-driven climate cycles. It would\n618 be very interesting to read the opinion of the authors on those findings, and how their\n619 own results fit into the already existing narratives.\n621 findings to previously published work, such as Berger (1978), Hagelberg et al (1991),\n623 the revised manuscript, we wrote as follows:\n625 power is most likely derived from precession, with a less likely an obliquity\n626 contribution. Nevertheless, given that the main mass of ice sheets remains in polar\n627 regions where daily insolation variation is dominated by obliquity variance, obliquity\n628 should persistently exert a strong influence on ice sheets (Huybers& Wunsch, 2005;\n631 from the intrinsic relationship between precession and eccentricity, as all the periods of\n632 eccentricity are related to precessional ones via the heterodyne mechanism (Berger,\n633 1978). The nonlinear resonance between the frequencies of two main precession cycles,\n634 producing a beat at ~100-kyr frequency, can amplify the ~100-kyr glacial power\n635 (Berger, 1977; Hagelberg et al., 1991; Berger and Loutre, 2010). However, it is difficult\n636 to see the physical significance of such resonance (Imbrie et al., 1993).\n639 the last part of the discussion).\n640\n641 3.3. A final minor point, as I am not a specialist on the field of spectral analyses, but on\n642 data collection and MPT mechanisms, I sometimes found it difficult to follow the\n643 terminology in the current manuscript. I therefore included a list of suggestions below\n644 that hopefully help to clarify some questions and make the text more accessible to a\n645 wider audience. This includes minor changes, such as including definitions for certain\n646 spectral analysis vocabulary, and adding mathematical expressions of these definitions\n647 in the supplement, as well as adapting the figures for easier access.\n649 estimating spectral strength within a specific frequency band. We have included this\n650 clarification in the last paragraph of the Introduction of the revised manuscript.\n652 between the timing of the local maximal change rate of the benthic δ18O stack within a\n653 given glacial cycle and the timing at the nearest precession minimum or obliquity\n654 maximum. This definition has now been supplemented in the last paragraph of the\n656 and obliquity phases are now expressed with two mathematical equations (see the\n658\n660 3.4. line 1: As I currently cannot support the final chapter of the manuscript (see\n661 comments below), I suggest to change the title to “CO2-forced change in glacial\n662 response to precession drives the Middle-Pleistocene Transition”\n2\n664 2700 ka, the title of the manuscript has been revised to ‘Cooling-induced shift in glacial\n665 response to precession drives the Middle-Pleistocene Transition’. This revision aims to\n666 mitigate unnecessary debate and better align with our results because our estimated\n667 precession phase displays a generally consistent trend with the stacked global mean\n668 surface temperature (Fig. 3C&D; Fig. 5B&F) (Clark et al., 2024).\n669\n670 3.5. line 38: are there any references that confirm the frequency of obliquity, precession,\n671 and eccentricity for a wider audience?\n673\n674 3.6. line 52: Other references state that the 100 kyr cycle is possibly not phase locked\n675 to eccentricity, but expresses multitides of precession and obliquity with internal\n676 feedback phase locking (e.g. Berger & Loutre, 2010, Global & Planetary change;\n680 validation of the obliquity influence on late Pleistocene glacial cycles through a\n681 statistical test analysis. They proposed that the climate state of skipping one or two\n682 obliquity beats before deglaciation can yield quantized glacial-cycle durations of either\n683 80 or 120 kyr through either a stochastic or deterministic process. Berger and Loutre\n684 (2010) proposed that the ~100-kyr glacial cycles are related to a non-linear\n685 amplification of the astronomical forcing involving a beat effect from the two main\n686 precession periodicities.\n688 cycles (Berger et al., 1999; Huybers and Wunsch, 2005; Drysdale et al., 2009; Cheng\n689 et al., 2009; Huybers, 2011; Cheng et al., 2016; Bajo et al., 2020; Liautaud et al., 2020;\n691 kyr glacial cycles and eccentricity (Imbrie et al., 1993; Tziperman et al., 2006; Lisiecki,\n692 2010), and the non-linear amplification of the astronomical forcing related to\n693 eccentricity and precession (Berger and Loutre, 2010) suggest that astronomically-\n694 forced variations in insolation must be indispensable for the glacial cycles\n695 (Milankovitch, 1941; Hays et al., 1976)”.\n696\n697 3.7. line 63: It would be helpful to include definitions of the main concepts that will be\n698 used throughout the manuscript. These include spectral power, spectral phase, spectral\n699 power ratio, power shift. It would be supportive to have mathematical expressions of\n700 those definitions in the supplement as well.\n702 obliquity/precession phase in the revised manuscript. For more details, please refer to\n703 our reply to comment 3.3 above.\n704\n705 3.8. line 66: Briefly describe what the 57 paleoclimate archives are for a wider audience.\n707 archives were chosen.\n709 following sentence ‘These datasets, globally derived from ocean sediments, loess\n710 deposits, ice cores, and other archives (Fig. 1), are all continuous, with well-constrained\n711 age models (Table S1) and possess a high temporal resolution with at least a 2-kyr\n712 sampling rate.’\n713\n714 3.9. line 82: The authors state in line 82 that the relationship between the 21kyr power\n715 ratio and the 100kyr power ratio holds during and before the MPT. However, the cross-\n716 plot in Figure 2 panel A has the highest r of 0.7. That surely shows that the relationship\n717 holds throughout? Line 83: Can the authors comment on why their results show no link\n718 with obliquity, against what published studies have found (e.g. Feng & Bailer-Jones,\n719 2015, QSR).\n721 robustness of the relationship between spectral power ratios of ~21-kyr, ~41-kyr, and\n722 ~100-kyr components (Fig. S10). These scenarios considered the potential influence of\n723 significance level of spectral power peaks, dataset length, age models used, and orbital\n724 tuning on the PDA results and their relationship. The results indicate that across all\n725 scenarios, the ~21-kyr power ratios negatively correlate with ~100-kyr power ratios,\n726 with no observed relationship between ~41-kyr and ~100-kyr power ratios (Fig. S10).\n728 kyr power ratios is robust.\n730 90% confidence level, we found that the resulting correlation between the ~21-kyr and\n731 ~100-kyr power ratios for datasets spanning the past 0–800 ka (r=0.75) remained\n732 consistent with those for datasets spanning the past 800–1200 ka (r=0.75) (Fig. 2D, F).\n734 their results show no link with obliquity, against what published studies have found’,\n735 we would like to clarify that we did not mean that obliquity has no influence on ice\n736 sheets. Our results show that the obliquity phase of glacial termination generally\n737 displays an opposite evolutionary trend compared to the precession phase (Fig. 3B, C),\n738 signifying that the obliquity’s role in modulating glacial cycles becomes less important\n739 over time. Moreover, the reference of Feng & Bailer-Jones (2015) has been cited in the\n740 revised manuscript.\n741\n742 3.10. line 94: can the authors confirm why the analysis was conducted on LR04, and\n743 not the 57 climate records which are plotted in Fig. 1?\n745 relative to the glacial terminations, which could be best determined using ice-volume-\n746 related proxies. Despite the benthic δ18O being influenced by local deep-sea\n747 temperature fluctuations, it remains a widely-used proxy for global ice volume\n748 (Huybers, 2006; Liautaud et al., 2020; Hobart et al., 2023). Moreover, Tzedakis et al\n749 (2017) provided a detrended version of the LR04 stack, removing the temperature\n750 contributions to the LR04 stack. Thus, in our manuscript, the definition of a glacial\n751 cycle uses the glacial roster derived from the detrended LR04 benthic δ18O stack\n752 (Tzedakis et al., 2017). Consequently, the estimated precession and obliquity phases are\n753 minimally affected by temperature impacts on benthic δ18O.\n755 from an ice sheet model driven by 65°N summer insolation (Lisiecki and Raymo, 2005),\n756 potentially involving circular reasoning when discussing orbital influences. Thus, in the\n757 revised manuscript, we employed the H07&NA stack as the reference curve to re\n758 calculated the precession and obliquity phases. The definition of the glacial-interglacial\n759 cycle still relies on the roster provided by Tzedakis et al (2017). This is because the\n761 640 kyr aligned to the North Atlantic benthic δ18O stack (NA stack) (Hobart et al., 2023)\n762 (Fig. S11; Table S2). The age model for the NA stack was established based on the\n763 speleothem chronology (Cheng et al., 2016), making the H07&NA stack independent\n764 of orbital tuning.\n765\n766 3.11. Figure 3: a) Can the authors include the mathematical expression in the\n767 supplement or reference the software, of how the derivative was calculated?\n768 b) Also I suggest to change the y-axes labels to something more intuitive for wider\n769 audience accessibility, especially for panels C, D, E, F, K. I suggest the authors to\n770 consider a re-organisation of the panels. The small cross-plots (H, I, J) on the right\n771 invoke connection to panels E, F, G, which I believe is not the intention of the\n772 authors?\n773 c) Considering B: can the authors comment on why this sea level record was chosen?\n775 et al., 2014, Sosdian et al. 2009, and others).\n776 d) Have the authors considered looking at other CO records, such as Yamamoto et al.,\n2\n777 2022? If so, why was the record in G chosen?\n778 e) drawing the vertical dotted lines through the entire graph might be more intuitive for\n779 the reader.\n782 its calculation in the Method section. It is presented as follows: “d(δ18O)/dt denotes the\n783 changing rate of benthic δ18O, calculated by dividing the difference between two\n784 adjacent δ18O values by the corresponding time difference.”\n786 sea level results (de Boer et al., 2014; Willeit et al., 2019; Berends et al., 2021-CP) were\n787 all incorporated in the revised manuscript (Fig. 3C).\n2\n789 reconstructed CO records, including δ11B-based reconstruction (Chalk et al., 2017;\n2\n791 2018; Bartoli et al., 2011), reconstructions derived from Allan Hills Blue Ice (Yan et al.,\n792 2019) and EPICA ice core (0-800ka; Bereiter et al., 2015), and a longer and continuous\n2\n794 (Yamamoto et al., 2022).\n795\n796 3.12. line 186: I would be interested to see the author’s comment on the following\n797 publication here: Huybers et al., 2006, Science.\n799 has been cited in the revised manuscript. Our study attributes the increasing influence\n800 of precession on ice sheets to the cooling-induced glacial dynamics. The enhanced\n801 stability of ice sheets, due to the cooling and the elevated deglaciation threshold (Fig.\n802 3D & Fig. 5E), promotes their long-term expansion and eventual merging into a\n803 continental-scale size (Bintanja and van de Wal, 2008). This process facilitates the\n804 gradual southward extension of the ice sheets (Barker et al., 2022) (Fig. 3E&F).\n806 responsive to middle-low latitude insolation, rendering them more sensitive to\n807 precession and less sensitive to obliquity over time (Fig. 3B, C), because precession\n808 variance dominates daily insolation at middle-low latitudes.\n810 sheets were believed to gradually transition to partly marine-based conditions after the\n812 ice volume on the precession scale and prevented the cancellation of precession effects\n813 on ice volume between hemispheres, thereby enhancing the precession expression in\n814 globally integrated proxies after the MPT (Lisiecki and Raymo, 2005; Raymo et al.,\n815 2006; Feng and Bailer-Jones, 2015; Liautaud et al., 2020).\n816\n817 3.13. line 200: It might be worth mentioning that ice sheets do not only drive internal\n818 variability, but can be affected by other internal feedback mechanisms themselves.\n820 have added the following narrative: Meanwhile, other climate subsystems may also\n821 influence ice sheets. Similarly, under the interplay of various subsystems, the spectral\n822 power ratio relationship within these subsystems will eventually propagate throughout\n823 the entire climate system through a series of positive and negative feedbacks (Clark et\n824 al., 1999; Denton et al., 2010).\n825\n826 3.14. Figure 4. In general, I did not find this figure very intuitive and it took me some\n827 time to understand what the authors are trying to express. I kindly listed my issues\n828 and questions below.\n829 a) Panel A: a wider audience will not be familiar with the term 1E1T-P. I assume this\n830 refers to the sum of all power spectra defined by the authors as Precession,\n832 suggest a small comment in the supplement, why this definition is chosen, and\n833 where it comes from. Maybe it is a common reference, but clarifying that would be\n834 useful.\n835 b) Since the term power ratio is not necessarily well known to a wider audience, it\n836 would be useful to either have a clear definition in the text, or re-word the definition\n837 very brief in the figure or y-axis figure caption. Sometimes annotations may also\n838 help.\n839 c) Is there a plot that shows the behavior of obliquity in comparison to the plot in panel\n841 correct?\n842 d) Is the threshold in panel D just an example of a threshold, given that panel H and I\n843 show multiple thresholds on the x-axis?\n844 e) I wonder whether it would be helpful to include a plot comparable to panel G, where\n845 each threshold’s power spectrum is plotted in the same plot, maybe using grey\n846 scales. Then the reader could see that changing the threshold will cause the\n847 appearance of the 95, 125, and 405 spectra.\n848 f) The jump from panel D to panel F is counter-intuitive, since the y-axis are flipped on\n849 panel F. This needs to be clearly stated, or highlighted on the y-axis.\n850 g) Using more descriptive titles on panels H, I, J would be useful.\n851 h) The figure caption must include a description of what the whisker plots represent. I\n852 assume they represent the amplitude of the spectrum at each threshold? What are\n853 the error bars on the whisker plots?\n854 i) In general, figure 4 is central to the hypothesis of the authors and shows the link\n855 between changes in climate sensitivity. Next to figure 2, it seems to be the most\n856 important figure. As such, I recommend the authors to consider re-organizing the\n857 panels, or re-writing the figure caption, to help the general reader follow this\n858 information-heavy figure more easily.\n860 (1) In terms of the 1E1T-P, we have supplemented a definition in the figure caption\n861 as follows: ‘1E1T-P curve over the past 2700 ka, generated by normalizing and\n862 averaging variations in eccentricity (E), tilt (T), and reversed precession (P) based on\n864 (2) Regarding the definition of the power ratio, in the revised manuscript, we have\n865 provided additional details in the Method section. ‘The spectral power ratios reflect the\n866 variance of precession, obliquity, and ~100-kyr cycles within datasets.’ This sentence\n867 has also been added at the end of the figure caption.\n868 (3) Is there a plot that shows the behavior of obliquity in comparison to the plot in\n869 panel C? For this question. Also, for clarification, I assume that panel C is a cross-plot\n870 of panel B, is that correct?\n872 precession, each with equal variance. Thus, the obliquity power ratio (variance) should\n873 show a complimentary correlation with the sum of the precession and ~100-kyr power\n874 ratios. In a linear climate system, the precession power ratio (variance) shows a similar\n875 evolutionary trend to the ~100-kyr power ratio (variance) (Fig. 4B-C). Thus, the\n876 obliquity power ratio (variance) is expected to show an inverse correlation with both\n877 the ~100-kyr and precession power ratios. The Fig.4 A-C were mainly used to\n878 demonstrate that in a linear climate system, the precession power ratio (variance) shows\n879 a similar evolutionary trend with the ~100-kyr power ratio; while in a nonlinear system,\n880 they should show a complementary correlation (Fig. 4H-J). Thus, to avoid confusion,\n881 we chose not to plot the correlation between the obliquity power ratio and ~100-kyr\n882 power ratio of the 1E1T-P in Fig. 4, because they are expected to display a negative\n883 correlation.\n885 shown in panel B. For clarity, the caption has been revised as follows: ‘C, correlation\n886 plot between the calculated precession and ~100-kyr spectral power ratios of the 1E1T-\n888 (4) Is the threshold in panel D just an example of a threshold, given that panel H\n889 and I show multiple thresholds on the x-axis?\n891 the rectified precession curve when the threshold is set to -0.025, alongside the\n892 corresponding power spectra.\n893 (5) I wonder whether it would be helpful to include a plot comparable to panel G,\n894 where each threshold’s power spectrum is plotted in the same plot, maybe using grey\n895 scales. Then the reader could see that changing the threshold will cause the appearance\n896 of the 95, 125, and 405 spectra.\n898 different thresholds within a single coordinate scheme would be impractical. This is\n899 because after rectifying the precession parameter, the unite/scale for the rectified\n900 precession parameter differs (please see the left Y axis in the Fig. 4D & Fig. 4F). Thus,\n901 the corresponding absolute power levels are not comparable. Hence, we used power\n902 ratio to quantitatively represent the variance of precession, obliquity, and ~100-kyr\n903 component within the rectified precession parameters. Nevertheless, we also included\n904 the MTM spectral plots for each threshold scenario. As illustrated in Fig. A4, the\n905 variance at ~1/100 kyr-1 increases with the increase in threshold effect (smaller Th\n906 values).\n908 components with the change in threshold, compared to Fig. 4E, where no threshold was\n909 applied to the precession parameter.\n910\n913 hemisphere summer insolation; thus, a more negative threshold value (Th values) indicates a\n914 stronger threshold effect.\n915\n916 (6) The jump from panel D to panel F is counter-intuitive, since the y-axis are\n917 flipped on panel F. This needs to be clearly stated, or highlighted on the y-axis.\n919 (7) Using more descriptive titles on panels H, I, J would be useful.\n921 used to rectify the precession parameter, plotted against the precession power ratios\n922 within the corresponding rectified precession parameters. The power ratios were\n923 displayed in a box whisker plot, where the horizontal line in the box whisker plot\n924 represents the median values of the calculated precession power ratios. I, same as H,\n925 but the y axis represents the ~100-kyr power ratio within the rectified precession\n926 parameter. J, correlation plot between the median of the power ratios at the precession\n927 and ~100-kyr bands in the rectified precession parameter generated in H and I.\n928 (8) The figure caption must include a description of what the whisker plots\n929 represent. I assume they represent the amplitude of the spectrum at each threshold?\n932 different threshold setting. This is similar to how the entire spectral power ratios of the\n933 collected datasets at specific frequency bands were displayed, as shown in Fig 2 A-C.\n935 the maximum and minimum values of the calculated power ratio. In the figure caption,\n936 we included this sentence ‘The spectral power ratios were displayed in a box whisker\n937 plot, where the horizontal line represents the median values of the calculated precession\n938 power ratios’.\n939\n940 3.15. line 233: Can the authors comment on why they do not find any influence of\n941 obliquity on climate in comparison to other studies? (e.g. Huybers & Wunsch, 2015;\n943 author’s results compare to other studies finding precession to be a dominant driver\n944 (e.g. Liautaud et al., 2020, EPSL)?\n946 ice sheets. We show that the obliquity phase of glacial termination displays a generally\n947 opposite evolutionary trend compared to the precession phase (Fig. 3B, C), signifying\n948 that the obliquity’s role in modulating glacial cycles becomes less important over time.\n950 power ratios with the results acquired by Liautaud et al. (2020) (Fig. S1). These spectral\n951 power ratios were calculated using the PDA method based on three different benthic\n952 δ18O records, including Prob_stack (Ahn et al., 2017), LR04 stack (Lisiecki and Raymo,\n953 2005), and the H07&NA stack (Huybers, 2007; Hobart et al., 2023) (Fig. S11). The\n954 resulting ~21-kyr and ~41-kyr spectral power ratios show similar long-term evolution\n955 trends with the envelops of precession and obliquity contributions to the derivative of\n957 findings of Liautaud et al. (2020). These envelopes represent the magnitude of the\n959 18O(Liautaud et al., 2020). The ENOF stands for the Empirical Nonlinear Orbital Fitting\n960 (Liautaud et al., 2020). As proposed by Liautaud et al. (2020), the ENOF method is\n961 complementary to spectral analysis. Thus, by adding a sliding time window, spectral\n962 analysis becomes equivalent to the ENOF method. The PDA methd is implemented\n963 with a sliding time window, which is why our PDA-based results are similar to ENOF-\n964 based results (Fig. S1).\n965\n966 3.16. line 234: I was unable to follow the statement made in this sentence. Can the\n967 authors clarify why they claim that “the nonlinear amplification of 100kyr power is an\n968 intrinsic characteristic of the entire climate system, not just of ice volume?\n970 ratios across 57 paleoclimate archives suggests that the amplification of ~100-kyr\n971 power occurs at the expense of precession power (Fig. 2A, B, D). This power shift\n972 persists during and before the MPT when glacial cycles predominantly exhibited ~41-\n973 kyr cyclicities (Fig. 2A, B, F, H). This indicates that even in the ‘41-kyr glacial world’,\n974 other climate subsystems could also show strong ~100-kyr cyclicities and display a\n975 reverse correlation between the ~100-kyr and ~21-kyr power ratios. Thus, the nonlinear\n976 amplification of ~100-kyr power, caused by energy transformation from precession, is\n977 not exclusive to the glacial subsystem but is an intrinsic feature of the entire climate\n978 system.\n979\n980 3.17. line 237: can the authors elaborate why multiple thresholds are needed to induce\n981 sensitivity change, and not just one threshold?\n983 of spectral power transferred from the precession band to the ~100-kyr modulation band\n984 become a constant. However, the PDA results of the benthic δ18O indicate that the ~100-\n985 kyr power ratio displays a gradual, rather than abrupt, evolution trend across the MPT\n986 (Fig. 5E; Fig. S1G), hinting at a gradual increase in the deglaciation threshold (Tzedakis\n987 et al., 2017; Legrain et al., 2023). It is estimated that the effective energy required to\n988 initiate deglaciation increased from an average of 6.025 GJ.m-2 to 6.142 GJ.m-2 before\n989 and after the MPT (Tzedakis et al., 2017) (Fig. 5E). Similarly, a conceptual model\n990 simulation suggests that reproducing the frequency shift over the MPT requires an\n991 increase in the deglaciation threshold over time (Leloup and Paillard, 2022).\n993 glacial cycles.\n995 parameter, and assuming a series of threshold insolation below which the climate\n996 system does not respond, result in more precession spectral power being transferred to\n997 the ~100-kyr band as the threshold increases (Fig. 4H–J), consistent with the PDA\n998 results from the collected climate archives (Fig. 2D, F, H).\n999\n1000 3.18. line 240: I believe this analysis was conducted on data from LR04 correct? As\n1001 such, why do the authors use LR04 as a proxy for global climate instead of ice volume,\n1002 even though they have a perfect collection of 57 cores that are better at representing\n1003 global climate than LR04?\n1006 datasets collected from 33 sites. The precession and obliquity phases of glacial\n1007 terminations are estimated using the benthic δ18O stack because it closely approximates\n1008 a proxy for ice volume, whereas the other collected datasets cannot be regarded as\n1009 proxis of ice volume.\n1011 remains a valid proxy for discussing global ice volume changes. This approach has been\n1012 extensively validate and continues to be used in numerous studies (e.g. Huybers and\n1015 multiple δ18O records from different oceans (e.g. the stacked δ18O) are combined like\n1016 in our study, the local deep-sea temperature effects are mitigated, providing a clearer\n1017 reflection of global ice volume changes.\n1019 derived from the detrended LR04 benthic δ18O stack and the two-threshold scheme\n1020 (Tzedakis et al., 2017). The purpose of detrending is to remove the temperature\n1021 contribution to benthic δ18O (Tzedakis et al., 2017). Thus, the estimated precession and\n1022 obliquity phases are minimally affected by temperature influences on benthic δ18O.\n1023\n1024 3. 19. line 244: I cannot quite follow the argument in the following lines. Are the authors\n1025 suggesting that the 100kyr cycle in LR04 occurs due to the 400 kyr power in eccentricity,\n1026 and therefore disagree with a collection of publications claiming that the 100kyr is a\n1027 “quasi 100kyr” cycle modulated by precession and obliquity due to internal feedbacks?\n1030 that the climate threshold effect not only enhances ~100-kyr power but also amplifies\n1031 ~400-kyr power as thresholds increase (Fig. 4F, G). Pleistocene ice-age records, such\n1032 as benthic δ18O records, have been noted for their lack of power at the ~1/400 kyr-1\n1033 frequency, a phenomenon known as the ‘400-kyr problem’ (Imbrie et al., 1993; Rial et\n1034 al., 2013). This absence could be explained by the frequency modulation of ~100-kyr\n1035 glacial cycles by the ~413-kyr eccentricity component (Rial, 1999; Rial et al., 2013).\n1037 the modulating signal as a spectral peak because its energy shifts to the ~100-kyr carrier\n1038 (Rial, 1999; Rial et al., 2013). One plausible interpretation is that enhancements in the\n1039 ~405-kyr components, caused by the precession threshold effect, are offset by internal\n1040 nonlinear synchronization or phase locking at the ~405-kyr band. This indicates that\n1041 the internal response of the climate system to external forcing is far more complicated\n1042 than a single, simple schematic model; multiple nonlinear mechanisms may be at work\n1043 simultaneously.\n1045 self-sustained oscillations, which are expected to synchronize with eccentricity-\n1046 modulated precession forcing (Rial et al., 2013; Mitsui et al., 2023).\n1047\n1048 3.20. line 262: At the moment, I cannot agree with this statement, since the authors\n1049 receive different timings depending on whether their analysis uses LR04, or the\n1050 probabilistic stack. Previously in the manuscript the authors argue that the use of the\n1051 probabilistic stack is better, due to a more robust age model. I suggest that the authors\n1052 stick with this logical approach, and present data from the probabilistic stack. The fact\n1053 that all the other analyses show no sensitivity to different age models allows me to\n1054 believe their results. In this case, I believe that the analysis shows sensitivity to age\n1055 modelling, which in itself is an interesting result, but should not be over interpreted in\n1056 terms of climate dynamics.\n1059 the revised manuscript by alternatively employing the H07&NA stack to estimate the\n1060 precession phase. As we elaborated above, the age model of this stack is independent\n1061 of orbital tuning. Based on these new results, we found a change in the correlation\n1062 patterns between the ~21-kyr and ~100-kyr power ratios of the H07&NA stack (Fig.\n1063 5E, G). Preceding ~800 ka, there is a positive correlation (r=0.63), suggesting a linear\n1064 mode; while after ~800 ka, a negative correlation (r=-0.7) between them indicates a\n1065 nonlinear threshold mode (Fig. S17H). Moreover, the correlation between the\n1066 precession phase and the ~21-kyr power ratio in the H07&NA stack changes from\n1067 negative to positive around ~800 ka (Fig. 5F, G; Fig. S17I). These contrasting\n1068 correlations before and after the MPT imply that the glacial response to precession has\n1069 shifted from a linear to a nonlinear mode.\n1071 occurrence of nonlinearity. Instead, we focus more on the reasons that cause the\n1072 transition in response mode and its implications for interpretating the glacial cycles.\n1074 the MPT) for further discussion.\n1075\n1076 3.21. line 259 and following section: In the following lines the authors present a\n1077 hypothesis for how precession may drive the MPT. At this point, I am unable to agree\n1078 with the hypothesis, due to (i) a lack of robust data which shows that the shift in\n1079 sensitivity to the precession cycle, and re-location of power to 100kyr cycle is due to\n1080 long-term reduction in pCO (see line 274). (ii) a mis-representation of reference 7 (see\n2\n1081 line 271). Reference 7 reviews processes that influence the MPT, but does not state\n1082 significant hypotheses on how the MPT was instigated. In contrast, ref 7 proposes a\n1083 collection of research questions which would be beneficial for the field to solve. (iii) a\n1084 misrepresentation of refs 16, 45 (see line 278) which consider the late Pleistocene, not\n1085 the early Pleistocene prior to the MPT. As such, their evidence needs to be used after\n1086 large-scale ice sheets are already occurring (see line 281). (iv) the statement that\n1087 “bedrock-mass-balance feedback and increased basal-sliding influences ice sheet\n1088 dynamics” would benefit from more than ref 9, as this process has been suggested to\n1089 influence ice sheets on multiple timescales, ie during deglaciations, or as an instigator\n1090 for the MPT (see Regolith hypothesis, Clark et al., 2006). (v) a misrepresentation of ref\n1091 11 which only tested the change in long-term CO coupled to regolith removal, vs. only\n2\n1092 regolith removal. It did not test whether abrupt or gradual changes in CO cause the\n2\n1095 the subtitle ‘Implications for the MPT’.\n1096 (1) Addressing the first concern regarding attributing the shift in glacial response\n1097 to precession to CO (i). In addition to comparing our estimated precession and\n2\n1098 obliquity phase to the simulated CO (Willeit et al., 2019), in the revised manuscript,\n2\n1099 we compare the precession phase to the stacked global mean surface temperature over\n1100 the last 2700 ka (Clark et al., 2024). This comparison reveals a covariance between\n1101 precession phase and temperature (Figs. 3&5), suggesting that the variation in the\n1102 glacial response to precession may be related to changes in background temperature.\n2\n1104 reconstructions (Foster, 2008; Hönisch et al., 2009; Seki et al., 2010; Bartoli et al., 2011;\n2\n1108 proxy spanning the last 1500 ka has been also collected (Yamamoto et al., 2022). These\n1109 proxy-based CO reconstructions do indicate a decline in glacial-stage CO across the\n2 2\n1111 glacial-stage pCO .\n2\n2\n1113 over the last 2700 ka, our revised manuscript we only established a direct link between\n1114 glacial response to precession and the cooling temperatures in the revised manuscript.\n1115 (2) In terms of the subsequent concerns about the mistaken citation of references\n1116 (ii- iv), for the citation of ref 7 (Berends et al., 2021-RG), we are simply referring to\n1117 their schematic two-threshold framework (in their Fig. 4) to demonstrate how the\n1118 insolation peaks were skipped. The difference is that they assumed that the skipped\n1119 insolation peaks are caused by obliquity rather precession. In our revised manuscript,\n1120 we no longer discuss the potential skipping of precession-induced insolation peaks.\n1123 (3) In terms of the misrepresentation of ref 11 (Legrain et al., 2023) (v). This article\n1124 argues that certain internal factors are necessary for explaining the MPT and suggests\n1125 that factors contributing to a gradual process and linear increase in deglacial threshold,\n1126 rather than an abrupt one, could better reproduce the MPT. In the revised manuscript,\n1127 we paid more attention to the role of background temperature in modulating the\n1128 variation in glacial response to precession and obliquity. This is because the global\n1129 mean surface temperature (Clark et al., 2024) displays a gradual declining trend over\n1130 the last 2700 ka. This gradual cooling could lead to a linear increase in deglaciation\n1131 threshold. Thus, in this discussion section, we cite this article (Legrain et al., 2023).\n1133 cited.\n1134 3.22. line 294: This statement is currently unsupported due to the above comments.\n1136 sentence has been deleted.\n1137\n1138 3.23. line 297: Given that the probabilistic stack analyses of the authors show an onset\n1139 at 1.2Ma, the literature seems to support the findings of the authors with respect to the\n1140 probabilistic stack. Although I would suggest not putting too much influence on this\n1141 analysis, as it seems highly variable to changes in age models, as the authors admit\n1142 themselves in line 300.\n1144 stack to recalculate the precession phase and power ratios at the ~21-kyr and ~100-kyr\n1145 bands, as it is based on an age model independent of orbital tuning, avoiding circular\n1146 reasoning. We focused on the implications of the shift in glacial response to precession\n1147 for interpretating the MPT and reduced the discussion on the exact timing of the shift.\n1149 more details.\n1150\n1151 3.24. line 301: This sentence might need reworking.\n1153\n1154 3.25. line 303: “High-amplitude” glacial cycles are generally believed to occur after\n1155 800kyr as they describe the 100kyr world glacial cycles, which occur after the MPT.\n1158 alternative idea regarding the timing of the occurrence of ‘high-amplitude’ glacial\n1159 cycles. In the revised manuscript, we have reorganized this section and deleted this\n1160 narrative to avoid unnecessary debate.\n1162 glacial response to precession and its implications for interpreting glacial cycles, rather\n1163 than delving into the detailed dynamic processes causing the skipping of insolation\n1164 peaks.\n1165\n1166 3.26. line 305: ref 47 needs to be discussed in line with more recent publications on this\n1167 topic (some already listed throughout the review).\n1169 emphasizing the occurrence of nonlinear responses to precession and their implications\n1170 for the interpretation of glacial cycles, rather than delving into the specific timing of\n1171 these occurrences. As a result, the discussion regarding the timing of the MPT\n1172 occurrence has been deleted.\n1173\n1174 3.27. line 308: At this point, I cannot agree with the conclusion that the MPT is a gradual\n1175 process. The authors do not provide robust evidence for that, and the list of arguments\n1176 is lacking. In fact, I am unsure whether the authors are able to supply this information\n1177 with the data they have at this stage. I therefore suggest to shorten the last section,\n1178 which should not take away from the novelty of the work presented before.\n1180 the MPT is a gradual or abrupt process, and the last section has been adjusted\n1181 accordingly.\n1182\n1183 3.28. Figure 5: a) The timing of 1452kyr seems dependent on the age model used which\n1184 makes this graph hard to accept as robust\n1185 b) Improvement of the y-axis label would help with the understanding of the graph.\n1186 c) To help a wider audience, it might be useful to guide the reader through the cross\n1187 plots and explain what the change in these patterns means for global climate. The\n1188 correlations are compelling and I believe the analysis should focus more on the change\n1189 in the correlation direction, which stays robust throughout different age models, rather\n1190 than focusing on the timing of the MPT.\n1191 d) In the cross plot in panel D, E, the yellow line is annotated with p=0.86. Should this\n1192 be r=0.86 or does it refer to a p-value? If so, maybe also report r.\n1195 the correlation between the astronomical forcings and factors internal to ice sheets .\n1197 because the calculated precession power ratios based on the LR04 stack, Prob_stack,\n1199 (Elderfield et al., 2012) (Fig. A2C) all show a similar evolution trend during the period\n1200 between 0–1200 ka, indicating a comparably low spectral power ratio. The estimated\n1201 ~100-kyr spectral power ratios based on these datasets also display similar evolution\n1202 trends, which are nearly the revere of the ~21-kyr power ratio (Fig. S1D, G; Fig. A2\n1204 longer emphasize the exact timing of the shift in glacial response to precession. Instead,\n1205 we focus more on the implication of the shift for the interpretation of glacial cycles.\n1207 grouped in Fig.S17, which is included in the supplementary materials.\n1209\n1210 3.29. line 439: Can this be expressed as a mathematical equation?\n1212 added a supplementary figure (Fig. S1). Moreover, we have rewritten the first part of\n1213 the Method. The revised description is as follows:\n1214 (1) Initially, the evolutionary power spectra of the datasets were calculated using\n1215 the 2π multi-taper method (MTM) (Thomson, 1982) with a sliding time window\n1216 (Meyers et al., 2012) (Fig. S1B). (2) Then, the time-dependent amplitude modulations\n1217 (power) at three frequency bands (~1/21 kyr-1, ~1/41 kyr-1, and ~1/100 kyr-1) were\n1218 quantitatively calculated by integrating the evolutionary power at the corresponding\n1219 frequency bands (Meyers et al., 2012; Li et al., 2019). (3) Thereafter, to facilitate\n1220 comparisons among different paleoclimate datasets, the integrated evolutionary spectral\n1221 power at each specific band was normalized by the total evolutionary power (integrated\n1222 from 0 to the Nyquist frequency) to obtain the evolution spectral power ratio. This\n1223 process is necessary because different paleoclimate proxies have varying units.\n1224\n1225 3.30. line 442: It would be useful to include a detailed and clear definition for the\n1226 different terms used throughout the manuscript, such as spectral power, spectral\n1227 strength, spectral power ratio, spectral phase etc.\n1229 definitions of the spectral power ratio, obliquity and precession phases. Additionally,\n1230 the definitions of the precession and obliquity phases are now expressed with two\n1231 mathematical equations (see the Method section).\n1232\n1233 3.31. line 444: I wonder how different the age models are between these different\n1234 approaches. Given the recent improvement in age modelling, the age models tested here\n1235 might not vary greatly which could lead to internal age model errors per age model\n1236 being larger than between- age model errors. Could the authors include a plot with one\n1237 record on all age models so the variability can be compared? Alternatively, the authors\n1238 could consider making a plot like S1 or S2, but using the maximum error values for one\n1239 analysis and comparing to the minimum error values from the same age model from a\n1240 second analysis. This should also show the influence of internal age model errors on\n1241 the analyses.\n1243 collected datasets. Most of them are established through orbital tuning process,\n1244 primarily constrained by the LR04 age model.\n1248 similar evolutionary trends (Fig. S1D, F, G;). These trends align with those acquired\n1249 through the Empirical Nonlinear Orbital Fitting (ENOF) method (Fig. S1C, E)\n1250 (Liautaud et al., 2020). Thus, the PDA results are insensitive to age model used.\n1252 age model usage, and orbital tuning on the PDA results and the correlation of spectral\n1253 power ratios. The results indicate that, across all scenarios, the ~21-kyr power ratios\n1254 negatively correlate with ~100-kyr power ratios, with no observed relationship between\n1255 ~41-kyr and ~100-kyr power ratios (Fig. S10). Thus, we can state that the negative\n1256 correlation between ~21-kyr and ~100-kyr spectral power ratios is robust.\n1257\n1258 3.32. line 458: What do the authors mean by “robust” age models?\n1260 collected datasets, most of which have been subject to orbital tuning. Please refer to\n1261 part one and Table S1 in the supplementary materials for details. The term ‘robust’\n1262 implies that their age models are correct and the temporal resolution of these datasets\n1263 is sufficient for detecting precession signals. Nevertheless, to avoid any potential\n1264 misunderstanding, we have removed the word ‘robust’ from the manuscript.\n1265\n1266 3.33. line 476: Can this be expressed as a mathematical equation?\n1268 between the timing of glacial terminations and the nearest precession minima/obliquity\n1269 maxima. They can be expressed as follows:\n1270 (1) 𝑃ℎ𝑎𝑠𝑒 = |𝑇 −𝑇 |\n𝑜𝑏𝑙𝑖𝑞𝑢𝑖𝑡𝑦 𝑚𝑎𝑥𝑖𝑚𝑢𝑚[𝑑𝛿18𝑂/𝑑𝑡] 𝑚𝑎𝑥𝑖𝑚𝑢𝑚[𝑜𝑏𝑙𝑖𝑞𝑢𝑖𝑡𝑦]\n1271 (2) 𝑃ℎ𝑎𝑠𝑒 = |𝑇 −𝑇 |\n𝑝𝑟𝑒𝑐𝑒𝑠𝑠𝑖𝑜𝑛 𝑚𝑎𝑥𝑖𝑚𝑢𝑚[𝑑𝛿18𝑂/𝑑𝑡] 𝑚𝑖𝑛𝑖𝑚𝑢𝑚[𝑝𝑟𝑒𝑐𝑒𝑠𝑠𝑖𝑜𝑛]\n1273 maximum/minimum of the nearest obliquity/precession parameter to the glacial\n1274 termination; d(δ18O)/dt denotes the changing rate of benthic δ18O, which is calculated\n1275 by dividing the difference between two adjacent δ18O values by the corresponding time\n1276 difference.\n1277\n1278 3.34. line 480 and following: Some of these definitions should be included in the main\n1279 text to help a wider audience understand the analyses more easily.\n1281 precession and obliquity phase expressed with two mathematical equations. Please refer\n1282 to our reply to comment 3.33 and the Method part in the manuscript for details.\n1283\n1285 1. Ahn, S., Khider, D., Lisiecki, L.E. & Lawrence, C.E. A probabilistic Pliocene–Pleistocene stack of benthic\n1286 δ18O using a profile hidden Markov model. Dyn. statist. clim. Syst 2, 1–16 (2017).\n1287 2. Bajo, P., Drysdale, R.N., Woodhead, J.D., Hellstrom, J.C., Hodell, D., Ferretti, P., et al. Persistent influence of\n1288 obliquity on ice age terminations since the Middle Pleistocene transition. Science 367, 1235–1239 (2020).\n1289 3. Barker, S., Starr, A., van der Lubbe, J. Doughty, A., Knorr, G. & Conn, S. et al. Persistent influence of\n1290 precession on northern ice sheet variability since the early Pleistocene. Science 376, 961–967 (2022).\n1291 4. Bartoli, G., Hönisch, B., & Zeebe, R. E. Atmospheric CO2 decline during the Pliocene intensification of\n1292 northern hemisphere glaciations. Paleoceanography 26, PA3206 (2011).\n1293 5. Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., et al. Revision of the\n1295 6. Berends, C.J., de Boer, B. & van de Wal, R.S.W. Reconstructing the evolution of ice sheets, sea level, and 385\n1296 atmospheric CO2 during the past 3.6 million years. Clim. Past 17, 361–377 (2021).\n1297 7. Berends, C.J., Köhler, P., Lourens, L.J. & van de Wal, R.S.W. On the cause of the Mid-Pleistocene Transition.\n1299 8. Berger A. Support for the astronomical theory of climatic change. Nature 268, 44–45 (1977).\n1300 9. Berger, A. & Loutre, M. F. Insolation values for the climate of the last 10 million years. Quat. Sci. Rev. 10,\n1301 297–317 (1991).\n1302 10. Berger, A. & Loutre, M. F. Modeling the 100-kyr glacial–interglacial cycles. Glob. Planet. Change 72, 275–\n1303 281 (2010).\n1304 11. Berger, A. Long-term variations of daily insolation and Quaternary Climatic Changes. J. Atmos. Sci. 35, 2362-\n1305 2367 (1978).\n1306 12. Berger, A., Li, X.S. & Loutre, M.F. Modelling northern hemisphere ice volume over the last 3 Ma. Quat. Sci.\n1308 13. Bintanja, R. & van de Wal, R.S.W. North American ice-sheet dynamics and the onset of 100,000-year glacial\n1309 cycles. Nature 454, 869–872 (2008).\n1310 14. Chalk, T. B., Hain, M. P., Foster, G. L., Rohling, E. J., Sexton, P. F., Badger, M. P. S., et al. Causes of ice age\n1311 intensification across the mid Pleistocene transition. PNAS 114, 13,114–13,119 (2017).\n1312 15. Cheng, H., Edwards, R.L., Broecker, W.S., Denton, G.H., Kong, X. & Wang, Y.J. et al. Ice age terminations.\n1314 16. Cheng, H., Edwards, R.L., Sinha, A., Spötl, C., Yi, L. & Chen, S. et al. The Asian monsoon over the past\n1315 640,000 years and ice age terminations. Nature 534, 640–646 (2016).\n1316 17. Clark, P. U., Alley, R. B. & Pollard, D. Northern Hemisphere ice-sheet influences on global climate change.\n1318 18. Clark, P. U., Shakun, J. D., Rosenthal, Y., Köhler, P. & Bartlein, P. J. Global and regional temperature change\n1319 over the past 4.5 million years. Science 383, 884–890 (2024).\n1320 19. Clark, P.U., Archer, D., Pollard, D., Blum, J.D., Rial, J.A., Brovkin, V., Mix, A.C., Pisias, N.G. & Roy, M. The\n1321 middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in\n1322 atmospheric pCO2. Quat. Sci. Rev. 25, 3150–3184 (2006).\n1323 20. de Boer, B., Lourens, L. & van de Wal, R. Persistent 400,000-year variability of Antarctic ice volume and the\n1324 carbon cycle is revealed throughout the Plio-Pleistocene. Nat. Commun. 5, 2999 (2014).\n1325 21. Denton, G.H., Anderson, R.F., Toggweiler, J.R., Edwards, R.L., Schaefer, J.M. & Putnam, A.E. The Last\n1327 22. Drysdale, R.N., Hellstrom, J.C., Zanchetta, G., Fallick, A.E., Sánchez Goñi, M.F. & Couchoud, I. et al.\n1329 23. Dyez, K. A., Hönisch, B. & Schmidt, G. A. Early Pleistocene obliquity-scale pCO2 variability at ~1.5 million\n1330 years ago. Paleoceanography and Paleoclimatology 33, 1270–1291 (2018).\n1331 24. Elderfield, H., Ferretti, P., Greaves, M., Crowhurst, S., McCave, I.N., Hodell., D. et al. Evolution of ocean\n1332 temperature and ice volume through the mid-Pleistocene transition. Science 337, 704–709 (2012).\n1333 25. Feng, F. & Bailer-Jones, C. A. L. Obliquity and precession as pacemakers of Pleistocene deglaciations. Quat.\n1335 26. Foster, G. L. Seawater pH, pCO2 and [CO2- 3] variations in the Caribbean Sea over the last 130 kyr: A boron\n1336 isotope and B/Ca study of planktic foraminifera. Earth Planet. Sci. Lett. 271, 254–266 (2008).\n1337 27. Hagelberg, T., Pisias, N. & Elgar, S. Linear and nonlinear couplings between orbital forcing and the marine\n1338 δ18O record during the Late Neocene. Paleoceanogr. Paleoclimatol. 6, 729–746 (1991).\n1339 28. Hays, J. D., Imbrie, J. & Shackleton, N. J. Variations in the Earth’s orbit: Pacemaker of the Ice Ages. Science\n1340 194, 1121–1132 (1976).\n1341 29. Henehan, M. J., Rae, J. W. B., Foster, G. L., Erez, J., Prentice, K. C., Kucera, M., et al. Calibration of the boron\n1342 isotope proxy in the planktonic foraminifera Globigerinoides ruber for use in palaeo-CO2 reconstruction. Earth\n1344 30. Hobart, B., Lisiecki, L.E., Rand, D. et al. Late Pleistocene 100-kyr glacial cycles paced by precession forcing\n1345 of summer insolation. Nat. Geosci. 16, 717–722 (2023).\n1346 31. Hodell, D., Crowhurst, S., Lourens, L. et al. A 1.5-million-year record of orbital and millennial climate\n1347 variability in the North Atlantic. Clim. Past 19, 607–636 (2023).\n1348 32. Hönisch, B., Hemming, N. G., Archer, D., Siddall, M. & McManus, J. F. Atmospheric carbon dioxide\n1349 concentration across the mid Pleistocene transition. Science 324, 1551–1554 (2009).\n1350 33. Huybers, P. & Wunsch, C. Obliquity pacing of the late Pleistocene glacial terminations. Nature 434, 491–494\n1351 (2005).\n1352 34. Huybers, P. Combined obliquity and precession pacing of late Pleistocene deglaciations. Nature 480, 229–232\n1353 (2011).\n1354 35. Huybers, P. Early Pleistocene glacial cycles and the integrated summer insolation forcing. Science 313, 508–\n1355 511 (2006).\n1356 36. Huybers, P. Glacial variability over the last two million years: an extended depth-derived age model,\n1357 continuous obliquity pacing, and the Pleistocene progression. Quat. Sci. Rev. 26, 37–55 (2007).\n1358 37. Imbrie, J., Berger, A., Boyle, E.A., Clemens, S.C., Duffy, A., Howard, W.R. & Kukla, G. et al. On the structure\n1359 and origin of major glaciation cycles 2. The 100,000-year cycle. Paleoceanogr. Paleoclimatol. 8, 699–735\n1360 (1993).\n1361 38. Köhler, P. & van de Wal, R. S. W. Interglacials of the Quaternary defined by northern hemispheric land ice\n1362 distribution outside of Greenland. Nat. Commun. 11, 5124 (2020).\n1363 39. Konijnendijk, T.Y.M., Ziegler, M. & Lourens, L.J. On the timing and forcing mechanisms of late Pleistocene\n1364 glacial terminations: Insights from a new high-resolution benthic stable oxygen isotope record of the eastern\n1366 40. Lambert, F., Bigler, M., Steffensen, J. P., Hutterli, M. & Fischer, H. Centennial mineral dust variability in high-\n1367 resolution ice core data from Dome C, Antarctica. Clim. Past 8, 609–623 (2012).\n1368 41. Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M. & Levrard, B. A long-term numerical\n1369 solution for the insolation quantities of the Earth. Astron. Astrophys. 428, 261–285 (2004).\n1370 42. Legrain, E., Parrenin, F. & Capron, E. A gradual change is more likely to have caused the Mid-Pleistocene\n1372 43. Leloup, G. & Paillard, D. Influence of the choice of insolation forcing on the results of a conceptual glacial\n1373 cycle model. Clim. Past 18, 547–558 (2022).\n1374 44. Li, M. S., Huang, C. J., Ogg, J., Zhang, Y., Hinnov, L., Wu, H. C. et al. Paleoclimate proxies for\n1375 cyclostratigraphy: Comparative analysis using a Lower Triassic marine section in South China. Earth-Sci. Rev.\n1376 189, 125–146 (2019).\n1377 45. Liautaud, P. R., Hodell, D. A. & Huybers, P. J. Detection of significant climatic precession variability in early\n1379 46. Lisiecki, L. Links between eccentricity forcing and the 100,000-year glacial cycle. Nature Geosci. 3, 349–352\n1380 (2010).\n1381 47. Lisiecki, L.E. & Raymo, M.E. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records.\n1383 48. Meyers, S. R., Sageman, B. B. & Arthur, M. A. Obliquity forcing of organic matter accumulation during\n1385 49. Milankovitch, M. Kanon der Erdbestrahlung und seine Anwendung auf das Eiszeitenproblem. Roy. Serb. Acad.\n1387 50. Mitsui, T., Willeit, M., Boers, N. Synchronization phenomena observed in glacial–interglacial cycles simulated\n1388 in an Earth system model of intermediate complexity. Earth Syst. Dynam. 14, 1277–1294 (2023).\n1389 51. Raymo, M. E., Lisiecki, L. E. & Nisancioglu, K. H. Plio-Pleistocene ice volume, Antarctic climate, and the\n1390 global δ18O record. Science 313, 492–495 (2006).\n1391 52. Rial, J. A. Pacemaking the Ice Ages by frequency modulation of Earth’s orbital eccentricity. Science 285, 564–\n1392 568 (1999).\n1393 53. Rial, J. A., Oh, J. & Reischmann, E. Synchronization of the climate system to eccentricity forcing and the\n1394 100,000-year problem. Nature Geosci. 6, 289–293 (2013).\n1395 54. Ridgwell, A.J., Watson, A.J. & Raymo, M.E. Is the spectral signature of the 100kyr glacial cycle consistent\n1396 with a Milankovitch origin? Paleoceanography 14, 437–440 (1999).\n1397 55. Seki, O., Foster, G. L., Schmidt, D. N., Mackensen, A., Kawamura, K. & Pancost, R. D. Alkenone and boron-\n1398 based Pliocene pCO2 records. Earth Planet. Sci. Lett. 292, 201–211 (2010).\n1399 56. Sun, Y.B, Clemens, S.C., Morrill, C., Lin, X.P., Wang, X.L. & An, Z.S. Influence of Atlantic meridional\n1400 overturning circulation on the East Asian winter monsoon. Nature Geosci. 5, 46–49 (2012).\n1401 57. Thomson, D. J. Spectrum estimation and harmonic analysis. Proc. IEEE 70,1055–1096 (1982).\n1402 58. Toggweiler, J.R., Russell, J.L. & Carson, S.R. Midlatitude westerlies, atmospheric CO2, and climate change\n1403 during the ice ages. Paleoceanogr. Paleoclimatol. 21, PA2005 (2006).\n1404 59. Tzedakis, P.C., Crucifix, M., Mitsui, T. & Wolff, E.W. A simple rule to determine which insolation cycles lead\n1405 to interglacials. Nature 542, 427–432 (2017).\n1406 60. Tziperman, E., Raymo, M.E., Huybers, P. & Wunsch, C. Consequences of pacing the Pleistocene 100 kyr ice\n1407 ages by nonlinear phase locking to Milankovitch forcing. Paleoceanogr. Paleoclimatol. 21, PA4206 (2006).\n1408 61. Willeit, M., Ganopolski, A., Calov, R. & Brovkin, V. Mid-Pleistocene transition in glacial cycles explained by\n1409 declining CO2 and regolith removal. Sci. Adv. 5, eaav7337 (2019).\n1410 62. Yamamoto, M., Clemens, S. C., Seki, O., Tsuchiya, Y., Huang, Y. S., O’ishi, R. et al. Increased interglacial\n1411 atmospheric CO2 levels followed the mid-Pleistocene Transition. Nat. Geosci. 15, 307–313 (2022).\n1412 63. Yan, Y. Z., Bender, M. L., Brook, E. J., Clifford, H. M., Kemeny, P. C., Kurbatov, A. V. et al. Two-million-\n1413 year-old snapshots of atmospheric gases from Antarctic ice. Nature 574, 663–666 (2019).\n1414\n2\n6 1.1. The first comes from the consequences of uncertainties in the age models used for frequency\n7 and phase analyses: errors of 5 kyr equate to a full ¼ wavelength of the precession cycle and\n8 certainly pose a challenge to the goals of the present manuscript of assessing the timing offset of\n9 maximum change in d18O relative to orbital components of obliquity and precession. A major\n10 problem here is that most age models (with the exception of the Huybers depth-based model) already\n11 assume a consistent response to the orbital components, so that timing to some degree is “baked\n12 into” the underlying age models. As I pointed out in my earlier review, there is evidence that signal\n13 alignment and age model methods do not capture the precessional component well (see Liautaud et\n14 al.).\n16 the middle-Pleistocene Transition. Consequently, calculating the precession phase relative to glacial\n17 terminations is unnecessary. Therefore, concerns about the robustness of the calculated precession\n18 phase due to uncertainties in existing age models are no longer relevant.\n19\n20 1.2. The second come from critical ambiguities in interpreting the primary signal analyzed here, the\n21 marine d18O record. As we have seen from the recent Clark et al. (2024) paper, whose interpretation\n23 the temperature versus ice volume components.\n25 ice volume. Thus, in the revised manuscript, when we estimated the evolution phase between\n26 benthic δ18O andδ13C, we also estimated the evolution phase between the sea-level equivalent and\n27 benthic δ13C. The sea-level equivalent was deconvolved from the stacked marine benthic δ18O of\n29 temperature (GMST) has not been used in the revised manuscript because our focus is not on\n30 addressing the mechanism behind the middle-Pleistocene Transition.\n31\n32 1.3.The last is distinguishing trends in orbital response (and variance) from long-term (e.g. tectonic\n33 scale) trends. In other words, the use and understanding of what is meant by “trend” is critical.\n2\n35 and temperature) and the influence of their evolutions on the MPT.\n36\n37 1.4. The PDA results on spectral power are robust, but and consistent with previous looks at the\n38 evolution of d18O time series over this span (e.g., Lisiecki and Raymo, 2007). The increase in the100\n39 kyr power band simply confirms rise of well-known ~100 kyr periodicity in the late Pleistocene\n40 (Figure 2A), so what’s interesting are the relations to the tilt and precession bands.\n42 ~100-kyr bands, represented by the horizontal lines in the box whisker plots, indicates consistent\n43 negative correlations between the 136 ~100-kyr and ~21-kyr power ratios (Fig. 2D, F, H), with no\n44 correlation between the ~100-kyr and ~41-kyr power ratios (Fig. 2E, G, I). Is largely consistent with\n45 the results of Lisiecki (2010) who wrote in her abstract: I find that the time-dependent 100,000-yr\n46 power of eccentricity has been anticorrelated with that of climate since 5 Myr ago, with strong\n47 eccentricity forcing associated with weaker power in the 100,000-yr glacial cycle. I propose that the\n48 anticorrelation arises from the strong precession forcing associated with strong eccentricity forcing,\n49 which disrupts the internal climate feedbacks that drive the 100,000-yr glacial cycle. The finding in\n50 the present manuscript is what one would empirically (e.g. consistent with Lisiecki’s observation)\n51 expect: a negative correlation of the 100 kyr power in d18O records relative to a linear precessional\n52 response, as precession is modulated by the eccentricity cycle. Put together a largely linear\n53 precession response and a (admittedly puzzling) anti-correlated 100 kyr envelope and you get the\n54 obvservations of Figure 2E, G, and I.\n56 ~21-kyr and ~100-kyr power ratios. In the revised manuscript, we discuss the connection between\n57 our observed reversed relationship (between ~21-kyr and ~100-kyr power ratios) and the\n58 anticorrelation between ~100-kyr power in climate and eccentricity. Nevertheless, we disagree with\n59 the interpretation proposed by Lisiecki (2010). Instead, we attributed this anticorrelation between\n60 ~100-kyr power in climate and eccentricity to the carbon cycles. This is because the strong ~100-\n61 kyr glacial cycles occurred as long as δ18O/ice volume and δ13C were phase-coupled. Moreover, this\n62 phase coupling occurred only at the ~2.4-myr eccentricity nodes over the past 4 Ma. Since\n63 eccentricity modulates precession, the ~21-kyr and ~100-kyr power ratios should exhibit an anti-\n64 phase relationship.\n65\n66 1.5. I am puzzled by this part of the response: Moreover, the extracted ~21-kyr spectral power ratio\n67 of δ18Ow shows a negative correlation with the ~100-kyr spectral power ratio after the MPT (r=-\n68 0.89), whereas their relationship is positive before the MPT (r=0.66) (Fig. A2E). Figure 2 D, F and\n70 am completely mis-reading the figures. The negative correlation of 100:21 kyr power persists in\n71 all 3 time periods- the weaker correlation in the earlier slices likely comes from poorer signal: noise\n72 in d18O and/or age model errors where the time scale are more uncertain than for the last 600 kyr.\n74 and ~21-kyr power ratios across 57 paleoclimate archives suggests that the amplification of ~100-\n75 kyr power occurs at the expense of precession power (Fig. 2A, B, D). This power shift persists both\n76 during and before the MPT when glacial cycles predominantly exhibited ~41-kyr cyclicities (Fig.\n77 2A, B, F, H). This indicates that even in the ‘41-kyr glacial world’, other climate subsystems could\n78 show strong ~100-kyr cyclicities with a reverse correlation between the ~100-kyr and ~21-kyr\n79 power ratios. Thus, the nonlinear amplification of ~100-kyr power, caused by energy transformation\n80 from precession, is not exclusive to the glacial subsystem but is an intrinsic feature of the entire\n81 climate system. The authors appear to conflate the change in their estimated phase shift (Figure A2)\n82 with the consistent negative correlation of 21 kyr to 100 kyr power. I fail to understand how Figure\n83 4B can support a change in the relation of precessional power to 100 kyr power…\n85 between ~21 kyr and ~100 kyr power ratios. Regarding the Figure. 4, we reorganize it in the revised\n86 manuscript, please see Fig. 3 in the revised manuscript.\n87\n88 1.6.The definition of a glacial cycle follows the glacial roster of Tzedakis et al.56, which is derived\n89 based on the detrended LR04 benthic δ18O stack and the two-threshold scheme. The purpose of\n90 detrending is to remove the temperature contribution to benthic δ18O56, ensuring that the estimated\n91 precession and obliquity phases are minimally affected by temperature influences on benthic δ18O.\n93 timescales- detrending in no way removes this ambiguity. The author’s additional response in the\n94 rebuttal similarly seems to underestimate the role of changes in deep ocean temperature: many δ18O\n95 records from different oceans are used like in our study: many deep ocean d18O records from\n96 different oceans are used, which is expected to mitigate the effect of local deep-sea temperature and\n97 provide a more accurate reflection of global ice volume changes. The issue isn’t local/regional deep\n98 ocean temperature changes; it’s whole ocean changes which will be registered everywhere. Likewise,\n99 a shift in the non-linearity of precession (response lines 74-77) in no way uniquely resolves the role\n100 of temperature to ice volume; it could be explained equally well by the evolution of a temperature\n101 as opposed to an ice volume response.\n103 δ18O and δ13C. To mitigate the influence of local temperature on the benthic δ18O, we also estimated\n104 the evolution phase between the sea-level equivalent and benthic δ13C. The sea-level equivalent was\n105 deconvolved from the stacked marine benthic δ18O of Westerhold et al (2020) (Rohling et al., 2021).\n106\n107 1.7. My issue here is how one defines trends. I would challenge the authors or a reader to deduce\n108 that the INTERGLACIAL state preceding glacials prior to or after the MPT differ in any meaningful\n109 way- see Figures 3A, 3D, 5A, 5B. Note that I do not include a number of the other curves, as they\n110 are really model results applied to the underlying isotope data (e.g. latitude of glaciation) and do not\n111 contain new information. One can of course fit smoothed trends to data; the question is whether\n112 those fits illuminate more than they confuse. The big conundrum is that the starting point for glacial\n113 cycles with regard to temperature, d18O, and CO2 has been approximately the same before and\n114 after the MPT. The smooth trends simply capture the asymmetry of the response pre and post-MPT,\n115 with much stronger glaciations corresponding the last 800 kyr. However the text asserts:\n117 believed to gradually develop into partly marine-based configurations after the MPT Just how robust\n118 is the evidence of a “cooling background”? It appears weak to non-existent to me. Thus, this\n119 statement in the conclusion is largely a description of the observed increase in the depth of\n120 glaciations, not in fact an interpretation of underlying climatic mechanisms: We posit that around\n121 the MPT, cooling of temperature due to the glacial-stage pCO drawdown reached a critical\n2\n122 threshold at which the glacial response to precession shifts from linear to nonlinear, transferring\n123 precession power into the ~100-kyr band.\n125 thus the datasets cited in previous Figure 3&5 are not used in the revised manuscript.\n126\n129 2.1. I thank the authors for responding to my comments thoroughly and thoughtfully. The revised\n130 manuscript is of high quality, and I look forward to its publication. I am especially encouraged by\n131 the updated Figure 2, where the blank spaces between panels A and B show clear violin shapes,\n132 indicating the inverse relation between precession and eccentricity. Such a relationship does not\n133 seem to exist between panels A and C, which is consistent with the authors' argument.\n135\n136 2.2. Line 272 - change"}]

Decision letter