Millennial-scale atmospheric CO 2 variations during the Marine Isotope Stage 6 period (190–135 ka)

. Using new and previously published CO 2 data from the EPICA Dome C ice core (EDC), we reconstruct a new high-resolution record of atmospheric CO 2 during Marine Isotope Stage (MIS) 6 (190 to 135 ka) the penultimate glacial period. Similar to the


Abstract.
Using new and previously published CO 2 data from the EPICA Dome C ice core (EDC), we reconstruct a new high-resolution record of atmospheric CO 2 during Marine Isotope Stage (MIS) 6 (190 to 135 ka) the penultimate glacial period. Similar to the last glacial cycle, where highresolution data already exists, our record shows that during longer North Atlantic (NA) stadials, millennial CO 2 variations during MIS 6 are clearly coincident with the bipolar seesaw signal in the Antarctic temperature record. However, during one short stadial in the NA, atmospheric CO 2 variation is small (∼ 5 ppm) and the relationship between temperature variations in EDC and atmospheric CO 2 is unclear. The magnitude of CO 2 increase during Carbon Dioxide Maxima (CDM) is closely related to the NA stadial duration in both MIS 6 and MIS 3 (60-27 ka). This observation implies that during the last two glacials the overall bipolar seesaw coupling of climate and atmospheric CO 2 operated similarly. In addition, similar to the last glacial period, CDM during the earliest MIS 6 show different lags with respect to the corresponding abrupt CH 4 rises, the latter reflecting rapid warming in the Northern Hemisphere (NH). During MIS 6i at around 181.5±0.3 ka, CDM 6i lags the abrupt warming in the NH by only 240±320 years. However, during CDM 6iv

Introduction
Ice core studies allow us to considerably extend our knowledge about natural climate-carbon cycle feedbacks by directly reconstructing atmospheric CO 2 from gas preserved in Antarctic ice sheets (Lüthi et al., 2008;Petit et al., 1999). Comparing atmospheric CO 2 records from Antarctic ice cores with proxies of paleoclimate helps us to understand how atmospheric CO 2 was controlled by carbon exchange with the ocean and land reservoirs on orbital to centennial timescales Brook, 2008, 2014;Bereiter et al., 2012;Higgins et al., 2015;Lüthi et al., 2008;Marcott et al., 2014;Petit et al., 1999).
Published by Copernicus Publications on behalf of the European Geosciences Union. source: https://doi.org/10.7892/boris.148172 | downloaded: 19.11.2020 Previous work on polar ice core records revealed that temperature variations in Greenland and Antarctica on millennial timescales appear to be a pervasive feature during the last glacial period. While Antarctic temperature varied gradually, Greenland temperature changes occurred abruptly. A phase difference can be observed between millennial-scale variations of temperature in the NH and SH (Northern Hemisphere and Southern Hemisphere, respectively), which is referred to as the bipolar seesaw phenomenon (Blunier and Brook, 2001;Pedro et al., 2018;Stocker and Johnsen, 2003). Potential triggers for this climatic variability on the millennial scale are fresh water perturbation in the North Atlantic (NA) or alterations of sea ice extent, surface temperature and salinity in the NA Broecker et al., 1992;Heinrich, 1988;McManus et al., 1998), which may reduce the strength of the Atlantic Meridional Overturning Circulation (AMOC). This would cause a reduction in heat transport from the SH to the NH, which leads to an abrupt cooling in the NA region and a gradual warming in the SH (Stocker and Johnsen, 2003) and the opposite behavior when AMOC is strengthened.
Existing CO 2 records show the presence of millennialscale oscillations on the order of 20 ppm over the last glacial period Brook, 2008, 2014;Bereiter et al., 2012), which generally co-vary with the major water isotope (δD) variations in Antarctic ice cores reflecting Antarctic temperature variations (Ahn and Brook, 2008;Bereiter et al., 2012) (Fig. 1). During cold periods in the NA, referred to as NA stadials, atmospheric CO 2 increased continuously and in parallel to Antarctic temperature increase. Again, at the onset of warming in Greenland, atmospheric CO 2 started to decrease (Ahn and Brook, 2008;Bereiter et al., 2012), generally in line with a co-occurring, slow Antarctic cooling. However, the CO 2 decrease did not always start at exactly the same time as the onset of the abrupt warming in the NA, and the lag itself varied. For example, during Marine Isotope Stage (MIS) 3, Carbon Dioxide Maxima (CDM) lagged behind abrupt temperature change in Greenland by 870±90 years. During MIS 5, the lag of CDM with respect to abrupt temperature warming in the NH was only about 250±190 years (Bereiter et al., 2012). Atmospheric CO 2 variations on millennial scales are thought to be related to the role of the Southern Ocean in carbon uptake and deep-ocean ventilation on millennial timescales (Fischer et al., 2010;Marcott et al., 2014;Schmitt et al., 2012;Sigman and Boyle, 2000;Toggweiler et al., 2006). In addition, atmospheric CO 2 can be affected by changes in the AMOC, which affects the ventilation of carbon from the deep ocean (Denton et al., 2010;Sigman et al., 2007). However, the mechanisms responsible for these oscillations are still a matter of debate (Bouttes et al., 2012;Bozbiyik et al., 2011;Gottschalk et al., 2019;Menviel et al., 2014).
Comparing CO 2 changes on millennial timescales during the past two glacial periods, MIS 3 (60-27 ka) and early MIS 6 (185-160 ka), can provide us with a better understand-ing of the carbon cycle, due to the similarities but also differences of climate conditions and events during the last two glacial periods (Fig. 1). Proxy evidence indicates that the states of several important components of the climate-carbon cycle were not the same between MIS 3 and MIS 6. Sea ice cover in the South Atlantic was more extensive in MIS 6, and sea surface temperature in the South Atlantic is thought to have been lower (Gottschalk et al., 2020). The bipolar seesaw phenomenon has also been observed during the early MIS 6 period (Cheng et al., 2016;Jouzel et al., 2007;Margari et al., 2010). However, the bipolar seesaw events during MIS 6 are longer than those found during MIS 3. Events of massive iceberg discharge into the NA, which are thought to have driven millennial-scale changes in the meridional overturning circulation during MIS 3 (de Abreu et al., 2003;Mc-Manus et al., 1999), appear to be much more frequent during MIS 3 than during MIS 6. During the early MIS 6, iceberg discharge was muted, and during the time period around 175 ka summer insolation levels in the NH approached interglacial values (Berger, 1978). Due to the stronger NH summer insolation, the Intertropical Convergence Zone (ITCZ) had shifted to the north, which intensified monsoon systems in low-latitude regions, such as in Asia, the Apennine Peninsula and the Levant (Ayalon et al., 2002;Bard et al., 2002;Cheng et al., 2016). This may have led to a weaker overturning circulation due to the reduction of the density of the North Atlantic surface water, making the AMOC cell shallower during MIS 6 than during MIS 3 (Gottschalk et al., 2020;Margari et al., 2010).
In order to investigate whether the different climate conditions between MIS 6 and MIS 3 could have impacted the relationship between atmospheric CO 2 and climate, we measured 177 new data points of atmospheric CO 2 concentrations from the EPICA Dome C (EDC) ice core (75 • 06 S, 123 • 24 E) spanning MIS 6 (190 to 135 ka). This new CO 2 data measured within the scope of this study is complemented with two published datasets (Lourantou et al., 2010;Schneider et al., 2013) and one unpublished CO 2 dataset measured in 2003. We use all of the available data from the EDC ice core to compile a composite dataset of atmospheric CO 2 covering MIS 6. In this work, we significantly improve existing records previously obtained from the Vostok ice core (Petit et al., 1999). The new composite dataset from the EDC core provides a temporal resolution of ∼ 230 years during MIS 6, as compared to 1000 years in the Vostok record. We also improved the relative age uncertainties between ice and gas in the EDC core using δ 15 N-based estimates of firn column thickness to better constrain leads and lags between peaks in δD (Antarctic Isotope Maxima, or AIM) in the EDC ice core and atmospheric CO 2 concentrations during the early MIS 6, thus establishing a new chronology using new δ 15 N data during the early MIS 6 in this study and published data from Landais et al. (2013). Finally, we improve the temporal resolution of existing CH 4 data from EDC (Loulergue et al., 2008) from 600 to ∼ 350 years to be  (de Abreu et al., 2003).
The 21 June insolation for 65 • N (Berger, 1978). (c) The δ 18 O calcite from Sanbao cave, indicative of the strength of the East Asian monsoon (Cheng et al., 2016). (d) Dust flux in the EDC ice core (Lambert et al., 2012). (e) Published atmospheric CH 4 in EDC (dark green dots) (Loulergue et al., 2008) and a composite atmospheric CH 4 record from EDC derived in this study (light yellow dots). (f) Composite CO 2 record from the EDC ice core derived in this study (orange dots) and a published composite CO 2 record from Antarctic ice cores (dark blue dots) (Bereiter et al., 2015). (g) The δD in the EDC ice core, Antarctica . Vertical grey bars indicate the timing of Heinrich events during the last glacial period. The event numbers are indicated at the bottom. The marine isotope substage numbers are written at the bottom (Railsback et al., 2015). able to more precisely calculate the shift of CDM relative to the rapid climate change in the NH. To avoid the age uncertainties between proxy data and atmospheric CO 2 data, CH 4 measurements are used as a time marker of rapid warming in the NH, as over the last glacial period CH 4 and Greenland temperatures are assumed to be synchronous, with a lag of CH 4 of less than ∼ 50 years on average (Baumgartner et al., 2014;Rosen et al., 2014).

