Precession-induced millennial climate cycles in greenhouse Cretaceous

Domain: earth_science
Published: 2025-11-27
Editorial outcome: minor_revision

Abstract
Millennial-scale climate cycles persist throughout greenhouse climates, yet the mechanisms remain unclear. Here, using proxy reconstructions, we present centennial-resolution geological records from early Campanian greenhouse deposits in both mid-latitude East Asia and low-latitude Southern Atlantic. These records document pronounced millennial (~1–6 kyr) wet-dry climate cycles. The amplitude modulation relationship between the most prominent ~4–5-kyr cycles (corresponding to the 1/4 precession cycle) and eccentricity aligns perfectly with the theoretically calculated equatorial insolation cycles, demonstrating the climate effect of this predicted insolation forcing on global climate. Other millennial cycles primarily emerge from these ~4–5-kyr cycles via nonlinear amplitude modulation and combination tones. Proxy reconstructions and theoretical calculation thus converge to demonstrate that during this warm greenhouse period, precession can directly and indirectly stimulate millennial climate cycles. The deterministic link between astronomical parameters and millennial climate cycles implies that high-frequency climate oscillations may be predictable in future greenhouse-like climates, particularly under anthropogenic warming.

Introduction
Earth’s climate has oscillated across multiple timescales1, from multimillion-year icehouse-greenhouse transitions to orbital-scale glacial-interglacial cycles. Superimposed on these gradual changes are abrupt millennial-scale climate cycles (1–10 kyr), demonstrating inherent climate instability. The most well-documented examples of millennial climate variabilities are the Pleistocene Dansgaard-Oeschger2 and Heinrich events3, manifested as abrupt temperature increases in Greenland and episodic iceberg discharges in the North Atlantic. These events were described within the context of the Quaternary icehouse world, especially during glacials, and they are thus mostly linked to cryospheric dynamics such as ice-sheet oscillations and concomitant ocean circulation shifts4. In the last decades, emerging evidence reveals that high-frequency climate oscillations persisted even during ancient greenhouse worlds5–9 (Fig. 1). In greenhouse climates, where polar ice sheets are minimal or absent, the generation of millennial climate cycles is hypothesized to arise from different mechanisms, including millennial shifts in monsoon intensity6,10 and oceanic redox conditions9, and nonlinear response to orbital forcing5,8. A key unresolved question is whether astronomical forcing can directly excite millennial climate cycles under such warm ice-free backgrounds. Theoretical models posit that the climate response to the largest amplitude of the seasonal cycle of insolation at the equator due to precession (equatorial insolation) could trigger millennial climate cycles, especially at ~5 kyr11. Yet empirical validation has been hindered by the paucity of well-dated, high-resolution paleoclimate datasets from deep-past greenhouse stages. Quaternary archives, while temporally precise, conflate millennial cycles caused by this theoretical insolation forcing with those triggered by ice dynamics (e.g., Heinrich events)3,11.

The early Late Cretaceous exemplifies a greenhouse climate12,13, offering an ideal candidate to study climate dynamics in warm, high-CO2 regimes. The Songliao Basin (SLB) in East Asia presents a unique opportunity for high-resolution analysis, with nearly the entire Late-Cretaceous terrestrial strata recovered during the International Continental Scientific Drilling Project14. Here, we target the early Campanian (83.428–83.538 Ma) deep-lacustrine rhythmic deposits of the lower Nenjiang Formation (Kn1+2; 1137.84–1145.2 m interval) from the SK2 borehole of the basin (46°14′26.89″N, 125°21′47.03″E) (Figs. 1, S1). The rhythmic deposition is primarily characterized by alternating black siliceous mudstone/shale and gray laminated calcareous mudstone (Fig. S2). This stratigraphic interval was selected because the strata were deposited in a geological setting (deep-lacustrine environment) where autogenic noise in the ‘source-to-sink system’ is expected to be minimal15, and its finer-grained mudstones better preserve allogenic signals than coarser deposits16. To quantify this rhythm, we present centennial-resolution (~1-cm stratigraphic spacing; ~0.14-kyr temporal spacing) grayscale and X-ray fluorescence (XRF)-derived element ratios, Rb/Sr and log (Ca/Ti) (Fig. 2, S1). Grayscale effectively captures lithological variation, with higher values indicating greater carbonate content and lower values reflecting enrichment in organic-rich siliceous sediments (Figs. S1, S2). Rb/Sr tracks chemical weathering intensity, since Rb is relatively immobile, whereas Sr is readily leached during weathering17; higher Rb/Sr values thus reflect stronger weathering18,19. Ca content primarily reflects carbonate deposition (Fig. S2), which tends to be chemically precipitated in shallow lakes under arid conditions, while Ti represents terrigenous clay input to deeper lakes and is resistant to diagenesis20. Accordingly, the Ca/Ti is a proxy for lacustrine salinity and lake levels, with higher values signifying a more arid climate and lower lake levels19–21.

Moreover, previous studies of coeval low-latitude South Atlantic marine successions recovered from the DSDP–516F borehole (30°16′35.4″S, 35°17′6″W; paleolatitude ~25°S) during Leg 72 (Fig. 1) revealed semi-precessional and millennial climate cycles7,22. However, these analyses were conducted on the stratigraphic domain with poor chronological constraints. Thus, for a comparable analysis from a marine setting, a high-quality and high-resolution color reflectance (L*) (~1.5-cm stratigraphic spacing; ~0.63-kyr temporal spacing) from this core (80.702–81.472 Ma)7 was reinvestigated after refining its chronology using magnetostratigraphy23 and cyclostratigraphy (Fig. 2, Figs. S3, S4; see “Methods”). By synthesizing these terrestrial-marine datasets, we aim to identify significant millennial climate cycles in different deposition settings of the greenhouse period, investigate their mechanisms, and explore their climate implications.

Results
Millennial cycles in Cretaceous East Asia
In the depth domain, the grayscale, Rb/Sr and log (Ca/Ti) ratios in SK2 core all display prominent cyclicities with ~0.12–0.52-m cycles. On this timescale, three datasets are phase-coupled with larger grayscale correlating with lower Rb/Sr and higher log (Ca/Ti) (Fig. S1D). An anchored astronomical timescale (ATS) for the Kn1+2 interval of the core, derived from radiometric zircon dating of four bentonite layers and cyclostratigraphic intercorrelations across three boreholes in the SLB24, constrains a sediment accumulation rate (SAR) of ~6.36 cm/kyr (Fig. S1). With this SAR constraint, the observed cycles correspond to periodicities of ~1.88–8.13 kyr. Using the anchored ATS24, we converted stratigraphic datasets to the time domain (83.428–83.538 Ma) (see Chronology in the “Methods”), revealing distinct cyclic patterns (Fig. 2A–C).

Spectral analysis indicates that both ~100-kyr and ~22-kyr power peaks are prominent in three datasets. The ~100-kyr power peaks in grayscale and the ~22-kyr peaks in log (Ca/Ti) display low confidence levels; nevertheless, both signals are significant in the Rb/Sr datasets (exceeding the 95% confidence level in ML96 robust AR1 model) (Figs. 3A, B, S5). Despite the total length of the datasets being ~110 kyr, the trend of ~100-kyr cycle remains evident (Figs. 2B, C, S6). Besides, the ~100-kyr and ~20-kyr filters in three datasets are phase-coupled (Fig. S6).

Higher-frequency millennial signals are also prominent in grayscale, with strong spectral peaks at ~5.8 kyr, ~4.5 kyr, ~3.7 kyr, moderate peaks at ~2.6 kyr, ~2.3 kyr, ~1.8 kyr, and weak signals at ~1.2–1.1 kyr (Fig. 3A). Similar millennial cycles occur in log (Ca/Ti), with strong peaks at ~5.7 kyr, ~4.4 kyr, ~3.6 kyr, and moderate peaks at ~2.5 kyr, ~2.0 kyr, ~1.8 kyr, ~1.7 kyr, and weak peaks at ~1.3–1.0 kyr (Fig. 3B). Rb/Sr retain significant millennial cycles only at ~5.8 kyr, ~4.5 kyr, ~3.7 kyr (Fig. S5). All of these peaks surpass 95% confidence levels in both conventional and ML96 robust AR1 models. Moreover, noise model tests using the False Discovery Rate (FDR) approach25,26 reveal test p-values below 0.05 for these peaks (see “Methods”) (Fig. 3A, B; Fig. S5), confirming their statistical significance.

Millennial cycles in Cretaceous South Atlantic
Magnetostratigraphic analysis indicates that the interval of DSDP–516F core between 1125.68 and 1198.44 m corresponds to magnetochron C33r23. This depth-age model yields an average SAR of 2.43 cm/kyr based on the geomagnetic polarity timescale 2020 (GPST 2020) chronology27 (Fig. S3). Under the first-order constraint of this SAR, a floating ATS was established for the L* dataset from 1163.42 to 1145.132 m of the core7 with the TimeOpt algorithm28 (see “Methods”). The acquired floating ATS returns an optimal SAR of 2.37 cm/kyr (Fig. S4) and was anchored at 80.702 Ma, the linearly-interpolated magnetochronology at 1145.132 m. The anchored ATS converts the stratigraphic L* into a continuous time interval of 81.473–80.702 Ma (see Chronology in the Methods for more details).

Spectral analysis of the tuned L* reveals significant power peaks at ~405 kyr, ~107 kyr, ~22.7 kyr, ~18.2 kyr, ~11.3 kyr, ~9.2 kyr, ~5.5 kyr, ~4.8 kyr, and ~4.1 kyr (Figs. 2D–F, 3C). All of these peaks exceed the 90% confidence level in the conventional AR1 model and the 95% confidence level in the ML96 robust AR1 model. Moreover, moderate power peaks at ~1.4 kyr and ~1.3 kyr are also significant (above 95% confidence level of both conventional and ML96 model, the multiple test p-values smaller than 0.1). The test p-values for the power peaks at ~405 kyr, ~22.7 kyr, ~18.2 kyr, ~4.1 kyr are below 0.1 (Fig. 3C). Although the p-values for the peaks at ~11.3 kyr, 9.2 kyr, 5.5 kyr, 4.8 kyr are above 0.1, the tuned L* data clearly show a semipartite and quadripartite structure of precession-related cycles, with each precession cycle comprising two semi-precessional and four quarter-precessional cycles (Fig. 2F). Moreover, the evolutionary correlation coefficients between L* and its ~10-kyr filter (0.095 ± 0.025 cycles/kyr) and between L* and its ~4–5-kyr filter (0.215 ± 0.055 cycles/kyr) display significant eccentricity signals that are phase-locked with the eccentricity in L* datasets (Fig. S7).

