CO 2 system calculations

For both the preindustrial ocean and down-core CO 2 system calculations, seawater carbonate system variables were calculated using the CO 2 sys.xls program28 with dissociation constants K 1 and K 2 according to Mehrbach et al.55 and \({K}_{{{\mathrm{SO}}_4}}\) according to Dickson56. Seawater total boron concentration was calculated from the boron–salinity relationship of Lee et al.57. For the GLODAP dataset, the anthropogenic CO 2 contribution was subtracted from the measured DIC to obtain preindustrial DIC values2.

Preindustrial Atlantic DIC as and [CO 3 2−] as

The GLODAP dataset2 is used to calculate preindustrial ocean CO 2 system variables. Following the established method of Broecker and Peng3, we account for DIC anomalies created by (1) freshwater addition or removal based on S, (2) soft-tissue carbon creation and respiration based on PO 4 , and (3) CaCO 3 formation and dissolution based on ALK and nitrate (NO 3 ). See Fig. 1 for the simplified concept. We adopt the term DIC as to represent net air–sea exchange component DIC signatures from:

$$ {\hskip-10pt}{{\mathrm{DIC}}_{\rm{as}}} = {{\mathrm{DIC}}_{\rm{s}}} - ({{\mathrm{PO}}_{\rm{4s}}} - {{{\mathrm{{PO}}}_4}^{\rm{mo}}}) \times {{\mathrm{C/PO}}_4}\\ - {\ }^{1}{\hskip-2.5pt} {\hbox{/}}{\ }_{{\hskip-5pt} 2}\times \left( {{\mathrm{ALK}}_{\mathrm{s}}-{\mathrm{ALK}}^{{\mathrm{mo}}} + {\mathrm{NO}}_{{\mathrm{3s}}}-{{\mathrm{{NO}}}_3}^{\rm{mo}}} \right)-{\mathrm{DIC}}_{{\mathrm{constant}}}$$ (1)

where the subscript “s” represents values normalized to S of 35 (e.g., DIC s = DIC × 35/S); the superscript “mo” denotes mean ocean values at S = 35 (PO 4 mo = 2.2 μmol/kg, ALKmo = 2383 μmol/kg, DICmo = 2267 μmol/kg, and NO 3 mo = 31 μmol/kg)29; C/PO 4 represents the soft-tissue stoichiometric Redfield ratio; and the arbitrary DIC constant ( = 2285 μmol/kg) is designed to bring zero DIC as close to the NADW–AABW boundary (Fig. 2). The term (PO 4s − PO 4 mo) × C/PO 4 corrects for DIC changes due to photosynthesis and soft-tissue degradation, and the term ½ × (ALK s − ALKmo + NO 3s − NO 3 mo) accounts for DIC changes caused by CaCO 3 formation and dissolution. To be consistent with previous work3,30, we used C/PO 4 = 127 to calculate DIC as and [CO 3 2−] as in Fig. 2. Using other C/PO 4 values10 does not significantly affect spatial DIC as and [CO 3 2−] as patterns (Supplementary Figs. 11 and 12). Neither are their patterns affected by using other PO 4 –ALK–NO 3 values to replace global mean values in Eq. (1) (Supplementary Figs. 13 and 14). Ideally, DIC constant would be the mean DIC value of an abiotic ocean (Fig. 1), but this value cannot be simply determined from modern observations. Because our interest lies in spatial DIC as contrasts instead of absolute values, the choice of DIC constant has no effect on our interpretation.

To obtain [CO 3 2−] as , we first calculate [CO 3 2−] PO4–T–S–P using (DIC as + DIC constant ), ALKmo, and PO 4 mo at T = 3 °C, S = 35, and P = 2500 dbar. [CO 3 2−] as is then calculated by [CO 3 2−] as = [CO 3 2−] PO4–T–S–P − [CO 3 2−] constant , where [CO 3 2−] constant ( = 78 μmol/kg, calculated using DIC constant and ALKmo) is designed to bring zero [CO 3 2−] as close to the NADW–AABW boundary. In essence, the [CO 3 2−] as distribution reflects the variation of [CO 3 2−] when normalized to the same PO 4 −T−S−P conditions.

CO 2 system sensitivities and calculation of [CO 3 2 −] Norm