CO 2 measurements
In this study, we measured atmospheric CO 2 from 177 depth intervals of the EDC ice core at the Institut des Géosciences de l'Environnement (IGE), France and Climate and Environmental Physics (CEP), Physics Institute, University of Bern, Switzerland. The majority of all atmospheric CO 2 samples (i.e. 150) were measured using the ball mill dry-extraction system coupled to a gas chromatograph at IGE (Schaefer et al., 2011). Each ball mill data point presented in this study corresponds to a single 40 g ice sample, which was measured five times by gas chromatography (five consecutive injections of the same extracted gas). Approximately 5 mm of ice was trimmed from the ice core surfaces before extraction in order to remove the external part that could be potentially contaminated with drilling fluid or might have been subject to gas loss during storage in the freezer (Bere-iter et al., 2009). The CO 2 measurements were referenced to a secondary gas standard (synthetic air from Air Liquide, Alphagaz 28416000) containing 233.7 ± 0.4 ppm of CO 2 in dry air, which was referenced to two primary standards (238.34±0.04 from NOAA, CB09707, and 260.26±0.2 from CSIRO, CSIRO1677).
Blank tests using 40 g of artificial bubble-free ice were conducted every 10 measurements to quantify the precision of the system and to correct for the CO 2 contamination caused by the crushing process. Blank tests were conducted in two steps: first, to validate the baseline of the system, a gas standard with 233.7 ± 0.4 ppm was injected over the bubblefree ice in the cell. The gas was then left to equilibrate in the cell for 10 min. The gas was analyzed twice by successive injections into the extraction line and sample loop. Afterwards, the bubble-free ice was crushed and the gas was analyzed five more times. The difference between the results before and after crushing was considered the contamination effect caused by the crushing process. These values were used to estimate the precision of the system. Measured CO 2 should be corrected for contamination caused by the analytical procedure by comparing measured CO 2 of the blank tests with the standard gas value. However, it is not feasible to correct CO 2 concentrations directly. The CO 2 mole fraction is calculated as the ratio between partial pressure of CO 2 and total pressure in the measurement line, which is a relative value. Thus, CO 2 concentration and the concentration of CO 2 con-tamination are dependent on the total pressure. To avoid this dependence, we corrected the absolute value (partial pressure) of CO 2 in the air by the expected partial pressure of CO 2 contamination, as estimated from the blank tests.
For this study, we used four different extraction chambers to hold the ice core samples during crushing and measurement. Each chamber showed different contamination levels. Therefore, blank tests were conducted on each chamber. The data were corrected by the average of each chamber. For these measurements, the precision of the system was estimated to be ∼ 1 ppm on average. On average, the extraction effect correction corresponds to a reduction of the measured CO 2 concentration by 1.7 ± 1.0 ppm (1σ ).
Additionally, we measured 65 samples from 27 depth intervals using the Centrifugal Ice Microtome (CIM) (Bereiter et al., 2013) at CEP. At each depth interval replicates were made and analyzed (2.4 on average). Approximately 5 mm of the ice surface were trimmed, and ∼ 8 g of samples was analyzed on average. The pooled standard deviation of replicates was ∼ 1 ppm. The WMOX2007 mole fraction scale (Tans et al., 2017) was used as a reference for the CO 2 measurements, via four different primary dry air standards for atmospheric CO 2 , which were calibrated at the NOAA Earth System Research Laboratory. The standard sample concentrations varied between 192.44 and 363.08 ppm .
Each CO 2 record was corrected for gravitational fractionation, using the δ 15 N isotope ratio (Craig et al., 1988). To this end, 88 new data points together with existing δ 15  covering the late MIS 6 (156.4-139.2 ka) were used. δ 15 N data were linearly interpolated in age to each corresponding CO 2 data point. On average, the correction corresponds to removing 1.2±0.1 ppm from the measured CO 2 . In total, 177 individual ice samples were measured for CO 2 in the depth range from 2036.7 to 1787.5 m along the EDC ice core, corresponding to the time period from 189.4 to 135.4 ka (Bazin et al., 2013).

