Age model

The logging, scanning and sampling of the sediment core (JM-FI-19PC) are described in Ezat et al.14 The sediment core JM-FI-19PC is aligned to the Greenland ice core NGRIP based on the identification of common tephra layers and by tuning increases in magnetic susceptibility and/or increases in benthic foraminiferal δ18O values to the onset of DO interstadials in the Greenland ice cores14,15,16 (Supplementary Fig. 1). In support of the reconstructed age model, eleven calibrated radiocarbon dates measured in N. pachyderma (with no attempt to correct for past changes in near-surface reservoir ages) show strong consistency with the tuned age model for the past 50 kyr (ref. 14).

Boron isotope and minor/trace element analyses

Only pristine N. pachyderma specimens with no visible signs of dissolution were picked from the 150 to 250 μm size fractions for boron isotope (200–450 specimens) and minor/trace element (70–160 specimens) analyses. For boron isotope measurements, the foraminifer shells were gently crushed, and cleaned following Barker et al.45 This cleaning protocol includes clay removal, oxidative and weak acid leaching steps. Thereafter, the samples were dried and weighed to determine the amount of acid required for dissolution. Immediately before loading, samples were dissolved in ultrapure 2N HCl, and then centrifuged to separate out any insoluble mineral grains. One μl of boron-free seawater followed by an aliquot of sample solution (containing 1–1.5 ng B per aliquot) were loaded onto outgassed Rhenium filaments (zone refined), then slowly evaporated at an ion current of 0.5A and finally mounted into the mass spectrometer. Depending on sample size, five to ten replicates were loaded per sample. Boron isotopes were measured as BO 2 - ions on masses 43 and 42 using a Thermo Triton thermal ionization mass spectrometer at the Lamont-Doherty Earth Observatory (LDEO) of Columbia University. Each sample aliquot was heated up slowly to 1,000±20 °C and then 320 boron isotope ratios were acquired over ∼40 min46. Boron isotope ratios are reported relative to the boron isotopic composition of SRM 951 boric acid standard, where δ11B (‰)=(43/42 sample /43/42 standard −1) × 1,000. Analyses that fractionated >1‰ over the data acquisition time were discarded. The analysis of multiple replicates allows us to minimize analytical uncertainty, which is reported as 2s.e.=2s.d./✓n, where n is the number of sample aliquots analysed. The analytical uncertainty in δ11B of each sample was then compared with the long-term reproducibility of an in-house vaterite standard (±0.34‰ for n=3 to ±0.19‰ for n=10) and the larger of the two uncertainties is reported (Supplementary Table 1). Two samples were repeated using the oxidative-reductive cleaning procedure from Pena et al.47 and yielded indistinguishable δ11B values (Supplementary Table 1).

Trace and minor element analytical procedures followed cleaning after Martin and Lea48 and included clay removal, reductive, oxidative, alkaline chelation (with DTPA solution) and weak acid leaching steps with slight modifications15 from Pena et al.47 and Lea and Boyle49. These modifications included rinsing samples with NH 4 OH (ref. 49) instead of using 0.01 N NaOH (ref. 48) as a first step to remove the DTPA solution, followed by rinsing the samples three times with cold (room temperature) MilliQ water, 5-min immersion in hot (∼80 °C) MilliQ water and two more rinses with cold MilliQ water47. After cleaning, the samples were dissolved in 2% HNO 3 and finally analysed by iCAPQ Inductively-Coupled Plasma Mass Spectrometry at LDEO. Based on repeated measurements of in-house standard solutions, the long-term precision is <1.4, 1.9 and 2.1% for Mg/Ca, B/Ca and Cd/Ca, respectively. Five samples were split after clay removal, reduction and oxidation steps; one half was cleaned by the full cleaning procedure, while the alkaline chelation step was omitted for the other half. This approach was applied to test the influence of the chelation step on Cd/Ca and B/Ca. The results with and without the alkaline chelation show an average difference of 0.0003 μmol mol−1 and 5 μmol mol−1 for Cd/Ca and B/Ca, respectively (Supplementary Table 2). The Mg/Ca values from the two cleaning methods are comparable, but two samples showed a significant decrease in Mg/Ca, Fe/Ca, Mn/Ca and Al/Ca values when the alkaline chelation step was applied (Supplementary Table 2). This might be due to a more efficient removal of contaminants that are rich in Mg, but not in Cd or B. All our Mn/Ca values from the full cleaning method are <105 μmol mol−1, indicating that our results are unlikely affected by diagenetic coatings50. Only minor/trace element results from the full cleaning method were used in this study. All cleaning and loading steps for boron isotope and minor/trace element analyses were done in boron-free filtered laminar flow benches and all used boron-free Milli-Q water.