Because the seawater CO 2 system is nonlinear, there is currently no simple way to derive these sensitivities based on CO 2 system equations16. We use GLODAP preindustrial data2 to calculate numerically [CO 3 2−] sensitivities to various physiochemical parameters. Use of LGM outputs from the LOVECLIM model58 yields comparable sensitivities. We first use hydrographic data, including T, S, P, DIC, ALK, PO 4 , and SiO 3 to calculate [CO 3 2−]. We then change S to 35‰ and other chemical concentrations proportionally. For example, ALK and DIC will change as follows:

$${\mathrm{ALK}}_{{\mathrm{s = 35}}}={\mathrm{ALK} \,\times 35/S,}\,{\mathrm{and}}$$ (2)

$${\mathrm{DIC}}_{{\mathrm{s = 35}}} = {\mathrm{DIC}} \; \times {{35/S}}{.}$$ (3)

We use S = 35‰, ALK S=35 , DIC S=35 , [PO 4 ] S=35 , and [SiO 3 ] S=35 along with hydrographic T and P to calculate [CO 3 2−] S=35 . The [CO 3 2−] to S sensitivity (Sen _S ) is calculated by:

$${\mathrm{Sen}}_{\_{\mathrm{S}}} = \left( {\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]-\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{S = 35}}}} \right)/\left( {{S}-{\mathrm{35}}} \right){.}$$ (4)

To estimate temperature effects, we calculate [CO 3 2−] S=35, T=3 °C using S = 35‰, ALK S=35 , DIC S=35 , [PO 4 ] S=35 , [SiO 3 ] S=35 , T = 3 °C, and hydrographic P. The sensitivity of [CO 3 2−] S=35 to temperature (Sen _T ) is defined by:

$${\mathrm{Sen}}_{\_{\mathrm{T}}}\,{\mathrm{ = }}\,\left( {\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{S = 35,}}\,{\mathrm{T}} = {\mathrm{3}}^ {\circ} {\mathrm{C}}}-\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{S}} = {\mathrm{35}}}} \right)/\left( {{\mathrm{3}}-{T}} \right){.}$$ (5)

Regarding pressure effects, we calculate [CO 3 2−] S=35, T=3 °C, P=2500 dbar using S = 35‰, ALK S=35 , DIC S=35 , [PO 4 ] S=35 , [SiO 3 ] S=35 , T = 3°C, and P = 2500 dbar. The sensitivity of [CO 3 2−] S=35, T=3°C to P (Sen _P ) is defined by:

$${\mathrm{Sen}}_{\_{\mathrm{P}}} = \,\left( {\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{S= 35, }}\,{\mathrm{T}} = {\mathrm{3}}^ {\circ} {\mathrm{C,}}\,{\mathrm{P}} = {\mathrm{2500}}\,{\mathrm{dbar}}}} \right. \\ \left. {-\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{S}} = {\mathrm{35,}}\,{\mathrm{T}} = {\mathrm{3}}^ {\circ} {\mathrm{C}}}} \right)/\left( {{\mathrm{2500}}-{P}} \right){\mathrm{ \times 100}}.$$ (6)

To estimate the influence on [CO 3 2−] from within-ocean ALK–DIC redistributions by biological processes, we assume a 0.1 μmol/kg increase in PO 4 (i.e., ΔPO 4 = 0.1 μmol/kg) due to biological respiration (photosynthesis has an opposite effect). The resultant ALK (ALK S=35+respiration ) and DIC (DIC S=35+respiration ) can then be calculated from:

$${\mathrm{ALK}}_{{\mathrm{S}} = {\mathrm{35}} + {\mathrm{respiration}}} = {\,} {\mathrm{ALK}}_{{\mathrm{s = 35}}} + {\mathrm{\Delta PO}}_{\mathrm{4}} \times {\mathrm{C/PO}}_{\mathrm{4}}\\ \div \, {R} \times {\mathrm{2}}-{\mathrm{\Delta PO}}_{\mathrm{4}} \times {\mathrm{N/PO}}_{\mathrm{4}}.$$ (7)

$${\mathrm{DIC}}_{{\mathrm{S}} = {\mathrm{35}} + {\mathrm{respiration}}} = {\mathrm{DIC}}_{{\mathrm{s}} = {\mathrm{35}}} + {\mathrm{\Delta PO}}_{{\mathrm{4}}} \times {\mathrm{C/PO}}_{\mathrm{4}} + {\mathrm{\Delta PO}}_{{\mathrm{4}}} \times {\mathrm{C/PO}}_{\mathrm{4}} \div {R}{.}$$ (8)