CH 4 measurements
We measured the atmospheric CH 4 content of 63 ice core samples, using the wet extraction method at IGE as described in detail in Spahni et al. (2005). This allowed us to improve the temporal resolution of existing CH 4 data (Loulergue et al., 2008) from ∼ 600 to 360 years on the AICC2012 chronology (see Fig. S1 in the Supplement). The precision of the system is ∼ 11 ppb on average, which was estimated by blank tests.
The previous CH 4 dataset (Loulergue et al., 2008) from EDC has been produced at both IGE and CEP. CH 4 measured at CEP used to be systematically higher than CH 4 measured at IGE by 6 ppb (Loulergue et al., 2008). The offsets are due to differences in corrections for contamination caused by the analytical procedure between the datasets. Previously, 6 ppb had been added to the data obtained at IGE (Loulergue et al., 2008).
In this study, a new data correction was applied to data measured at IGE, and it is no longer necessary to add 6 ppb to the new data. The systematic offset between the new data and the previously published data points (Loulergue et al., 2008) is only 1.7 ± 2.4 ppb (n = 63, the standard error of the mean) during MIS 6, (the datasets agree with each other within the measurement uncertainty).

Nitrogen isotopes
Isotopes of molecular nitrogen in air bubbles were measured by a melting technique at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE), France. The gas was extracted from the ice by a wet extraction technique, and the released air was analyzed by a dual-inlet mass spectrometer (Delta V Plus; Thermo Scientific). The analytical method and data correction are described in detail in Bréant et al. (2019). In total, 151 samples from 88 depth intervals (63 duplicates) between the depths of 2124.7 and 1875.0 m below the surface were measured, corresponding to the time interval from 205 to 154 ka ( Fig. S2 in the Supplement). The average resolution on the AICC2012 chronology is ∼ 580 years.

Gas age revision by estimating ∆depth from δ 15 N
The water isotopic signature (δD), is, unlike CO 2 , measured on the ice matrix. Air enclosed in an ice core moves through the porous firn layers at the top of the ice sheet by molecular diffusion and becomes gradually trapped in the ice at the socalled Lock-In Depth (LID, around 100 m below the surface in the case of the EDC ice core). An age difference, thus, exists between the ice and the air at a given depth. For the conditions of the EDC ice core, the age difference can reach up to 5 kyr during glacial maxima and is associated with a large uncertainty (several hundred years; Loulergue et al., 2007). The δ 15 N of molecular nitrogen in air bubbles can be used to determine the LID and to calculate the depth difference between synchronous events in the ice matrix and air bubbles, called delta depth ( depth), thus creating a more precise relative chronology of the gas -with respect to the ice-phase . We use the depth calculation to adjust the gas chronology of EDC, while the ice chronology from AICC2012 remained unchanged.
We calculate the height of the firn column, h, from δ 15 N of N 2 measurements (Craig et al., 1988;Dreyfus et al., 2010;Sowers et al., 1989) using the following equation: In this equation, T diff is the temperature difference between the top and the bottom of the diffusive zone as estimated by the Goujon-Arnaud model, where surface temperature and accumulation are estimated from the stable water isotope record . (T ) is the thermal diffusion sensitivity, which has been estimated from laboratory measurements by Grachev and Severinghaus (2003).
m is the mass difference between 14 N and 15 N (kg mol −1 ), g is the gravitational acceleration (9825 m s −2 for Antarctica) , and R is the universal gas constant (8314 J mol −1 K −1 ). Finally, h conv is the height of the convective zone at the top of the firn column, which is considered negligible at EDC according to current observations . The variation of the convective height is related to changes in wind stress. According to Krinner et al. (2000), wind on the East Antarctic plateau varied little during the LGM (Last Glacial Maximum), and we assume that this is also the case for late MIS 6. Moreover, the convective zone was confirmed to be very thin during the last deglaciation by Parrenin et al. (2012). Thus, we assume that h conv is negligible during MIS 6. depth is calculated from the height of the air column using a constant average firn density and a modeled vertical thinning function suggested by Parrenin et al. (2013).
Raw δ 15 N data cannot be used directly to calculate depth because bubbles at a given depth in the ice sheet can be trapped at slightly different times, due to small physical variations that affect the gas trapping process. This can lead to age scale inversions on depth intervals up to 1 cm (Etheridge et al., 1992;Fourteau et al., 2017). These layer inversions should not strongly affect our record, as its resolution is larger than the scale of age inversion events (Fourteau et al., 2017). However, the change in depth between two different depths (z) in the ice core, denoted ∂ depth/∂z, deduced from the raw δ 15 N data, shows some values higher than 1, that could correspond to small stratigraphic inversions in the gas phase ( Fig. 2). Therefore, a three-point moving average weighted by the time difference between a point and its two neighbors was applied to the δ 15 N dataset. The weights for the three points are equal if the time difference is less than or equal to 500 years, which is close to the average sampling resolution of the δ 15 N dataset. Otherwise, the neighboring data points are weighted by 500/ T , where T is the time difference. This avoids assigning too much weight to neighboring points where the resolution in the record is lower, which would lead to local over-smoothing.
depth as estimated using the three-point moving average weighted by the time difference is shown in Fig. 2. The difference of the Lock-In Depth in Ice Equivalent (LI-DIE) calculated on the AICC2012 age scale (Bazin et al., 2013) and deduced from δ 15 N in this study is 0.5 ± 3.0 m, or 0.7 ± 5.6 % on average (see Fig. S3 in the Supplement). The AICC2012 LIDIE was calculated using a background scenario derived from δD, using a linear relationship between δD and δ 15 N. Our results show this relationship to be relatively unbiased but not entirely exact (Fig. S3 in the Supplement). The depth values are used to update the EDC gas chronology from the original AICC2012 age scale, while the ice chronology remains unchanged. For MIS 6, this method significantly reduces the relative age uncertainty between air and ice to 900 years on average with respect to the original AICC2012 chronology (see Fig. S4 in the Supplement).

Age scale of the MD01-2444 marine sediment record
The MD01-2444 marine sediment core from the Iberian Margin (Margari et al., 2010) represents an important archive for the interpretation of our CO 2 record, as it is well resolved during MIS 6, and it provides high-resolution benthic and planktonic foraminiferal records. The original MD01-2444 age scale was built by matching the benthic δ 18 O record (Shackleton et al., 2000) to the δD record from EDC on the EDC03 ice age scale  using 11 tie points (Margari et al., 2010). In order to make optimal use of the record, its relative chronology with respect to EDC requires additional improvement using the latest version of the EDC age scale. Thus, in this study, we recalculated the age of the tie points using the AICC2012 chronology (Bazin et al., 2013), and the age of the sediment record was linearly interpolated between the tie points (Table S1 in the Supplement).

