Synchronizing climate-carbon cycle heartbeats in the Phanerozoic vegetated icehouses

Domain: earth_science
Published: 2025-10-16
Editorial outcome: major_revision

Abstract
Earth experienced state-specific climate-carbon cycle feedbacks during the Late Cenozoic Ice Age (LCIA). Whether similar feedbacks existed in the penultimate icehouse, the Late Paleozoic Ice Age (LPIA), remains uncertain. Here, we present phase relationships between eccentricity-paced climate cycles and carbonate carbon isotope across ~337–300 Ma. Up to 307 Ma, low-latitude continental carbon reservoirs expanded during eccentricity-forced coolings, resembling the Oligocene and Miocene climate-carbon cycle dynamics. After 307 Ma, this relationship reversed, analogous to the Plio-Pleistocene dynamics. We attribute this reversal to the increasing importance of high-latitude biome dynamics, comparable to what occurred at 6 Ma in the LCIA. Paralleling LPIA (335–301 Ma) and LCIA (past 34 Myr) records using this event reveals quasi-synchronization in the interaction of astronomical forcing, carbon cycling and glacial events from onset to apex of two icehouses. We propose that, despite different boundary conditions, extraterrestrial forcing shaped the evolutionary trajectory of Phanerozoic vegetated icehouses.

Introduction
The Cenozoic Era underwent a significant transition from a greenhouse to an icehouse state beginning ~34 million years ago (Ma)1, after which global climate followed a long-term cooling trend, leading to our current icehouse with bipolar ice sheets2,3. On an astronomical (Milanković) timescale (20–400 kyr), the Late Cenozoic Ice Age (LCIA) exhibited state-dependent patterns of climate-carbon cycle interactions4,5. Between 34 and 6 Ma, low-latitude continental carbon reservoirs expanded mostly during the astronomically-forced cooling intervals, which changed at 6 Ma5, along with a decrease in atmospheric CO2 concentrations during late Miocene cooling6. Conversely, after 6 Ma to the present, global carbon cycle dynamics were likely influenced by the spatial redistribution of high-latitude biomes on orbital timescales due to their increasingly important dynamics5. As the tundra ecosystems have much smaller carbon storage capacities than the boreal forests7, the continental carbon reservoirs shrunk in the colder periods on eccentricity timescales5.

Earth’s penultimate icehouse, the Late Paleozoic Ice Age (LPIA, ~340–280 Ma), is viewed as a deep-time analog for the current glacial world8. Despite very different boundary conditions determined by continental configurations (Fig. 1a, b) and oceanic circulation patterns9, the potential for parallels is suggested by similarities between the two icehouse periods. First, the LPIA archives a dynamic glacial evolution under low atmospheric CO2 concentrations comparable to those of the Pleistocene glacial cycles to present-day10, and secondly the paleo-continents were also covered by complex vegetation ecosystems11. Rich palaeobotanical evidence from the Pangaean tropics reveals repeated shifts in vegetation from glacial wetland forests (e.g., arborescent lycopsids and medullosan dominated) to interglacial seasonally-dry flora (e.g., cordaitalean, conifer, and marattialean tree fern dominated) throughout the late Carboniferous12, in response to eccentricity-paced variations in solar insolation and pCO2 10,13,14. This was a finding consistent with simulations made using the Earth System15 and ecologic models16. Net positive terrestrial carbon sequestration occurred during late Carboniferous glacials driven by the predominance of wetland coal forests10. At 307 Ma, arborescent lycopsids abruptly disappeared from Euramerican biomes17,18, coincident with the onset of long-term tropical aridification19, remaining only in Cathaysia (China)20,21 that held ever wet conditions throughout the late Paleozoic19. The continental-scale aridification likely intensified the repeated vegetational turnovers that characterized the remainder of the LPIA12. Thus, it is possible that the widespread loss of tropical coal forests may have led to a change in the response of the global carbon cycle to astronomical forcing. This hypothesis, however, remains untested.

Here, we report a high-temporal resolution (~25 kyr) δ13C time series developed using ~1500 bulk carbonate samples collected from a 250-m deep-marine succession (Narao section) in South China (Fig. 1a and Supplementary Fig. 1), which we calibrated to an astronomical time scale constructed by using the cyclostratigraphic analysis on 5000 gamma ray (GR) data (Fig. 1f; see “Methods” section). This new δ13Ccarb, integrated with a GR time series, spans the interval from the onset (~335–330 Ma8) to the apex (300 ± 5 Ma22; referred to as 300 ± 1 Ma23 herein) of the LPIA, providing new insight into the response of global carbon cycling to astronomical forcing during the penultimate icehouse.

Results
Cyclostratigraphic analysis
The Viséan to Gzhelian sediments of the Narao section were deposited on an open marine deep-water (carbonate) slope24,25 and are characterized by an average sedimentation rate of ~0.7 cm/kyr. This estimate is derived from dividing the 187.81 m of stratigraphy between the basal Serpukhovian (41.8 m24) and basal Gzhelian (229.61 m26) (Fig. 1c), by a total duration of ~26.7 Myr27. Increased terrestrial input during sea-level falls would have provided more clay to the studied region28, resulting in higher GR values. Time series analysis of the new GR data reveals significant wavelengths of 4.5–2.8 m, 0.82 m, 0.63 m, 0.32–0.27 m, 0.22–0.20 m, 0.18–0.16 m, and 0.12–0.11 m, above a 99% Chi-squared confidence level (Fig. 2a; see “Methods” section), after removing a >8 m trend (red line in Fig. 1f). These cycles can be clearly observed in the raw data series and data after applying different ways of detrending (Supplementary Fig. 2). Given the estimated sedimentation rate, the 4.5–2.8 m cycle most likely corresponds to the 405-kyr eccentricity cycle. Although the stability of the 405-kyr eccentricity cycle in deep time remains a topic of ongoing debate29,30, geochronological evidence supports its relatively stable period during much of the Carboniferous31–33, which is further substantiated by results derived from 64 solutions in the ZB23 astronomical solution34 that currently are available (Supplementary Fig. 3).

Ninety 4.5–2.8 m cycles were recognized through bandpass filtering of the GR time series (Fig. 1d). Calibrating these cycles to the 405-kyr eccentricity period reveals dominant spectral peaks at 405 kyr (long eccentricity), 134.5 kyr (short eccentricity), 32.5–29.8 kyr (obliquity), and 20.5 kyr, 18.6 kyr and 15.6 kyr (precession; Fig. 2c, d), corresponding to average periods of the Carboniferous Milanković cycles (i.e., ~405 kyr, ~132 kyr, ~31.5 kyr, ~20.2 kyr, ~19.3 kyr and ~16.7 kyr)34 (Supplementary Fig. 3). The calculated change in sedimentation rates based on astronomical calibration (Fig. 1e) tracks with those obtained using an objective statistical methodology, evolutionary time scale optimization35 (Supplementary Fig. 4; see “Methods” section), providing confidence in our astronomical calibration. The Kasimovian/Gzhelian boundary at Narao succession was assigned a numerical age of 303.7 ± 0.10 Ma27, and an uncertainty of ± 0.20 Myr (half a 405 kyr cycle) was given for the astronomical calibration. This yields a floating astronomical timescale (i.e., ~337–300 Ma) for the interval spanning the stepwise onset to peak of the LPIA, with a combined uncertainty of ± 0.3 Myr (Fig. 3a; Supplementary Text 1). An abrupt increase in the sedimentation rate (dark blue arrow in Fig. 1e), accompanied by the enhanced terrestrial inputs (Fig. 1f), is observed at a stratigraphic depth of ~196 m, which has been astronomically calibrated to an age of ~306.9 Ma. We suggest that it may represent a regional response to glacio-eustatic sea-level lowering during the early Kasimovian33, occurring within the context of secular cooling towards the Carboniferous–Permian transition22,23.

Astronomical cycles of million-year-scale can be determined by the amplitude modulation assessment using Hilbert transforms36 (see “Methods” section). Amplitude modulation analysis was conducted on the bandpass-filtered ~32-kyr obliquity (Fig. 3h) and ~100-kyr short eccentricity signal (Fig. 3g). Results show main periods of ~1.17 Myr and ~2.36 Myr for the long-period obliquity and eccentricity amplitude modulation cycles, respectively (Supplementary Fig. 5c, d). The periods of these million-year-scale astronomical cycles during ~337–300 Ma are indistinguishable from those in the astronomical solution (e.g., ZB18a) for the past 36 Myr37 (Supplementary Fig. 6).