Resultant [CO 3 2−] ([CO 3 2−] Norm+respiration ) values are calculated using DIC S=35+respiration , ALK S=35+respiration , and ([PO 4 ] S=35 + ΔPO 4 ) at constant physical conditions of T = 3 °C, S = 35, and P = 2500 dbar. The sensitivity of [CO 3 2−] Norm to PO 4 is defined by:

$$ \left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}}}{\mathrm{/PO}}_{\mathrm{4}}\,{\mathrm{sensitivity}} = \\ \left( {\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}} + {\mathrm{respiration}}}{\mathrm{- }}\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}}}} \right)/{\mathrm{\Delta PO}}_{\mathrm{4}}.$$ (9)

We consider four Redfield stoichiometric scenarios: C/PO 4 = 127, R = 4 (the reference composition; Fig. 4d); C/PO 4 = 140, R = 4; C/PO 4 = 127, R = 8; and C/PO 4 = 140, R = 8 (Supplementary Fig. 15). In all cases, strong exponential correlations exist between [CO 3 2−] Norm /PO 4 sensitivities and [CO 3 2−] Norm . The correlations may reflect the buffering effect of the seawater CO 2 system: for seawater with high DIC (low [CO 3 2−] and high buffering capability), [CO 3 2−] would be relatively less sensitive to biological DIC and ALK disturbances. All of the above sensitivity calculations assume no net air–sea CO 2 change.

To calculate air–sea exchange sensitivities, we assume a 10 μmol/kg increase in DIC S=35 due to atmospheric CO 2 invasion (i.e., ΔDIC as = 10 μmol/kg). We calculate [CO 3 2−] Norm+as using S = 35‰, ALK S=35 , DIC S=35+as ( = DIC S=35 + ΔDIC as ), [PO 4 ] S=35 , [SiO 3 ] S=35 , T = 3 °C, and P = 2500 dbar. The sensitivity of [CO 3 2−] as to DIC as is defined by:

$$ \left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{as}}}{\mathrm{/DIC}}_{{\mathrm{as}}}\,{\mathrm{sensitivity}} = \\ \left( {\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}} + {\mathrm{as}}}-\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}}}} \right)/\Delta {\mathrm{DIC}}_{{\mathrm{as}}}.$$ (10)

Using sensitivities shown in Fig. 4, [CO 3 2−] Norm can be calculated by:

$$\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right]_{{\mathrm{Norm}}} \,= \,\left[ {{{\mathrm{CO}}_{\mathrm{3}}}^{{\mathrm{2 - }}}} \right] + \left( {{\mathrm{35}}-{S}} \right) \times {{\mathrm{Sen}}_{\_{\mathrm{S}}}} + \left( {{\mathrm{3}}-{T}} \right)\\ \times \, {{\mathrm{Sen}}_{\_{\mathrm{T}}}} + \left( {{\mathrm{2500}}-{P}} \right){\mathrm{/100}} {\hskip1pt}\times {{\mathrm{Sen}}_{\_{\mathrm{P}}}}.$$ (11)

Excel spreadsheets are provided in Supplementary Data 7–8 to calculate [CO 3 2−] Norm and the biological curves shown in Fig. 5.

LGM–Holocene North Atlantic carbon budget