Definition of NA stadial duration
Due to the absence of a Greenland temperature record for MIS 6, the durations of the six NA stadials must be identified using other proxies. Here, we used two methods using both ice core and sediment core proxies. First, we estimated the durations of the six NA stadials using δ 18 O of planktonic foraminifera and tree pollen in MD01-2444, which reflect temperature variability in the NH (Margari et al., 2010) (Fig. 3). The midpoint of the stadial transitions in both δ 18 O of planktonic foraminifera and tree pollen in MD01-2444 were used to identify the NH stadial transitions. The time interval between two stadial transition points were defined as the NA stadial duration. Using this approach the duration of stadial at 6ii in the NH could not unambiguously be determined due to small variations in the foraminifera and the pollen record. However, the average age difference between the durations identified using the two records is only 205 years, which is less than the sampling resolution of MD01-2444 during MIS 6. The uncertainty of the timing of each stadial transition was estimated as half of the temporal difference between maxima and minima of δ 18 O of planktonic foraminifera. The mean values for stadial onsets and terminations are shown as dashed lines in Fig. 3. However, not all of the stadial durations during MIS 6c to 6d are entirely unambiguous using this method. Second, we recalculated the durations of the six NA stadials during MIS 6 period using a method developed by Margari et al. (2010). Figure 3 also shows a synthetic Greenland δ 18 O ice record (Barker et al., 2011) and Antarctic δD   , (b) synthetic Greenland δ 18 O ice record (Barker et al., 2011), (c) tree pollen percentage in the MD01-2444 core (Margari et al., 2010) and (d) δ 18 O of planktonic foraminifera in the MD01-2444 core (Margari et al., 2010). Proxy data shown here are given on the AICC2012 age scale. Red lines indicate the midpoints of the stadial transition of both δ 18 O of planktonic foraminifera and tree pollen in MD01-2444. Light green bars indicate the uncertainty of the duration of each stadial transition, estimated as half the temporal difference between maxima and minima of δ 18 O of planktonic foraminifera before and after the transition. Red dots indicate minima and maxima of δD in the EDC ice core as selected in this study. The event numbers are indicated at the top.
variations on the AICC2012 age scale. The interval between the maximum and the preceding minimum of δD in the EDC record can also be used to estimate the duration of the stadial transitions (Gottschalk et al., 2020;Margari et al., 2010). In most cases, the synthetic Greenland δ 18 O ice record and the interval between the maximum and the preceding minimum of δD in the EDC record confirm the definition of NA stadials selected by δ 18 O of planktonic foraminifera and tree pollen in MD01-2444. However, the duration of the NA stadial in MIS 6iii is not clearly confirmed by synthetic Greenland δ 18 O ice and δD in the EDC (Fig. 3).
The interval between each maximum and the preceding minimum of δD in the EDC record was calculated to estimate the stadial durations using the second method (Gottschalk et al., 2020;Kawamura et al., 2017;Margari et al., 2010). Maxima and minima in the δD record were selected by finding zero values in the second (Savitzky-Golay filtered) derivative of the data (the same method we used to pick minima and maxima of atmospheric CO 2 and temperature; Figs. S5-S6, Tables S2-S3 in the Supplement). In one case (the minimum before MIS 6i) two potential minima were selected by this method; these values were then averaged.
The red dots and error bars in the EDC δD record in Fig. 3 show the estimated minima and maxima of temperature corresponding to stadial transitions using this method, along with their uncertainties. However, using this tool, durations of MIS 6ii and 6i are apparently overestimated due to ambiguities concerning the definition of the maximum in MIS 6i and minimum in 6ii.
Neither our method nor that of Margari et al. (2010) can be considered absolutely correct. To account for the differences between the two methods, we took the stadial duration to be the mean of the durations estimated by the two methods ( Table S4 in the Supplement).

Data compilation
We measured 177 samples of atmospheric CO 2 from EDC during MIS 6 using two extraction systems (the ball mill at IGE and CIM system at CEP). To improve the resolution of our new dataset even further, we made a composite dataset by aligning previous sets of CO 2 measurements made over the MIS 6 period on the EDC ice core to our new data. First, we compared the two existing CO 2 datasets and the new CO 2 dataset from EDC ( Fig. 4 and Table S5 in the Supplement). There are two published CO 2 datasets for EDC during MIS6 -the first measured using the ball mill system at IGE (Lourantou et al., 2010) and the second using a sublimation extraction system at CEP (Schneider et al., 2013). We also added unpublished data measured in 2003 using the CEO needle cracker measurement system in 2003 (Monnin et al., 2004;Siegenthaler et al., 2005) (see Supplement for details). All records are on the AICC2012 air age scale (Bazin et al., 2013). All datasets are corrected for the gravitational fractionation effect using the new δ 15 N data in our study. In total, the datasets contain 237 CO 2 measurement points. Two samples were excluded because of system operator error. Figure 4 shows CO 2 concentrations measured by the ball mill system, sublimation (Schneider et al., 2013), the CIM and the needle cracker. Concentrations from the CIM, sublimation and the needle cracker are systematically higher than CO 2 concentrations measured using the ball mill system (Table S5 in the Supplement).
The offsets between each additional dataset and our data were calculated and corrected using a Monte Carlo procedure to account for uncertainty and a Savitzky-Golay filter used to remove noise from our dataset (Supplement). The offsets between the multiple datasets may be at least in part linked to differences in extraction efficiency between the measurement methods. The sublimation and CIM systems have high extraction efficiency on clathrates and should therefore present more unbiased baseline CO 2 values. However, since these datasets are of much lower resolution, we have aligned all datasets to the baseline value of our ball mill dataset.
Despite the systematic offsets, the other datasets confirm the millennial-scale variations shown in the data from this study ( Fig. 4 and Fig. S7 in the Supplement). Although none of the individual additional datasets is of high enough resolution to show millennial-scale variations in detail, when aligned to our data the new data follow the millennial-scale variations with very few outliers. The new composite atmospheric CO 2 record from EDC is shown in Fig. 5. 3.2 Comparison with Vostok CO 2 record for MIS 6 The new composite atmospheric CO 2 record from EDC was compared to the existing CO 2 data from the Vostok ice core, measured using the ball mill system (Fig. S7 in the Supplement), where the CO 2 record from Vostok was aligned to the AICC2012 gas age scale (Bazin et al., 2013). Atmospheric CO 2 data from Vostok are corrected for the gravitational fractionation effect using the existing δ 15 N data (Bender, 2002).
Atmospheric CO 2 variability in the EDC ice core shows similar general patterns as those retrieved from the Vostok ice core, although previously unidentified features are resolved due to the improved temporal resolution. CO 2 concentrations from Vostok appear systematically higher than those from EDC by 4.6±3.0 ppm on average and the difference increases to more than 10 ppm at the beginning of Termination 2 at around 135 ka.
Three main mechanisms may explain the offsets between the EDC and Vostok measurements: firstly, contamination due to the extraction system was assumed to be negligible at the time of the Vostok sample measurements at IGE (Petit et al., 1999). Based on recent ball mill measurements conducted within the scope of this study, we presume that the neglected correction concerning the data published by Petit et al. (1999) amounts to approximately 1.7 ppm, reducing the average apparent offset to about 3 ppm. Secondly, part of this remaining offset between the new EDC and previously published Vostok record may be related to age scale uncertainties due to a limited number of stratigraphic tie points between the two cores (Bazin et al., 2013). In particular, this effect may explain part of the larger offset during the glacial termination in air younger than 140 ka (Fig. S7 in the Supplement). Thirdly, the ball mill system has a different extraction efficiency depending on the presence of bubbles and/or clathrates in the ice sample, which may affect the accuracy of the reconstructed absolute mean CO 2 level. When the air is extracted from an ice core sample where bubbles and clathrates coexist, different proportions of bubbly and clathrate ice may lead to biased CO 2 concentrations Schaefer et al., 2011). The Vostok measurements were made on recently drilled ice, in which clathrates had less time to transform into secondary bubbles (Vladimir Lipenkov, 2015 personal communication). Accordingly, at least part of the approximately 3 ppm offset may be due to these systematic extraction differences (see Supplement). . Atmospheric CO 2 from EDC and Vostok ice cores, compared to δD of the ice in the EDC core (temperature proxy) during 190-135 ka. Blue dots represent atmospheric CO 2 from EDC as measured by the ball mill system (this study). The error bars of the new CO 2 data using the ball mill extraction system indicate the standard deviation of five consecutive injections of the gas extracted from each sample into the gas chromatograph added to the precision of the measurement estimated by the reproducibility of the control measurement (∼ 0.8 ppm) using a quadratic sum. Yellow dots represent atmospheric CO 2 from EDC measured using the ball mill system (Lourantou et al., 2010). The error bars indicate the standard deviation of five consecutive injections of the gas extracted from each sample into the gas chromatograph. Red triangles represent atmospheric CO 2 from EDC measured using the needle cracker (this study). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval. Black inverted triangles represent atmospheric CO 2 from EDC measured using the CIM (this study). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval. Green diamonds represent atmospheric CO 2 from EDC measured using sublimation extraction (Schneider et al., 2013). The error bars of data points with replicates indicate the standard deviation of the mean of replicates from the same depth interval (Schneider et al., 2013). Grey dots represent atmospheric CO 2 from the Vostok ice core (Petit et al., 1999). The error bars of the CO 2 data from Vostok (Petit et al., 1999) show the estimated overall accuracy for CO 2 measurements. The grey line represents δD of ice in the EDC core .