Astrochronologic constraint on the Carboniferous climate-carbon cycle history
The astronomically calibrated GR time series for the Narao succession provides a high-resolution relative sea-level history (Fig. 3f), that aligns on the million-year time scale with inferred sea-level reconstructed for the late Viséan (~333 Ma) through the latest Carboniferous (~300 Ma) from other low-latitude records (i.e., Donets Basin, Ukraine33 and North America38; Fig. 3i) (Supplementary Text 2). On astronomical time scales, the onlap (transgression)–offlap (regression) record from the Donets Basin and the GR time series from the Narao section display a constant anti-phase relationship within the 1.17 Myr (–160° ± 18°) and 400-kyr (–173° ± 23°) bands (Supplementary Fig. 7a), indicating the potential synchronicity of the astronomically-paced eustasy at these low latitude locations. On the ~1.2 Myr band, the long-period obliquity amplitude modulation cycle shows a nearly consistent anti-phase relation with the Narao GR time series (–135° ± 18°; Supplementary Fig. 7b), and an in-phase relationship with the onlap–offlap record (40° ± 11°; Supplementary Fig. 7c). These phase relationships support a linkage between sea-level falls and obliquity nodes (minima of long-period obliquity amplitude modulation cycles) during the Carboniferous, this is congruent with findings for the LCIA39.

The astronomically calibrated marine δ13Ccarb time series from the Narao section is globally correlated, and it provides high temporal resolution insight into the nature of carbon cycle perturbations. Three positive δ13Ccarb excursions during the late Visean–earliest Serpukhovian interval (red bars in Fig. 3e) overlap with the proposed timings of the stepwise onset of the LPIA8 (Fig. 3b). Astrochronologic constraints for England shallow marine carbonates from the U.K. indicate a basal late Asbian age of 334.48 ± 0.35 Ma40. Thus, the δ13Ccarb positive excursion at ~335–334.6 Ma is likely time-equivalent with a cooling pulse inferred from an earliest Asbian positive shift in biogenic apatite δ18O values and ice buildup41, suggesting a possible mechanistic link between the onset of the LPIA and the enhanced carbon burial. Three subsequent and prominent negative δ13Ccarb excursions (black bars in Fig. 3e) at the end-Serpukhovian (excursion of ~ –1.2‰; ~324.9–324.4 Ma), the end-Moscovian (~ –2‰; ~307.3–307.0 Ma) and the end-Kasimovian (~ –2‰; ~303.8–303.7 Ma) were observed in the Narao record, which can be correlated to the excursions archived in western paleotropical Pangaea successions42,43 (Supplementary Fig. 8). Collectively, our new record provides an astronomically-calibrated, high-resolution stratigraphic reference for the evolution of low-latitude glacio-eustasy and changing carbon cycle dynamics for the time interval from ~337–300 Ma.

Eccentricity-scale climate-carbon cycle
A time-evolutive phase analysis between the astronomically-calibrated GR and δ13Ccarb time-series was used to differentiate in-phase and anti-phase behavior in the 405-kyr and 100-kyr eccentricity bands (Fig. 3j). For the interval of 337 to 307 Ma and on the 405-kyr eccentricity time scale, GR and δ13Ccarb exhibit oscillations around an in-phase relationship, indicating that maxima in δ13Ccarb are associated with the eccentricity-paced glacial-eustatic sea-level falls. The in-phase behaviour inverts to anti-phase at ~307 Ma (purple shading in Fig. 3j), along with an abrupt increase in sedimentation rate and GR values (Fig. 1e, f). It suggests a major shift in the response of the continental carbon sink(s) to the glacial advances switched from expansion to contraction. This shifted phase relation occurred almost simultaneously in the short eccentricity band (blue shading in Fig. 3j and Supplementary Fig. 9). A similar shift in phasing is observed in another carbonate slope succession in South China (Naqing section; Supplementary Fig. 10). The shifted phase relation can also be observed via the stationary phase relation analysis between δ13Ccarb from Narao and onlap–offlap record from the Donets Basin (Supplementary Fig. 7f, g), likely reflecting a global signal for the reorganization of climate-carbon cycle system on eccentricity timescale.

Based on these observed phase relationships and robust paleobotanical evidence12,14, we propose the climate-carbon cycle dynamics during the period of 335–307 Ma. The dynamics were mechanistically linked by the expansion of wetland coal forests and increased carbon sequestration and storage of isotopically light carbon (12C) through enhanced peat burial12 during eccentricity minima. Earth system44 and coupled climate and ecosystem modeling45,46 indicate that the cooler and wetter conditions of Carboniferous eccentricity minima promoted the widespread expansion of wetland forests. In turn, increased burial of 12C-enriched organic matter on land during eccentricity-forced glacial periods would have contributed to an increase in the δ13C of marine dissolved inorganic carbon, as is observed in the increases in the δ13Ccarb at maxima of the eccentricity-paced GR cycle.

In contrast, we propose that for the subsequent interval between ~307–300 Ma, eccentricity minima led to cold and dry climate conditions, resulting in a shrinkage of the 12C-enriched continental carbon reservoirs at high latitudes. This shift in carbon cycle dynamics coincides with a global cooling (~307–304.5 Ma), characterized by ~50% drop in pCO2 10,47 (Fig. 3d), a short-lived but intense glaciation33 (Fig. 3f, i) (Fig. 3c). The subsequent increase of two to three times in pCO2 during the Kasimovian–Gzhelian transition (Fig. 3d) appears to be insufficient to reverse this dynamic (Fig. 3j). We speculate that the center of terrestrial carbon burial may have shifted from the low to high latitudes in response to the onset of long-term tropical aridification at the close of the Middle Pennsylvanian. Key evidence supporting this hypothesis includes precipitous loss of coal forests and peat burial by the end of the Middle Pennsylvanian14,48, which served as the primary source of terrestrial carbon sequestration prior to 307 Ma20. Furthermore, ecosystem modeling of late Carboniferous terrestrial ecosystems indicates major changes in high-latitude forest cover during glacial and interglacial extremes in particular the higher tolerance of drought-adapted vascular plants to freezing conditions16, highlighting the potential for the increasing importance of different high-latitude carbon reservoirs, e.g., tundra and boreal forests, in response to climate deterioration and amelioration. We suggest that global carbon cycle dynamics could have been substantially influenced by such spatial redistribution in high-latitude biomes. Glacial advances during attenuated eccentricity cycles would have caused an upsurge in tundra areal extent at high latitudes, accounting for the anti-phase relationship observed between δ13Ccarb and GR after ~307 Ma.

Although we tend to support the biome areal competition hypothesis, other studies have stressed the relevance of organic carbon burial in the ocean during the LPIA49–51. Iron-rich dust flux likely increased substantially during late Carboniferous glacials, which may have stimulated marine productivity that would have enhanced marine carbon burial49,51. This hypothesis was proposed to explain the secular decoupling of continental weathering and pCO2 after ~303 Ma50. Despite more vigorous ocean circulation at glacial phases during late Carboniferous52, it remains uncertain as to which boundary conditions would have driven reversals in the response of organic carbon burial in the ocean to astronomical forcing. Given the great capacity of the terrestrial biosphere to influence the LPIA carbon cycling10, we view the areal competition between low- and high-latitude biomes as a more plausible explanation for the observed shift in phase relationships of the eccentricity-scale climate-carbon cycling at ~307 Ma. Future investigations will be necessary to determine whether climate-carbon cycle feedbacks shifted along major deglaciation events and vegetation turnovers during the early Permian period53.

Climate-and carbon-cycle evolution in icehouse conditions
The shift in the carbon-climate cycle dynamics was linked to the growing significance of competing effects between high-latitude biomes during the LCIA5. This shift occurred at 6 Ma5, coinciding with maxima in ~1.2-Myr and ~2.4-Myr cycles (Fig. 4a2–4) under a declining atmospheric pCO2 concentration6 (Fig. 4a5). A comparable scenario is revealed at ~307 Ma for the LPIA (Fig. 4b2–5). Thus, this event occurred in high-amplitude insolation variations superimposed on secular cooling trend in both icehouses. This similarity indicates that the increased importance of high-latitude vegetation dynamics may be a critical threshold for the secular evolution of the vegetated icehouses. Since the late Paleozoic Era, continental areas in the Northern Hemisphere have consistently increased in extent relative to those in the Southern Hemisphere9. Specifically, the Pangaean supercontinent had large areas of southern high-latitude land and the Central Pangaean Mountains extended across the tropical continent (Fig. 1a), which is substantially different with the late Cenozoic paleogeography9 (Fig. 1b). The discovery in this study of analogous thresholds, despite the different paleogeographic background, enhances confidence in a more comprehensive comparison of the evolutionary patterns between the two icehouses, e.g., climate rhythmicity and transient events.

Our study documents the climate- and carbon-cycle evolution of the two icehouses follows similar rhythms from their onset to apex (Fig. 4). Astronomically tuned deep-sea benthic foraminifer δ18O and δ13C records are utilized to track glacial variation (Fig. 4a6) and the global carbon cycle (Supplementary Fig. 11) during the LCIA3. Paralleling the δ18O data (36–0 Ma) and GR series (337–301 Ma) shows a nearly in-phase relation on the pronounced “heartbeats” at ~405-kyr, ~1.2-Myr and ~2.4-Myr periods (Supplementary Fig. 7h), pacing with the multiple astronomical cycles (Fig. 4). The δ13C time series for the two icehouses show an in-phase relation for 405-kyr eccentricity cycles (Supplementary Fig. 7i), likely revealing the quasi-synchronicity of the climate-carbon cycle. It strengthens our hypothesis that astronomical forcing mediated the climate-carbon cycle dynamics on kilo-year time scale in the two icehouse periods. No stable relationship is observed between the δ13C time series of the two icehouses on million-year time scales, probably indicating a more complex carbon cycle dynamic on longer time scales.

