Coupling of oceanic carbon and nitrogen facilitates spatially resolved quantitative reconstruction of nitrate inventories

Anthropogenic impacts are perturbing the global nitrogen cycle via warming effects and pollutant sources such as chemical fertilizers and burning of fossil fuels. Understanding controls on past nitrogen inventories might improve predictions for future global biogeochemical cycling. Here we show the quantitative reconstruction of deglacial bottom water nitrate concentrations from intermediate depths of the Peruvian upwelling region, using foraminiferal pore density. Deglacial nitrate concentrations correlate strongly with downcore δ13C, consistent with modern water column observations in the intermediate Pacific, facilitating the use of δ13C records as a paleo-nitrate-proxy at intermediate depths and suggesting that the carbon and nitrogen cycles were closely coupled throughout the last deglaciation in the Peruvian upwelling region. Combining the pore density and intermediate Pacific δ13C records shows an elevated nitrate inventory of >10% during the Last Glacial Maximum relative to the Holocene, consistent with a δ13C-based and δ15N-based 3D ocean biogeochemical model and previous box modeling studies.


Supplementary Note 1: The influence of microhabitat preferences on foraminiferal δ 13 C
Since U. peregrina lives shallow infaunal, an offset would be expected between δ 13 CFORAM in their tests and the δ 13 CDIC signature in the bottom water. This has been shown for U. peregrina specimens from the North Pacific and the Atlantic. The offset is strongly controlled by the gradient of δ 13 CDIC between bottom and pore waters 1 . The gradient of δ 13 CDIC between bottom and pore waters is controlled by bottom water oxygenation and diminishes with decreasing oxygen penetration depth into the sediments 2,3 . Since the Peruvian upwelling system is one of the most oxygen depleted regions in today´s ocean this gradient is supposed to be rather low at our sampling location. This might explain why there is no significant difference in the intercept between the downcore correlation of [NO3 -]BW and δ 13 CFORAM and [NO3 -] and δ 13 CDIC in the recent water masses caused by vital effects. Nevertheless, Virgulinella fragilis from the anoxic Namibian diatomaceous mud belt show an extreme carbon fractionation of -12.5‰ probably derived from a mixture of carbon sources by anaerobic methane oxidation and enhanced incorporation of metabolic CO2 derived from organic matter decomposition under oxygen depleted conditions 4 . For a comprehensive approach to reconstruct intermediate [NO3 -] using benthic δ 13 CFORAM, it should be considered to focus on epifaunal species or to do a local species specific calibration to account for possible regional offsets between δ 13 CFORAM and bottom water δ 13 CDIC due to vital effects or changes within the microhabitat.

Supplementary Note 2: Local fluctuations of [NO3 -]BW
Local [O2] fluctuations are directly influencing NO3loss processes at the Peruvian ODZ. As discussed in the main text, elevated [O2] during the LGM probably contributed to an increase in global [NO3 -] during the LGM. Since sediment core M77/2 52-2 is located in intermediate water depths well below the most oxygen depleted center of the ODZ a change in water column denitrification probably did not directly influence [NO3 -]BW at this site. Benthic denitrification at depths below the Peruvian ODZ is negligible due to the lack of bacterial activity and low abundances of denitrifying foraminifera [5][6][7] , probably related to the reduced Corg supply compared to the shelf sediments. The decrease in water column denitrification might locally lead to an increase in [NO3 -] but cannot explain a decrease in δ 13 CDIC which should be decoupled from denitrification. Although [O2] was probably increased during the LGM at the M77/2 52-2 site 8,9 , a local decrease in denitrification cannot alone explain the tendencies observed in our sediment record, which follow to a major part the global changes in the oceanic [NO3 -] inventory. Not all [NO3 -]BW fluctuations in the record of core M77/2 52-2 can be explained by changes in the global reactive N reservoirs, though. The most well-marked offset to the global model predictions appears during H1, when [NO3 -]BW depletes distinctively for ~4 kyrs ( fig. 1B). The most important sources of nutrients in the Eastern Tropical Pacific are high latitude intermediate and deep waters. In the modern North Pacific deep waters are highly enriched in NO3 -(see fig. 2A). During H1 a breakdown in stratification in the shallower North Pacific water masses has been observed 10 . This probably led to mixing of nutrient depleted δ 13  Supplementary Note 3: The influence of anthropogenic CO2 on δ 13 C Penetration of anthropogenic CO2 into the modern Pacific is already deeper than 2000 m 12,13 . The δ 13 C values of anthropogenic CO2 from the burning of fossil fuels, coal or natural gas are very low with a range from about -23.8‰ to -44.2‰ 14 . At some stations in the modern Pacific the δ 13 CDIC signature is relatively low and shows a distinctive offset from the correlation between δ 13 C and [NO3 -] which has been found in the modern Pacific as well as downcore over the last deglaciation. The deglacial correlation between δ 13 CFORAM and reconstructed [NO3 -]BW is not influenced by anthropogenic CO2. Since both linear regressions, for the modern Pacific and the downcore data (see Figure 2C and Supplementary Figure 4B), do not differ significantly from each other, the influence of anthropogenic CO2 onto the correlation in eq.1 seems to be only of minor influence, yet. Offsets from the downcore correlation (Eq. S1), e.g. the distribution of stations with very low δ 13 CDIC, might indeed be used to trace anthropogenic CO2 in the modern Pacific. which has been found by previous modeling studies 12,13 . This indicates that Supplementary Eq. 1 might indeed be used to trace anthropogenic CO2 and shows that a correlation which was even stable on glacial/interglacial timescales, although still consistent in most parts of the modern Pacific, already begins to break down at several locations which are strongly influenced by anthropogenic CO2.

