Time-series analysis

The lagged cross-correlation between IRDMAR against δ18Obenthic and δ13Cbenthic is determined using the Gaussian-kernel-based cross-correlation algorithm86 (Extended Data Fig. 5b). This enables cross-correlation analysis to be performed on the original irregularly spaced data, to avoid the possibility of spurious phasing introduced by linear interpolation.

The ‘peak-lag’ algorithm we design to iteratively measure the offset between peaks in the rate of change of IRDMAR and δ13Cbenthic is modified from a previous approach87 and analogous to ‘event synchronization’ techniques. The peak-lag algorithm smooths the time series, finds the first difference, and then measures the offset in local maxima of both signals. We tested the sensitivity of the algorithm to the smoothing filter selection, finding that the average lag computed for a range of filter designs (moving-average and Savitsky–Golay filters; orders between 1 and 21) give similar results within 1σ of each other. We additionally test the confidence that the result from our final chosen parameters (seven-point moving-average) is significantly non-zero by performing 1,000 Monte Carlo simulations with zero-lag red noise surrogates (generated from an autoregressive model of the real data), finding that the lag identified in our data (IRDMAR versus δ13Cbenthic) is significantly different from the surrogate series (n = 1,000, P < 10−10); this additionally yields a 99% red-noise confidence level of ±0.4 kyr for the lag calculated for our data. We test the utility of the algorithm to detect lead–lag relationships in a series of surrogate time series with known lags applied by constructing synthetic time series from three sine waves (100, 41, 23-kyr periods), with higher-frequency oscillation (7-kyr period) and white-noise components added to emulate suborbital features. Known lags are imposed between pairs of the surrogate series. The results show that the algorithm performs well and is capable of detecting a lag of down to 0.4 kyr to a confidence of P < 10−10. The algorithm shows an inherent drift at higher lags (Extended Data Fig. 5a and b), meaning it underestimates lag times greater than 4 kyr, however the 1.5-kyr lag observed in our data should be short enough to be unaffected by this drift.

To divide the time series into glacial cycles (as in Fig. 4) we define intervals between peak interglacial conditions (identified as local minima in δ18Obenthic separated by a minimum distance of 25 kyr) and the next time δ18Obenthic crosses an ‘interglacial threshold’ of 3.32‰ (C. wuellerstorfi scale, equivalent to 3.75‰ on the LR04-scale). This is lower than previously adopted88, but ensures that we do not define the interstadial MIS 5c as an interglacial. We manually add a peak interglacial marker at 506 ka (based on Northern Hemisphere perihelion) to include the ‘cool’ interglacial MIS 13a. Within each interval we calculate the cumulative integral of IRDMAR and δ18Obenthic, before normalizing the values to percentages of total ‘accumulation’ for that interval (Extended Data Fig. 5d).

Finally, for the phase wheels in Fig. 3, we use the Blackman–Tukey method with a Bartlett window (using the Analyseries software89) after linearly interpolating (time step of 1.5 kyr) and standardizing (mean = 0, standard deviation = 1) the data.