We also observe that the timing of certain transient climatic events in the Carboniferous and late Cenozoic eras was approximately comparable, following the transition of Earth’s climate system into an icehouse state (Fig. 4a6, b6). Since the onset of the LCIA at ~34 Ma during the Eocene–Oligocene climatic transition, global climate was characterized by stepwise cooling occurring around ~23 Ma (Mi-1 Glaciation at Oligocene/Miocene boundary), ~14 Ma (onset of the middle Miocene climate transition), ~7 Ma (initiation of late Miocene cooling) and ~3.3 Ma (Marine Isotope Stage M2 in the late Pliocene, indicating the onset of the Pleistocene Ice Ages)2,3. These cooling trends were punctuated by distinct warming intervals or events, such as the late Oligocene warming (~26.3–23.7 Ma), the Miocene Climatic Optimum (~17–14 Ma) and the Middle Pliocene Warm Period (~3.3–3.0 Ma)2,3 (Fig. 4a6). The climate evolution of the LPIA (Fig. 4b6) exhibits a similar pattern to that observed in the LCIA. Stepwise expansion of continental ice sheets during the LPIA was associated with notable cooling events in the late Visean (~335 Ma), mid-Carboniferous (~324 Ma), early Moscovian (~315 Ma), early Kasimovian and latest Carboniferous (~302 Ma) glaciations. Overall cooling was interrupted by the warming intervals in Serpukhovian, Bashkirian and Kasimovian–Gzhelian transition8,33,38 (Fig. 4b6).

With a comprehensive understanding of the Cenozoic climate, there is a tendency to attribute the Cenozoic cooling/warming events to multiple contributing factors. Revisiting the Eocene–Oligocene climatic transition as an example, its primary drivers include the tectonic-oceanic hypothesis and the CO2 hypothesis. The first hypothesis posits that the enhanced thermal isolation of Antarctica, due to the widening of oceanic passages (such as the Drake and Tasman seaways), facilitated the expansion of Antarctic glaciation54. The second suggests that the dominating factor in the latest Eocene ice expansion is the decline in atmospheric CO2 levels down to a specific threshold55. However, the predominant forcing of climate changes remains a subject of debate, due to the intricate feedback mechanisms among the ocean, atmosphere, and Earth’s surface. Classically, climate regulation on Earth is predominantly associated with atmospheric pCO2 levels, which are influenced by the variations in carbon sources (e.g., volcanism) and sinks (e.g., continental weathering, organic carbon burial on land or in the ocean)56. There is no clear consensus regarding the factors responsible for long-term climatic shifts; debates often focus on whether carbon sources or sinks have been more influential. For instance, studies have connected both carbon sources or sinks to major climate changes of the LCIA57–59, and this is also the case for the LPIA47,53,60. Both factors were largely controlled by plate tectonics, resulting in gradual changes that predominantly occurred unidirectionally over multi-million-year timescales.

Beyond atmospheric CO2 levels, extraterrestrial factors may play a significant role in maintaining quasi-synchronization of these transient climatic events between two icehouses. We propose that tectonically controlled CO2 variations brought the climate to threshold conditions, and then orbital forcing may have caused these transient climatic events (thereby determining their timing), by pushing the system beyond the critical thresholds. The shifts in other boundary conditions—such as paleogeography and topography, oceanic gateway locations and bathymetry—cannot fully account for the climatic quasi-synchronization due to their independent changes within two icehouses. Efforts have been made to link astronomical forcing and cryospheric responses in the LCIA2,61–63, which would shed light on the climatic counterparts or analogs that occurred in the LPIA. The rapid cooling and glacial building maybe associated with periods of low eccentricity and low-amplitude change in obliquity (see blue vertical lines in Figs. 4a2–3, 4b2–3), permitting the accumulation of summer snow becoming glaciation. Such orbital configurations maybe one of the important forcings to tip the balance towards glaciation during intervals of declining CO2 levels below the threshold, initiating the glaciation (see the glaciation/cooling events in Figs. 4a6, 4b6). On the contrary, the secular cooling trends would have been interrupted by high eccentricity and 1.2-Myr obliquity-forced deglaciations (see red vertical lines in Figs. 4a2–3, 4b2–3). During eccentricity maxima, highly seasonal precipitation would have led to prolonged dry seasons, when the terrestrial organic carbon may have undergone oxidization due to seasonally drier climates64. Seasonal extremes would have been further amplified during orbital configurations of maxima in all 405-kyr eccentricity, 1.2-Myr obliquity, and 2.4-Myr eccentricity cycles, likely acting as positive feedbacks for the transient warming events, e.g., a thermal maximum at ~15.6 Ma during the Middle Miocene climatic optimum65 and the Tortonian thermal maximum at ~10.8 Ma66 (Fig. 4a7) in the LCIA, and short-term warming in middle Moscovian (~311.8 Ma33) in the LPIA.

Implications for the deep-time climate change in icehouses
Our proxy-based study demonstrates that the development of the main phases of the LCIA and LPIA follows a comparable climate trajectory, even though significantly different boundary conditions exist in two icehouses. This trajectory involves secular trends superimposed with multiple astronomically forced climate-carbon cycles, transient climatic events and certain critical threshold. The LPIA should be regarded as a climate analog for the LCIA. Trends, rhythms, and aberrations in climate changes in one icehouse may provide critical information for their counterpart, thereby elucidating the dynamics, such as forcings, feedbacks, responses and climate thresholds, underlying the climate-carbon cycle linkages. The magnitude of environmental perturbations during the LPIA may have been amplified due to the super-continentality at the time, providing important insight into the functioning of physical, geochemical, and biological processes under a more dynamic glacial/interglacial variation. Moreover, the demise of the LPIA and turnover to greenhouse condition holds critical information regarding processes driving ice sheet destabilization, large-scale hydrological changes, and key feedback mechanisms influencing biological responses.

Methods
Sampling and proxy data measurements
GR intensity was measured directly on the succession using an RS-230 SuperSpec portable spectrometer, with a 150-s recording time for one data. A total of 5,000 data points were taken from the Narao section, measuring in total 250 m (13–263 m) with a sample distance of 0.05 m. Rock samples were consistently collected at ~10-cm intervals along the horizons where GR measurements were conducted. No lithological bias occurred during the field sampling process. Based on the established astronomical timescale, the temporal resolution of the δ13Ccarb measurements is ~25-kyr. The Narao δ13C record was generated using samples of micritic limestone matrix drilled from fresh rock surfaces. Carbonate powders were reacted with supersaturated 100% phosphoric acid at 70°C using a Kiel IV Carbonate Device coupled with a Finnigan 253 plus mass spectrometer. The maximum number of samples measured each time is 88, of which 18 are standard samples. All isotope values are reported relative to Pee Dee Belemnite (PDB) using standard delta notation; analytical precision for both δ13C and δ18O is better than 0.1‰. The weak relation between δ18O and δ13Ccarb (Pearson’s r = 0.108) indicates that the influence of diagenesis on the δ13C signal appears to be negligible. 1307 samples from the Narao section were analyzed for δ13Ccarb, and integrated with published δ13Ccarb data (n = 193)24,25,67, which were used for further study. Bulk δ13Ccarb from the adjacent Naqing succession was also conducted to further validate the findings obtained from the Narao succession. For the Naqing section, 475 samples were analyzed for δ13Ccarb and integrated with 245 published data points24,67. The δ13Ccarb time series were calibrated using a published astronomical timescale28 for cyclostratigraphy analysis (Supplementary Fig. 10).

Cyclostratigraphic analysis
A trend (>8 m) in the GR data was calculated using the LOWESS method (locally weighted scatterplot smoothing method) and subtracted (Supplementary Fig. 4a). It is essential to identify regular cyclicity in successions that lack high-resolution dating. We employed smoothed window averaging (SWA)68 to determine the spectral background. The SWA for red-noise spectra represents an empirical, non-parametric approach that makes no assumption about the underlying cause of the sloping spectral background68. The final results of the SWA background fitting for the time series power spectrum include the standard Chi-squared 99% confidence level. In cyclostratigraphic studies, confidence levels of >99.9% may effectively minimize the risk of Type I (false positive) errors; however, such excessive stringent confidence levels may increase the likelihood of Type II (false negative) errors69. In this study, the confidence levels of 99–99.9% were employed to validate the presence of significant cycles. These levels could balance the probability of incorrectly rejecting and incorrectly accepting a null hypothesis of no astronomical forcing69. Sliding-window spectral analysis was conducted using evolutive Fast Fourier Transform spectrograms (eFFT)70. These analyses were conducted in Acycle71.