Figure 5.
Composite atmospheric CO 2 (left axis) from the EDC ice core (this study) compared to the EDC water isotopic record (right axis) . The blue line indicates a Savitzky-Golay filter using a 150-year cut-off period (dotted blue line). Vertical dotted blue lines indicate the six CDM events that we identify during the early MIS 6. The numbers of the CDM events are indicated at the top of the figure.
In summary, the individual datasets used for our data compilation and the Vostok ice core record all show the same millennial CO 2 variability over MIS 6. Analytical and icerelaxation-related offsets could be removed to synthesize all records. This allows us to compile a precise and high-resolution record; however, the true absolute mean level of atmospheric CO 2 is only known to be better than 5 ppm, due to this offset correction. However, this has no impact on the conclusions drawn from the millennial CO 2 variability, which is the focus of this study.  (Berger, 1978), (b) ice-rafted debris (IRD) input in the Iberian Margin core MD95-2040 (de Abreu et al., 2003), (c) atmospheric CH 4 in the EDC ice core (Loulergue et al., 2008; this study), (d) δ 18 O of planktonic foraminifera in the Iberian Margin marine core MD01-2444 (Margari et al., 2010), (e) δ 18 O of benthic foraminifera in the Iberian Margin marine core MD01-2444 (Margari et al., 2010), (f) δD of the EDC ice core  and (g) our new composite CO 2 record during the MIS 6 period. The numbers of CDM events are indicated at the top. Margari et al. (2010) suggested that MIS 6 can be divided into three sections depending on the degree of climatic variability observed in δD (indicative of Antarctic climate variability) and CH 4 (reflecting NA climate variability) in EDC: early (185.2-157.7 ka), transition (157.7-151 ka) and late MIS 6 (151-135 ka) (Figs. 5 and 6). Climatic oscillations on millennial timescales are pervasive during the early MIS 6 period (185-160 ka) (Barker et al., 2011;Cheng et al., 2016;Jouzel et al., 2007;Margari et al., 2010Margari et al., , 2014, which is similar to MIS 3 ( Figs. 1 and 6). However, during the late MIS 6 period, i.e., the Penultimate Glacial Maximum (PGM), millennial variability is subdued and more resembles the climate variability on millennial timescales during MIS 2 ( Figs. 1 and  6). During the transitional period from 157 to 151 ka, δD in EDC slowly increased . Like δD in EDC, CO 2 variations on millennial timescales are pervasive during the early MIS 6 period (185-157.7 ka). During the transitional period from 157 to 151 ka, atmospheric CO 2 increased slowly, while during the late MIS 6 period CO 2 variation is subdued (Fig. 5).

Atmospheric CO 2 variability on millennial timescale
While there was an indication of millennial CO 2 variability in the Vostok record, the resolution of that record did not allow us to constrain this variability. We now present clear evidence of millennial variability of CO 2 concentrations during MIS 6 that are associated with Antarctic Isotope Maxima (AIM) events (EPICA Community Members, 2006), thanks to the improved time resolution and precision of the obtained CO 2 data (Fig. 5). To better discuss millennial variability in MIS 6, a Savitzky-Golay filter with a 150-year cut-off period was used to filter out centennial-scale variability and noise (see Fig. 5 and Fig. S5 in the Supplement). Five prominent and one subdued CO 2 variations were detected in atmospheric CO 2 during early MIS 6 ( Fig. 5 and Fig. S5 and Table S2 in the Supplement), where we name CDM according to the numbering by Margari et al. (2010) and Gottschalk et al. (2020) (Fig. 5). The five prominent peaks are observed at 160.7 ± 0.3 (CDM 6vi), 164.2 ± 0.3 (CDM 6v), 169.6 ± 0.2 (CDM 6iv), 174.3±0.2 (CDM 6iii) and 181.3±0.1 ka (CDM 6i). Note that the provided uncertainty was calculated with respect to the position of each maximum and does not include the absolute age uncertainty of the ice core (in each case around 3.0 kyr, 1σ ) (Bazin et al., 2013). We can also identify one low-amplitude CO 2 peak at around 150 ka, representing another potential candidate for a CDM (Fig. 5). This atmospheric CO 2 variation is of triangular shape and follows the δD pattern. The change of direction is also associated with a CH 4 peak. This variation has analogues in MIS 4 and MIS 10 .
Each CDM is associated with an AIM event. The short AIM 6ii event corresponds to CDM 6ii at around 178 ka. CDM 6ii has an amplitude of only ∼ 4 ppm, and the CO 2 variation becomes even less pronounced after filtering. Due to the smoothing of the CO 2 variation at CDM 6ii, atmospheric CO 2 and δD composition in EDC appear decoupled. This observation seems confirmed when considering the relationship between atmospheric CO 2 change and the duration of NA stadials calculated using tree pollen and the δ 18 O composition of planktic foraminifera in an Iberian Margin core (Margari et al., 2010) for MIS 6 ( Fig. 3 and Table S4 in the Supplement), or using isotopic records from Greenland ice cores (Rasmussen et al., 2014) for MIS 3 as shown in Fig. 7. Here, the magnitude of atmospheric CO 2 change is generally correlated with the NA stadial duration (r = 0.87, n = 6) during the early MIS 6 period, and CDM 6ii is by far the shortest of all six events detected in MIS6.
We note a similar correlation between the NA stadial duration and atmospheric CO 2 change during MIS 3 (r = 0.85, n = 14). When the NA stadial duration was shorter than 1500 years, atmospheric CO 2 varied less than 5 ppm Bereiter et al., 2012) as is the case for CDM 6ii. Margari et al. (2010) note one exception, AIM 14, during which ice discharge may have led to a stronger perturbation to AMOC.
Both Bereiter et al. (2012) and Ahn and Brook (2014) observed that during short NA stadials which last less than 1300 years, the CO 2 maxima do not appear to have a consistent phase relationship with AIM and CO 2 and δD anomalies are not correlated . We observe the same to be the case when MIS 6 is included (r = 0.07, n = 11), however, we cannot exclude that the firn column-induced smoothing at EDC obliterated small atmospheric responses with magnitudes smaller than 5 ppm for such short stadials. In both MIS 3 and 6, CO 2 is highly correlated with δD anomalies during the longer stadials (r = 0.78, n = 9) with a clear increase in CO 2 Bereiter et al., 2012).
We observe that during the last two glacial periods, the amplitude of CO 2 is strongly related to the NA stadial duration (r = 0.88, n = 20). This observation implies that despite the different climate boundary conditions during the last two glacials, the behavior of atmospheric CO 2 was similar (see Fig. 1) and the overall bipolar seesaw coupling of climate and atmospheric CO 2 acted the same way.