The total extra carbon increase (Δ∑C LGM−Holocene ) in Fig. 6 is calculated by Δ∑C LGM−Holocene = V × density × %GNAIW × ([CO 3 2−] as_ODP999-BOFS LGM/0.61) × 12 − V × density × %NADW × ([CO 3 2−] as_ODP999-BOFS Holocene/0.59) × 12, where V is the global deep ocean volume (>1 km water depth) at 100.8 × 1016 m3, density = 1027.8 kg/m3 (ref. 29), %GNAIW and %NADW, respectively, represent their volume fractions in the deep ocean, [CO 3 2−] as_ODP999-BOFS Holocene = 56 μmol/kg, [CO 3 2−] as_ODP999-BOFS LGM = 114 μmol/kg (Fig. 5), terms 0.61 and 0.58, respectively, represent the absolute LGM and Holocene [CO 3 2−] as /DIC as sensitivities (Fig. 4e) used to transfer [CO 3 2−] as_ODP999-BOFS into ODP999–BOFS DIC as contrasts (LGM: 186 μmol/kg; Holocene: 95 μmol/kg), and the number 12 converts C from moles into weight. Based on previous estimates, %NADW is thought to be ~50% (refs. 38,39), while %GNAIW remained roughly similar to %NADW or shrank (refs. 35,36). These estimates are debated and have large uncertainties, and we thus calculate Δ∑C LGM−Holocene for a range of %NADW and %GNAIW values (Fig. 6). Any influence from AAIW is ignored because of its similar [CO 3 2−] as signals to Gulf Stream during the Holocene (Supplementary Fig. 3) and much reduced northward advection during the LGM23,31,32,33. We tentatively treat Δ∑C LGM-Holocene of ~100 PgC using %NADW = 50% and %GNAIW = 30% as our best estimate. Assuming no Holocene–LGM DIC as gradient change (i.e., the same CO 2 uptake efficiency) and everything else being equal, Δ∑C LGM–Holocene would be −240 PgC at %NADW = 50% and %GNAIW = 30%.

Cores, age models, samples, and analytical methods

We used ODP Site 999 for Gulf Stream surface-water reconstructions (Fig. 2). The age model is from Schmidt et al.59. Planktonic foraminiferal Globigerinoides ruber (sensu stricto, white variety) δ18O, Mg/Ca, and δ11B data are from refs. 21,22,59. Briefly, about 25 and 55 shells from the 250–350 μm size fraction were used for δ18O and Mg/Ca analyses, respectively. Samples for δ18O analyses were sonicated in methanol for 5–10 s, roasted under vacuum at 375 oC for 30 min, and analyzed on a Fisons Optima IRMS with a precision of <0.06‰. Shells for Mg/Ca were cleaned following the reductive cleaning procedure60 and measured on an inductively-coupled plasma mass spectrometer (ICP–MS) with a precision of ~1.7%. For δ11B analyses, about 100–120 G. ruber (w) shells from the 300–355 μm size fraction were cleaned following the “Mg-cleaning” procedure61, to minimize material loss during cleaning62. G. ruber (w) δ11B was measured on a Neptune multicollector (MC)–ICP–MS with an analytical error in δ11B of about ±0.25‰ (ref. 21).

Three cores (BOFS 17, BOFS 11, and BOFS 14 K) from the polar North Atlantic Ocean are used for deep-water reconstructions (Fig. 3). Their age models are based on published chronologies24,63,64,65. For each sample (~2 cm thickness), ~10–20 cm3 of sediment was disaggregated in de-ionized water and was wet sieved through 63 μm sieves. To facilitate analyses, we picked the most abundant species for measurements. For each B/Ca analysis, ~10–20 monospecific shells of the benthic foraminifera C. mundulus (BOFS 17 K) and C. wuellerstorfi (BOFS 14, 11 K) were obtained from 250 to 500 μm size fraction. The shells were double checked under a microscope before crushing to ensure that consistent morphologies were used throughout the core. On average, following this careful screening the starting material for each sample was ~8–12 shells, which is equivalent to ~300–600 μg of carbonate. For benthic B/Ca analyses, foraminiferal shells were cleaned with either the “Mg-cleaning” method61 or the “Cd-cleaning” protocol61, to investigate cleaning effects on trace element/Ca in foraminiferal shells62,66. No discernable B/Ca difference is observed between the two cleaning methods25,62. Benthic B/Ca ratios were measured on an ICP–MS using procedures outlined in ref. 67, with an analytical error better than ~5%.

For each benthic Cd/Ca analysis, ~10–20 shells of the benthic foraminiferal taxa C. mundulus (BOFS 17 K), C. wuellerstorfi (BOFS 14 K, 11 K), and Uvigerina spp. (BOFS 17 K) were picked from the 250–500 μm size fraction. Previous studies26,27,68 showed similar Cd/Ca ratios between infaunal Uvigerina spp. and epifaunal Cibicidoides, and we thus combined Cd/Ca data from these taxa to obtain continuous downcore PO 4 records. We used the “Cd-cleaning” method60,69 to clean benthic shells for Cd/Ca measurements. Cd/Ca ratios were measured on an ICP–MS with an analytical error better than ~5% (ref. 67)

