A potential feedback loop underlying glacial-interglacial cycles

The sawtooth-patterned glacial-interglacial cycles in the Earth’s atmospheric temperature are a well-known, though poorly understood phenomenon. Pinpointing the relevant mechanisms behind these cycles will not only provide insights into past climate dynamics, but also help predict possible future responses of the Earth system to changing CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} levels. Previous work on this phenomenon suggests that the most important underlying mechanisms are interactions between marine biological production, ocean circulation, temperature and dust. So far, interaction directions (i.e., what causes what) have remained elusive. In this paper, we apply Convergent Cross-Mapping (CCM) to analyze paleoclimatic and paleoceanographic records to elucidate which mechanisms proposed in the literature play an important role in glacial-interglacial cycles, and to test the directionality of interactions. We find causal links between ocean ventilation, biological productivity, benthic δ18\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta ^{18}$$\end{document}O and dust, consistent with some but not all of the mechanisms proposed in the literature. Most importantly, we find evidence for a potential feedback loop from ocean ventilation to biological productivity to climate back to ocean ventilation. Here, we propose the hypothesis that this feedback loop of connected mechanisms could be the main driver for the glacial-interglacial cycles.


Introduction
The past 0.9 million years have been characterized by large cycles in global Earth temperature with a periodicity of about 100 ky. During warm (interglacial) periods, ice volume was small and CO 2 concentrations were high. During cold (glacial) time periods, ice volume was large and CO 2 concentrations were on average 90 ppm lower (Ruddiman 2001). Understanding the processes involved in these cycles can help to pinpoint the relevant processes in the carbon cycle. Given the rapid change that our climate system is undergoing today, knowledge on the existing feedbacks is of particular importance.
There is an ongoing debate about which mechanisms drive the reduction of CO 2 during glacials and which mechanisms drive the rapid increases in CO 2 during deglaciations (Menviel et al. 2018). Interactions with the ocean's carbon reservoir are an obvious first candidate for two main reasons: 1. The ocean is by far the largest carbon reservoir interacting with the atmosphere on the relevant timescale (Fasham 2003). 2. The 13 C of the ocean-atmosphere system was isotopically lighter at the Last Glacial Maximum than during the Holocene (especially in the deep ocean) (Eggleston et al. 2016;Peterson and Lisiecki 2018). Therefore, it is unlikely that the excess carbon was stored in the terrestrial biosphere or methane hydrates, which are both isotopically light carbon reservoirs (Zeebe and Wolf-Gladrow 2001).
It has been hypothesized that the excess carbon was stored in an isolated abyssal reservoir (Broecker and Barker 2007). At the glacial-interglacial transition, carbon from this reservoir would have been released to the atmosphere through upwelling (Broecker and Barker 2007;Marchitto et al. 2007). However, it has turned out to be difficult to locate a stagnant reservoir of sufficient size to account for the excess carbon (Broecker and Barker 2007;Broecker and Clark 2010). Another proposed mechanism that has received considerable attention is that changes in ocean circulation caused the increase in CO 2 (Siegenthaler and Wenk 1984) by bringing carbon from the deep to the sea surface (Anderson et al. 2009). Ocean circulation or ocean ventilation could have caused the rise of CO 2 in other ways as well. For example, meltwater pulses in the North Atlantic during the last deglaciation may have led to a temporary shutdown of the Atlantic Meridional Overturning Circulation (AMOC) and an associated increase in Antarctic Bottom Water (AABW) (Broecker 1998). According to a hypothesis by Toggweiler (1999), such a shift would in turn have led to outgassing of CO 2 from the ocean to the atmosphere due to the poor nutrient utilization subpolar Southern Ocean, where AABW is formed. Even so, modeling studies have been ambiguous about the impact of meltwater pulses on the ocean carbon cycle. In some simulations, the addition of meltwater in the North Atlantic led to carbon release from the ocean to the atmosphere (Schmittner et al. 2007;Bouttes et al. 2012;Matsumoto and Yokoyama 2013;Schmittner and Lund 2015), whereas it led to a net uptake of carbon by the ocean from the atmosphere in others (Obata 2007;Bozbiyik et al. 2011;Chikamoto et al. 2012). According to a set of model experiments by Menviel et al. (2014), the net effect depends strongly on the detailed salt budget. Other potential mechanisms for CO 2 fluctuations during glacial-interglacial cycles are based on biological activity. For a long time, the main hypothesis was that enhanced productivity in polar regions during glacial times drove down carbon from the atmosphere into the ocean (Sarmiento and Toggweiler 1984;Sigman and Boyle 2000). An early hypothesis by Martin (1990) states that enhanced biological productivity during glacial times was caused by iron deposition (or dust deposition). This idea is supported by experiments that show that current densities of phytoplankton (biological productivity) are limited by iron, for example within the Southern Ocean Iron RElease Experiment (SOIREE) project (Boyd and Law 2001). It has been suggested that changes in iron concentration can be amplified by local feedbacks, for example by increased ocean surface temperatures. These feedbacks could further increase the effect of iron fertilization on climate (Ridgwell 2002). However, according to the iron fertilization hypothesis, primary productivity in the subpolar Southern Ocean should be higher during glacial times, which is not supported by proxy records (Kohfeld and Chase 2011;Kohfeld et al. 2005Kohfeld et al. , 2013. Furthermore, general circulation models have shown that atmospheric CO 2 did not respond as much to increased productivity as suggested by earlier box models (Archer et al. 2000). A hypothesis by Omta et al. (2013) describes another role for biological productivity: it suggests that spikes in the densities of marine calcifiers may lead to a quick reduction in sea surface alkalinity which explains the observed rapid increase in temperature and CO 2 during deglaciations. Instead of extracting CO 2 from the sea surface and thus reducing atmospheric CO 2 and temperature, marine productivity increases atmospheric CO 2 according to this explanation.
It is generally assumed that the different mechanisms are linked to each other (Crucifix et al. 2017), and many of them are not mutually exclusive. Therefore, most modeling studies use a combination of proposed mechanisms (Brovkin et al. 2007). Furthermore, many geological variables are synchronized with the Milankovitch forcing (the Summer insolation at subpolar latitudes in the Northern Hemisphere), which makes it hard to distinguish between links between the variables themselves or common links with the forcing (Daruka and Ditlevsen 2016). What makes the topic even more complicated is the cyclic behaviour that is observed, which may be linked to positive or negative feedback loops involving different elements (Lenton et al. 2008). One example is the ice-albedo effect, where higher temperatures lead to reduced ice caps which leads to a reduction in the albedo, which in turn increases the temperature. Thus, two processes reinforcing each other can give a strong positive feedback loop. The addition of more processes may weaken the feedback, particularly if different scales are involved. This suggests that a potential feedback loop underlying glacial-interglacial cycles would likely be dominated by a few key causal links.
It is still unclear what these few dominant causal links would be, although many different potential causal interactions have been described (see Table 1 and Fig. 1). Here we take a statistical approach to investigate whether existing data suggest any particularly strong causal links. The recent increase in available records that go back hundreds of thousands of years and that have improved temporal resolution and quality allows for such an approach. A suitable method for unravelling causal links from nonlinear time series is Convergent Cross-Mapping (CCM) (Sugihara et al. 2012). It has been applied in ecological systems to find causal links between temperature and anchovy and sardine abundance (Sugihara et al. 2012), in physiology to distinguish between normal (healthy) and impaired cerebral autoregulation (Heskamp et al. 2014), in social systems for predicting the behaviour of users of social media (Luo et al. 2014) and in the climate system for detecting directionality in the relation between temperature and CO 2 (Van Nes et al. 2015). In this study, we apply CCM to published marine and ice-core records spanning the past 800,000 years in order to shed light on the (directions of) interactions between climate, ocean circulation, biological productivity and dust deposition (Fig. 2).

Data
We use benthic 18 O from the LR04 stack of 57 globally distributed sites (Lisiecki and Raymo 2005). Benthic 18 O is a function of the temperature and the 18 O of the surrounding water, which is influenced by the global ice volume. Therefore,

18
O is considered a proxy for temperature and ice volume. In the Supplementary Materials, we also analyze the role of atmospheric CO 2 (as opposed to 18 O) to see if a climate proxy that is more related to alkalinity gives similar results. We use atmospheric CO 2 data from the European Project for Ice Coring in Antarctica [EPICA (Lüthi et al. 2008)]. This CO 2 record can likely also be used as a proxy for ocean alkalinity (ALK) for the following reason. [HCO − 3 ] and [CO 2− 3 ] account for more than 95% of ocean alkalinity (ALK) and 99% of dissolved inorganic carbon (DIC) (Sarmiento and Gruber 2006; Williams and Follows 2011) and therefore: Thus, ALK ≈ DIC + [CO 2− 3 ]. Whole-ocean DIC is likely higher during glacials than during interglacials, as excess atmospheric carbon is stored in the ocean. Furthermore, both the lysocline depth and the B/Ca proxy indicate relatively small glacial-interglacial changes in deep-ocean [CO 2− 3 ] (Catubig et al. 1998;Yu et al. 2010;Raitzsch et al. 2011;Yu et al. 2014). Hence, ocean alkalinity probably increases from an interglacial to a glacial period and decreases from a glacial to an interglacial.
Ocean circulation cannot be measured directly, but the 13 C of benthic foraminifera is a proxy for ocean ventilation and ocean circulation (Curry and Oppo 2005), because it decreases in "older" water (measured from the time it was last in contact with the atmosphere) due to the remineralization of organic material. We use the Lisiecki (2010a) benthic foraminifera 13 C data compilation.
We use a Ba/Fe record from the Southern Ocean as a proxy for biological productivity . It has been shown that the vertical flux of marine Ba is a viable proxy for productivity (Paytan and Griffith 2007) and that the Ba/Fe in turn reflects the vertical flux of marine Ba (Studer et al. 2015).
Lastly, we use Antarctic dust data from the EPICA dome C ice core (Lambert et al. 2008).

Fig. 1
A schematic depiction of the links described in the literature between 13 C , biological productivity, 18 O and dust, indicating that all links are potentially possible. Mechanisms are described in Table 1 1 3 The data are available at resolutions of 0.2 ky (Ba/Fe), 1 ky ( 18 O), 0.015 ky (CO 2 ), 2 ky ( 13 C), and varying steps with a mean of 0.6942 ky (dust). The data dates back to the last 1557 ky (Ba/Fe), 5320 ky ( 18 O), 798 ky (CO 2 ), 802 ky ( 13 C), and 800 ky (dust). Therefore, we only use data up to 798 ky before present and we use the data at intervals of 2 ky (linear interpolation if necessary) for all interactions with 13 C and with intervals of 1 ky (linear interpolation if necessary) for all other interactions.

CCM
For data analysis we use convergent cross-mapping (CCM) (Sugihara et al. 2012). The method is based on Takens' Theorem (Takens 1981) that shows that dynamical systems involving more than one variable can be reconstructed from time-lagged series of only one of those involved variables, as long as sufficient time-lagged states are included: the socalled embedding dimension. In CCM, time-lagged time series of two variables are analysed to see if time-lagged states of those two variables can predict each other's current state. If time series of X can be predicted based on the timelagged time series Y, it is concluded that X causes Y. The prediction is tested by calculating the CCM skill, which is the correlation between predicted values and true values. The method also tests for convergence, which means that the predictions improve with the length L of the time series that is used to predict X from Y. We calculate the 'CCM skill' for X → Y (X causes Y) with a slowly increasing length of data. We start with a library size (the number of data points used for the prediction) of 15, and then increase with steps of 15, to end with a library size of 270. More information about the CCM method can be found in the Supplementary Materials or in Sugihara et al. (2012).
The significance of causal links is tested by generating 100 surrogate data sets by cutting the time series at a random point and then exchanging the first and second part. In that way, the surrogate data have exactly the same characteristics as the original data in terms of power spectrum, autocorrelation and variance, but the linkage between the two data sets is broken. This method has been used previously by Van Nes et al. (2015). To test the robustness of our results, we repeat the analysis with an alternative way of surrogate data generation by Ebisuzaki (1997).
If causal links are found in two directions, one cannot distinguish between the following options of having 1/ truly bidirectional causal links, 2/ synchronization with a third variable, or 3/ a one-way strong interaction where one variable behaves as a 'follower' of the other. One possible method to distinguish between these three options is to perform CCM analysis on time lags of the variables by shifting the two time series relative to each other . If one link has a clear optimum at a negative time lag, this means that past values from one variable can better predict future values from another variable, which is what you expect from a true causal link. If however the optimum is found at a positive time lag (in the future), this suggests that the link might not be causal but a mere artifact from a strong coupling between the variables . If one link has a negative time lag and one link has a zero time lag, we cannot distinguish between a bidirectional link and a strong unidirectional link. To be conservative and because we are not sure about the link with the time lag of zero, we will call these unidirectional links where we exclude the link with a time lag higher than or equal to zero. Lastly, it is important to note that the method cannot distinguish between positive and negative interactions.

Choice of parameters
Using the method by Cao (1997) (for details, see Supplementary Materials) for finding an optimal embedding dimension, we find a dimension of 10. The choice of parameters specific to a method is often a topic for debate. Therefore, we also calculate the sensitivity to the choices we made, and we find that a deviation from E = 10 does not significantly affect our results (Supplementary Materials Fig. 6). Subsequently, we tried an alternative method from Sugihara (1994) for finding E and that also yields a maximum embedding dimension of E = 10 (Supplementary Materials Fig. 2-4).
A time lag = 1 is used. Since climate time series are not over-sampled, it seems reasonable to take the time lag as small as possible. Increasing can cause a loss of the signal, which is to be expected for the data resolution we have ( Supplementary Materials Fig. 7).

Results
The causal links found in our analysis are summarized in Fig. 3. Importantly, these links include a causal loop, consisting of the link from 13 C to Ba/Fe to 18 O back to 13 C in both the Pacific and East-Atlantic ocean. Within this loop, bidirectional significant links are found between Ba/Fe (proxy for productivity) and 18 O (climate proxy) ( Fig. 4a). As these interactions are strong, this could either indicate truly bidirectional causal links, or a coupling that is so strong that one variable becomes a 'slave' of the other. The time-delayed results from convergent crossmapping shows that the link from Ba to 18 O has a clear negative time lag (Fig. 5b) whereas the other link has a zero time lag (Fig. 5a). This suggests that productivity (Ba/Fe) affects climate ( 18 O) with a time lag of 4 ky whereas 18 O either causes Ba/Fe as well, or behaves as a "slave" of Ba/Fe.
As discussed in the Introduction, glacial-interglacial changes in ocean productivity are thought to be caused by (i) changes in the ocean circulation (for which deep-ocean 13 C is a proxy) or (ii) variations in iron input (for which dust from Antarctic ice cores is a proxy).
In case (i), we expect a causal link from deep Pacific and/ or deep Atlantic 13 C to Ba/Fe. The results for the analysis with 13 C show that indeed there are links from deep Pacific and deep East Atlantic 13 C to and from 18 O (Fig. 4c, g). Furthermore, there are links from deep Pacific and deep East Atlantic 13 C to productivity (Fig. 4b, f). The time-lag analysis shows that all the links have an optimum for a negative time lag (Fig. 5d-f, l-n), indicating causal links with a time lag. These results are consistent with the hypothesis that changes in ocean circulation could have had an effect on glacial-interglacial oscillations. Furthermore, the closed loop suggests a potential mechanism for nonlinear dynamics, which could explain the observed cycles (see Discussion). Note that in the North-East Atlantic, the pattern is different from that in the Pacific and East Atlantic: no significant causal links are found between the deep North-East Atlantic 13 C, Ba/Fe and CO 2 (Fig. 4d, e). In case (ii), we expect a causal link from dust to Ba/Fe. Our results are consistent with this hypothesis too, but only East-Atlantic 13 C (EA), biological productivity (Ba), 18 O and dust, according to our CCM analysis. Arrows indicate causal links, i.e., X → Y means that X causes Y via a confounded route. Even though the link from dust to Ba/Fe is significant (Fig. 4h), our time lags analysis shows that the optimum time lag for this interaction is positive, suggesting that this link is not a causal link but an artifact of the coupling between Ba/Fe and dust. However, we find causal links between dust and 18 O (Fig. 4i) and as mentioned before from 18 O to 13 C and from 13 C to Ba/Fe. Time lags analysis suggest that these are truly causal links (Fig. 5q, r, e, m, d, l).
An alternative way of surrogate data generation to test for significance does not affect our results (Supplementary Materials Figs. 8-9). Furthermore, replacing 18 O for CO 2 , as an alternative climate and temperature proxy that is closely linked to ocean alkalinity, only alters one link in the results: the link from CO 2 to Ba also has an optimum for a negative time lag, suggesting a causal link with a time lag of 4 kyr (Supplementary Materials Figs. 10-11).

Discussion
For the first time, we have applied CCM to identify causal relationships between ice-core and paleoceanographic proxy records. Our results suggest the existence of several closed feedback loops. These feedback loops could provide   (Fig. 1) that we used for our analysis does not suggest particularly high productivity during times of low atmospheric CO 2 . Rather, the highest productivity appears to be associated with the glacial-interglacial transitions. Therefore, we think that a causal link between productivity and atmospheric CO 2 as described in Omta et al. (2018) would be more consistent with this record. That is, a spike in burial of organic matter leads to a decrease in deep-ocean alkalinity and the eventual outgassing of CO 2 from the ocean to the atmosphere due to carbonate compensation. This appears consistent with the recently reported finding that the variation in productivity in the Southern Ocean became larger as the amplitude of the glacialinterglacial cycles increased across the Mid-Pleistocene Transition ( Fig. 3E in Hasenfratz et al. 2019). In other words, small-amplitude glacial cycles are associated with relatively constant levels of productivity, whereas large-amplitude glacial cycles are associated with large periodic spikes in productivity. 3. 18 O and CO 2 may impact the overturning circulation in various ways. For example, it has been suggested that a more northerly position of the Southern Hemisphere westerly winds during glacial periods would lead to weaker upwelling around Antarctica, which would in turn weaken the Southern cell of the overturning circulation (Toggweiler et al. 2006). Alternatively, a northward shift of the regions of upwelling and deep-water formation could be buoyancy-driven through changes in the Southern Ocean sea-ice extent (Ferrari et al. 2014;Jansen 2017). Regardless of the precise mechanism, it appears that the structure of the overturning circulation varies significantly between glacial and interglacial times (Curry and Oppo 2005;Lynch-Stieglitz et al. 2007;Burke et al. 2015).
Beside this three-variable loop, CCM detects two loops consisting of two variables, namely between 18 O and Pacific and East-Atlantic 13 C and between 18 O and dust. However, we want to interpret these loops a bit more cautiously than the others, because CCM cannot always distinguish between a strong unidirectional link and a bidirectional link. Even so, these two loops can also be considered candidate potential drivers for the nonlinear behaviour observed in the glacialinterglacial dynamics. In particular, lower sea levels during glacial times led to larger areas of the continental shelf off Patagonia being exposed (Iriondo 2000). Since most Antarctic dust seems to have a Patagonian origin (Fischer et al. 2007), changes in climate cause changes in dust deposition. Furthermore, it has been suggested that the larger equatorto-pole surface temperature gradient during glacial times led to more storms in the Southern Hemisphere, which in turn led to more dust transport (McGee et al. 2010). The link from dust to 18 O may be explained from the radiative impacts of airborne dust (Tegen 2003;Yoon et al. 2005;Rosenfeld et al. 2006) and the impact of dust deposited in Antarctica on the snow albedo (Bar-Or et al. 2008).
In addition to our standard CCM analysis, we also performed the same analysis for different shifts in the data to see if the optimum causal link was pointing in the right direction: from the past to the future. That is, a true causal link implies a negative optimum time lag. For the link from 18 O to Ba/Fe, we found an optimum at a zero time lag and  Morris and Farrell (1971) At lower temperatures assimilation of CO 2 happens more rapidly. Temperature has a non-linear effect on most enzyme activity Productivity → Dust Martin (1990) Uptake of iron by plankton decreases iron availability in sea surface Dust → Productivity Martin (1990) Increased iron in ocean leads to higher biological productivity which is iron-limited 18 O → Dust Ridgwell (2002) A colder dryer glacial climate has a less vigorous hydrological cycle leading to higher dust deposition rates Tegen (2003) Absorption or reflection of incoming solar radiation affects temperature which is linked to 18 O the link from dust to Ba/Fe had a positive time lag. To be conservative, we therefore interpret these interactions as unidirectional links from Ba/Fe to 18 O and to dust. We want to point out that with the current quality of the data, the time lags analysis should not be used to interpret the precise time scale of the interactions we find. However, it can be used as an additional test to see if a link could be causal (i.e., past values can only cause the future and not the other way around). It should be noted that, like any statistical method, the absence of a significant causal link in our analysis does not prove that this link is absent in the real world. It is possible that, as longer and more accurate data become available, links will be revealed that currently remain undetected. The unidirectional link from Ba/Fe to dust suggested by our analysis may not reflect a direct causal link but may rather be the consequence of the transportation route of the dust. From Patagonia, the dust first travels over the Southern Ocean where it may enhance productivity. Subsequently, some of the remaining dust is transported further and ultimately deposited on the Antarctic ice sheet (See Supplementary Figure 1 for a location of these sites). As a result of this chronological order, the signal measured by the productivity proxy Ba/Fe will be present in the Antarctic dust signal (and therefore we find that productivity causes dust). However, elements of the Antarctic dust signal that developed in the later stages of the dust's journey will not be visible in the Ba/ Fe data (and therefore we find a positive optimal time lag for the causal link from dust to Ba/Fe).
An interesting implication of the existence of feedback loops is the possibility for nonlinear dynamics in the oceanclimate system. Feedback loops can be negative, in which case they work stabilizing or cause oscillatory behaviour if there are long time lags in the interactions (Levins 1974). Feedback loops can also be positive, in which case they can cause self-amplifying dynamics, but time lags in positive feedback loops can either weaken or amplify the destabilizing effect. Based on our results, one cannot determine whether the identified dominant feedback loop is stabilizing or not, and under which conditions. There is no reason to believe that the interactions we investigate here are linear interactions, and therefore the loops can vary in their response under different circumstances. Furthermore, the present study cannot prove the existence of positive feedback loops, both because CCM does not provide information about the way the effect-variable responds to the cause variable (i.e., positive or negative) and because of large uncertainties in the precise time lags of the analysis. However, as one of our objectives of this study is to find an explanation for the nonlinear glacial-interglacial cycles, we believe that these loops are worth investigating.
Positive feedback in a dynamic system can generate alternative system states that reinforce themselves. Under certain conditions the system may 'flip' from one state to another -the conditions under which this occurs are generally referred to as 'tipping points'. The model hypothesized by Omta et al. (2013) involves calcifiers in a feedback loop, and in a numerical study it was shown that 1/ this model can generate the glacial-interglacial sawtooth pattern under Milankovitch forcing of 20 ky, and 2/ that the period of the sawtooth can vary under minor perturbations with the same forcing (Omta et al. 2016), effectively suggesting the possibility of a tipping point. This can explain the glacial interglacial cycles of the past, but it can also guide the predictions for current climate dynamics. Currently, there is a significant perturbation in the form of a considerable increase of atmospheric CO 2 . The existence of a feedback loop as found in the data could indicate a potential change in the future from the sawtooth pattern observed in the data over the last 800,000 years to some unknown behaviour.
As we mentioned in the Introduction, various models have been formulated to explain glacial-interglacial dynamics. Crucifix (2012) reviewed many such models, based on a range of different mechanisms. It turned out to be difficult to decide between these models, since they all generate a sawtooth oscillation in temperature/ice volume and/or CO 2 . Therefore, Crucifix (2012) argued, we need more stringent criteria to operate our model selection. In Omta et al. (2016), one such criterion was formulated: models will need to reproduce the observed linear proportionality between the length and amplitude of the glacial-interglacial cycles. The causal relationships that emerged from our analysis provide further criteria. More specifically, model output could be analysed with CCM to test whether the same causal relationships emerge as from the analysis of the proxy records. However, it should be kept in mind that interpreting these links requires knowledge of the system under study (see for example for our link from productivity to dust).
Quite remarkably, our analysis suggests that mid-depth Northeast Atlantic 13 C (Lisiecki 2010a) has no significant causal relations with any of the other measured records. The locations of the cores that form the basis of this record are likely within the Northern overturning cell (North Atlantic Deep Water and Glacial North Atlantic Intermediate Water) during both glacial and interglacial times. Therefore, one possible interpretation would be that the strength of the Northern cell has no impact on, and is not impacted by, the strength of the Southern overturning and climate in general. However, ocean models generally indicate that a stronger and deeper Northern cell is associated with a weaker and shallower Southern cell and vice versa (Shin et al. 2003;Liu et al. 2005;Jansen 2017). Another possible interpretation would be that variations in the mid-depth Northeast Atlantic 13 C primarily reflect atmospheric changes, since the water at these sites is so 'young'. Thus, our analysis would suggest that atmospheric 13 C does not have a strong causal relationship either way with the overturning circulation, climate or atmospheric CO 2 . This appears in contradiction with the generally accepted view that the release of isotopically light 'old' carbon from the abyss causes both the observed minimum in atmospheric 13 C and the rise in atmospheric CO 2 during the deglaciation (Schmitt et al. 2012;Broecker and McGee 2013). Lastly, we show that the deep Northeast Atlantic 13 C data is not clearly distinguishable from noise according to one of the tests we did ( Supplementary Fig. 5), making it hard to pick up signals with CCM. It is possible that if data of higher resolution becomes available in the future, signals might become clear that are not significant when using the currently available data sets.
The Ba/Fe record is from a single location in the Southern Ocean, whereas the other records represent either global signals or averages for oceanic regions. Nevertheless, the deglacial marine productivity maxima appear to be a global phenomenon. Deglacial maxima in productivity proxies such as opal accumulation and biogenic Ba have been found at various sites throughout the Atlantic (Kasten et al. 2001;Romero et al. 2008;Gil et al. 2009;Meckler et al. 2013), the Pacific (Crusius et al. 2004;Jaccard et al. 2005;Galbraith et al. 2007;Jaccard et al. 2010;Kohfeld and Chase 2011;Hayes et al. 2011), and the Southern Ocean (Anderson et al. 2009). A global compilation of marine sediment proxies reveals an increase in suboxic conditions in Oxygen Minimum Zones (OMZ) during the deglaciation (Jaccard and Galbraith 2012) which could result from enhanced productivity. Furthermore, particularly strong deglacial maxima in 15 N have been recorded in OMZ (Pride et al. 1999;Thunell and Kepple 2004;Deutsch et al. 2004), a possible indication of enhanced denitrification due to the low oxygen conditions. Paleoceanographic records are often dated by overlapping them with other records, for which the dating ultimately relies on orbital tuning. For example, the 13 C records Lisiecki (2010a) were dated by aligning their 18 O records to the LR04 benthic 18 O stack (Lisiecki and Raymo 2005), which was in turn orbitally tuned. It has been pointed out that this may cause spurious correlations with orbital cycles (Huybers and Wunsch 2004). However, the Fourier spectra of records using a depth-derived age model and an orbitally tuned age model appear rather similar ( Fig. 9 in Huybers and Wunsch (2004)). Perhaps more importantly, any shifting of the time series due to tuning should only affect the estimated lag between records and not the direction of the causality. As CCM does not primarily establish causality based on time lags but rather based on the structure of the data, we therefore believe the causalities found with the method are true causalities and not the result of artifacts.

Conclusion
Our CCM analysis on paleoceanographic and ice-core records suggests a dominant loop from ocean ventilation to biological productivity to climate and then back to ocean ventilation. While all separate mechanisms have been suggested earlier as explanations of the sawtooth-patterned glacial-interglacial dynamics, the existence of the dominant loop as a whole has to our knowledge not been identified from any data before. These loops provide possible mechanisms that explain why the climate system dynamics might not react in a linear way to perturbations in elements of this loop. The current considerable increase in atmospheric carbon dioxide could present a perturbation that results in a shift of the climate system dynamics from the glacial-interglacial cycle to new (unpredictable) behaviour.
To elucidate the ramifications of feedback mechanisms in the climate system further, numerical models could be used to quantify these mechanisms and simulate possible glacial-interglacial dynamics as well as possible transitions between dynamical regimes. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.