Leads and lags between CO 2 and the abrupt warming in NH
To better understand the link between the bipolar seesaw mechanism and atmospheric CO 2 variability on millennial timescales, we calculated the varying time lag for each CDM following abrupt warming events in the NH (see Figs. S8-S11 in the Supplement) (Bereiter et al., 2012). Due to the lack of an MIS 6 temperature proxy in Greenland, and due to the difficulty of placing marine temperature proxies (Shackleton et al., 2000) on a precise common chronology with the EDC ice core, in this work CH 4 measurements performed on the EDC ice core were used as a time marker of rapid warming in the NH (Baumgartner et al., 2014;Brook et al., 1996;Huber et al., 2006). Because CH 4 and CO 2 signals are both imprinted in the air bubbles, there is no chronological uncertainty when comparing the timing of changes of those two signals. The only remaining uncertainty is related to analytical uncertainties and to the temporal resolution of the two records. We pick intervals when CH 4 increases rapidly by at least 50 ppb over a time period of less than 1 kyr that correspond to AIM (Buizert et al., 2015;Loulergue et al., 2008). The timing of abrupt CH 4 increases was defined as the midpoint between the beginning of the increase of CH 4 and its peak (see Supplement for details). The age uncertainty of the midpoint is defined by half the time difference between the two endpoints. Figure 8 shows the shifts of CDM with respect to the onset of the abrupt warming in the NH. During the MIS 6 period, three abrupt NH warmings (as inferred from the CH 4 signal) at 181.5 ± 0.3, 175.4 ± 0.40 and 171.1 ± 0.2 ka (2σ ) were found. These events correspond to CDM 6i, CDM 6iii and CDM 6iv, respectively. CDM 6vi, CDM 6v and CDM 6ii do not have corresponding rapid changes in the methane record that fulfill our detection criterion of a 50 ppb increase; this may be due to slow gas trapping as compared to interglacial periods, which could smooth out smaller changes. A synthetic Greenland temperature record (Barker et al., 2011) shows abrupt temperature jumps at CDM 6vi, CDM 6v and CDM 6ii as well. However, this record is calculated using EDC δD, and the large relative chronological uncertainty (∼ 900 years on average) between ice and gas phases does not allow us to make any conclusions about leads and lags using this record. We therefore exclude these events from our lag analysis.
From MIS 6i to MIS 6iv, the lag of CO 2 with respect to abrupt warmings in the NH, which were identified from this chronological comparison between EDC CH 4 and CO 2 , becomes significantly larger. During the earliest MIS 6, atmospheric CO 2 increases rapidly (by ∼4.2 ppm in 240 ± 320 years) following the abrupt CH 4 increase at 181.5 ± 0.3 ka. The peak of CDM 6i is nearly synchronous with the onset of the NH abrupt warming (nonsignificant lag of 240±320 years, Fig. 8). During CDM 6iv and 6iii, CO 2 concentrations show a much slower increase over a duration of ∼ 3.3 kyr. Here, the CO 2 maximum lags significantly behind the onset of the NH abrupt warming by 1460 ± 270 years and 1110 ± 460 years, respectively (1290 ± 540 years on average, with the error calculated by propagation of the uncertainties in the individual events). Interestingly, these two CDM events occurred during MIS 6d (Fig. 1), when iceberg discharge was muted and the ITCZ is thought to have shifted northward, intensifying monsoon systems in lowlatitude Northern Hemisphere regions, such as in Asia, the Apennine Peninsula and the Levant (Ayalon et al., 2002;Bard et al., 2002;Cheng et al., 2016). This may have led to a weaker overturning circulation due to the reduction of the density of the NA surface water, making the AMOC cell shallower with a smaller threshold in NA during MIS 6 than during MIS 3 (Margari et al., 2010). Therefore, the two different CO 2 lag timescales with respect to abrupt warming in NH during MIS 6 might be explained by this difference in background climate conditions. Indeed, these features during MIS 6 appear fully compatible with those observed during the last glacial period (Bereiter et al., 2012). To compare CO 2 variations with respect to abrupt warming in the NH during the last two glacial periods, in this study we also re-estimated the abrupt CH 4 rise during the last glacial period with the same methodologies. Figure 8 also shows the CO 2 evolution during the onset of abrupt warming in the NH during MIS 3 and MIS 5 Bereiter et al., 2012). Atmospheric CO 2 during MIS 3, as shown in Fig. 8, was reconstructed from the Talos Dome ice core (TALDICE). For MIS 5 it was obtained from Byrd and the EPICA Dronning Maud Land (EDML) ice core (Ahn and Brook, 2008;Bereiter et al., 2012). In Bereiter et al. (2012), both TALDICE and EDML records were used during MIS 3 and compared to the onset of abrupt warming in the NH as indicated by the co-occurring CH 4 rise. However, here we only use data from TALDICE, which are more accurate due to the narrower age distribution of the gas trapped in the LID (Bereiter et al., 2012). Using the same method as above, the average value of CDM lag with respect to the abrupt warming in NH was calculated. The average CDM lags with respect to the abrupt warming in the NH for the MIS 3 and 5 periods are 770±180 and 280 ± 240 years. Thus, over the course of the last glacia-tion, the lag of CO 2 maxima with respect to the abrupt NH warming events significantly increased. We observe the same trend through the millennial events depicted during MIS 6, albeit with different absolute lags.