For δ11B measurements, about 20 benthic shells from the 250–500 μm size fraction were picked for each sample. Shells used for δ11B analyses were cleaned using the “Mg-cleaning” method, to minimize loss of shell material61. After cleaning, shells were dissolved and pure boron was extracted using column chemistry as described by Foster21. Benthic δ11B was measured on a Neptune multi-collector (MC)–ICP–MS following ref. 21. The analytical error in δ11B is about ± 0.25‰. Due to the relatively large sample size requirement, shell availability, and lengthy chemical treatments for δ11B, we present low-resolution δ11B for C. mundulus from BOFS 17 K and for C. wuellerstorfi from BOFS 11 K. Note that consistent [CO 3 2−] results from B/Ca and δ11B strengthen the reliability of our reconstructions (Fig. 3).

Published benthic Cd/Ca and B/Ca results are included in Fig. 3. Altogether, we generated 180 new measurements of benthic δ11B, B/Ca, and Cd/Ca. All data are listed in Supplementary Data 1–9.

ODP 999 reconstructions

ODP Site 999 was used to constrain past physical conditions and carbonate chemistry of the Gulf Stream (Supplementary Fig. 4). Following previous approaches21,22, surface water temperature (T surface ) and salinity (S surface ) were estimated based on G. ruber Mg/Ca (ref. 59) and sea level changes21,22,59, respectively. We first convert G. ruber δ11B to borate δ11B (δ11B borate ), following the conversion method of ref. 22. Surface water pH (pH surface ) was calculated from seawater δ11B borate along with T surface and S surface . To constrain the CO 2 system, two CO 2 system variables are necessary16. In addition to δ11B-derived pH, literature studies21,22,41 generally estimate past surface-water ALK (ALK surface ) changes. Following refs. 21,22, we estimate ALK surface from S surface using the modern S surface –ALK surface relationship (ALK surface = 59.19 × S surface + 229.08, R2 = 0.99)21. Together with T surface and S surface , pH surface , and ALK surface were used to calculate other CO 2 system variables including surface-water [CO 3 2−] ([CO 3 2−] surface ) and DIC (DIC surface ) using the CO 2 sys program28. Surface-water PO 4 concentration at ODP 999 is assumed to be zero over the last 27 ka.

Following refs. 21,22,59, errors are estimated to be 1 °C, 1‰, 100 μmol/kg, and ~0.43‰ for T surface , S surface , ALK surface , and δ11B borate , respectively. Integrated average uncertainties in [CO 3 2−] surface and DIC surface for a single reconstruction are, respectively, ~20 (Holocene: ~18, LGM: ~24) and ~90 μmol/kg, based on quadratic addition of all individual errors sourced from T surface ([CO 3 2−] surface : 2 μmol/kg, DIC surface : 3 μmol/kg), S surface ([CO 3 2−] surface : 2 μmol/kg, DIC surface : 5 μmol/kg), ALK surface ([CO 3 2−] surface : 14 μmol/kg, DIC surface : 86 μmol/kg), and δ11B borate ([CO 3 2−] surface : 16 μmol/kg, DIC surface : 24 μmol/kg; note that δ11B borate leads to an error in [CO 3 2−] via pH). Uncertainties for calculated CO 2 system variables at ODP 999 are tabulated in Supplementary Data 1. Use of other methods to estimate ALK would have little impact on our conclusions (Supplementary Figs. 16 and 17).

From pH to [CO 3 2−]

For palaeo-studies, surface-water pH is generally obtained from planktonic foraminiferal δ11B. To calculate [CO 3 2−], a second CO 2 system variable is needed16. Following the previous approach21,22, past ALK surface at ODP 999 have been estimated from S using the S surface –ALK surface relationship. Due to limited knowledge about the past S surface –ALK surface relationship, a generous uncertainty has been assigned to ALK surface at ±100 μmol/kg (ref. 21,22), which is about half of the entire ALK range in the present global ocean2. Using ALK surface and pH surface along with T surface and S surface , [CO 3 2−] surface and DIC surface can be calculated using the CO 2 sys program28. Because of the large uncertainty in ALK surface , large errors in DIC surface might be expected (Supplementary Fig. 4). However, given the constraint from pH surface , seawater ALK surface and DIC surface variations are not random but must vary systematically within ALK–DIC space (Supplementary Fig. 5). Because of the close relationship between pH and [CO 3 2−] (i.e., roughly parallel patterns of pH and [CO 3 2−] within ALK−DIC space; Supplementary Fig. 5), this systematic ALK−DIC variation allows us to confine [CO 3 2−] with acceptable uncertainty. For a given pH at ODP 999, an error of 100 μmol/kg in ALK only leads to an error of about ±14 μmol/kg in [CO 3 2−] (Supplementary Fig. 5).