Stable isotope analyses

Pristine specimens of the benthic foraminifera Melonis barleeanus (∼30 specimens, size fraction 150–315 μm) and the planktic foraminifera N. pachyderma (∼50 specimens, size fraction 150–250 μm) were picked for stable isotope analyses. The stable oxygen and carbon isotope analyses were performed using a Finnigan MAT 251 mass spectrometer with an automated carbonate preparation device at MARUM, University of Bremen. The external standard errors for the oxygen and carbon isotope analyses are ±0.07‰ and ±0.05‰, respectively. Values are reported relative to the Vienna Pee Dee Belemnite (VPDB), calibrated by using the National Bureau of Standards (NBS) 18, 19 and 20. The oxygen isotope data were previously presented14,15,16, while the carbon isotope results are presented here for the first time (Supplementary Data 1).

Salinity and temperature reconstructions

We used the calcification temperature and δ18O SW values from Ezat et al.15 based on parallel δ18O and Mg/Ca measurements in N. pachyderma (Supplementary Data 1). Previous studies suggested that carbonate chemistry may exert a significant secondary effect on Mg/Ca in N. pachyderma20. The possible influence of secondary factors on temperature reconstructions are discussed in detail in Ezat et al.15 In brief, the main effect of the secondary factors appears to be the elevated pH and carbonate ion concentration during the LGM; a correction for this effect may lower the temperatures by 0–2 °C. However, the exact effect remains uncertain15. Here we used the temperature and δ18O SW reconstructions with no correction for non-temperature factors on Mg/Ca (see section ‘Propagation of error’ below).

In the absence of a direct proxy for salinity, we estimated the salinity from our reconstructed δ18O SW . There is a quasi-linear regional relationship between salinity and δ18O SW in the modern ocean, as both parameters co-vary because of addition/removal of freshwater51. However, temporal changes in the δ18O SW composition of freshwater sources and/or their relative contribution to a specific region, as well as changes in ocean circulation complicate using a local modern δ18O SW -salinity relationship to infer past changes in salinity. We therefore estimate salinity using the δ18O SW -salinity mixing line from the Norwegian Sea51 for the Holocene and the Eemian, when the hydrological cycle and ocean circulation were likely similar to modern. For the deglacial and last glacial periods, we use the δ18O SW -salinity mixing line52 based on data from the Kangerdlugssuaq Fjord, East Greenland, where the dominant source of freshwater is glacial meltwater from tidewater glaciers with δ18O SW values ranging from −30 to −20‰. These conditions are probably more representative of the sources of glacial meltwater during deglacial and glacial times53. Our salinity estimates during the deglacial and last glacial periods would have been ∼1.5‰ lower if we had used the modern δ18O SW -salinity mixing line from the Norwegian Sea. Although this salinity difference may appear large, it has little consequence for our pH and pCO 2 reconstructions and our conclusions (see ‘Sensitivity tests’ below).

pH and pCO 2 estimations

The boron isotopic composition of biogenic carbonate is sensitive to seawater-pH (ref. 17), because the relative abundance and isotopic composition of the two dominant dissolved boron species in seawater, boric acid [B(OH) 3 ] and borate [B(OH) 4 −] changes with pH (ref. 54), and borate is the species predominantly incorporated into marine carbonates. Culture experiments with planktic foraminifera provide empirical support for using their boron isotopic composition as a pH proxy30,55,56, but species-specific δ11B offsets are also observed, which are widely ascribed to ‘vital effects’57.

Linear regressions of δ11B CaCO3 versus δ11B borate relationships allow to infer δ11B borate from δ11B CaCO3 (ref. 30) as follows:

where ‘c’ is the intercept and ‘m’ is the slope of the regression. pH can then be estimated from foraminiferal δ11B-based δ11B borate using the following equation17:

where pK B is the equilibrium constant for the dissociation of boric acid for a given temperature and salinity58, δ11B SW is the δ11B of seawater (modern δ11B SW =39.61‰; ref. 59), and α (B3-B4) is the fractionation factor for aqueous boron isotope exchange between boric acid and borate. Klochko et al.54 determined the boron isotope fractionation factor in seawater α (B3-B4) = 1.0272±0.0006.

Because δ11B in the symbiont-barren N. pachyderma has so far only been calibrated from core top sediments, with large uncertainties and over a very limited natural pH range32, the pH sensitivity of this species is uncertain. However, we can use evidence from other calibrated symbiont-barren planktic foraminifera species to further constrain the pH sensitivity of this species. Martínez-Botí et al.60 suggested a pH sensitivity for the symbiont-barren planktic foraminifera G. bulloides similar to values predicted from aqueous boron isotope fractionation (that is, slope m in eq. (1) =1.074). We therefore used a slope value of 1.074 in equation (1). In addition, we calculated the intercept c=2.053‰ in equation (1) for N . pachyderma by calibrating our core top foraminiferal δ11B to a calculated pre-industrial pH (that is, δ11B borate ). Pre-industrial pH was estimated from modern hydrographic carbonate data (total Dissolved Inorganic Carbon ‘DIC’, total alkalinity, phosphate, silicate, temperature, salinity; ref. 21) from the southern Norwegian Sea (Fig. 2a, Supplementary Fig. 4), and subtracting 50 μmol kg−1 from DIC (ref. 61) to correct for the anthropogenic CO 2 effect. We used the hydrographic data collected during June 2002 and from the 22nd of September to the 13th of October 2003 (that is, within the assumed calcification season of N. pachyderma; refs 19, 20) and at our assumed calcification depth (that is, 40–120 m). This approach allows us to determine δ11B borate from δ11B CaCO3 (equation 1), which can then be used to calculate pH based on equation 2.

Although the slope determined for G. bulloides60 is similar to the coretop calibration of N. pachyderma32, neither calibration encompasses a wide pH range, and the uncertainty of the slopes is therefore large. In contrast, laboratory culture experiments with (symbiont-bearing) planktic foraminifera cover a much wider pH-range but display a lesser pH sensitivity (slope in equation (1)=∼0.7) than predicted from aqueous boron isotope fractionation30,55,56. However, this difference in slope has little consequence for our pH and pCO 2 reconstructions. A sensitivity test using slopes m=1.074 (ref. 60) and m=0.7 (refs 30, 55, 56) shows little difference between the two estimates (see section ‘Sensitivity tests’ below).

If two of the six carbonate parameters (total Dissolved Inorganic Carbon (DIC), total alkalinity, carbonate ion concentration, bicarbonate ion concentration, pH and CO 2 ), are known in addition to temperature, pressure and salinity, the other parameters can be calculated62. We used the modern local salinity-total alkalinity relationship (Alkalinity=69.127 × Salinity−116.42, R2=0.76, ref. 21) to estimate total alkalinity. Because weathering processes are slow and alkalinity is relatively high in the ocean, alkalinity can be considered a quasi-conservative tracer on these time scales, and we do not consider potential past changes in the salinity-total alkalinity relationship. Nonetheless, if we use the modern alkalinity-salinity relationship from the polar region as a possible analogue for our area during the last glacial, this would decrease the error in total alkalinity (because of the uncertainty in salinity) by up to 65 μmol kg−1 (Supplementary Fig. 5). Aqueous pCO 2 is then calculated using CO 2 sys.xls (ref. 63), with the equilibrium constants K 1 and K 2 from Millero et al.64, K SO4 is from Dickson59 and the seawater boron concentration from Lee et al.65

Sensitivity tests of pCO 2 reconstructions