Discussion
4.1 Atmospheric CO 2 variability on millennial timescales Similar to the AIM amplitude EPICA Community Members, 2006), we found that the amplitude of atmospheric CO 2 variations is well correlated to the NA stadial duration during MIS 6 and MIS 3, which implies that the amplitude of CO 2 variations might also be affected by the duration of AMOC disruption during the early MIS 6 period (Margari et al., 2010). This hypothesis is also supported by a recent study using oceanic sediment cores from the Southern Ocean (Gottschalk et al., 2020). The authors report that respired carbon levels in the deep South Atlantic decrease when AMOC is weakened during both glacial periods, and the amount of carbon loss in the deep South Atlantic is highly correlated with the duration of NA stadials.
As mentioned above, atmospheric CO 2 on millennial timescales can be controlled by CO 2 exchange between the ocean and the atmosphere, as well as changes of terrestrial carbon stocks. Coupled climate carbon cycle models reported that the variations of atmospheric CO 2 concentration on millennial timescales are mainly dominated by the deep ocean inventory, requiring a few millennia to re-equilibrate to climate change (Schmittner and Galbraith, 2008). On the other hand, the response of the terrestrial biosphere is usually fast (decadal to centennial timescale) (Bouttes et al., 2012;Menviel et al., 2014;Schmittner and Galbraith, 2008). Although different models differ significantly in the CO 2 response to AMOC changes, the initial CO 2 evolution of the terrestrial biosphere and deep ocean to AMOC perturbations are opposite in model simulations (Gottschalk et al., 2019). Thus, due to the opposite direction of CO 2 change of ocean and terrestrial reservoirs, atmospheric CO 2 variations might be muted if the NH duration is short (Bouttes et al., 2012;Menviel et al., 2014;Schmittner and Galbraith, 2008), while during long stadials the carbon release from the ocean dominates. There is, on the other hand, evidence that not all of the processes of CO 2 exchange follow these general trends. For example, atmospheric CO 2 might be changed on centennial timescales by carbon exchange between the deep and surface ocean (Rae et al., 2018) or atmospheric CO 2 might be influenced slowly by soil decomposition (Köhler et al., 2005), and it is important to note that modeling studies do not agree in the amplitude or even the direction of the modeled net CO 2 exchange.
The more prominent CO 2 changes during stadials involved with the Heinrich events may be related to a stronger reduction of the North Atlantic Deep Water (NADW) forma- Figure 8. CDM lags relative to abrupt temperature increases in the NH. Dotted grey lines indicate when climate changes abruptly in the NH as indicated by the CH 4 jumps. During the last glacial period, the AIM number corresponds to the DO number for corresponding DO and AIM events, and for MIS6 numbers correspond to event numbers. (a) Atmospheric CO 2 as measured using the TALDICE ice core during MIS 3, (b) atmospheric CO 2 as measured using the Byrd ice core during MIS 5, (c) atmospheric CO 2 as measured using EDML ice core during MIS 5, (d) new CO 2 composite as derived in this study for MIS 6. For MIS 6, we selected the 3 CDM that correspond to an abrupt methane increase; the other CDM do not correspond to an abrupt change that fulfill our 50 ppb detection criterion (see main text for details). Note that the scale of the y axis is not the same for the four panels. tion during Heinrich events (Henry et al., 2016;Margari et al., 2010), which would cause a stronger upwelling of deep water in the Southern Ocean (Menviel et al., 2008;Schmittner et al., 2007). These events may reduce stratification in the Southern Ocean due to an increase in salinity of the surface waters and a relative freshening of the deep water (Schmittner et al., 2007). As a result, atmospheric CO 2 can be increased due to upwelling and outgassing of CO 2 in the Southern Ocean (Schmittner et al., 2007;Schmitt et al., 2012). The co-occurring upwelling in the SO during AIM for the last termination has been examined (Anderson et al., 2009), but due to the lack of proxy data with precise age scale for upwelling in the Southern Ocean, this hypothesis cannot be confirmed during MIS 6.
According to the results of Margari et al. (2010), the six AIM events of MIS 6 were likely affected by AMOC perturbations of similar strength. During these events, there is no clear evidence for freshwater perturbation in the NA (Fig. 1), and the strength of the associated AMOC perturbations is estimated to be similar to that during non Heinrich stadials or non-classic Heinrich (different ice-rafted detritus source signatures from Heinrich events, occurring at transitions in ice volume) AIM events of MIS 3. The durations of NA stadials during early MIS 6 (except for AIM 6ii) appear to be longer than those during non-Heinrich AIM events of MIS 3, which might be caused by the different climate boundary conditions during MIS 3 and the increase of the hydrological cycle during the early MIS 6. A longer duration of the AMOC disruption apparently impacts the amplitude of CO 2 variations (Bouttes et al., 2012;Menviel et al., 2008). Considering that the duration of the AMOC disruption may be related with climate background conditions (Bouttes et al., 2012;Menviel et al., 2008), this observation suggests how different climate background conditions may impact atmospheric CO 2 . The strength of AMOC perturbations also appears to be an important factor in determining the amplitude of CO 2 variations. For example, the duration of the Heinrich events in MIS 3 (AIM 8, 12 and 14) is shorter than any of the MIS 6 events except for 6ii, but atmospheric CO 2 varied significantly in all three.
The relationship between the amplitude of atmospheric CO 2 variations and the NA stadial duration is explained by the duration of AMOC disruption during the early MIS 6 period. However, the temporal resolution of δ 18 O composition of planktonic foraminifera in MD01-2444 and the precision of the age scale are too low to precisely define the duration of stadials during MIS 6. Additional proxy data providing information about climate change in the NH are needed to confirm the relationship between atmospheric CO 2 variations and the NA stadial duration (Fig. 3). The limited available proxy data permit only to formulate a hypothesis for the mechanisms responsible for CO 2 variability during MIS 6 but not to rigorously test it. To compare the behavior of the bipolar seesaw with atmospheric CO 2 variations, additional investigations about AMOC disturbances and their associated climate responses are needed.
4.2 Why did CO 2 lag the abrupt warming in the NH during MIS 6d?
Two different lags of the CO 2 maxima with respect to NH warming are present in the MIS 6 period (Fig. 8). CDM 6i is nearly synchronous with the abrupt warming in the NH (no significant lag of 240 ± 320 years), while the lags for CDM 6iii (1110 ± 460 years) and CDM 6iv (1460 ± 270 years) are much longer. Two modes of CO 2 variations are also observed during the last glacial period. As the last glaciation progressed from MIS 5 to MIS 3 ( Figs. 1 and 8), the lag of CO 2 maxima with respect to NH millennial-scale warming significantly increased. This observation may be explained by the different AMOC settings in MIS 5 and MIS 3 (Bereiter et al., 2012). We speculate that, as observed during the last glacial period, the configuration of oceanic circulation during MIS 6d might also be the cause of the change in the time lags between NH abrupt warming events and CO 2 variations during the early MIS 6. Bereiter et al. (2012) explained that during MIS 3 the oceanic circulation in the Atlantic was in a more "glacial" state, with shallower NADW and carbon-rich Antarctic Bottom Water (AABW) extended to the north, while during MIS 5 ocean circulation was similar to the present, in what can be referred to as a "modern-like" state. At the onset of abrupt warming in the NH, AMOC is thought to accelerate rapidly, delivering heat to the north and resuming the formation of NADW. When the NADW cell expands, AABW is withdrawn and the upwelling of carbon-rich deep water in the Southern Ocean is enhanced. Essentially, over the timescale of ocean overturn, part of the previously expanded carbon-rich southern-sourced water is converted to carbonpoor northern-sourced water.
Thus, during MIS 3 CO 2 continues to be released into the atmosphere for another 500 to 1000 years after the NADW resumption (Bereiter et al., 2012). However, during MIS 5 where AABW was not expanded in the NA, less CO 2 can be released from the ocean to the atmosphere after NADW resumption.
The temporal resolution of proxy data related to oceanic circulation during MIS 3 and 5 is unfortunately not sufficient to validate from the marine realm itself whether the two different modes of CO 2 variations reflect the hypothesized mechanism described above. Modeling studies of the carbon stock in AABW and NADW during MIS 5 and 3 have been attempted. However, dependent on the chosen model, the modes of atmospheric CO 2 variation are different (Gottschalk et al., 2019). Some studies (for example, Menviel et al., 2008) mention an atmospheric CO 2 decrease when the AMOC is resumed. Others, for example Bouttes et al. (2012), confirm an atmospheric CO 2 release from the ocean when AMOC resumed.
In spite of the inconclusive modeling studies (Gottschalk et al., 2019), limited proxy evidence does not exclude the possibility that the configuration of AMOC and its changes over MIS 6 may explain the presence of two different CDM lags. We find this hypothesis to be worth at least a speculative discussion. According to the benthic δ 13 C record in the MD01-2444 core (Margari et al., 2010), the value of benthic δ 13 C during 180-168 ka was lower than during MIS3, which indicates that the NA overturning cell during MIS 6 was likely even shallower than that during MIS 3 (Margari et al., 2010). This implies southern-sourced water masses were more expanded to the north, and the density difference between the northern-sourced water masses and southern-sourced water masses increased. This shallower oceanic circulation during MIS 6 (Margari et al., 2010) could have caused the millennial-scale delays with respect to abrupt NH warming events. It could also explain the longer lag between the abrupt warming in NH and CDM during MIS 6d (1290 ± 540 years on average) when compared to the lags of CDM (770 ± 180 years on average) during MIS 3. However, because of the low accumulation at EDC and its wider age distribution, the estimation of the exact timing of CDM from the EDC ice core might be less accurate compared to that from the TALDICE ice core with a narrower gas age distribution (Bereiter et al., 2012). The remaining uncertainty is related to analytical uncertainties and to the temporal resolution of the two records. To further investigate the exact relationship between CDM and abrupt warming in the NH, additional CO 2 measurements from a higher accumulation site would be helpful.