For clarity, Supplementary Fig. 5a, b only consider the effect of ALK errors on [CO 3 2−] estimates assuming constant pH and T–S–P conditions. To fully propagate errors from various sources including T surface , S surface , ALK surface , and pH surface , we use a Monte Carlo approach (n = 10,000) to calculate the integrated error in [CO 3 2−] (ref. 70). As can be seen from Supplementary Fig. 5c–f, the final errors (~20–25 μmol/kg) in an individual [CO 3 2−] reconstruction based on the Monte-Carlo are similar to those (~18–24 μmol/kg) based on quadratic addition of individual errors, justifying our major error estimation approach (i.e., quadratic addition).

Subtropical western North Atlantic surface [CO 3 2−]

Because most of North Atlantic subtropical gyre waters circulate through the Caribbean Sea before being transported to the subpolar North Atlantic via the Gulf Stream, ODP 999 from Caribbean Sea is used to constrain past Gulf Stream carbonate chemistry20. To further test the feasibility of using ODP 999 to represent the first-order Gulf Stream [CO 3 2−] changes during the Holocene and LGM, we have estimated surface-water [CO 3 2−] for four sites from the wider subtropical western Atlantic region (latitude: 12–33°N, longitude: 61–91°W). Among these sites, KNR140–51GGC (33°N, 76°W) is located within the Gulf Stream today71. Because subtropical surface waters cycle multiple times through the upper ocean gyre circulations, it is possible that surface waters have been close to equilibrium with past atmospheric pCO 2 (refs. 21,22). Therefore, we assume surface-water pCO 2 of 270 and 194 ppm for the Holocene and LGM, respectively72. We assign a ±15 ppm error to surface-water pCO 2 to account for any potential air–sea CO 2 disequilibrium. For these sites, we use surface temperature and salinity reconstructions from previous publications71,73,74,75. ALK is calculated based on the same approach for ODP 999. The reconstructed in situ [CO 3 2−] values show some differences between cores, due to local T–S conditions. Since we are interested in air–sea CO 2 exchange signals, we convert reconstructed in situ [CO 3 2−] into [CO 3 2−] Norm using Eq. (11). As can be seen from Supplementary Fig. 6 and Supplementary Data 2, these cores show similar [CO 3 2−] Norm values for the Holocene (~260 μmol/kg) and LGM (~300 μmol/kg) as ODP 999. Therefore, we argue that ODP 999 sufficiently records first-order Gulf Stream air–sea exchange carbonate chemistry for the Holocene and LGM. Because we aim to obtain a proxy-based estimates, we use ODP 999 data for calculations in the main text.

Benthic B/Ca and δ11B to deep-water [CO 3 2−]

Most deep-water [CO 3 2−] values are reconstructed using benthic B/Ca (refs. 25,47) from [CO 3 2−] downcore = [CO 3 2−] PI + ΔB/Ca downcore–coretop /k, where [CO 3 2−] PI is the preindustrial (PI) deep-water [CO 3 2−] value estimated from the GLODAP dataset2, ΔB/Ca downcore–coretop represents the deviation of B/Ca of down-core samples from the core-top value, and k is the B/Ca–[CO 3 2−] sensitivity of C. wuellerstorfi (1.14 μmol/mol per μmol/kg) or C. mundulus (0.69 μmol/mol per μmol/kg)25. We use a reconstruction uncertainty of ±10 μmol/kg in [CO 3 2−] based on global core-top calibration samples25,76.