Methods
Materials from the East Asia (SK2 borehole)
The Songliao Basin (SLB) is one of the largest Mesozoic terrestrial inland basins, with its paleogeographical location being located at ~45°N during the Cretaceous14. Under the scheme of the International Continental Scientific Drilling Project of Cretaceous Songliao Basin, the SK2 borehole (125°21′47.03″E, 46°14′26.89″N) was drilled in the basin to recover the successions, including early Mesozoic to Paleozoic basement, Early Cretaceous strata, and part of Late Cretaceous strata. Of these strata, the Santonian-Campanian lower Nenjiang Formation, including the first and second Member of the Nenjiang Formation (Kn1+2) (correlating to the interval of 1249.32–1145.23 m and 1145.23–1086 m in the SK2 core, respectively), was a set of deep-lacustrine successions, characterized by black mudstones with intercalations of thin carbonate layers, black shales, oil shales, fine sandstones, and bentonites54 (Fig. S1A). These successions were proven to preserve the paleoclimate information, including the contemporary humid/arid climate perturbations and astronomical signals19,24. Especially, the lowest part of the second member of Nenjiang Formation (Kn2) (1145.2–1137.84 m), deposited in the early Campanian24,55 when polar areas were devoid of major glaciation13, consists largely of alternating gray calcareous and black siliceous mudstones (Fig. S1C; Fig. S2). Thus, here we selected this interval as researching materials.

Materials from the South Atlantic (DSDP–516F borehole)
Borehole DSDP–516F, drilled on the Rio Grande Rise in the South Atlantic (30°16.59′S, 35°17.1′W) during Leg 72, was designed to acquire the Neogene-Late Cretaceous successions with the purpose of decoding the coeval paleoecologic characteristics of the southwestern Atlantic56. The Rio Grande Rise was located at ~25°S during 84 Ma57. Part of the early-Campanian interval of DSDP–516F has been subjected to high-resolution climate analysis. Park et al.22. recognized semi-precession cycles in the interval between 1145 m and 1166 m using magnetic susceptibility (MS) datasets. Later, de Winter et al.7. provided high-resolution stratigraphic color reflection proxy (L*) for the interval between 1145.132 m and 1163.42 m, and unveiled millennial-scale cycles. However, their analyses are all implemented on the depth domain, which may induce a large bias when interpreting these cycles, especially for the millennial-scale cycles. This high-resolution L* was reevaluated here using the TimeOpt method28, and periodogram spectral approach, as well as bicoherence spectra58,59.

Chronology
The age model of the Kn1+2 in the SK2 core was established based on the U-Pb chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS) zircon dating and cyclostratigraphic analysis24. The zircon U-Pb CA-ID-TIMS dating of four bentonite layers in the depositions of Kn1+2 of SK2 core (at depths of 1139.39 m, 1142.49 m, 1218.61 m, and 1236.81 m) yielded ages of 83.35 ± 0.11 Ma, 83.498 ± 0.052 Ma, 84.45 ± 0.14 Ma, and 84.709 ± 0.064 Ma, respectively (Fig. S1). Based on this depth-age model, the average sedimentary accumulation rate (SAR) for the entire Kn1+2 in the SK2 core is ~7.45 cm/kyr. Under this SAR constraint, long-eccentricity signals were recognized and filtered from the Gamma Ray (GR) logging datasets of SK2. Because long-eccentricity cycles have remained stable over the past 250 Ma30, the filtered and interpreted long eccentricity signals in GR were tied to the 405-kyr periodicity to establish a floating astronomical timescale (ATS) for the Kn1+2 in the SK2 core24. This floating ATS was then anchored to the U-Pb age at 1142.49 m, yielding a SAR of ~6.36 cm/kyr for the studied interval of this study (1137.84–1145.2 m). This tuned SAR agrees with Monte Carlo numerical simulations based on GR datasets from the same core, which estimated an optimal SAR of ~6.3–7.1 cm/kyr24. In addition, optimal SARs estimated for the Kn1+2 in adjacent cores are ~6.79 cm/kyr in the SK1S55 and ~6.6–8.3 cm/kyr in SK1N24. Using this anchored ATS, the stratigraphic grayscale, log (Ca/Ti), and Rb/Sr datasets were transformed into the time domain.

The age model of the DSDP–516F core was established based on magnetostratigraphic and cyclostratigraphic analyses. Magnetostratigraphy places the 1198.44–1125.68 m interval within magnetochron C33r23 (Fig. S3A). According to the geomagnetic polarity time scale 2020 (GPTS 2020)27, the C33n/C33r and C33r/C34n boundaries occur at 79.9 Ma and 82.9 Ma, respectively. This depth-age model yields an average SAR of 2.43 cm/kyr for 1198.44–1125.68 m interval, consistent with Park et al.22. who estimated a SAR of 2.24 cm/kyr. Based on this magnetochronology and linear interpolation, the stratigraphic interval covered by the L* datasets (1163.42–1145.132 m) is approximately dated to 81.456–80.702 Ma (Fig. S3B). Because the timespan of this interval is shorter than two long eccentricity cycles (~754 kyr), tunning the interpreted long-eccentricity cycle in L* to a ~405-kyr periodicity is unreliable. However, previous cyclostratigraphy identified significant eccentricity and precession cycles, particularity with prominent double peaks at precession band in both magnetic susceptibility (MS) and L* datasets7,22. In the theoretical astronomical solution, precession amplitude is modulated by eccentricity. The TimeOpt algorithm28 in the R package Astrochron58 estimates the optimal SAR and constructs the ATS from this modulation relationship. We therefore applied TimeOpt to transform the stratigraphic L* datasets into the time domain. With a Monte Carlo procedure, TimeOpt identifies an optimal SAR within a given range by simultaneously optimizing: (1) eccentricity amplitude modulations within the precession band, and (2) the concentration of spectral power at specified target astronomical periods28. The optimization is quantified by a value, r2opt (Fig. S4F), which is the product of two separated correlations: (1) r2envelope, the correlation between eccentricity cycles and the precession envelope in the reconstructed time-domain datasets (Fig. S4C–E); (2) r2spectral, the correlation between the spectra of datasets and the coeval theoretical astronomical solutions (gray line in Fig. S4B, E). Here, under the constraint of magnetochronology, the tested SAR range for TimeOpt algorithm was set to 0.1–5 cm/kyr, and the target/theoretical eccentricity and precession periodicity were set at 409.6, 132.129, 126.0308, 99.9024, 96.3765, 22.6925, 21.5013, 18.4921, 18.3677 kyr, based on the La04 astronomical solution at 83 Ma30. The resulting floating ATS returns an optimal SAR of 2.37 cm/kyr, and was then anchored at the magnetochronological age of 80.702 Ma at 1145.132 m.

Grayscale extraction
The recovered cores of lower Nenjiang Formation from SK2 core were segmented with 1 m spacing. These segments were cut into two halves, with one of them being polished and covered with transparent crystal resins. These working halves were photographed. After being transformed into 8-bit format of these core photographs, the grayscale data for the lowest part of the second member of Nenjiang Formation (1145.23–1137.84 m) were extracted with 1 cm spacing using the ImageJ software. The scale for grayscale is pixel. The cracks and bentonite layers were ignored.

XRF-core scanning
Based on these working halves, the bulk elemental composition of the lower Nenjiang Formation in the SK2 core was measured at the State Key Laboratory of Biological Geology and Environmental Geology of the China University of Geosciences using an Itrax X-ray Fluorescence Core Scanner (Cox Analytical Instruments). The scanner was operated at a stratigraphic resolution of 3 mm with a Cr tube as the X-ray source. Each scan lasted for 40 s with current and voltage setting of 30 mA and 30 kV. The acquired elements are scaled with counts per second (cps), representing a semi-quantitative result. Data points associated with crocks and bentonites were removed, and the acquired datasets were resampled at 1 cm spacing to align with the resolution of grayscale. Here, the log (Ca/Ti) and Rb/Sr ratios were used due to their high quality.

Thin-section petrography, quantitative XRF, LA-ICP-MS, and XRD analyses
To calibrate the semi-quantitative XRF-core scanning results, and to confirm the petrographic and mineralogical characteristics of the lower Nenjiang Formation of SK2 core in the SLB, four thin-section and twenty-five bulk samples were taken from the box 59 of the core, between 1141.0 and 1141.7 m (Fig. S2). The core samples are preserved in the Cores and Samples Center of Nature Resources of China. Thin-section examination: The four thin-section samples were cut, mounted on glass slides, and ground and polished to a thickness of ~30 μm. Thin sections were examined under a Leica DM2500 polarizing microscope using transmitted light at the Chinese University of Geosciences, Beijing, to identify mineral phases, microstructures, and sedimentary textures. XRF: Bulk-rock major elements (Si, Ca, Ti) were analyzed by X-ray fluorescence spectrometry (XRF) using a Bruker S8 TIGER wavelength-dispersive spectrometer. The analytical procedure followed the Chinese industry standard JY/T 016–1996 General Rules for Wavelength Dispersive X-ray Fluorescence Spectrometry. The analytical process generally involved sample powdering by agate abrasion, pellet preparation, and quantitative measurement under standard conditions. Certified reference materials were analyzed together with the samples for calibration and quality control. Analytical precision is better than 5% relative standard deviation (RSD). Detection limits are typically <0.01 wt% for these major elements. LA-ICP-MS: Bulk-rock trace-elements (Rb and Sr) were analyzed by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) using a Thermo Scientific iCAP RQplus ICP-MS instrument. Analytical procedures followed the Chinese industry standard JY/T 015-1996 General Rules for Inductively Coupled Plasma-atomic Emission Spectrometry. In general, powdered rock samples were digested and introduced into the plasma for quantitative analysis. The minimum detection limit for Rb and Sr is 0.01 ppm. Certified reference materials were analyzed together with the samples to monitor instrumental drift and ensure analytical accuracy and reproducibility. XRD: The mineral composition (e.g., calcite and quartz) of bulk mudstones was determined by X-ray diffractometer (XRD) using a Bruker D8 ADVANCE diffractometer, following the Chinese industry standard SY/T 5163-2010 Analysis Method for Clay Minerals and Ordinary Non-clay Minerals in Sedimentary Rocks by the X-ray Diffraction. Powdered samples were analyzed under standard conditions, and relative mineral abundances were determined by Rietveld refinement. Analytical precision for major mineral phases was generally better than ±2 wt%. XRF, LA-ICP-MS, and XRD analyses were implemented at the Beijing Beida Zhihui Microstructure Analysis & Testing Center. Prior to analysis, the bulk samples were ground to powders (<200 mesh) by agate abrasion.