Supplementary Fig. 6 shows that pH and pCO 2 reconstructions based on very different temperature, salinity and total alkalinity scenarios are very similar and do not significantly affect the large pCO 2 increases during HS1, HS4 and HS11. Because the intercept ‘c’ in the \({\rm \delta }^{{\rm 11}} {\rm B}_{{\rm CaCO}_{\rm 3} } \) versus δ11B borate calibrations (see Methods) is dependent on our choice of calcification depth for N. pachyderma, and corresponding selection of depths of hydrographic data to calculate the pre-industrial pH (after removing the anthropogenic carbon effect), we alternatively calculated the pre-industrial pH and the intercept ‘c’ based on hydrographic data from both 50 and 200 m water depths. This sensitivity test shows that the uncertainty in the calcification depth of N. pachyderma has insignificant effect on the amplitude of our down core pCO 2 variations (Supplementary Fig. 7).

In addition, to assess the uncertainty in our pH and pCO 2 estimations because of the uncertainty in the \({\rm \delta }^{{\rm 11}} {\rm B}_{{\rm CaCO}_{\rm 3} } \) versus pH sensitivity in N. pachyderma, we recalculated the δ11B borate using slope value of m=0.7 instead of m=1.074 in equation (1) as suggested for some symbiont-bearing planktic foraminifera species30,55,56, and re-adjusted the intercept ‘c’ accordingly (=−4.2‰). This test shows that the uncertainty in species-specific pH-sensitivity has no effect on our pCO 2 reconstructions for the Heinrich stadial events, while the main difference is an increase in the glacial/interglacial pCO 2 by ∼30 μatm, when a slope value of m=0.7 is used (Supplementary Fig. 8). This brings ΔpCO 2cal-air for the LGM to values of −30 μatm (and ΔpCO 2sea-air =−70 μatm), strengthening our conclusion about enhanced oceanic CO 2 uptake in our area during the LGM.

Finally, because our ΔpCO 2cal-air record can be biased because of errors in the age model especially for the Heinrich stadials (times with increasing atmospheric pCO 2 ), we performed a sensitivity study, in which 500 and 1,000 years were both added and subtracted from our age model (Supplementary Fig. 9). This arbitrary sensitivity study shows that such errors in the age model do not significantly affect the large increases in ΔpCO 2cal-air during HS1, HS4 and HS11 (Supplementary Fig. 9).

Error propagation in pCO 2 reconstructions

The uncertainty of each pCO 2 value in our record (Fig. 3c) is based on the propagated error of the effect of individual uncertainties in δ11B, calcification depth of N. pachyderma, temperature, salinity and total alkalinity on the pH and pCO 2 calculations. The error propagation (2σ) was calculated as the square root of the sum of the squared individual uncertainties. Note that total alkalinity has no effect on the pH estimations; it only affects the pCO 2 calculations.

The analytical uncertainty in δ11B ranges from ±0.22 to ±0.43‰, which translates to ∼±10 to ±40 μatm in pCO 2 . The error in pCO 2 due to the uncertainty in the calcification depth of N. pachyderma is equal to ±11 μatm on average (see previous Section and Supplementary Fig. 6). The uncertainty in salinity due to the choice of different salinity-δ18O SW mixing models for the last glacial period and the deglaciation is ∼±1.5‰, which translates to ∼±4 μatm pCO 2 . The error in total alkalinity due to the uncertainty in salinity estimations is up to ±100 μmol kg−1, which is equivalent to ∼±9 μatm pCO 2 .

For the assessment of uncertainty in our temperature estimates, one should ideally consider uncertainties associated with empirical calibrations and other non-temperature factors that affect Mg/Ca in N. pachyderma. Because the sensitivity of Mg/Ca in N. pachyderma to factors other than temperature (for example, carbonate chemistry) is not known20, we only include an error of ±0.7 °C, based on the calibration and analytical uncertainties of Mg/Ca (see ref. 15). This uncertainty translates to ±7 μatm pCO 2 on average. Ezat et al.15 discussed that the correction for elevated carbonate ion concentration during the LGM on Mg/Ca may lower the LGM temperature by 0–2 °C; however, the exact effect is very uncertain. A decrease in LGM temperatures would decrease our reconstructed pCO 2 values (∼−10 μatm decrease per 1 °C decrease), strengthening our conclusion that our study region was an intense area for CO 2 uptake at that time.

Data availability

The data generated and analysed during the current study are available along the online version of this article at the publisher’s web-site.