For cores BOFS 17 K and BOFS 11 K, new monospecific epifaunal benthic δ11B values were converted into deep-water [CO 3 2−] following the approach detailed in ref. 77. Briefly, benthic δ11B is assumed to directly reflect deep-water borate δ11B, as suggested by previous core-top calibration work78. Deep-water pH is calculated using benthic δ11B along with T deep and S deep , similar to the approach to calculate surface-water pH at ODP 999 (refs. 21,22). We assume constant ALK at the studied sites (2313 μmol/kg at BOFS 17 K and 2310 μmol/kg at BOFS 11 K) in the past. Following ref. 77, a generous error of 100 μmol/kg is assigned to ALK estimates. We then calculate deep-water [CO 3 2−] from pH and ALK using the CO 2 sys program28. The integrated average uncertainty in deep-water [CO 3 2−] is ~±10 μmol/kg, based on quadratic addition of individual errors of ~±2 μmol/kg sourced from T deep (±1 °C), ~±2 μmol/kg from S deep (±1‰), ~±5 μmol/kg from ALK (±100 μmol/kg), and ~±8 μmol/kg from δ11B borate (~±0.25‰). As demonstrated by Supplementary Fig. 5, the large ALK error only contributes a small uncertainty to the final [CO 3 2−] estimate. As shown in Fig. 3, benthic B/Ca and δ11B yield consistent deep-water [CO 3 2−] reconstructions for the Holocene and LGM.

Benthic Cd/Ca to deep-water PO 4

We follow the established approach26,46,79 to convert benthic (C. wuellerstorfi, C. mundulus, and Uvigerina spp.) foraminiferal Cd/Ca into deep-water Cd concentrations. Partition coefficients (D Cd ) are used to calculate deep water Cd from: Cd (nmol/kg) = [(Cd/Ca) foram /D Cd ] × 10. Bertram et al.65 used empirical D Cd values of 2.3, 2.2, and 2.7 for BOFS 17, 14, and 11 K, respectively. However, these D Cd values would result in Holocene Cd of 0.3–0.4 nmol/kg, higher than the observed value of ~0.25 nmol/kg from modern hydrographic measurements (Supplementary Fig. 7)80. This offset may suggest higher D Cd values for the North Atlantic Ocean, which has been acknowledged recently81. We thus adjust D Cd (~25% increase) so that the calculated Holocene deep-water Cd concentrations match modern measurements. This adjustment is supported by consistent Cd reconstructions from this study and previous reconstructions based on Cd/Ca measurements for Hoeglundina elegans. Compared to Cibicidoides, D Cd into H. elegans is far less variable79. As can be seen from Supplementary Fig. 8, for cores with similar benthic δ13C from similar water depths (i.e., bathed in similar water masses), our Cd reconstructions match favorably with those based on H. elegans measurements82. Deep water Cd is converted into PO 4 using the relationship based on the latest North Atlantic Ocean measurements (Supplementary Fig. 7)80. Using older published Cd–PO 4 relationships26,83 only marginally affects our PO 4 estimates.

Uncertainties associated with Cd and PO 4 reconstructions are estimated as follows. Error for Cd is estimated using 2σ Cd = \(\sqrt{(2\sigma _{D_{\mathrm{Cd}}})^2 + (2\sigma _{{\mathrm{Cd}}/{\mathrm{Ca}}})^2}\), where \(2\sigma _{D_{\mathrm{Cd}}}\) and \(2\sigma _{\mathrm{Cd}/\mathrm{Ca}}\) (=5%) are errors for D Cd and Cd/Ca, respectively. Due to poorly defined uncertainty for D Cd from the literature, we assume an error of 50%, and then compare our final errors with literature estimates to assess the appropriateness of our calculations. Seawater PO 4 is calculated from Cd using: PO 4 = \(\frac{{\mathrm{Cd} - b \pm (2\sigma _b)}}{{a \pm (2\sigma _a)}}\), where 2σ a and 2σ b , respectively, represent 95% confidence errors associated with a and b (Supplementary Fig. 7b). The PO 4 uncertainty was calculated from: \(2\sigma _{{{\rm{PO}}_4}} = \sqrt {\left( {\partial _{{{\rm{PO}}_4}}{\mathrm{/}}\partial _a \cdot 2\sigma _a} \right)^2 + \left( {\partial _{{{\rm{PO}}_4}}{\mathrm{/}}\partial _b \cdot 2\sigma _b} \right)^2 + \left( {\partial _{{{\rm{PO}}_4}}/\partial _{\mathrm{Cd}} \cdot 2\sigma _{\mathrm{Cd}}} \right)^2}\), where \(\partial _{{{\rm{PO}}_4}}/\partial _a\) = \(\frac{{ - (\mathrm{Cd} - b)}}{{a^2}}\), \(\partial _{{{\rm{PO}}_4}}/\partial _b\) = \(\frac{{ - 1}}{a}\), and \(\partial_{{{\rm{PO}}_4}}/\partial _{\mathrm{Cd}}\) = \(\frac{1}{a}\). Our final errors on individual Cd and PO 4 are ~0.12 nmol/kg (~55%) and ~0.5 μmol/kg (~50%), respectively. When compared with previously published uncertainties (~0.08 nmol/kg for Cd and ~0.17 μmol/kg for PO 4 )46,68, our error estimates are possibly too generous. Here we use ~50% error to be conservative. We encourage future work to improve uncertainty estimates for the benthic Cd/Ca proxy.