Timeseries analysis
The precise recognition of millennial-scale cycles needs high spectral resolution. Thus, we used the periodogram method to detect the periodical components within the datasets. Both conventional autoregressive-1 (AR-1)60 and robust AR-1 noise model testing61 were employed to evaluate spectra against the stochastic null model. To avoid the incorrect rejection of the null hypothesis of stochastic variability during the noise model test (false positives), multiple-testing corrections are necessary. Here, the False Discovery Rate (FDR) approach of Benjamini and Hochberg25,26 was employed to correct the null hypothesis for the periodogram spectra. The evaluative power spectrum of L* datasets was realized using Evolutive Harmonic Analysis (EHA)62 with eha function in the R package Astrochron58, the parameter setting refers to R codes in the Supplementary Information. Periodogram, AR-1 noise model testing, the multiple-testing corrections, and the TimeOpt analysis were realized in the Astrochron package58 based on R codes (Supplementary Information). A band pass filter63 was employed to isolate the presumed eccentricity, precession, and millennial-scale components from the datasets. The filtered millennial-scale components were then subjected to a Hilbert transform to estimate the amplitude modulation (AM) of millennial cycles, allowing evaluation of the relationship between millennial and astronomical cycles. The Multi-Taper Method (MTM)62 was also used as an alternative to the periodogram spectral method to detect the periodicity in the time series. Filter, Hilbert transform, and MTM analyses were implemented in Acycle software64.

Bicoherence spectra
The high-order bicoherence spectrum analysis was utilized to detect the potential nonlinear intermodulation/frequency mixing of different frequencies65,66. The coherence here refers to the coupling of two fundamental frequencies, f1 and f2, to generate a third frequency, f3, where f3=f1+f2, with their phases being coupled8,59. To reduce the spurious bicoherence peaks, the bicoherence method based on WOSA spectra67,68 was employed. The bicoherence spectra were realized using the ‘bicoherence’ function in the R-package Astrchron58,59.

Bioturbation simulation experiment
Bioturbation can destroy primary climate signals and even obscure the high-frequency cycles. To assess its impact on detecting millennial-scale components in the SK2 and DSDP–516F records under their respective sampling and sedimentation rates, we employed the bioturb function in the R package Astrochron42,58. The bioturbation algorithm is a function of mix layer depth (ML), diffusion/bioturbation intensity (G), and sedimentation accumulation rate (v)42. Two scenarios were tested (Fig. S8), representing typical marine mixing depths (5 and 10 cm), with diffusion/bioturbation intensity fixed at G=1.

Latitudinal evolution of the quarter-precessional variance in the equatorial insolation
The equatorial insolation over the last 400 kyr at latitudes between 20°S and 20°N were initially calculated following the procedure of Berger et al.11. with a 0.1° step. Then, the power ratios at the 5.5-kyr band in these calculated special insolation curves were extracted using the integratePower function in the R package Astrochron58,69. The frequency bandpass was set as 0.18 ± 0.02 cycles/kyr. See the R codes in the Supplementary Information for details.

Data availability
All relevant proxy datasets of SK2 core that support the findings of this research are available in the Supplementary Datasets 1–3. The L* datasets of DSDP–516F borehole are from de Winter et al.7. and can be downloaded at https://doi.org/10.1594/PANGAEA.828372.