Astronomical calibration was used to transform data from depth to time domains by fixing cycles to a constant duration. Calibration was conducted using the inferred 405-kyr eccentricity cycles. We conducted evolutionary time scale optimization (eTimeOpt)35 using Astrochron package72 to track the variable sediment accumulation rates throughout the Narao (Supplementary Fig. 4b–d). Long-term cycles modulate short-term ones according to astronomical theory29, and TimeOpt is a statistical method to test each sedimentation accumulation rate using this well-known phenomenon. Short eccentricity amplitude modulation was used here, with a frequency of 1/405 for the long eccentricity band, and frequencies of 1/132, 1/124, 1/100 and 1/95 for the short eccentricity bands (ZB23 solution34; Supplementary Fig. 3). The main procedure of eTimeOpt includes: (1) the detrended Narao GR data are filtered using Taner bandpass (low- and high-frequency cut-off was here used as 1/160 and 1/80, respectively), and Hilbert transformed to isolate the short eccentricity index amplitude envelope, which is then linearly regressed on a synthetic dataset with long eccentricity. The regression result is shown with the squared Pearson correlation coefficient r2envelope (Supplementary Fig. 4b); and (2) the data are linearly regressed on a synthetic dataset with long and short eccentricity index frequencies, and the regression result is shown together with the squared Pearson correlation coefficient r2power (Supplementary Fig. 4c). The product (r2opt = r2envelope × r2power) can combine information from both r2envelope and r2power (Supplementary Fig. 4d). The test sedimentation rates (n = 200) are from 0 to 2 cm/kyr, with sliding windows of 8 m.

Amplitude modulation (AM) analysis was conducted as follows36: (1) the 405-kyr-calibrated GR time series was bandpass-filtered to isolate the interpreted short eccentricity or obliquity signal, and (2) the amplitude modulations were obtained from the filtered signal by application of the Hilbert transform.