The oceanic residence time of PO 4 is ~100,000 years84. The LGM deep ocean was possibly more reducing85, which might have facilitated sediment organic matter preservation, and, thus, PO 4 removal from the ocean. However, this effect might have been compensated by decreased organic burial on continental slopes due to shallower LGM sea levels86,87. Considering the short (~10,000 years) last deglacial84, we assume that global PO 4 and Cd reservoirs remained constant between the Holocene and LGM. Our reconstructions (Fig. 3) are consistent with high benthic δ13C and low benthic Cd/Ca at numerous glacial North Atlantic mid-depth sites23,31,46,65,88,89.

Deep-water temperature and salinity estimates

Deep-water temperature (T deep ) is estimated from the ice volume corrected benthic δ18O (δ18O IVC ) and the δ18O-temperature equation of Marchitto et al.90 from T deep = 2.5 − (δ18O IVC − 2.8)/0.224, where δ18O IVC = δ18O benthic − δ 18 O global_sealevel . δ18O global_sealevel was estimated from sea level curves86,87 with a global δ18O seawater −sea level scaling of 0.0085‰/m (ref. 91). Deep-water salinity (S deep ) is calculated by: S deep = S core_top + 1.11 × δ18O global_sealevel , where S core_top is the modern S deep (35.06, 34.926, and 34.893 at BOFS 17 , 11 , and 14 K, respectively2) and the term 1.11 is the scaling term for a global S−δ18O global_sealevel relationship29,91. We assume ±1 °C and ±1‰ uncertainties in T deep and S deep , respectively. Use of other methods to estimate T deep and S deep negligibly affects our conclusions, due to relatively weak sensitivities of [CO 3 2−] Norm to T and S changes (Fig. 4).

Uncertainties and statistical analyses

Uncertainties associated with [CO 3 2−] and PO 4 were evaluated using a Monte-Carlo approach92,93. Errors associated with the chronology (x-axis) and [CO 3 2−] and PO 4 reconstructions (y-axis) are considered during error propagation. Age errors are assumed to be ±3000 years for the three BOFS cores. Methods to calculate errors associated with individual [CO 3 2−] and PO 4 reconstructions (y-axis) are given above. All data points were sampled separately and randomly 5000 times within their chronological and [CO 3 2−] or PO 4 uncertainties and each iteration was then interpolated linearly. At each time step, the probability maximum and data distribution uncertainties of the 5000 iterations were assessed. Figure 3 shows probability maxima (bold curves) and ±95% (light gray; 2.5−97.5th percentile) probability intervals for the data distributions, including chronological and proxy uncertainties. For details, see refs. 92,93.

For a time period (e.g., Holocene) where multiple analyses are available, uncertainties are calculated following the method from ref. 94 by 2σ = \(\sqrt {[\mathop {\sum }

olimits_{i = 1}^n (2\sigma _i)^2]/n}\), where n is the number of reconstructions and 2σ i is the error associated with individual reconstruction. For [CO 3 2−] or [CO 3 2−] Norm offsets between the Holocene and LGM, 2σ = \(\sqrt {(2\sigma _{\mathrm{Holocene}})^2 + (2\sigma _{\mathrm{LGM}})^2}\), where \(2{\sigma} _{\mathrm{Holocene}}\) and \(2{\sigma} _{\mathrm{LGM}}\) are 2σ of Holocene and LGM values, respectively. Other methods (e.g., weighted mean)95 would give similar results.