Code availability
The grayscale extraction was implemented using ImageJ software, downloaded at https://imagej.net. The periodogram spectra, multitest, bicoherence spectra, and bioturbation simulation were analyzed using R package Astrochron, which can be accessed at https://cran.r-project.org/package=astrochron58. The corresponding R codes please refer to the Supplementary Information. The Bandpass filter, Hilbert transform, and MTM spectrum were implemented using the Acycle software64, which can be downloaded at https://acycle.org/.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The ms. by Zhang et al. provides evidence for the persistent occurrence of millenial climate cycles under the greenhouse state of the Campanian. Evidence is derived from two expanded archives that offer an optimal resolution for the observation of sub-milankovian cycles.\nThe paper is well-written and organized. However, I found the use of specific terminology a little confusing (see comments below).\nThe topic is of interest to a wider Earth Sciences community and beyond. Most of the literature on millennial cycles in stratigraphic archives older than the Quaternary does not attempt giving an explanation to their occurrence. In this respect, the ms. is innovative even if based on previous work done by one of the authors. Because of this, its publication would offer some new discussion points to the community.\nThere are, however, some points of concern that should be addressed.\nI report them below relative to specific sections of the ms.\nRESULTS (MILLENNIAL CYCLES IN CRETACEOUS EAST ASIA)\nWhy the precession signal only emerges in the grayscale proxy for the SK2 borehole record? This discrepancy should be addressed, considering the following results of yours: (1) the three records are phase-coupled at the 4-5 kyr timescale, (2) the 4-5 kyr signal is itself modulated by the precession signal. Moreover, while the precession signal does not appear in the Rb/Sr proxy data here, Yang et al. (2020) - also cited in this paper - identified a clear precession influence using the same proxy, from the same section, with identical laboratory conditions and resolution.\nOn the other hand, short eccentricity is not detected in the grayscale record but appears in the log (Ca/Ti) and Rb/Sr records from SK2. This also, requires some discussion/explanation.\nRESULTS (MILLENNIAL CYCLES IN CRETACEOUS SOUTH ATLANTIC)\nFor what mechanistic reason should the 405 kyr eccentricity component (influence) be detectable through spectral analysis or EHA, while the 100 kyr component (influence) remains undetected? As mentioned in the text (line 116), the 405-kyr component, and not the 100 kyr component, seems to more strongly modulate the quarter-precessional cycle.\nThe claim that the AM of the ~4–5-kyr cycles is \"phase-locked with eccentricity cycles (main text) and the subsequent explanation in Extended Data Fig. 5 caption represent a circular reasoning. The caption states that the ~95.2-kyr, ~120-kyr, and ~405-kyr filters (representing eccentricity cycles) are derived \"from the amplitude modulation curves\" themselves.\nTherefore, it becomes circular to subsequently claim a phase-lock between AM and those eccentricity signals, and it does not provide independent evidence of the authenticity of those eccentricity signals.\nTo provide robust evidence for this claim, the authors should compare the derived AM of the ~4–5-kyr cycles (or its filtered “eccentricity” components) directly against independently obtained eccentricity signals such as those obtainable from the original data (whether it is merely a trend believed to reflect short eccentricity or a real signal identified through spectral analysis).\nThe same validation method should be applied to the precession component to demonstrate its authenticity after filtering it from spectra of the 4–5-kyr cycles’ AM. This can be easily done for the grayscale and the L* records - which show a clear precession signal.\nAt line ine 134 you state that “Although the pvalues for the peaks at 9.2 kyr, 5.5 kyr, 4.8 kyr are above 0.1, the tuned L* data clearly show a semipartite and quadripartite structure of precession-related cycles, with each precession cycle comprising two semi-precessional and four quarter-precessional cycles (Fig. 2F).”\nFair claim, even though non statistically significant. For this reason, to further strengthen this finding, a rolling correlation analysis could be implemented - alongside the visual inspection method adopted by the authors and shown in Fig. 2F - to quantify the degree of agreement between the millennial-scale signal and the L* data. This is especially needed as the authors only show the visual agreement for a small window of the entire L* record.\nDISCUSSIONS (ORBITAL AND SEMI- PRECESSIONAL CYCLES)\nAt line 145 you state that: “Besides, a distinct ~9.2-kyr signals were recognized in the L* record of DSDP–516F core, aligning with semi-precessional cycles regarding the periodicities.”\nWhy should only the semi-precession component be resulting from half the 18.2kyr precession component manifest? What about the ~22.7/2 kyr component?\nHow is then demonstrated that this signal is a real semi-precession signal and not an artifact produced by a precession clipping as also discussed by the authors?\nDISCUSSION (INDIRECT PRECESSIONAL ORIGIN OF MILLENNIAL CYCLES)\nThe usage of a precession frequency of 1/18 kyr-1 to demonstrate that the 4.5 kyr cycles could be the fourth order harmonic of precession is in contradiction with the 22 kyr periodicity of precession clearly emerging from the spectral analysis of the grayscale record from SK2. On the other hand, at Site 516 where the two main components of precession (22.7 kyr and 18.2 kyr) are present, the 5.5 kyr could well be the fourth harmonic of the 22.7 kyr component, instead of a combination tone between the 4.5 kyr and 22 kyr components.\nAlso, note that the statement that “the harmonics of precession were generated based on the detected precession frequency in grayscale” (lines482-483) does not stand. Again, the detected precession at SK2 has period of 22 kyr NOT 18 kyr. Using a synthetic 18 kyr signal to verify that its fourth harmonic (4.5 kyr) is modulate by the signal itself does not prove anything or, at least, redundant. You should filter the precession from the grayscale and compare it with the AM of the 4.5 kyr component, as already mentioned.\nMINOR COMMENTS\nAt line 199 you state that “Equatorial insolation is calculated from the differences between the two insolation maxima at spring and fall equinoxes and the two insolation minima at summer and winter solstices”\nThis is potentially misleading, to clarify. While it is clear that the authors intend to define a parameter representing the seasonal amplitude of insolation (as introduced for the first time by Berger et al. 2006 to explain the ~5 kyr cycles), the phrasing suggests this is how the instantaneous or mean equatorial insolation itself is computed. “The seasonal contrast of equatorial insolation is calculated ..” would be then more appropriate.\nThe terms bi-hemisphere or bi-hemispheric are potentially confusing. As an example, the following sentence (lines 150-152) “Bi-hemisphere maximum insolation corresponds to climate responses to summer solstice insolation peaks in both hemispheres.” could probably be written as “Alternate insolation maxima occur on summer solstices in the two hemispheres”.\nFIG. 4D the curve does not represent tropical insolation", "reviewer_2": "This manuscript is on an interesting topic, but I feel that large portions of the science in the paper is not presented clearly enough in order to be published in a journal with a wide readership like Nature Communications.\nThis study seeks to show millennial-timescale variability during a greenhouse period, without ice sheets, during the Cretaceous. They do this by comparing records from an ocean core and from sedimentary strata on land in China. In order to make the argument that they are resolving millennial-scale variability, they need to convert the depth domain of the sedimentary records into time. They do this by using various techniques including paleomagnetics, fossil stratigraphy, and most importantly astrochronology.\nOne way in which their paper falls short is that the chronological model is not presented in a clear way, and are not explained thoroughly – the text simply refers to SOM, the SOM refers to additional figures, and the additional figures are simply plots of their age modelling, but there is no thorough explanation of the errors, assumptions, and limitations of the techniques used. I found the additional figure captions to be of little help, especially for a non-specialist to understand, but even for someone who works on paleoclimate. For example Extended Data Fig. 3 is an important figure in linking up their chronological controls, but the figure caption simply refers to various techniques. All of the figures about age modelling are equally obscure. They may have done a lot of work developing their age model, but there is no explanation for the reader on what the exactly the constraints are other than trying to reconstruct the authors thought-process from the many figures with many panels. I know that Andre Berger is a well-respected researcher in this field, and does excellent work on this topic, so it may be that the age model presented is well supported, but the presentation feels rushed. I think that this type of research would be better presented in a specialist journal, where they have more text space to fully explain and substantiate how they developed and their age model. As it stands, with little explanation in the text the reader (and this reviewer) would be hard-pressed to make a connection between the data in the figures.\nThe reason why it is really important that their methodology is substantiated with a clear explanation, is that they are claiming to distinguish millennial cycles and use them to fingerprint climate dynamics. Given that the core has an “optimal SAR of 2.37 cm/kyr”, which seems quite low, what is the actual sample spacing or integrated beam size of the measurement relative to the cyclicity, and how does this relate to bioturbation?\nI also have questions about what the chemical measurements are caused by. In the SOM, the authors suggest that their proxies (Rb/Sr and Ca/Ti) are controlled by calcite vs detrital silicate sediment. I understand that they used an XRF to scan the core, but have they calibrated to measurements to actual wet-chemical measurements? Do they know what sedimentary phases or minerals are giving them the high Ca, Sr, Ti signal?\nAlso can we be sure that this chemical variation is actually a climate proxy? Calcite need not be primary, in fact many terrestrial sedimentary strata have diagenetic calcite. Is there any evidence that the SLB Ca and Sr is reflecting the amount of primary calcite, rather than secondary or diagenetic calcite? Rocks of this age have often undergone secondary calcite precipitation if there are abundant porosity, which could mean that the apparent cyclicity is due to changes in grain size or porosity followed by post-depositional fluid flow and cementation, rather than a primary climate signal.\nIt is also possible that these chemical cycles are due to autogenic sedimentary processes, which are internal to a sedimentary basin. Examples of this include river channel migration or delta lobe switching controlling sediment deposition within the system itself, rather than allogenic processes which are external to the depositional system, such as changes in climate. Even if this apparent cyclicity is similar to that occurring in other basins or sediment records, given the large range of cycle periodicity (see lines 96-101), it would be surprising not to find some kind of match. I would like to see some substantiation that what they are measuring are primary sediment phases, and that the geological and sedimentary setting precludes autogenic sedimentary processes."}, "author_response": "Response to Reviewer #1:\n(General comments) The ms. by Zhang et al. provides evidence for the persistent occurrence of millenial climate cycles under the greenhouse state of the Campanian. Evidence is derived from two expanded archives that offer an optimal resolution for the observation of sub-milankovian cycles.\nThe paper is well-written and organized. However, I found the use of specific terminology a little confusing (see comments below).\nThe topic is of interest to a wider Earth Sciences community and beyond. Most of the literature on millennial cycles in stratigraphic archives older than the Quaternary does not attempt giving an explanation to their occurrence. In this respect, the ms. is innovative even if based on previous work done by one of the authors. Because of this, its publication would offer some new discussion points to the community.\nThere are, however, some points of concern that should be addressed.\nI report them below relative to specific sections of the ms.\nOur reply: Thank you for your generous comments and for your appreciation of our work. We also sincerely appreciate your constructive concerns, which will greatly help enhance the overall quality of the manuscript.\nMajor points:\nRESULTS (MILLENNIAL CYCLES IN CRETACEOUS EAST ASIA)\nComment 1.1 Why the precession signal only emerges in the grayscale proxy for the SK2 borehole record? This discrepancy should be addressed, considering the following results of yours: (1) the three records are phase-coupled at the 4-5 kyr timescale, (2) the 4-5 kyr signal is itself modulated by the precession signal. Moreover, while the precession signal does not appear in the Rb/Sr proxy data here, Yang et al. (2020) also cited in this paper - identified a clear precession influence using the same proxy, from the same section, with identical laboratory conditions and resolution. On the other hand, short eccentricity is not detected in the grayscale record but appears in the log (Ca/Ti) and Rb/Sr records from SK2. This also, requires some discussion/explanation.\nOur reply: We acknowledge that the ~22-kyr precession signal is strongest in the grayscale record. The 22-kyr signal is also present in the Rb/Sr datasets, and the corresponding power peak exceeds the ML96 95% robust confidence level (Extended Data Fig. 5), so this result is consistent with previous findings of Yang et al. (2020). Moreover, a power peak at 1/22 kyr-1 is observed in the log (Ca/Ti), but it falls below the 90% confidence level in both conventional and robust AR1 estimations (Fig. 3B). Nevertheless, the ~22-kyr filter in log (Ca/Ti) is phase-locked with those in grayscale and Rb/Sr (Extended Data Fig. 6). In the revised Fig. 3B and Extended Data Fig. 5, we labeled this ~22-kyr power peak in the corresponding spectra.\nIn terms of the ~100-kyr signal, both log (Ca/Ti) and Rb/Sr display strong and significant power at this frequency, whereas grayscale shows no significant power peak at this frequency. Nevertheless, the detrended grayscale (robust LOWESS detrending with a 55% window) displays a prominent ~100-kyr trend (Fig. R1), and the ~100-kyr filter in the detrended grayscale is phase-locked with the same filters in both log (Ca/Ti) and Rb/Sr (Extended Data Fig. 6A, D, G). Moreover, the power spectra of detrended grayscale show weak power in the ~100-kyr band, although it falls below the 90% confidence level (Fig. 3A; Fig. R1). We propose that the monotonically increasing long-term trend in the grayscale record suppresses the manifestation of the ~100-kyr cycle. In the revised manuscript, we have replaced the spectra of the raw grayscale dataset with those of the detrended grayscale dataset, and the ~100-kyr power peak has been labelled in Fig. 3A.\nFig. R1 | The ~100-kyr cycle in grayscale of SK2 core. A, the raw grayscale dataset (black) and its long-term trend (robust LOWESS detrending with 55% window; blue); B, the detrended grayscale dataset (black) and its ~100-kyr filter (0.0105±0.0035 cycle/kyr; blue).\nFig. R2 | A&C, the summer insolation at 45°N over the past 4000 ka, and its periodogram spectra. B&D, the equatorial insolation over the past 4000 ka, and its periodogram spectra.\nOverall, the spectral characteristics of grayscale, Rb/Sr, and log (Ca/Ti) indicate that both ~100-kyr and ~22-kyr signals are prominent and significant in the Rb/Sr dataset. In contrast, the ~100-kyr signal is weak in the grayscale but strong in log (Ca/Ti), whereas the ~22-kyr signal is weak in log (Ca/Ti) yet most pronounced in grayscale. This contrasting expression of the ~100-kyr and ~22-kyr signals in log (Ca/Ti) and grayscale likely reflects differences in the proxies’ sensitivity to local versus equatorial insolation forcing. Equatorial insolation exhibits dominant ~100-kyr power but no power at the precession band, as precession-band power is redistributed to its harmonics and amplitude modulation bands (Berger et al., 2006) (Fig. R2). In contrast, local insolation forcing is dominated by precession power, with minimal variance at the ~100-kyr eccentricity band (Fig. R2).\nIn the first paragraph of the discussion of the revised manuscript, we added an extra explanation and wrote “Although the ~100-kyr and ~22-kyr signals in grayscale, log (Ca/Ti), and Rb/Sr show different confidence level, the ~100-kyr and ~20-kyr filters across these datasets are phase-coupled with each other (Extended Data Fig. 6A, D, G), substantiating their authenticity. This differential expression likely reflects differences in proxy sensitivity to different insolation forcing; for example, local insolation is dominated by precession peaks, whereas tropical insolation is dominated by eccentricity cycles (Berger et al., 2006).”\nRESULTS (MILLENNIAL CYCLES IN CRETACEOUS SOUTH ATLANTIC)\nComment 1.2 For what mechanistic reason should the 405 kyr eccentricity component (influence) be detectable through spectral analysis or EHA, while the 100 kyr component (influence) remains undetected? As mentioned in the text (line 116), the 405-kyr component, and not the 100 kyr component, seems to more strongly modulate the quarter-precessional cycle.\nOur reply: Thank you for your comments. The ~100-kyr components in the L* datasets are also prominent and significant (above the 95% ML96 robust confidence level). In Fig. 3C, we did not label them previously due to layout constraints.\nThe presence of short-eccentricity signals in the L* dataset was also validated by Winter et al (2014). Their wavelet analysis of L∗ in the stratigraphic domain revealed a significant ∼2.5 m periodicity. Based on the interpretation of Park et al. (1993), Winter et al. (2014) inferred an average sedimentation accumulation rate (SAR) of ∼ 2.3 cm/ kyr. In our study, the TimeOpt method yields an optimal SAR of ~2.37 cm/kyr. Under this constraint, the ~2.5 m wavelength corresponds to ~105-kyr cycles. Accordingly, we now label the ~100-kyr component in the revised Fig. 3C.\nTheoretically, equatorial insolation contains 405-kyr components that are as strong as, or even stronger than, the ~100-kyr components (Berger et al., 2006) (Fig. R2D). This makes it reasonable that the ~405-kyr signal, rather than the ~100-kyr signal, exerts the dominant modulation on the quarter-precessional cycle. Another possible explanation for the stronger manifestation of the ~405-kyr component is that longer-period signals are generally more readily preserved in the climate system.\nComment 1.3 The claim that the AM of the ~4–5-kyr cycles is \"phase-locked with eccentricity cycles (main text) and the subsequent explanation in Extended Data Fig. 5 caption represent a circular reasoning. The caption states that the ~95.2-kyr, ~120-kyr, and ~405-kyr filters (representing eccentricity cycles) are derived \"from the amplitude modulation curves\" themselves. Therefore, it becomes circular to subsequently claim a phase-lock between AM and those eccentricity signals, and it does not provide independent evidence of the authenticity of those eccentricity signals. To provide robust evidence for this claim, the authors should compare the derived AM of the ~4–5-kyr cycles (or its filtered “eccentricity” components) directly against independently obtained eccentricity signals such as those obtainable from the original data (whether it is merely a trend believed to reflect short eccentricity or a real signal identified through spectral analysis).\nThe same validation method should be applied to the precession component to demonstrate its authenticity after filtering it from spectra of the 4–5-kyr cycles’ AM. This can be easily done for the grayscale and the L* records - which show a clear precession signal.\nOur reply: Thank you for your comment and suggestions. Owing to the chaotic dynamical nature of the Solar System, precise astronomical solutions are unavailable beyond ~50 Ma (Laskar et al., 2004; Ma et al., 2017). Therefore, we compared the ~100-kyr and ~20-kyr components in the amplitude modulation (AM) of the ~4–5-kyr filters with those in the raw datasets. As shown in Extended Data Fig. 6, the amplitude modulation of the ~4–5-kyr cycles is phase-coupled with eccentricity.\nComment 1.4 At line 134 you state that “Although the p-values for the peaks at 9.2 kyr, 5.5 kyr, 4.8 kyr are above 0.1, the tuned L* data clearly show a semipartite and quadripartite structure of precession-related cycles, with each precession cycle comprising two semi-precessional and four quarter-precessional cycles (Fig. 2F).”\nFair claim, even though non statistically significant. For this reason, to further strengthen this finding, a rolling correlation analysis could be implemented - alongside the visual inspection method adopted by the authors and shown in Fig. 2F - to quantify the degree of agreement between the millennial-scale signal and the L* data. This is especially needed as the authors only show the visual agreement for a small window of the entire L* record.\nOur reply: Thank you for your suggestions. We calculated the evolutionary correlation coefficient between L* and its ~4–5-kyr filter (0.215±0.055 cycles/kyr), as well as between L* and its ~8–14-kyr filter (0.095±0.025 cycles/kyr). Both evolutionary correlation coefficients display significant eccentricity signals that are phase-locked with eccentricity in the raw L* datasets (Extended Data Fig.7). This indicates that the strength of the semi- and quarter- processional cycles is modulated by eccentricity, consistent with the theoretical curve of tropical (bi-hemisphere maximum insolation and equatorial maximum insolation) and equatorial insolation. The reason is that the semi- and quarter-processional cycles are most prominent above the average amplitude modulation of precession cycles (or eccentricity) in these theoretical insolation curves (Berger et al., 2006; Hagelberg et al., 1994). Thus, we conclude that the semi- and quarter- processional cycles in the L* datasets are significant and authentic, even though the p-values for the peaks at 9.2 kyr, 5.5 kyr, 4.8 kyr are above 0.1.\nAt the end of the Results of the revised main file, we supplemented that ‘Moreover, the evolutionary correlation coefficients between L* and its ~8–14-kyr filter (0.095±0.025 cycles/kyr) and between L* and its ~4–5-kyr filter (0.215±0.055 cycles/kyr) display significant eccentricity signals that are phase-locked with the eccentricity in L* datasets (Extended Data Fig.7).’\nDISCUSSIONS (ORBITAL AND SEMI- PRECESSIONAL CYCLES)\nComment 1.5 At line 145 you state that: “Besides, a distinct ~9.2-kyr signals were recognized in the L* record of DSDP–516F core, aligning with semi-precessional cycles regarding the periodicities.” Why should only the semi-precession component be resulting from half the 18.2kyr precession component manifest? What about the ~22.7/2 kyr component?\nHow is then demonstrated that this signal is a real semi-precession signal and not an artifact produced by a precession clipping as also discussed by the authors\nOur reply: Thank you for your careful review. The double peaks in the precessional band (~22.4 kyr and 19.0 kyr) were previously reported by Park et al., (1993) based on magnetic susceptibility (MS) from the same interval of the same core. In our study, the astronomically-tuned L* datasets further unveiled these characteristics (~22.7 kyr and ~18.2 kyr) in DSDP–516F. Accordingly, we acknowledge that there should be two semi-precessional bands, centered at ~2/22.7 and ~2/18.2 kyr-1 bands, respectively.\nIn the L* record of DSDP–516F core, in addition to the 9.2 kyr power peaks, other power peaks also occur near the ~10-kyr band, including ~10.2 kyr, ~10.8 kyr, and especially the ~11.3 kyr peak (above 95% significant level of robust AR1 model). These were not labeled previously due to figure layout constraints. In the revised Fig. 3C, we have now included labels for these power peaks.\nFig. R3 | The clipping effect causing artifacts. A, the clipped spring insolation at equator (in blue) and its semi-precessional filter (in brown). B, the periodogram spectra of the clipped spring insolation based on ML96 robust AR-1 estimations. C&E, the maximum equatorial insolation and the equatorial insolation (in blue), and their semi-precessional filters (in brown). SE, FE, SS, WS denote the spring equinox, fall equinox, summer solstice, winter solstice, respectively. D&F, the periodogram spectra of the maximum equatorial insolation and the equatorial insolation. Despite the spectra of the clipped insolation displays prominent and significant power peaks at precession harmonic bands (e.g., ~11.5, ~10.4, ~5.8 and ~5.5kyr), the comparison between clipped insolation and its semi-precessional filter shows a notable mismatch, especially where the former exhibits a basal plateau and the latter an arching trend. This spectral ‘false negative’ is called as an artifact.\nIn terms of your subsequent concern about how it is demonstrated that this signal represents a real semi-precession signal rather than an artifact produced by a precession clipping as also discussed by the authors?\nThe clipping effect on precession can redistribute power from the fundamental precessional frequency into frequency bands associated with its amplitude modulation (such as ~100-kyr and ~405-kyr eccentricity) and its harmonics (such as semi-precessional cycles) (Fig. R3B). Consequently, the spectral characteristics of clipped precession are same as those of equatorial insolation and maximum equatorial insolation (Fig. R3B, D, F). However, a comparison between the clipped precession signal and its semi-precessional filter shows a distinct mismatch, especially where the former exhibits a basal plateau, whereas the latter shows an arching trend (Fig. R3A). This spectral ‘false negative’ is therefore identified as an artifact. Please also refer to the discussion by Hagelberg et al (1994).\nOn the contrary, in Fig 2F, we compared the raw L* and its ~10-kyr filter, the semipartite structure of precession-related cycles is evident (each precession cycles comprise two semi-precessional cycles). Besides, our calculated evolutionary correlation coefficient between the L* and its ~8–14-kyr filter (0.095±0.025 cycles/kyr) also corroborates the significance and authenticity of the semi-precessional cycles in the L* datasets (Extended Data Fig. 7). Please also see our reply to the comment 1.4 above. Thus, these results confirm that the semi-precessional signal in the L* datasets is not an artifact of precession clipping.\nIn the last paragraph of the ‘Orbital and semi- precessional cycles’ of the revised manuscript, we added that “Notably, simple clipping of precession can also produce prominent semi-precesional power peaks in the corresponding spectra; however, the resulting semi-precesional cycle is not a true single and is identified as a spectral artifact (Hagelberg et al., 1994). In contrary, the comparison of L* with its ~10-kyr filter reveals a clear semipartite structure, with each precession cycle comprising two semi-precessional cycles (Fig. 2F). Moreover, the ~10-kyr variance in L* is high (>50%) during intervals when eccentricity is above its average value (Extended Data Fig. 7). These findings demonstrate the significance of the semi-precessional cycles in L* datasets. Thus, the observed semi-precessional cycles cannot be attributed to precession clipping.”\nDISCUSSION (INDIRECT PRECESSIONAL ORIGIN OF MILLENNIAL CYCLES)\nComment 1.6 The usage of a precession frequency of 1/18 kyr-1 to demonstrate that the 4.5 kyr cycles could be the fourth order harmonic of precession is in contradiction with the 22 kyr periodicity of precession clearly emerging from the spectral analysis of the grayscale record from SK2. On the other hand, at Site 516 where the two main components of precession (22.7 kyr and 18.2 kyr) are present, the 5.5 kyr could well be the fourth harmonic of the 22.7 kyr component, instead of a combination (difference) tone between the 4.5 kyr and 22 kyr components.\nOur reply: Thank you for your comments. The amplitude modulation of the ~ 4–5-kyr filters in all datasets from both the SK2 and DSDP-516F displays distinct eccentricity cycles, perfectly aligning with the theoretically-calculated equatorial insolation curve (Berger et al., 2006) (Fig. 4). Consequently, we conclude that the ~4–5-kyr signals in the datasets reflect the 4th harmonics of precession cycles arising from equatorial insolation (despite we also proposed other interpretations, such as combined effects of equatorial maximum insolation and bi-hemisphere maximum insolation).\nIn the equatorial insolation, precessional components are absent because precession power is redistributed into its amplitude modulation and harmonic components (Berger et al., 2006) (Fig. R2B). Thus, the ~ 4–5-kyr cycles in grayscale, log (Ca/Ti), Rb/Sr, and L* are independent of the precession component detected in these datasets.\nAs the ~4.4–4.5 kyr cycles are most prominent, we infer that the corresponding precession periodicity generating these harmonics is ~18 kyr (since 4.5×4 = 18 kyr). Thus, using a precession frequency of 1/18 kyr-1 to demonstrate that the 4.5 kyr cycles represent the fourth-order harmonic of precession does not contradict the ~22 kyr periodicity of precession detected in these datasets.\nIn terms of the ~5.5-kyr cycles, we proposed that the ~5.5-kyr cycles in datasets from the Songliao Basin are most likely derive from the difference tone (4/18 - 1/22 = 1/5.7) between the modulator (1/22) and carrier (4/18). This is because the ~4–5-kyr cycles are amplitude-modulated by the 22-kyr cycles. According to the amplitude modulate mechanism, two sidebands should form around the ~4.5-kyr peaks: the combination tone (4/18 + 1/22 = 1/3.7) and difference tone (4/18 - 1/22 = 1/5.7) between the modulator (1/22) and carrier (4/18). This interpretation is underpinned by the bicoherence analysis of the grayscale data, which reveals strong coherence at the triad point (1/5.8, 1/22, 1/4.5) (Fig. 3A).\nHowever, in terms of the ~5.5-kyr signals in the L* datasets, it may reflect the third harmonic of the ~18.2-kyr precessional cycles, as indicated by the triad point (1/9.2, 1/18.2, 1/6.1) in the bicoherence plot (Fig. 3C). Alternatively, they could arise from the fourth harmonic of the ~22.7-kyr precessional cycles, since 4/22.7 = 1/5.7.\nIn the first paragraph of ‘Indirect precessional origin of millennial cycles’ of the revised manuscript, we wrote “Besides the most prominent ~4–5-kyr cycles, there are two millennial-scale peaks at 1/3.7 (1/3.6) kyr-1, and 1/5.7 (1/5.8) kyr-1 in datasets of the SLB, which bracket the quarter-precessional cycles at 1/4.5 kyr-1 (Fig. 3A&B; Extended Data Fig. 5). These two power peaks can be interpreted as the consequence of amplitude modulation between precession and quarter-precessional cycles. This is because the amplitudes of the ~4–5-kyr cycles are also modulated by the ~22-kyr precession cycle (Extended Data Fig. 6). According to the amplitude modulation mechanism, two sidebands should form around the ~4.5-kyr peaks: the combination tone (4/18 + 1/22 = 1/3.7) and difference tone (4/18 - 1/22 = 1/5.7) between the modulator (1/22) and carrier (4/18). This interpretation is underpinned by the bicoherence analysis of the grayscale data, which reveals strong coherence at the triad point (1/5.8, 1/22, 1/4.5) (Fig. 3A). However, the ~5.5-kyr signals in the L* may instead reflect the third harmonic of the ~18.2-kyr precessional cycles, as indicated by the triad point (1/9.2, 1/18.2, 1/6.1) in the bicoherence plot (Fig. 3C). Alternatively, they could arise from the fourth harmonic of the ~22.7-kyr precessional cycles, since 4/22.7 = 1/5.7.”\nComment 1.7 Also, note that the statement that “the harmonics of precession were generated based on the detected precession frequency in grayscale” (lines482-483) does not stand. Again, the detected precession at SK2 has period of 22 kyr NOT 18 kyr. Using a synthetic 18 kyr signal to verify that its fourth harmonic (4.5 kyr) is modulate by the signal itself does not prove anything or, at least, redundant. You should filter the precession from the grayscale and compare it with the AM of the 4.5 kyr component, as already mentioned.\nOur reply: Thank you for your suggestion. As we explained in Comment 1.6, equatorial insolation contains no precession components; thus, the ~4–5-kyr cycles in grayscale, log (Ca/Ti), Rb/Sr, and L* are independent of the precession component detected in these datasets. Nevertheless, to avoid of confusion, we have remove the generation of precession harmonics, as it is unnecessary and, as you noted, redundant.\nFollowing your advice, it is clear that the ~4–5-kyr cycles in the datasets are amplitude-modulated by the ~22-kyr precession cycles. The amplitude modulation of ~4–5-kyr cycles displays ~22-kyr precession cycle, which is phase-locked to the ~22-kyr filter in raw datasets (Extended Data Fig. 6). Based on this modulation relationship, two sidebands should emerge at ~1/3.7 (1/4.5+1/22=1/3.7) and ~1/5.7 kyr-1 (1/4.5-1/22=1/5.7), consistent with the observed frequencies.\nMINOR COMMENTS\nComment 1.8 At line 199 you state that “Equatorial insolation is calculated from the differences between the two insolation maxima at spring and fall equinoxes and the two insolation minima at summer and winter solstices”. This is potentially misleading, to clarify. While it is clear that the authors intend to define a parameter representing the seasonal amplitude of insolation (as introduced for the first time by Berger et al. 2006 to explain the ~5 kyr cycles), the phrasing suggests this is how the instantaneous or mean equatorial insolation itself is computed. “The seasonal contrast of equatorial insolation is calculated.” would be then more appropriate.\nOur reply: Thank you for your very careful check. In the revised manuscript, we revised it following your suggestion, “The equatorial insolation, defined as the largest amplitude of the seasonal cycle of mean irradiance at the equator, is estimated from the difference between mean irradiance at the equinoxes (spring and fall) and at the solstices (summer and winter) (Berger et al., 2006) (Fig. 4A).”\nComment 1.9 The terms bi-hemisphere or bi-hemispheric are potentially confusing. As an example, the following sentence (lines 150-152) “Bi-hemisphere maximum insolation corresponds to climate responses to summer solstice insolation peaks in both hemispheres.” could probably be written as “Alternate insolation maxima occur on summer solstices in the two hemispheres”.\nOur reply: Thank for your suggestions. We revise the definition of Bi-hemisphere maximum insolation as follows: Bi-hemisphere maximum insolation corresponds to climate responses to insolation maxima at both hemispheric summer solstices.\nComment 1.10 FIG. 4D the curve does not represent tropical insolation.\nOur reply: Thank you for your reminder. This was a mistake: Fig. 4A&D represent the equatorial insolation, not the tropical insolation. We have corrected this error in the revised manuscript.\nReviewer #2 (Remarks to the Author):\nThis manuscript is on an interesting topic, but I feel that large portions of the science in the paper is not presented clearly enough in order to be published in a journal with a wide readership like Nature Communications.\nOur replay: Thank you for your general comments. Your concerns would intensively help us improve the quality of the manuscript.\nComment 2.1 This study seeks to show millennial-timescale variability during a greenhouse period, without ice sheets, during the Cretaceous. They do this by comparing records from an ocean core and from sedimentary strata on land in China. In order to make the argument that they are resolving millennial-scale variability, they need to convert the depth domain of the sedimentary records into time. They do this by using various techniques including paleomagnetics, fossil stratigraphy, and most importantly astrochronology. One way in which their paper falls short is that the chronological model is not presented in a clear way, and are not explained thoroughly – the text simply refers to SOM, the SOM refers to additional figures, and the additional figures are simply plots of their age modelling, but there is no thorough explanation of the errors, assumptions, and limitations of the techniques used. I found the additional figure captions to be of little help, especially for a non-specialist to understand, but even for someone who works on paleoclimate. For example, Extended Data Fig. 3 is an important figure in linking up their chronological controls, but the figure caption simply refers to various techniques. All of the figures about age modelling are equally obscure. They may have done a lot of work developing their age model, but there is no explanation for the reader on what the exactly the constraints are other than trying to reconstruct the authors thought-process from the many figures with many panels. I know that Andre Berger is a well-respected researcher in this field, and does excellent work on this topic, so it may be that the age model presented is well supported, but the presentation feels rushed. I think that this type of research would be better presented in a specialist journal, where they have more text space to fully explain and substantiate how they developed and their age model. As it stands, with little explanation in the text the reader (and this reviewer) would be hard-pressed to make a connection between the data in the figures. The reason why it is really important that their methodology is substantiated with a clear explanation, is that they are claiming to distinguish millennial cycles and use them to fingerprint climate dynamics.\nOur replay: Thank you for your concerns and helpful suggestions. We apologize that the chronology of the DSDP–516F datasets was not clearly presented, which may have caused confusion.\nIn the revised manuscript, we have added a ‘Chronology’ section in the ‘Materials and Method’ to describe in detail how the age models of the two cores were established.\nSpecifically, the age model of the SK2 core was constructed based on U-Pb CA-ID-TIMS zircon dating of bentonite layers, combined with cyclostratigraphic analysis of gamma-ray logging data. The chronology of DSDP–516F was established based on magnetostratigraphy and cyclostratigraphic analyses of L* datasets.\nThe following are our ‘Chronology’ part supplemented in the ‘Materials and Method’ section: The age model of the K n1+2 in the SK2 core was established based on the U-Pb 2 chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS) zircon dating and cyclostratigraphic analysis (Zhang et al., 2022). The zircon U-Pb CA-ID-TIMS dating of four bentonite layers in the depositions of K n1+2 of SK2 core (at depth of 2 1139.39 m, 1142.49 m, 1218.61 m, and 1236.81 m) yielded ages of 83.35 ± 0.11 Ma, 83.498 ± 0.052 Ma, 84.45 ± 0.14 Ma, and 84.709 ± 0.064 Ma, respectively (Extended Data Fig. 1). Based on this depth-age model, the average sedimentary accumulation rate (SAR) for the entire K n1+2 2 in the SK2 core is ~7.45 cm/kyr. Under this SAR constraint, long-eccentricity signals were recognized and filtered from the Gamma Ray (GR) logging datasets of SK2. Because long-eccentricity cycles have remained stable over the past 250 Ma (Laskar et al., 2004), the filtered and interpreted long eccentricity signals in GR were tied to the 405-kyr periodicity to establish a floating astronomical time scale (ATS) for the K n1+2 in the SK2 core (Zhang et al., 2022). 2 This floating ATS was then anchored to the U-Pb age at 1142.49 m, yielding a SAR of ∼6.36 cm/kyr for the studied interval of this study (1137.84–1145.2 m). This tuned SAR agrees with Monte Carlo numerical simulations based on GR datasets from the same core, which estimated an optimal SAR of ∼6.3–7.1 cm/kyr (Zhang et al., 2022). In addition, optimal SARs estimated for the K n1+2 in adjacent cores are ∼6.79 cm/kyr in the SK1S (Wu et al., 2022) and ∼6.6–8.3 2 cm/kyr in SK1N (Zhang et al., 2022). Using this anchored ATS, the stratigraphic grayscale, log (Ca/Ti), and Rb/Sr datasets were transformed into the time domain.\nThe age model of the DSDP–516F core was established based on magnetostratigraphic and cyclostratigraphic analyses. Magnetostratigraphy places the 1198.44–1125.68 m interval within magnetochron C33r (Hamilton and Suzyumov, 1983) (Extended Data Fig. 3A). According to the geomagnetic polarity time scale 2020 (GPTS2020; Ogg, 2020), the C33n/C33r and C33r/C34n boundaries occur at 79.9 Ma and 82.9 Ma, respectively. This depth-age model yields an average SAR of 2.43 cm/kyr for 1198.44–1125.68 m interval, consistent with Park et al (1993), who estimated a SAR of 2.24 cm/kyr. Based on this magnetochronology and linear interpolation, the stratigraphic interval covered by the L* datasets (1163.42–1145.132 m) is approximately dated to 81.456–80.702 Ma (Extended Data Fig. 3B). Because the timespan of this interval is shorter than two long eccentricity cycles (~754 kyr), tunning the interpreted long-eccentricity cycle in L* to a 405-kyr periodicity is unreliable. However, previous cyclostratigraphy identified significant eccentricity and precession cycles, particularity with prominent double peaks at precession band in both magnetic susceptibility (MS) and L* datasets (Park et al., 1993; de Winter et al., 2014). In the theoretical astronomical solution, precession amplitude is modulated by eccentricity. The TimeOpt algorithm (Meyers, 2015) in the R package ‘Astrochron’ (Meyers, 2014) estimates the optimal SAR and constructs the ATS from this modulation relationship. We therefore applied TimeOpt to transform the stratigraphic L* datasets into the time domain. With a Monte Carlo procedure, TimeOpt identifies an optimal SAR within a given range by simultaneously optimizing: (1) eccentricity amplitude modulations within the precession band, and (2) the concentration of spectral power at specified target astronomical periods (Meyers, 2015). The ‘optimization’ is quantified by a value, r2 (Extended Data Fig. 4F), which is the product of two separated correlations: (1) opt r2 , the correlation between eccentricity cycles and the precession envelope in the envelopment reconstructed time-domain datasets (Extended Data Fig. 4C–E); (2) r2 , the correlation spectral between the spectra of datasets and the coeval theoretical astronomical solutions (gray line in Extended Data Fig. 4B&E). Here, under the constraint of magnetochronology, the tested SAR range for TimeOpt algorithm was set to 0.1–5 cm/kyr, and the target/theoretical eccentricity and precession periodicity were set at 409.6, 132.129, 126.0308, 99.9024, 96.3765, 22.6925, 21.5013, 18.4921, 18.3677 kyr, based on the La04 astronomical solution at 83 Ma (Laskar et al., 2004). The resulting floating ATS returns an optimal SAR of 2.37 cm/kyr, and was then anchored at the magnetochronological age of 80.702 Ma at 1145.132 m.\nComment 2.2 Given that the core has an “optimal SAR of 2.37 cm/kyr”, which seems quite low, what is the actual sample spacing or integrated beam size of the measurement relative to the cyclicity, and how does this relate to bioturbation?\nOur replay: Thanks for your concerns. As we elaborated in the ‘Chronology’ section, the magnetostratigraphic analysis (Hamilton and Suzyumov, 1983) yields an average SAR of ~2.43 cm/kyr for the 1198.44–1125.68 m interval. By comparison, Park et al. (1993) estimated a SAR of 2.24 cm/kyr for the similar interval by correlating the 1.9-m, 2.35-m, and 2.9-m cycles in the magnetic susceptibility (MS) datasets with the three major precession periods of 22.4 kyr, 19 kyr, and 16.9 kyr. Here, using the TimeOpt algorithm, we obtained an optimal SAR of ~2.37 cm/kyr (Extended Data Fig. 4). The close agreement among these three independent estimates supports the robustness of the SAR.\nFor the L* datasets in the DSDP-516F, the average sampling rate is 0.015 m. Under the SAR constraint of 2.37 cm/kyr, the average temporal sampling rate is 0.63 kyr (Fig. R4), which is sufficient to recognize the millennial-scale cycles.\nFig. R4 | The sampling rate of the L* datasets of DSDP–516F core. A, the stratigraphic sampling rate of the raw L* datasets, with an average spacing at 0.015 m. B, the sampling rate of the astronomically-tuned L* datasets, with a temporal spacing rate of 0.63 kyr.\nRegarding bioturbation, it preferentially attenuates the amplitude of higher-frequency millennial-scale cycles relative to lower-frequency orbital cycles (Liu et al., 2021). This indicates that the observed magnitude of millennial-scale variability is strongly underestimated, and the true oscillations are likely much stronger. To evaluate this effect, we applied the bioturbation mixing model (function “bioturb” in Astrochron) to simulate the impact of bioturbation on the L* datasets of DSDP-516F (Extended Data Fig. 8). The bioturbation mixing model is a function of mix layer depth (ML), sedimentation rate (v), and diffusion/bioturbation intensity (G) (Liu et al., 2021). We performed simulations with ML values of 5 cm and 10 cm at G = 1, using a sedimentation rate of 2.37 cm/kyr (the SAR of DSDP–516F). The simulated timeseries were generated at the sampling resolution of the L* datasets (0.63 kyr).\nIn the case of a 10-cm mixing depth (blue in Extended Data Fig. 8), the simulation shows that the amplitudes of ~100-kyr, ~20-kyr cycles decreased by 5.18% and 56.28%, respectively. The ~10-kyr and ~5-kyr components remain detectable, retaining 15.18% and 3.38% of their original amplitudes, respectively. In contrast, the ~2.5-kyr and ~1.25-kyr components are almost entirely smoothed out.\nSimilarly, we simulated the effect of bioturbation on datasets from the Songliao Basin (Extended Data Fig. 8). Under both scenarios of mixing depth, millennial-scale cycles, even at ~1.25 kyr, can be clearly identified.\nThese simulations indicate that caution is warranted when interpreting much higher-frequency (~1–3 kyr) cycles in the L* of DSDP–516F, as bioturbation may completely smooth climate oscillations at this timescale under the deposition condition of the core. However, the ~1–3 kyr cycles recognized in SK2 are significant and are most likely derived from the external forcing (Extended Data Fig. 8).\nIn this study, our primary goal is to provide robust geological evidence for the theoretically proposed equatorial insolation forcing (Berger et al., 2006). Accordingly, we focus on the quarter-prcessional cycles at ~5-kyr bands. Extended Data Fig. 8 demonstrates that even under strong bioturbation and deep mixing, the ~5-kyr signals remain detectable at the relatively low SAR of DSDP–516F. We thus conclude that the current sampling resolution is sufficient to resolve the significant millennial-scale cycles in the ~5-kyr band, and that our results, interpretation, and main conclusion remain unchanged.\nIn the revised manuscript, when discussing the possible origin of the ~1–1.7-kyr signals, we have revised the previous content as follows: The millennial-scale power observed between 1/1 kyr-1 and 1/1.7 kyr-1 does not correspond to known astronomical cycles or nonlinear combination tones (Fig. 3). Such high-frequency signals are often affected by bioturbation (Liu et al., 2021). Given the low SAR at DSDP–516F in the North Atlantic, our bioturbation simulations, based on typical sediment mixing depths of 5–10 cm in the ocean, suggest that these cycles in L* datasets may be entirely smoothed (Extended Data Fig. 8A). Thus, caution is warranted when interpreting these millennial cycles in the L* record of DSDP–516F. By contrast, the higher SAR and finer sampling resolution of the SLB records allow millennial cycles, including those at ~1.25 kyr to be clearly identified (Extended Data Fig. 8B). Furthermore, the presence of laminations through the study interval of SK2 (Extended Data Fig. 1&2) indicates minimal mixing by bioturbation, supporting the preservation of such cycles in the SLB, which may reflect a climate response to solar activities (Braun et al., 2005). Nevertheless, the spectral variances in the ~1–1.7 kyr bands is negligible in both SK2 and DSDP–516F, with most millennial-scale variance resided in the ~1.8–6-kyr band and attributable to the climate response to precession.\n2.3 I also have questions about what the chemical measurements are caused by. In the SOM, the authors suggest that their proxies (Rb/Sr and Ca/Ti) are controlled by calcite vs detrital silicate sediment. I understand that they used an XRF to scan the core, but have they calibrated to measurements to actual wet-chemical measurements? Do they know what sedimentary phases or minerals are giving them the high Ca, Sr, Ti signal?\nOur replay: Thank you for your comments. In the revised manuscript, we select box 59 from the SK2 core (1141.0–1141.7 m) and collected bulk samples for thin-section petrography, quantitative XRF (X-ray fluorescence spectrometry), LA-CP-MS (laser ablation inductively coupled plasma mass spectrometry), and XRD (X-ray diffraction) wet-chemical analyses. As shown in Extended Data Fig. 2, the semi-quantitative Ca, Sr, Si, Ti, Rb, Ca/Ti and Rb/Sr ratios obtained from XRF core scanning closely match the quantitative values derived from XRF and LA-CP-MS wet-chemical analysis.\nMineralogical evidence further supports our interpretation. The black mudstones intervals are characterized by high quartz content (confirmed by XRD and thin-section observations), elevated Si, Ti, Rb and Rb/Sr values, but negligible calcite, and correspondingly low Ca, Sr, and Ca/Ti values. In contrast, the gray laminated mudstones intervals are dominated by calcite, showing high Ca, Sr, and Ca/Ti values but low quartz, Si, Ti, Rb, and Rb/Sr values (Extended Data Fig. 2).\nThin-section observations under cross-polarized light also reveal that the calcite occurs in microcrystalline and micritic forms.\nThese results collectively demonstrate that Ca and Sr are mainly hosted in primary calcite, whereas Ti, Si, and Rb are associated with detrital silicates, thereby validating the use of Ca/Ti and Rb/Sr ratios as reliable proxies for carbonate versus siliciclastic input.\n2.4 Also can we be sure that this chemical variation is actually a climate proxy? Calcite need not be primary, in fact many terrestrial sedimentary strata have diagenetic calcite. Is there any evidence that the SLB Ca and Sr is reflecting the amount of primary calcite, rather than secondary or diagenetic calcite? Rocks of this age have often undergone secondary calcite precipitation if there are abundant porosity, which could mean that the apparent cyclicity is due to changes in grain size or porosity followed by post-depositional fluid flow and cementation, rather than a primary climate signal.\nOur replay: Thank you for your comments. Multiple independent lines of evidence support our interpretation that Ca and Sr in the SLB reflect the primary calcite abundance rather than a late diagenetic overprint. Specifically:\n(1) XRF results from both core scanning and wet-chemical measurements, together with XRD analyses, indicate that calcium is mainly hosted in calcite (Extended Data Fig. 2E);\n(2) Sr covaries with Ca and is elevated in calcite-rich intervals (Extended Data Fig. 2E–F), whereas recrystallized/secondary carbonates typically exhibit reduced Sr concentrations due to leaching during diagenesis;\n(3) Petrographic thin sections reveal that calcite occurs predominantly in microcrystalline and micritic forms (Extended Data Fig. 2B), characteristic of primary depositional textures rather than secondary cementation;\n(4) In the gray calcareous (calcite-rich) mudstone intervals, laminations are well developed (Extended Data Fig. 2A&B), with calcite occurring in parallel with the laminations as their primary constituent.\nTaken together, these observations indicate that the chemical variations, reflected by Ca/Ti and Rb/Sr ratios, record original depositional/climatic signals rather than being generated by later porosity-controlled cementation and fluid-flow cementation.\n2.5 It is also possible that these chemical cycles are due to autogenic sedimentary processes, which are internal to a sedimentary basin. Examples of this include river channel migration or delta lobe switching controlling sediment deposition within the system itself, rather than allogenic processes which are external to the depositional system, such as changes in climate. Even if this apparent cyclicity is similar to that occurring in other basins or sediment records, given the large range of cycle periodicity (see lines 96-101), it would be surprising not to find some kind of match. I would like to see some substantiation that what they are measuring are primary sediment phases, and that the geological and sedimentary setting precludes autogenic sedimentary processes.\nOur replay: Thank you for this this insightful comment. Several lines of evidence indicate that the millennial-scale cycles observed in the Songliao Basin are unlikely to have resulted from autogenic sedimentary processes:\n(1) Stable deep-lacustrine setting. The K n1+2 Member strata in the SK2 core were deposited 2 in a deep-lacustrine environment (Feng et al., 2010), where autogenic noise in the ‘source-to-sink system’ is expected to be minimal (Hajek and Straub, 2017). Besides, the K n1+2 Member consists 2 mainly of fine-grained mudstones, which are more capable of preserving allogenic signals than coarser deposits (Yang et al., 2024). Sedimentary facies analyses further indicate that the depositional environments were stable, the K n1+2 Member strata are laterally 2 continuous in the Songliao Basin, covering nearly the entire basin and forming a regional barrier to fluid flow (Feng et al., 2010), with no evidence of large-scale lateral facies shifts that would be expected from autogenic reorganization.\nFig. R5 | Lithology of the Second Member of the Nenjiang Formation in the SK1N borehole, and XRF-scanning ratios of log (Ca/Ti) and Rb/Sr. The passbands for the ~10.8-kyr and ~5.3-kyr filters are 0.0925±0.0175 cycles/kyr and 0.19±0.02 cycles/kyr, respectively. The vertical bars indicate correlations among the two datasets at ~10.8- or ~5.3-kyr timescales. It is worth noting that these datasets have not been published elsewhere and are being prepared for inclusion in a new manuscript.\n(2) Primary sedimentary signals. Petrographic (thin-section), mineralogical (XRD), and geochemical (XRF and LA-ICP-MS wet-chemical) analyses show that geochemical variations primarily reflect changes in primary sedimentary signals, reflecting alternating deposits between gray laminated calcareous mudstones and black siliceous mudstones (Extended Data Fig. 2; see also our response to Comments 2.3 and 2.4). The regular rhythmic alternation of these lithologies supports an allogenic, climatic-driven origin rather than an autogenic origin.\n(3) Phase-locking across proxies. Observed cycles in different proxies display phase-locking relationships on orbital and millennial timescales (Extended Data Fig. 7), whereas cycles caused by autogenic processes are generally irregular (quasi-periodic) and lack consistent phase relationships across proxies.\n(4) Spatial reproducibility. Similar 10 kyr, 5 kyr, and 2.5 kyr cycles are observed in the same stratigraphic interval of the Songliao Basin from the adjacent SK1 core (Fig. R5), ~145 km away from the SK2 core, indicating that these cyclicities are regionally reproducible rather than local autogenic features.\nTaken together, these lines of evidence support our interpretation that the measured cycles primarily reflect allogenic, climate-driven processes rather than autogenic sedimentary dynamics."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have carefully reviewed the revised version of the manuscript. The authors have clearly made a substantial effort to address all the points raised in the previous round of review. The changes implemented have improved the clarity and overall quality of the paper. I have no additional comments or concerns, and I consider the manuscript suitable for publication in its present form.", "reviewer_2": "I feel that the authors have addressed my concerns.\nI feel that the text they added clarifies my concerns (text in italics in Response to Reviewer).\nIt would be good if a the text that they wrote in reply within the Response to Reviewer concerning autogenic sedimentary processes and phase control of Rb/Sr and Ca/Ti by calcite vs detrital silicate could also appear in the SOM, just in case readers have similar questions."}, "author_response": "Response to Reviewer #1:\n(General comments) I have carefully reviewed the revised version of the manuscript. The authors have clearly made a substantial effort to address all the points raised in the previous round of review. The changes implemented have improved the clarity and overall quality of the paper. I have no additional comments or concerns, and I consider the manuscript suitable for publication in its present form.\nOur reply: Thank you for your generous comments and for your appreciation of our work.\nResponse to Reviewer #2:\n(General comments) I feel that the authors have addressed my concerns.\nI feel that the text they added clarifies my concerns (text in italics in Response to Reviewer).\nIt would be good if a the text that they wrote in reply within the Response to Reviewer concerning autogenic sedimentary processes and phase control of Rb/Sr and Ca/Ti by calcite vs detrital silicate could also appear in the SOM, just in case readers have similar questions.\nOur reply: Thank you for your positive comments and suggestions. In the revised manuscript, we added a Supplementary Note to the Supplementary Information (subheading: Allogenic origin of millennial-scale cycles in the Songliao Basin) to discussion the influence of autogenic sedimentary processes on our interpretation of millennial-scale cycles, as follows:\nAutogenic sedimentary processes, such as river channel migration or delta lobe switching, can produce quasi-periodical cycles (Li et al., 2016; Burgess et al., 2019). These autogenic cycles are internal to a sedimentary basin, thus do not contain climatic implications. However, several lines of evidence indicate that the millennial-scale cycles observed in the Songliao Basin are unlikely to have resulted from autogenic sedimentary processes:\n(1) Stable deep-lacustrine setting. The K n1+2 Member strata in the SK2 core were deposited 2 in a deep-lacustrine environment (Feng et al., 2010), where autogenic noise in the ‘source-to-sink system’ is expected to be minimal (Hajek and Straub, 2017). Besides, the K n1+2 Member consists 2 mainly of fine-grained mudstones, which are more capable of preserving allogenic signals than coarser deposits (Yang et al., 2024). Sedimentary facies analyses further indicate that the depositional environments were stable, the K n1+2 Member strata are laterally 2 continuous in the Songliao Basin, covering nearly the entire basin and forming a regional barrier to fluid flow (Feng et al., 2010), with no evidence of large-scale lateral facies shifts that would be expected from autogenic reorganization.\n(2) Primary sedimentary signals. Petrographic (thin-section), mineralogical (XRD), and geochemical (XRF and LA-ICP-MS wet-chemical) analyses show that geochemical variations primarily reflect changes in primary sedimentary signals, reflecting alternating deposits between gray laminated calcareous mudstones and black siliceous mudstones (Fig. S2). Moreover, the calcite in the intervals of calcareous mudstone occurs predominantly in microcrystalline and micritic forms (Fig. S2B), characteristic of primary depositional textures rather than secondary cementation. Thus, the regular rhythmic alternation of these lithologies, as reflected by the geochemical variations at millennial timescale, supports an allogenic, climatic-driven origin rather than an autogenic origin.\n(3) Phase-locking across proxies. Observed cycles in different proxies display phase-locking relationships on orbital and millennial timescales (Fig. S7), whereas cycles caused by autogenic processes are generally irregular (quasi-periodic) and lack consistent phase relationships across proxies.\nTaken together, these lines of evidence support our interpretation that the measured cycles primarily reflect allogenic, climate-driven processes rather than autogenic sedimentary dynamics."}]

Decision letter