A relationship (e.g., cause and effect) could be evaluated using the multivariate time series analysis70. Cross-spectral analysis to determine the coherency and phasing between the two separated records was performed using the Blackman–Tukey cross-spectral analysis in Analyseries 2.0.873. To identify the exact age for the changed phase relation, we conducted the cross-spectral analysis between GR and δ13Ccarb using a Bartlett window, with sliding windows of 4 Myr and moving steps of 0.25 Myr. For the results with error bars >95% (or >90% and >85%), non-zero coherence is higher than 0.549641 (or 0.476565 and 0.425021), and the error estimation on the power spectrum is between 0.489097 and 3.06605 (or between 0.54708 and 2.52863, and between 0.590019 and 2.24219).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In their manuscript (ms), Fang et al. present new gamma ray (GR) and d13C bulk-rock carbonate measurements of the upper Carboniferous Narao section (China) to investigate astronomical climate forcing and carbon cycle dynamics during the Late Paleozoic Ice Age. Using time-series analysis, the authors identify a shift from in-phase to anti-phase behavior between variations in GR (a proxy for glacio-eustatic sea level (SL) change) and d13C on eccentricity timescales around the Moscovian-Kasimovian boundary (~307 Ma), that is, from coupled low SL-high d13C to low SL-low d13C values. This shift is interpreted as a reorganization of Earth’s climate-carbon cycle system in response to orbital forcing – from a growing terrestrial biosphere during glacial intervals to a shrinking terrestrial biosphere during glacial intervals (?) – analogous to the shift from in-phase to antiphase behavior between benthic foraminifera d18O and d13C during the Late Cenozoic.\nThe proposed phase shift and analogy to the Late Cenozoic Ice Age is fascinating, if true. However, I see a couple of potential issues that need to be addressed first and/or require further clarification before I can recommend publication:\n1. Insufficient info on d13C. At the moment, I cannot properly assess the robustness of the proposed phase relation shift between GR and d13C on eccentricity timescales, because no results of cyclostratigraphic analysis (Milanković hypothesis testing) are presented for the Narao d13C record (at least as far as I am aware). If the goal is to investigate phase relationships between changes in global SL/ice volume and the carbon cycle in response to orbital forcing, one also needs to show (convincingly) that both components are indeed controlled by orbital forcing in the first place. In addition, the authors need to provide more information about the sampling strategy – in particular, how are d13C sample depths related to GR depths? Currently, the Methods section mentions that 5000 GR data points from the Narao section were taken “..with an average sample distance of 0.05 m”, and a total of 1307 samples for d13C (lines. 473-477). In the Excel data file, however, I find d13C data interpolated also at 0.05 m steps, so no raw depths are given. What was the original sampling density of d13C, and how is it related to GR? If, for example, d13C data are offset from GR by 1-2 m, this could artificially introduce a change from in-phase to anti-phase behavior. Also, I would recommend to include a figure in the ms that shows to readers a zoom-in of a typical/well-expressed eccentricity cycle in GR and d13C from where these proxies are in phase, and one from where they are out of phase. (So in addition to the phase-relationship plot in Fig. 3j).\n2. Uncertainty in the astronomical timescale –> this should definitely get more attention in the ms (currently uncertainty is not discussed at all). Specifically, what would be the consequences of errors in the astrochronologic framework for the inferred phase relationships? First, tuning was done based on the long (~405 kyr) eccentricity cycle, which the authors refer to as “.. a stable astronomical component throughout the geological history.” However, see the paper by Zeebe and Lantink (2024), which shows that the long eccentricity cycle can become unstable on 10^8-10^9 year scales. Second, I find it very hard to believe that a complete record of Milanković forcing of ~37 Myrs long is preserved in a single section, especially considering the nature of the sedimentary environment (“deep-marine slope deposits”, see l. 602). I think the longest deep-sea (composite) record used in Westerhold et al. (2020) is about ~20 Myr. The Newark-Hartford record of Olsen et al (2019) is also on that order; astronomical signal quality and the density of independent chronostratigraphic constraints is much higher for these records than is the case for the Narao data (which has zero independent U-Pb constraints…). What about the possibility of turbidites or hiatuses associated with sea level change? In Extended Data Fig. 1c, I can clearly see irregular bedding above cycle E7 and E32. In Fig. 2, the EHA plot shows a disruption of the (inferred) long eccentricity-scale signal around 0-30 m and ~165-185 height (just before the inferred phase shift occurs).\n3. A global signal. It’s a big leap from GR and d13C bulk-rock carbonate data of one section (Narao) to interpreting changes in the global climate-carbon cycle. The GR and d13C curves from the Narao section have been correlated to time-equivalent records from Donets and the GTS2020 site, indicating they may be representative for global-scale processes (Fig. 3i, Extended Data Fig. 7). However, the Donets and GTS datasets are of a lower resolution and have not been analyzed at (and correlated on) Milanković timescales. So while the general trend in the Narao series may indeed reflect a global signal, I wouldn’t dare to say that the correlation to Donets and GTS provides direct support for global-scale climate-carbon cycle processes on orbital timescales. (Cyclostratigraphic study has supposedly been conducted also one data from the lateral Naqing section – see l. 136-137 and Extended Data Fig. 8 - but the actual analyses are not shown in the current ms and a precise cyclostratigraphic correlation to Narao is also lacking).\nMinor points:\nl. 42-43: “low-latitude continental carbon reservoirs expanded mostly during the astronomically-forced cool intervals (4) -> Is the reference to Kirkland-Turner (2014) appropriate here? Because from my understanding, that paper talks about a “gradual expansion of vegetated terrestrial biosphere during warm intervals”.\nl. 65 “to test this hypothesis” -> elaborate.\nl. 74-75: how is it possible that no significant short eccentricity-scale peaks are identified in the raw depth series, only after tuning to a 405 kyr cycle?\nl.83-86: “The calculated change in sedimentation rate…tracks with those obtained from two objective statistical methodologies, i.e.,…. and evolutionaly correlation coefficient.” I’ve attached a document to this review from prof. Stephen Meyers (personal comm.). In this document Meyers explains why he thinks that the eCOCO method of Li et al. (2018) is unfortunately flawed. Meyers has given me permission to share this document with the authors of the present ms. I haven’t discussed the current ms (Fang et al.) with Meyers.\nl. 93-96: the ~1.2 and ~2.4 Myr Grand cycles are not stable in deep time.\nl. 108-111: I don’t understand how the long-term obliquity AM cycle can show an anti-phase relation with the Narao GR series and in-phase relation with the onlap-offlap record of the Donets, if Donets and Narao GR record display a constant antiphase relationship (l. 105-106).\nl. 134-135: “This shifted phase relation is delayed by ~0.5 Myr for the 95-132 kyr eccentricity band..” -> in Extended Data Fig. 2, I see no concentration of spectral lower in the short eccentricity range for Kasimovian-Gzehlian. Are the result of phase analysis then reliable for this interval?\nReferences:\nOlsen, Paul E., Jacques Laskar, Dennis V. Kent, Sean T. Kinney, David J. Reynolds, Jingeng Sha, and Jessica H. Whiteside. \"Mapping solar system chaos with the Geological Orrery.\" Proceedings of the National Academy of Sciences 116, no. 22 (2019): 10664-10673.\nZeebe, R. E., & Lantink, M. L. (2024). A Secular Solar System Resonance that Disrupts the Dominant Cycle in Earth’s Orbital Eccentricity (g 2− g 5): Implications for Astrochronology. The Astronomical Journal, 167(5), 204.", "reviewer_2": "This study presents an impressive record of novel bulk carbonate δ18O measurements, with a new astronomical time scale based on gamma-ray data. The authors argue that the Late Paleozoic Ice Age shows interesting parallels to our current glacial world, and that this record can be used to test this hypothesis. I think the work is significant and warrants publication. However, there are some major issues with the astronomical solutions and data analysis that should be addressed, requiring major revisions.\nSee attached PDF for the Major and Minor comments. I have also included a quick reproduction of the spectral analysis.", "reviewer_3": "The manuscript uses a new dataset of densely sampled marine carbon isotopes from the Carboniferous of China to argue for the cyclicity of orbital forcing to be driving the patterns and then draws comparison between Paleozoic and Cenozoic ice ages to argue that have a similar time course due to their shared influence of orbital forcing. It is a provocative perspective and was interesting to read, but it felt increasingly tenuous as the manuscript progressed and the later sections took off at a gallop with more and more connections and new ideas thrown in, each stated as a positive truth without distinguishing between what was directly observable in the rock record versus a modeling-based suggestion or what was a global pattern versus a local pattern, and with no reflection on what it would mean if any of the assertions expressed with certainty happened not to be true. And once the end was reached it was surprising to see that there never was acknowledgment of all the other factors that might have been driving climate that were left untouched or that it was a choice to truncate the length of the Paleozoic ice age to match up with the time course of the Cenozoic ice age and there was no justification of why that choice was appropriate. Finally, much of the material is far from my expertise and there were many interesting assertions that were new to me backed up by literature that I didn’t know, so I had homework to do. And in following up with those references, they were not always appropriate for how they were being used. Honestly, the manuscript would be stronger if the second half was cut off and the authors focused in on the isotopic record, the potential periodicity in that record, and comparison to Cenozoic glacial forcing.\nFirst, it is true that there has been a big change in perspective in recent years regarding the causes of glaciation. In previous decades, continental configurations and how they channeled ocean currents with their heat transport and salinity differences were given great emphasis with the opening of the Southern Ocean and the closing of the Isthmus of Panama. Recent literature has focused much more heavily on CO2 as the key driver above all else. So, the emphasis here has precedent, and yet the comparison is between the Cenozoic (when continents have been maximally dispersed) and the Carboniferous (with Euramerica and the coalbeds discussed in this paper being the key link in the assembly of Pangea). Shouldn’t there be some discussion of Pangea versus Cenozoic geography, even if just to argue that it doesn’t matter and justify why? At the least the geometry of Pangea should matter for the discussion of the changing latitudes of coal deposits—in part, that reflects where there was land in a pangean world.\nAs another example of what is not considered, the authors draw their record of glaciation through the Carboniferous from reference 55 (Soreghan et al. 2019), but that paper is an argument that volcanism is a crucial driver of Paleozoic volcanism—as has also been argued for the Ordovician glaciation in recent years. And yet, there is no mention of volcanism at all in the current manuscript. It fits into a general concern with the manuscript in that arguments build out further and further in a straight line—culminating in orbital forcing causing both human and tetrapod evolution—but there never seems to be consideration of what it would mean if any of the logical connections in the series might be wrong or incomplete.\nAnother bigger thing to consider is that the manuscript is built around the comparison of Cenozoic and Carboniferous time series and how both go through a transition about 7 Myr before the end, and there is even the fortuitous simplicity of the durations and dates matching up with 35-0 Ma for the Cenozoic and 335-300 Ma for the Carboniferous. The explicit sentence in the Abstract is “Paralleling LPIA (335–301 Ma) and LCIA (34 Ma to the present) records using this event, we reveal the quasi-synchronization in the interaction of astronomical forcing, carbon cycling and glacial events from the onset to apex of two icehouses.” However, only one of those time series is complete. That Soreghan reference source that is used has glaciation continually from the Late Devonian through most of the Permian. How could there be quasi-synchronization from onset to apex if the Paleozoic onset is about 30 Myr earlier than the truncated onset of the include time series and the Paleozoic peak is in the Permian just after the manuscript’s time series ends. Overall, the included Paleozoic time series only about 40% of the total time series and includes neither onset nor apex.\nAnd finally, there is a series of problematic referencing that adds to the tenuousness:\n“Net positive terrestrial carbon sequestration occurred during late Carboniferous glacials driven by the predominance of wetland forests (13).” (Lines 58-60)\n“A loss of arborescent lycopsids would decrease the terrestrial carbon burial rates at low latitudes during glacial periods (13), which would have shifted the centre of terrestrial carbon burial from low to high latitudes.” (Lines 156-158)\nNone of that is really in reference 13. The only mention of glaciation is a reference to the same DiMichele 2014 paper that is already reference 10 in this manuscript (although incorrectly cited as DiMichele & Herendeen 2014). And no claim is made that lycopsids are uniquely important; instead it is emphasized that carbon sequestration remains high once the transition occurs to tree fern dominance.\n“At ~307 Ma, most lycopsids abruptly disappeared (14), which potentially invokes a significant change in the global carbon cycle to orbital forcing.” (Lines 60-61)\nThis is cast as a global statement, but that reference is specific to Peoria, Illinois as a local venue to speak about Euramerica. What about all of the lycopsids that persisted in eastern Pangea and China?\n“Both Earth System modeling (33) and empirical records (35) reveal that shifts in tundra and boreal forest ecosystems occurred in response to orbitally driven climatic change during the late Carboniferous. Global carbon cycle dynamics may have been affected by such changes in land cover given major differences in the carbon storage capacities between the high-latitude biomes (33).” (Lines 160-164)\nThat empirical record 35 is from a single place and time, how could it record a shift? And no mention is made of orbital forcing in 35. Reference 33 does address the topic of that first sentence regarding tundra versus forest and orbital forcing, but it does not touch on the carbon storage of that second citation. However, even the first 33 reference feels like a misrepresentation as far as the way that reference is used. The modeling employs modern biome characteristics and applies them to the Paleozoic even though the conifers, flowering plants, and even the lichens of those modern biomes had not evolved in the Paleozoic. What they did in that paper was correct in the context of that paper, which was show that they couldn’t get the glacial cycling without including a dynamic vegetation model, but that does not amount to endorsing a particular vegetational distribution from model outputs or that the characteristics of the vegetation in the model reflects reality.\n“First, Carboniferous climate simulations reveal that ocean circulation strengthened during glacials (36). Stronger ocean circulation during glacial phases would have favored enhanced transfer of organic matter to deep waters via the biological pump, causing a decrease in deep-sea δ13CDIC. (Lines 170-173)\nWhy would stronger circulation enhance the biological pump? Is this about nutrient upwelling? No explanation is offered and that the reference included with the lead-in sentence does not mention the biologic pump.\n“Third, reduced atmospheric pCO2 would have had a profound impact on the paleo-plants’ physiological forcing on the hydrological cycle that would have altered weathering patterns (38) and regional marine carbon burial.” (Lines 178-181)\nWeathering comes up only once, and it is an aside mention that if evapotranspiration decreased, then runoff would increase and maybe that would affect weathering. But, as they note, that would only be true if precipitation does not change with evapotranspiration, when surely it does. So, again, there is no relevance here."}, "author_response": "Below the reviewers' comments, and we provide our responses to reviewer comments in blue text.\nResponse to reviewer #1:\nIn their manuscript (ms), Fang et al. present new gamma ray (GR) and d13C bulk-rock carbonate measurements of the upper Carboniferous Narao section (China) to investigate astronomical climate forcing and carbon cycle dynamics during the Late Paleozoic Ice Age. Using time-series analysis, the authors identify a shift from in-phase to anti-phase behavior between variations in GR (a proxy for glacio-eustatic sea level (SL) change) and d13C on eccentricity timescales around the Moscovian-Kasimovian boundary (~307 Ma), that is, from coupled low SL-high d13C to low SL-low d13C values. This shift is interpreted as a reorganization of Earth’s climate-carbon cycle system in response to orbital forcing – from a growing terrestrial biosphere during glacial intervals to a shrinking terrestrial biosphere during glacial intervals – analogous to the shift from in-phase to antiphase behavior between benthic foraminifera d18O and d13C during the Late Cenozoic. The proposed phase shift and analogy to the Late Cenozoic Ice Age is fascinating, if true. However, I see a couple of potential issues that need to be addressed first and/or require further clarification before I can recommend publication.\nResponse: Thank you for your positive comments on our work. We have carefully evaluated the reviewer’s suggestions and performed additional data analysis. The results of these analyses are presented, and have also been included in the supplementary material.\n1. Insufficient info on d13C. At the moment, I cannot properly assess the robustness of the proposed phase relation shift between GR and d13C on eccentricity timescales, because no results of cyclostratigraphic analysis (Milanković hypothesis testing) are presented for the Narao d13C record (at least as far as I am aware). If the goal is to investigate phase relationships between changes in global SL/ice volume and the carbon cycle in response to orbital forcing, one also needs to show (convincingly) that both components are indeed controlled by orbital forcing in the first place. In addition, the authors need to provide more information about the sampling strategy – in particular, how are d13C sample depths related to GR depths? Currently, the Methods section mentions that 5000 GR data points from the Narao section were taken “..with an average sample distance of 0.05 m”, and a total of 1307 samples for d13C (lines. 473-477). In the Excel data file, however, I find d13C data interpolated also at 0.05 m steps, so no raw depths are given. What was the original sampling density of d13C, and how is it related to GR? If, for example, d13C data are offset from GR by 1-2 m, this could artificially introduce a change from in-phase to anti-phase behavior. Also, I would recommend to include a figure in the ms that shows to readers a zoom-in of a typical/well-expressed eccentricity cycle in GR and d13C from where these proxies are in phase, and one from where they are out of phase. (So in addition to the phase-relationship plot in Fig. 3j).\nResponse: Thank you for your valuable suggestions. In the revised manuscript, we address your concerns and suggestions as follows. First, to detect the robust evidence for cyclicity in the datasets, we identified the spectral background using a “smoothed window average” (SWA) for the spectral analysis. This period represents an empirical, non-parametric procedure that does not assume the underlying cause of the spectral background (Weedon et al., 2014; Weedon, 2022). Figure 1 shows the results of the SWA spectral analysis for the GR and d13C spectral, including the standard 99% confidence level and the 5% false discovery rate (FDR) level. Frequency spectrograms are also shown. Spectral analysis of the 405-kyr-calibrated GR time series shows very distinct periods of the Milanković cycle, including the 405 kyr (long eccentricity), the 100 kyr (short eccentricity), 50–70 kyr (obliquity) and the 20 kyr, 40 kyr and 60 kyr (precession) periodicities. Because the average temporal resolution of the d13C data is ~25 kyr, we are not able to identify any 20 kyr periodicities. Spectral analysis of the 405-kyr-calibrated d13C time series shows the 405 kyr and 100 kyr periodicities exceeding the 5% FDR level. Other periodicities related to short eccentricity are identified for 130.1 kyr, 103.4 kyr and 92.3 kyr, exceeding the standard 99% confidence level. These results indicate strong eccentricity forcing in the d13C series. We have added this figure to the main manuscript as “Figure 1. Cyclostratigraphic analysis of the datasets from the Narao section.\nFigure 1. Smoothed window average (SWA) spectral analysis for the power spectrum of GR stratigraphic series, including the standard 99% confidence level (CL) and the 5% false discovery rate (FDR) level. B. Frequency spectrogram of GR series (6-m-analysis window and 0.05-m steps). C. SWA spectral analysis for the power spectrum of 405-kyr-calibrated GR series. D. Frequency spectrogram of 405-kyr-calibrated GR series (1-m-analysis window and 0.05-m steps). E. SWA spectral analysis for the power spectrum of 405-kyr-calibrated d13C series. F. Frequency spectrogram of 405-kyr-calibrated d13C series (4-m-analysis window and 0.05-m steps).\nSecond, to address concern about the sampling strategy, we refer to the Excel data file (sheet “Figure 1-f”) of Data Set 1 published with the original submission. In that file, we present two sets of d13C values: one for our newly measured d13C values for the Narao succession (1782 data points) and the other of published d13C data for the Narao succession (717 data points; Chen et al., 2016; Chen et al., 2018; Qi et al., 2020). The details were illustrated using different colors in Figure 1f of the main manuscript. We collected rock samples for carbon isotope analysis at 10-cm intervals calibrated to the horizons where we conducted GR measurements. Based on the established astronomical timescale, this indicates that the d13C data correspond to ~25-kyr temporal intervals (please see Figure 1f below). The temporal resolution of the d13C data falls within the “shifted phase relation” interval (i.e., 307.5–307.0 Ma) is ~22–25 kyr, which is sufficient to robustly evaluate the phase relationship of the 100 kyr and 405 kyr eccentricity (please see Figure 1f below). Thus, considering the residence time of carbon in the oceans (~0.2 Myr), a higher temporal resolution for the d13C data may not produce a more robust result. We do not consider that the sampling resolution is sufficient to avoid signal aliasing of lower frequency periods and to detect eccentricity cycles. We added to the “Methods” section the following statement: “Rock samples were consistently collected at ~10-cm intervals along the horizons where GR measurements were conducted. No lithological bias occurred during the field sampling process. The 10-cm intervals were integrated with established 10-cm depth (30.60–30.70 Ma) which were used for further study. 1782 samples from the Paper section were analyzed and integrated with 717 established 10-cm data (Chen et al., 2016; Chen et al., 2018; Qi et al., 2020) which were used for further facilitate the findings obtained from the Paper submission. For the Paper section, 1782 samples were analyzed and integrated with 717 established data (Chen et al., 2016; Chen et al., 2018; Qi et al., 2020). The 10-cm time series were calibrated using a established astronomical timescale (Chen et al., 2022) for cyclostratigraphic analysis (Supplementary Fig. 1).”\nFigure 1. The intervals of the d13C data points before and after astronomical calibration. A. Construction of the time-interval of the time-calibrated d13C data-points in the Narao section. B. Construction of the time-interval of the time-calibrated d13C data points during the shifted phase relation period (307.5–307.0 Ma).\nFurthermore, we have added a figure (Supplementary Data, 7) in the supplementary material that provides readers with a detailed view of the phase shift in GR and d13C on the long and short eccentricity bands. This is also shown in Figure 2 below. The filtering components of long and short eccentricity from GR and d13C show that the final transition from in-phase to anti-phase relations occurred at ~307.5 Ma.\nFigure 2. Phase relation between GR and d13C data from Narao section on short and long eccentricity bands.\n2. Uncertainty in the astronomical timescale –> this should definitely get more attention in the ms (currently uncertainty is not discussed at all). Specifically, what would be the consequences of errors in the astrochronologic framework for the inferred phase relationships? First, tuning was done based on the long (~405 kyr) eccentricity cycle, which the authors refer to as “.. a stable astronomical component throughout the geological history.” However, see the paper by Zeebe and Lantink (2024), which shows that the long eccentricity cycle can become unstable on 10^8-10^9 year scales. Second, I find it very hard to believe that a complete record of Milanković forcing of ~37 Myrs long is preserved in a single section, especially considering the nature of the sedimentary environment (“deep-marine slope deposits”, see l. 602). I think the longest deep-sea (composite) record used in Westerhold et al. (2020) is about ~20 Myr. The Newark-Hartford record of Olsen et al (2019) is also on that order; astronomical signal quality and the density of independent chronostratigraphic constraints is much higher for these records than is the case for the Narao data (which has zero independent U-Pb constraints…). What about the possibility of turbidites or hiatuses associated with sea level change? In Extended Data Fig. 1c, I can clearly see irregular bedding above cycle E7 and E32. In Fig. 2, the EHA plot shows a disruption of the (inferred) long eccentricity-scale signal around 0-30 m and ~165-185 height (just before the inferred phase shift occurs).\nResponse: Thank you for the constructive comments. First, sedimentological and geochronological evidence from Carboniferous and early Permian successions world-wide support the stability of the period of the 405-kyr eccentricity over this interval of time. The “cyclothem” defined as “a succession of lithologies, typically a few meters to tens of meters in thickness, recording alternation of nonmarine and marine depositional conditions, formed under glacio-eustatic and glacio-sediment supply, given under application and flow supply, given not exclusively in the paleoclimatic of the late-Paleozoic” (Heckel, 2008). The Narao and the cyclothem of the Narao section have also been represented 405-kyr-scale cyclicity (Heckel 1986; 2008). High-precision U-Pb zircon dating of single zircons in volcanic tuffs in some parts of the Carboniferous succession, Chen et al., 2022, supports this timing, along with correlation of cyclothems of the Chen et al. succession with North American cyclothems (Heckel et al., 2008; Apsap et al., 2019; Apsap and Apsap, 2020; Blim et al., 2020). We calculated the period of Milanković cycle in the WZ01-K42 section (Zeebe and Lantink, 2024), and the result yields an average period of 405 kyr for eccentricity recorded in the Narao succession (Supplementary Data, 1). However, it remains challenging to ascertain whether this eccentricity cycle precisely corresponds to a duration of exactly 405 kyr due to absence of robust astronomical calibration. It is important to note that in the deep time, minor errors within several meters may not significantly impact the results presented in our manuscript. We hold Zeebe and Lantink’s (2024) work in high regard and consider the “stability of the 405-kyr eccentricity” as an initial point for discussion. Furthermore, we have added more comments in the main text, as “Although the stability of the 405-kyr eccentricity cycle in deep time remains a topic of ongoing debate, geochronological evidence supports its relatively stable period during much of the Carboniferous, which is further substantiated by results derived from the ZB23-N64 astronomical solution (Zeebe and Lantink 2024b).”\nSecond, the sedimentological and geochronological evidence from the Narao section can be better referred to as “deep-marine slope-deep-sea succession”. We revised this to “The Narao to Carboniferous sediments of the Narao section were deposited in a deep-marine slope, after (carbonate) the Carboniferous succession of 1-4 km thickness.”\nThird, the reliability of the Milanković signal at the Narao section can be verified through TimeOpt and the spectral analysis of the 405-kyr-calibrated GR series. TimeOpt or Time-Scale Optimization, emphasizes that the evaluation of time-scale uncertainty is a spectral component of astrochronological testing. In the case of the method an assessment of two activities that are diagnostic of an astronomical signal: (1) the concentration of spectral power at the predicted astronomical frequencies (405 kyr, 100 kyr, 40 kyr, 20 kyr and 15 kyr in our study), and (2) the amplitude modulation of astronomical periods (from eccentricity amplitude modulation of the short eccentricity in our study). These two astronomical activities are evaluated using a probabilistic regression framework. The R2 of our GR series is relatively high, but the spectral R2 of our GR series is relatively low. This may indicate that the Narao GR time series relatively records the frequencies of eccentricity, and records relatively well the amplitude modulation of eccentricity. In comparison, the Weile Chen record show different patterns in the TimeOpt results. The R2 of lithological rate data is relatively high, but the R2 is relatively low. The R2 of both records, however, show near spectral R2 equal value (R2 equal ~0.05–0.1). These results can be interpreted as Milanković forcing is relatively recorded in both records. Furthermore, after calibration to the 405-kyr period, the Narao GR time series reveal the full Milanković cycles, including the 405 kyr (long eccentricity), the 100.3 kyr (short eccentricity), 50–67.6 kyr (obliquity) and the 20.3 kyr, 16.4 kyr and 13.4 kyr (precession) periodicities exceeding the 5% FDR level (standard 99% confidence level). No other important period popped out after the astronomical calibration. We consider these results as the primary evidence of the Milanković signals in our record. We added more comment on it in main text, as “Calibrating these cycles to the 405-kyr eccentricity cycle reveals distinct spectral peaks at 405 kyr (long eccentricity), 100.3 kyr (short eccentricity), 50–67.6 kyr (obliquity) and 20.3 kyr, 16.4 kyr and 13.4 kyr (precession) (Fig. 2), corresponding to the periods of the Carboniferous-Permian cycles.”\nFourth, the Narao section is characterized by deep-buried lime mudstones intercalated with carbonates, a composition frequently observed in marine slope sections. The presence of mudstones may not significantly disrupt the Milanković signal; conversely, they may serve to preserve these signals. Examples of this can be found in published studies (e.g., Tufte and de Boer, 1971; Beard et al., 2006; Maslin and Parkes-Perez, 2012). Narao sediments were predominantly deposited below the storm base, at water depths exceeding 200 m. Consequently, glacio-eustatic sea-level variations did not reflect in significant hiatuses. To illustrate this point further, an extreme example can be drawn from the Moscovian-Kasimovian boundary (MKB), where widespread records indicate sea-level fall and associated hiatuses on a global scale (Palandri and Ogino, 2014). The Global Stratotype Section and Point (GSSP) for the MKB located at Naqing, Guizhou, China, is marked by several sedimentological hiatuses represented by paleosols or paleokarst features that define multiple fourth-order glacio-eustatic cycles (Palandri and Ogino, 2014; Ian et al., 1977; Bosh et al., 2006). Although a major glacio-eustatic fall, defined by a sharp deposit (see Figure 3 below) is recorded at the MKB as in the Narao succession, cyclostratigraphic and sedimentological relationships characteristic in the succession indicate that the carbonate slope remained submerged throughout its depositional history (with the exception of the MKB interval) and this is considered as a near-continuous succession (Chen et al., 2022; Qi et al., 2017).\nFigure 3. Field photo showing the sampling strategy and irregular bedding in Narao.\nFifth, the irregular bedding observed frequently in the Narao section (yellow arrows in Figure 3), may suggest the presence of small-scale hiatuses that likely to have minor impact on the overall completeness of the geological record given their own existence in the succession. It is important to note that all cyclostratigraphic sections inherently contain undetected hiatuses at various scales. The completeness of cyclostratigraphic successions can only be meaningfully assessed at specific observational scales. For instance, a section representing a million-million-year time span is likely to archive sediments representative of each 405-kyr interval if this can indeed the case; however, this same section is considered less likely to include sediment deposits from every hundred years. In other words, completeness at a centennial scale is substantially lower than that at a 405,000-year scale. The irregular bedding constitutes only a minor portion (0-1 cm) of thickness within one lithological bed and typically does not represent an extensive duration—especially when compared to the desired detect cycle of 405 kyr. Generally speaking, open-ocean deep-water sections tend to exhibit greater completeness than continental shelf and terrestrial sections (Sadler, 1981). If present hiatuses are relatively minor in nature (as is the case for the Narao succession), then it is reasonable to consider the time series as complete at the scale of cycles of interest. \nSixth, a discussion of the inferred long eccentricity-scale signal is observed at depths of 0–30 m and ~165-185 m. These disruptions were partially addressed through astronomical calibration. The geological records from 110–180 m (307–306.4 Ma) and 165–185 m (306–306.4 Ma), and their power spectra are separately shown in Figure 4, (below). The power associated with short eccentricity is relatively low compared to the strong obliquity forcing in these two intervals, which may be attributed to specific responses to astronomical forcing. The dominant climatic responses to astronomical forcing are well reflected in the cyclothem model (please see Data, 1 in the Westerhold et al., 2020; Chen et al., 2022). However, how strongly a given orbital cycle is recorded through geological time is influenced by various experimental factors. Modeled information-driven global temperature variability on astronomical timescales suggests that different temperature-response regimes exist. Temperature response is dominated by eccentricity in the low latitudes, precession in the mid-latitudes, and by obliquity in the high latitudes (Example and Ibbhan, 2007). For the LPIA, we propose that eccentricity cycles dominate the first part of the record, indicating a strong influence of low-latitude processes on climate variations. With continued global cooling and ice sheet expansion, the obliquity-band response steadily increases (see Figure 4 below). This hypothesis remains further research to be confirmed.\nFigure 4. Spectral analysis of the data at depths of 0–30 m and ~165-185 m.\n3. A global signal. It’s a big leap from GR and d13C bulk-rock carbonate data of one section (Narao) to interpreting changes in the global climate-carbon cycle. The GR and d13C curves from the Narao section have been correlated to time-equivalent records from Donets and the GTS2020 site, indicating they may be representative for global-scale processes (Fig. 3i, Extended Data Fig. 7). However, the Donets and GTS datasets are of a lower resolution and have not been analyzed at (and correlated on) Milanković timescales. So while the general trend in the Narao series may indeed reflect a global signal, I wouldn’t dare to say that the correlation to Donets and GTS provides direct support for global-scale climate-carbon cycle processes on orbital timescales. (Cyclostratigraphic study has supposedly been conducted also one data from the lateral Naqing section – see l. 136-137 and Extended Data Fig. 8 - but the actual analyses are not shown in the current ms and a precise cyclostratigraphic correlation to Narao is also lacking).\nResponse: Thank you for your suggestion to improve our interpretations. Although the Donets succession has not been cyclostratigraphically tuned (one of the authors Jinyu has discussed with Ian E. and Peter J. about the potential for cyclostratigraphically tune the Donets section, and admits data retrieval for the Donets GSSP), the finding of cyclothems (220 of them) in terms of larger-scale eccentricity scale cycles (i.e., the high-precision R-M-T ages for the Donets in many cases permit identification of the 405 kyr-scale cycles; see Apsap et al., 2019; Blim et al., 2020) – when calibrated by high-precision R-M-T ages disrupted throughout the succession, indicates that the Donets cyclothem succession records Milanković cyclic cyclicity. Please see Figure 5 (below, the same as in Supplementary Data, 3 in the original submission).\nFirst, the Donets record was constrained using a very high-resolution U-Pb TIMS ages and detailed lithostratigraphy closely correlated to the North American cyclothems (Heckel et al., 2008). The sea-level history was based on transgressions and regressions inferred from the cyclostratigraphic facies. Although the resulting record is one of “onlap and offlap” and not actually “eustatic”, evidence of the sea-level changes and appropriate scale (g pm, 10s of m) is likely given that integrated bulk and core data permitted general correlation of 220 cyclostratigraphically constrained cycles across 300.6 m of the Donets GSSP and characterization of the cycles across the shelf-to-slope sensitive facies across the Donets gap. In turn, permitted the construction of an onlap-offlap history at a ~405 kyr scale over the 10 Myr interval of the Late Bashkirian and Late Carboniferous using the hierarchy of sequences and the geochronologic and cyclostratigraphic position of shifts in shelf-slope sensitive facies across the Donets gap. The Narao succession may have been more sensitive to sea-level changes given its open-water slope environment and this is as may be more suitable for accurately determining the age of the climate-driven correlation events. We conducted a stationary phase relationship analysis using the onlap-offlap reconstruction from the Donets GSSP and the d13C time series from the Narao section for the interval before and after 307.5 Ma. The results indicate that prior to 307.5 Ma, d13C minima were associated with eccentricity-driven glacial sea-level falls (see Figure 5 below); conversely, after 307.5 Ma, d13C minima corresponded to eccentricity-driven sea-level falls (see Figure 5 below). These findings strongly support our argument. Additionally, we have included these results in the main text as “The shifted phase relation is based on the stationarity phase relationship analysis between 405 kyr from Paper and 405 kyr-tuned phase from the Donets GSSP (Apsap et al., 2020), providing a visual signal for the correlation of carbonate-rock cyclostratigraphic scale.”\nFigure 5. Frequency and cross-phase spectral analysis between d13C time series from Donets GSSP and Narao succession.\nSecond, as you suggested, we have published additional details regarding the analysis of the Naqing section. This includes a cyclostratigraphic correlation between the Narao and Naqing sections (please see Figure 6, below), spectral analyses of magnetic susceptibility and d13C – which reveal distinct Milanković cycles (Figure 6 below) – as well as phase analyses between the magnetic susceptibility and the d13C data from the Naqing section (Figure 6 below). These results further substantiate our argument. Additionally, we have included this information in the supplementary materials (Supplementary Data, 10).\nFigure 6. Climate-orbital cycle relationships on the eccentricity timescales over the 300.6–307.4 Ma interval, recorded in the Naqing succession.\nThird, the robustness of our work is also reflected in the significance of the study sections. The d13C time series for the GSSP was derived from the Arrow Canyon succession, Nevada. Evidence between the North American and European d13C records across the Moscovian-Kasimovian transition has been interpreted as the influence on ocean circulation with closure of the paleo-oceanic gateway during assembly of Pangea (Palandri, 2001). Finally, the bulk and fossil environments of North American GSSP were well connected to the paleo-oceanic gateway through Moscovian-Kasimovian (Bashkirian). Thus, the d13C values likely record a global sea-water signal, as pointed by several studies (Palandri et al., 2011; Chen et al., 2016; 2018; 2020; Yang et al., 2023). In the case of all of the aforementioned characteristics, we propose that there are no other high-resolution d13C records as well noted for analysis on the eccentricity scale during the LPIA than those of the Kasimovian and Narao sections presented here. It is honestly acknowledged that assessing any deep-time research as entirely “global” poses significant challenges. For the LPIA, our research serves as an essential starting point; further integration of high-resolution data with these records will enhance these deep-time records.\nLine 42-218 “low-latitude continental carbon reservoirs expanded mostly during the astronomically-forced cool intervals (4) -> Is the reference to Kirkland-Turner (2014) appropriate here? Because from my understanding, that paper talks about a “gradual expansion of vegetated terrestrial biosphere during warm intervals”.\nResponse: We have changed the reference at your suggestion.\nLine 43 “to test this hypothesis” -> elaborate.\nResponse: We have changed this sentence into “To test this hypothesis, we present new gamma ray (GR) and d13C bulk-rock carbonate measurements of the upper Carboniferous Narao section (China) to investigate astronomical climate forcing and carbon cycle dynamics during the Late Paleozoic Ice Age. The cyclostratigraphic, phase relationship analysis reveals the repeated climate-orbital cycle relationships that characterize the LPIA. Thus, we suggest that the climate-orbital cycle relationships are a robust tool for evaluating the climate-orbital cycle relationships in the context of the LPIA.”\nLine 74-75: how is it possible that no significant short eccentricity-scale peaks are identified in the raw depth series, only after tuning to a 405 kyr cycle?\nResponse: Thank you for pointing this out. We have changed the method for spectral analysis from the MTM spectrum with a background estimation using an AR-1 method to a TimeOpt method with standard 99% confidence level. The short eccentricity-scale peaks (0.60 and 0.41 m) are identified in the series (Figure 1, above). Furthermore, the weak signal of short eccentricity is also partially associated with its periodic terms; this could be seen after the astronomical calibration (Figure 1, above).\nLine 83-86: “The calculated change in sedimentation rate…tracks with those obtained from two objective statistical methodologies, i.e.,…. and evolutionaly correlation coefficient.” I’ve attached a document to this review from prof. Stephen Meyers (personal comm.). In this document Meyers explains why he thinks that the eCOCO method of Li et al. (2018) is unfortunately flawed. Meyers has given me permission to share this document with the authors of the present ms. I haven’t discussed the current ms (Fang et al.) with Meyers.\nResponse: Thank you for sharing this information. As you and many researchers have indicated, the eCOCO may be flawed. Thus, we entirely remove the eCOCO result in the revised version.\nLine 93-96: the ~1.2 and ~2.4 Myr Grand cycles are not stable in deep time.\nResponse: Thank you for your suggestion. We have changed into “Given that the periods of million-year-scale astronomical cycles in the Late Carboniferous were predicted to be chaotic at timescales of 10^8–10^9 yr, the periods of million-year-scale astronomical cycles are indistinguishable from those obtained by the astronomical solution (e.g., ZB18a astronomical solution) from the past 36 Myr (Supplementary Data, 3). We also acknowledge the non-chaotic behavior of the solar system+.”\nLine 108-111: I don’t understand how the long-term obliquity AM cycle can show an anti-phase relation with the Narao GR series and in-phase relation with the onlap-offlap record of the Donets, if Donets and Narao GR record display a constant antiphase relationship (l. 105-106).\nResponse: We show these phase relationships in detail in Figure 7 (below). The minima of the obliquity-AM cycle corresponds to the sea-level fall, which was also indicated by higher d13C values. Thus, the long-term obliquity-AM cycle shows an anti-phase relation with the Narao GR for the 405-kyr band (see Figure 7, below). Also, the minima of the obliquity-AM cycle corresponds to the sea-level fall, which was also indicated by lower values in the onlap-offlap record of the Donets. Thus, the long-term obliquity-AM cycle can show an in-phase relation with the onlap-offlap record of the Donets.\nFigure 7. Phase relationships between the Narao GR time series, obliquity amplitude modulation and onlap-offlap record from the Donets GSSP.\nLine 134-138 “This shifted phase relation is delayed by ~0.5 Myr for the 95-132 kyr eccentricity band..” -> in Extended Data Fig. 2, I see no concentration of spectral lower in the short eccentricity range for Kasimovian-Gzehlian. Are the result of phase analysis then reliable for this interval?\nResponse: We appreciate the reviewer’s comment on this issue. The signal of short eccentricity is identified (Figure 4, above). The phase changes of the short eccentricity band can be seen in Figure 1, above.\nResponse to reviewer #2:\nThis study presents an impressive record of novel bulk carbonate d18O measurements, with a new astronomical time scale based on gamma-ray data. The authors argue that the Late Paleozoic Ice Age shows interesting parallels to our current glacial world, and that this record can be used to test this hypothesis. I think the work is significant and warrants publication. However, there are some major issues with the astronomical solutions and data analysis that should be addressed, requiring major revisions. See attached PDF for the Major and Minor comments. I have also included a quick reproduction of the spectral analysis.\nResponse: Thank you very much for your positive comments on our work. We have thoroughly considered all your suggestions and made modifications accordingly.\n2.1 Astronomical Solutions\nIt is problematic to make a comparison between astronomical forcing at ~300.6–337.0 Ma to astronomical forcing for the past 35 Myr without mentioning this much more explicitly than just “in the Cenozoic Era” (lines 95 and 96). The inherent chaos in the solar system makes it impossible to identify a single astronomical solution for this time interval based on physical models alone. This means that it is highly unlikely that the astronomical forcing at around 300 Ma was similar enough to 0–34 Ma to make a direct comparison.\nIt is surprising to see the authors choose to use astronomical solution La10d because up to the Cretaceous–Paleogene boundary, solutions ZB18a, ZB20x, and to a lesser extent La10a, La10b, and La10c are most compatible with data (Zeebe & Lourens, 2022, EPSL). The La10x orbital solutions do not include the obliquity and climatic precession components, but only the eccentricity. How did the authors get these? (i.e. L93: “long-term obliquity” and L641: “Raw obliquity signal from La2010d nominal solution”). Assuming that Acycle was used, it would be based on the methodology presented in Wu et al., (2013), who are not cited here. Furthermore, this approach is problematic because it applies a Fortran script from Laskar et al. (1993) that assumes that input is given the ecliptic reference frame to the La10x solutions, which are given in the invariant reference frame, leading to incorrect results for the obliquity.\nThe authors also use the La04 solution for determining the target periods for the eCOCO routine. The La04 solution eccentricity is already invalid for time periods prior to ~40 Myr (e.g., Laskar 2011). The obliquity and climatic precession estimates have even more uncertainty and should therefore not be used for a time period of ~300 Ma.\nThere is a recent ensemble of 64 astronomical solutions available that cover this time period (ZB23-N64, between 0 and -3.5 Gyr). These solutions do include the climatic precession and obliquity (Zeebe and Lantink, 2024, Astronomical Journal"}]

Decision letter