Conclusions
Using new and existing CO 2 data from the EPICA Dome C ice, we reconstruct a high temporal resolution record of atmospheric CO 2 during the MIS 6 period (189-135 ka). In this study, we investigate how different climate background conditions during the last two glacial periods may have impacted atmospheric CO 2 . Millennial-scale atmospheric CO 2 changes are revealed during the last two glacial periods, with amplitudes ranging between 15 to 25 ppm, mimicking similar patterns in Antarctic δD variations Bereiter et al., 2012). On the other hand, during short NA stadials which last less than 1300 years, atmospheric CO 2 variations are negligible and decoupled from δD in EDC. This finding suggests that during the last two glacial periods the amplitude of millennial CO 2 variations is strongly influenced by the NA stadial duration (r = 0.88, n = 20).
In the earliest MIS 6 (MIS 6i and 6iv, corresponding to 189 to 169 ka), a change of CO 2 lags with respect to the onset of rapid NH -warming as deduced from atmospheric CH 4 changes -is revealed. CDM 6i (at ∼ 182 ka) is nearly synchronous with the abrupt warming in the NH (nonsignificant lag of 240 ± 320 years), while the lags during MIS 6d corresponding to CDM 6iv and 6iii (at ∼ 171 and ∼ 175 ka, respectively) are much longer, 1290 ± 540 years on average. Similar observations are drawn for the time period covered by our study as for previous studies on MIS 3 and MIS 5 periods, although the lag of CO 2 with respect to NH warming reaches even larger values during MIS 6d. We tentatively attribute this to a generally weaker and shallower AMOC during MIS 6 compared to MIS 3 as suggested by the results from Margari et al. (2010). However, the limited avail-able proxy data from the marine realm only permits an exploratory discussion of the mechanisms responsible for CO 2 variability during MIS 6. Because the boundary conditions of the last glacial period cannot be applied to MIS 6, additional proxy data and multiple modeling studies conducted during the MIS 6 period are needed.
Data availability. Data available in the Supplement. All data will be available on the NOAA (National Oceanic and Atmospheric Administration) and PANGAEA (Paleoclimatology database websites) shortly.
Author contributions. The research was designed by JS, RG, FP and JC. The CO 2 measurements were performed by JS with contributions from GT, LS and BB. The data analyses were led by JS with contributions from JCB, RG, FP and JC and HF. The methane data was provided by GT and JC. The nitrogen isotopes data was provided by AL. JS wrote the manuscript with inputs from all authors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This work is a contribution to the "European Project for Ice Coring in Antarctica" (EPICA), a joint European Science Foundation/European Commission scientific program, funded by the European Union and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The main logistical support was provided by IPEV and PNRA. This is EPICA publication no. 316. It also received funding from the European Community's Seventh Framework Programmes ERC-2011-AdG under grant agreement no. 291062 (ERC ICE&LASERS). As it is part of the PhD work of Jinhwa Shin, it was also supported by the LabEX OSUGat2020 project of the Grenoble Observatory of Sciences of the Universe (OSUG). The Swiss authors also acknowledge long-term financial support for ice core research at the University of Bern by the Swiss National Science Foundation under grants 200020_159563, 200020_172745, 200020_172506 and 20FI21_189533. The authors would like to thank Grégoire Aufresne for providing assistance with the additional methane measurements, Xavier Faïn and Kévin Fourteau for their help with the CH 4 analytical system, and Dominique Raynaud for the valuable discussions. We thank Eric Monnin and Urs Siegenthaler for providing additional CO 2 data. We would like to thank Julia Gottschalk for the discussions about CO 2 variability and climate change during the last two glacial periods.
Financial support. This research has been supported by the ERC-2011-AdG (grant no. 291062, ERC ICE&LASERS), the Swiss Na-