Supplementary Note 4: Offsets between measurements and modeling of the modern δ 13 CDIC-[NO3 -]coupling
The model reproduces the modern observations in the intermediate Pacific Ocean reasonably well (R 2 =0.638 , see Supplementary Figure 7). However, the model predicted δ 13 CDIC-[NO3 -] slope (-0.066) is slightly lower than observed (-0.093) (Supplementary Figure 2). This underestimated slope is due to the too high δ 13 CDIC values in the deep North Pacific, the most notable bias in the modern model-data comparison (Supplementary Figure 7). This is caused by underestimated organic matter remineralization in the deep North Pacific that would further decrease δ 13 CDIC in the deep ocean. This slight systematic model bias likely applies to all time periods, which explains why the model predicts virtually no change to the δ 13 CDIC-[NO3 -] slope between the modern and glacial ocean, particularly in the [NO3 -] range >40µM where our core location exists (Supplementary Figure 2b). This model result supports that the different environmental condition of the glacial ocean did not significantly change the δ 13 CDIC -[NO3 -] slope and that the observed modern ocean relationship is applicable in the glacial ocean.
Supplementary Table 3 shows the model NO3sensitivity simulations when changing the amount of iron deposition to the Southern Ocean (LGM_lowSOFe, LGM_highSOFe).
LGM_control simulates the Southern Ocean δ 15 Nsed.org changes the best, suggesting it simulates the most realistic changes of enhanced remineralized [NO3 -] in the deep ocean and reduced preformed [NO3 -] in AAIW that affect our core location bottom water. The low and high Southern Ocean Fe simulation perform worse than the control simulation with respect to observations, but are within a reasonable uncertainty level considering the complexity of processes not directly tested here (see reference 15 for full discussion). These sensitivity simulations altered bottom water [NO3 -] at our core location by an additional ±0.7 µM in addition to the change to global [NO3 -], which can be considered as the uncertainty range of changes to our bottom water [NO3 -] related to Southern Ocean iron fertilization that may not be reflecting global [NO3 -] changes.

Supplementary Figures
Supplementary Figure 1 fig. 2C but not cut at a δ 13 C of -1‰.
Supplementary Figure 5: Offsets (Δδ 13 C) between δ 13 C from DIC in the modern Pacific and δ 13 C which would be predicted by Supplementary Eq. 1 for the measured [NO3 -] at the same station. The distribution and quantity of negative Δδ 13 C is directly reflecting the distribution of anthropogenic CO2 in the modern Pacific which has been found in previous studies 12,13 . The Ocean Data View software has been used to compile this plot 19 . Supplementary

Supplementary Tables
Supplementary Table 1: Downcore data for core M77/2 52-2. PD is the mean pore density of Bolivina spissa for this depth. N is the number of specimens used for the pore density analysis in this sample.
[NO3 -]BW was calculated from the pore density of Bolivina spissa using eq. 3. δ 18   LGM_control (2.44‰) simulates the most realistic iron fertilization effect there, which is the simulation our main analysis is based upon, with LGM_lowSOFe (1.58‰) and LGM_highSOFe (3.31‰) representing our low-and high-end uncertainty sensitivity simulations. Supplementary