Coupled climate and subarctic Pacific nutrient upwelling over the last 850,000 years

High latitude deep water upwelling has the potential to control global climate over glacial timescales through the biological pump and ocean-atmosphere CO2 exchange. However, there is currently a lack of continuous long nutrient upwelling records with which to assess this mechanism. Here we present geochemical proxy records for nutrient upwelling and glacial North Pacific Intermediate Water (GNPIW) formation in the Bering Sea over the past 850 kyr, which demonstrates that glacial periods were characterised by reduced nutrient upwelling, when global atmospheric CO2 and temperature were also lowered. We suggest that glacial expansion of sea ice in the Bering Sea, and the simultaneous expansion of low nutrient GNPIW, inhibited vertical mixing and nutrient supply across the subarctic Pacific Ocean. Our findings lend support to the suggestion that high latitude sea ice and the resultant intermediate water formation, modulated deep water upwelling and ocean-atmosphere CO2 exchange on glacial-interglacial timescales.

section_content: Recent studies have shown a sharp increase in both surface ocean subarctic Pacific productivity and CO2 ventilation during the last deglaciation (Galbraith et al., 2007; Gray et al., 2018), and suggested that the equivalent of ~30 ppm atmospheric CO2 was released into the atmosphere due to increased overturning (Rae et al., 2014). This was proposed to coincide with a collapse of glacial North Pacific Intermediate Water (GNPIW), allowing overturning to bring nutrient and CO2-rich deep water to the surface (Gray et al., 2018). Although a pervasive link between subarctic Pacific export productivity and Pleistocene glacial-interglacial cycles has been demonstrated (Jaccard et al., 2010 (Jaccard et al.,, 2005)), it is not yet known if this was primarily controlled by changes in iron supply (Galbraith et al., 2008; Praetorius et al., 2015), water column stratification and/or light limitation (Jaccard et al., 2010; Lam et al., 2013) and/or the extent of GNPIW (Knudson and Ravelo, 2015a; Li et al., 2017) and nutrient upwelling (Gray et al., 2018). If the upwelling of nutrient and CO2-rich water was the predominant control on subarctic Pacific productivity variations (Kender et al., 2018), this would imply the region has been an active source of pCO2 throughout the Pleistocene and not just the last deglacial, because they are coupled with ice core temperature records (Jaccard et al., 2010 (Jaccard et al.,, 2005)). One way to test for GNPIW and its effect on upwelling nutrients is from the Bering Sea, which today contains an active upwelling zone. 

Surface waters in the Bering Sea are fed by the subarctic Pacific Alaskan Stream, which enters through passes and straits in the Aleutian arc and circulates northward around the basin, forming the shelf-adjacent Bering Slope Current (BSC) (Stabeno et al., 1999) (Fig. 1 ). 

Turbulent eddies in the BSC extend below the surface waters at ~300 m water depth, causing upwelling of North Pacific Deep Water (NPDW), which contains some of the highest nutrient concentrations in the marine system (Stabeno et al., 1999). This seasonally raises photic zone CO2 and fuels a high productivity region at the Bering shelf break, known as the "green belt" (Springer et al., 1996) (Fig. 1 ). Due to its terminal position in the ocean conveyor and deep bathymetric eddies, any glacial reduction in NPDW upwelling across the subarctic Pacific should also be apparent as reduced nutrient availability in the green belt. 

Modern North Pacific Intermediate Water (NPIW) is sourced from the Okhotsk Sea and is widely distributed across the North Pacific ocean at a water depth between ~300 -800 m characterized by salinity minima with density centred at 26. 8 σ θ (Talley, 1993; Yasuda, 1997). 

During the last glacial maximum, intensified NPIW (termed GNPIW) expanded down to 2000 m (Matsumoto et al., 2002), with previous studies suggesting that the Bering Sea was a key source of GNPIW (Horikawa et al., 2010; Ohkushi et al., 2003) as a result of enhanced brine rejection on the Beringian shelf, following increased sea ice growth since ~900 kyr (Jang et al., 2017; Kender et al., 2018; Knudson and Ravelo, 2015b). Given its density, the formation of GNPIW from the Bering Sea has the potential to suppress upwelling by preventing eddies in the BSC from reaching NPDW. Expansion of this low nutrient GNPIW from the Bering Sea would crucially act to enhance stratification across the subarctic regionnot just the smaller Bering Seacausing region-wide impact on vertical mixing, nutrient supply and the biological pump. 

Here, we present a new high resolution long sedimentary nitrogen isotope record (δ 15 N), in combination with previously published lower resolution results (Kim et al., 2017), from IODP Site U1343 in the Bering Sea green belt (Fig. 1 ) to assess variation in nitrate utilisation between Marine Isotope Stages (MIS) 1 and 21, where we assume that the sedimentary bulk δ 15 N is a representation of total macronutrient utilisation. To account for changes in export production, we normalise the opal mass accumulation rate (MAR) and produce the first high latitude record of nutrient delivery from deep to surface water (termed the "upwelling index"; see Methods). 

By semi-quantitatively assessing the offset between the opal MAR and δ 15 N in this way, we constrain the impact of local growth rate and iron supply, to reveal long-term variability in macronutrient supply from deep water upwelling. We also present a new high-resolution multispecies benthic foraminiferal oxygen isotope record for an improved age model and produce the first assessment of the influence of GNPIW at 2 km water depth in the Bering Sea. From this, we discuss the extent to which oceanographic changes in Bering Sea upwelling and GNPIW production may have contributed towards NPDW upwelling and glacial-interglacial pCO2 changes in the subarctic Pacific over the last 850 kyr.
section_name: Introduction

section_name: Materials and Methods

section_content: IODP Site U1343 (54°33. 4'N, 176°49. 0'E, water depth 1950 m) is situated on a topographic high in the Bering Sea green belt adjacent to the northern continental shelf, proximal to the modern winter sea ice edge and at a depth just below the oxygen minimum zone (Fig. 1 ). 

Marine sediments are composed primarily of fine clays and biogenic material, and are characteristically distinct from shelf-transported materials (Takahashi et al., 2011).
section_name: Core materials

section_content: Foraminiferal preservation is notably good at this site given the shallow North Pacific-wide lysocline which limits carbonate preservation at other Bering Sea sites (Takahashi et al., 2011), providing an excellent opportunity to create a robust high resolution δ 18 OU1343 age model (Asahi et al., 2016; Kender et al., 2018). In the absence of a single species consistently occurring in all samples, 262 samples of ~100 μg of foraminiferal calcite were measured for δ 18 O using ten different species of benthic foraminifera (Cassidulina parkerianus, Cassidulina teretis, Elphidium batialis, Globobulimina auriculata, Globobulimina pacifica, Islandiella norcrossi, Nonionella labradorica, Uvigerina bifurcata, Uvigerina senticosa and Valvulineria spp. ) at an average sampling resolution of 0. 84 m from 2 -222 m CCSF-A. Measurements were made using an IsoPrime© mass spectrometer with a Multicarb device at the British Geological Survey, where the analytical precision of duplicates was <0. 05‰. The stable isotope data were calibrated to the VPDB scale through the international NBS standards. Species-specific offsets from Kender et al. (2018) were applied to fit the data to the most commonly occurring species, E. batialis, and are presented in Figure 2 and Supplementary Table 2. Given that benthic foraminiferal δ 13 C is sensitive to diagenetic alteration, as calcite dissolution, re-precipitation and overgrowth can alter the isotopic composition of their tests (Cook et al., 2016), poor linear regressions between foraminiferal δ 13 C and both raw δ 18 O (R 2 = 0. 09, p < 0. 001) and δ 18 O adjusted to E. batialis (R 2 = 0. 00, p > 0. 05) confirms that diagenetic alteration of foraminiferal shells cannot explain the glacial-interglacial variability in the δ 18 O (Fig. S1 ).
section_name: Foraminiferal isotopes

section_content: Benthic foraminiferal isotope data presented here were combined with previously published data from Site U1343 (Asahi et al., 2016), to create a total of 751 δ 18 OU1343 data points with an average time step of 1. 1 kyr at an average resolution of 0. 30 m from 0 -222 m CCSF-A. The age model, was defined by correlating to the LR04 global composite stack (Lisiecki and Raymo, 2005), and constructed by selecting 17 age-depth tie points during periods of rapid isotopic change primarily at deglacial intervals (Table S2 ). Although there is a possible time lag, where Pacific benthic δ 18 O records may lag the Atlantic Ocean (Asahi et al., 2016; Lisiecki and Raymo, 2005), the high sample resolution and strong correlation between δ 18 OU1343 and LR04 (r = 0. 82, p = <0. 001) suggests that these are within the expected uncertainty of the age model. 

The high sedimentation rate at IODP Site U1343 and a lack of hiatuses in core material provides confidence in the continuity of its record (Asahi et al., 2016; Takahashi et al., 2011). Based on the new composite age model, sedimentation rates were calculated between each age-depth tie point. The average sedimentation rate was 29. 15 cm kyr -1, with higher rates during interglacials and a notable sedimentation rate crash at ~630 kyr (likely due to low biogenic opal; Fig S2 ). These new sedimentation rates were applied to existing opal data (Kim et al., 2014) to calculate an opal MAR record using linear interpolation of the ship dry bulk density data and a ten-point smoothing average (Takahashi et al., 2011).
section_name: Age Model, Sedimentation rates and Opal MAR calculations

section_content: Bulk samples at an average resolution of 0. 81 m CCSF-A, totalling 276 samples between 0 -222 m CCSF-A were analysed for δ 15 N with 500 mg of raw material run on a Carlo Erba 1108 elemental analyzer, interfaced to a Thermo Finnigan Delta Plus XP IRMS at the University of California, Santa Cruz, with a precision of 0. 15‰ based on both sediment standards and replicates. Results, relative to atmospheric nitrogen (AIR), were combined with previously published data from Site U1343 (Kim et al., 2017), to produce a composite record of 373 data points (δ 15 NU1343), with an average resolution of 0. 60 m between 0 and 222 m CCSF-A; average resolution of 2. 3 kyr on the new age model (Fig. S2 ).
section_name: Bulk sedimentary δ 15 N

section_name: Oceanographic indices

section_content: To assess regional variation and discrepancies between Site U1343 and the deep North Pacific Ocean, a Pacific end-member record from ODP Site 849 (Mix et al., 1995)

3851 m water depth) was used to remove high frequency and glacial-interglacial variability. A consistent offset of ~0. 5 ‰ was found between δ 18 OU1343 and δ 18 O849, and so this was added to δ 18 OU1343 first, to ensure offset values reflected only differences in δ 18 O amplitude. After being linearly interpolated at an interval of 2 kyr (Fig. 2 ), a three-point smoothing spline was applied to the δ 18 OU1343, and was subtracted from the similarly processed δ 18 O849 record (Δδ 18 O849-U1343), following a similar methodology to Knudson & Ravelo (2015a). Lower δ 18 OU1343 become more similar to those of shallower water Bering Sea Site U1342 at ~800 m water depth (Fig. 1, 3a), which is suggested to be influenced by GNPIW (Knudson and Ravelo, 2015a). 

Therefore, where Δδ 18 O849-U1343 is approximately 0, Site U1343 is assumed to reflect the same water mass as the deep North Pacific Ocean, which upwells along the Bering Sea shelf edge today (Stabeno et al., 1999). Whereas more positive Δδ 18 O849 -U1343 values occur when δ 18 OU1343 is lower than δ 18 O849, suggesting a colder and/or more saline intermediate water mass is present at Site U1343, similar to the GNPIW signature at Site U1342.
section_name: Δδ 18 O849 -U1343

section_content: To account for the changes in the isotope values of deep North Pacific source water (δ 15 Nsource) and gain insights into regional nitrate utilisation, we follow previous work (Galbraith et al., 2008; Knudson and Ravelo, 2015b) thought to be a site of near complete nitrate utilisation throughout the Pleistocene (Galbraith et al., 2008). δ 15 N1012 was linearly interpolated at a time interval of 2 kyr (the average resolution of δ 15 N1012) and subtracted from the similarly interpolated δ 15 NU1343, to produce a record driven by nutrient utilisation for Site U1343 (Δδ 15 NU1343-1012) (Fig S3 ). 

Given that high Δδ 15 NU1343-1012 values can either reflect higher export production (more rapid removal of nitrate from the euphotic zone) or a slower nutrient resupply (slower replenishment of the surface nitrate pool), and as diatoms are the dominant contributors to export production, the opal MAR can be used to constrain the impact of variable productivity on Δδ 15 NU1343-1012. 

Comparison of the opal MAR with the opal (%) confirms that the MAR record is not an artefact of the age model and calculated sedimentation rates (Fig S2). The opal MAR record was also linearly interpolated at an interval of 2 kyr and then normalised so that all values ranged between 0 and 1. Then the similarly normalised Δδ 15 NU1343-1012 was subtracted to create a semiquantitative proxy termed the 'upwelling index', with low (higher) values indicating decreased (increased) upwelling nutrient supply on the Bering shelf slope (Equation 1 ; Fig. 4b ). 

= - 15 1343-1012 (1)

In doing this we assume that the rate of export production is a reflection on total primary productivity, and that the majority of macronutrient delivery to the photic zone at IODP Site U1343 was from upwelling of NPDW and that rates of nutrient utilisation were controlled by both nutrient upwelling and delivery of iron (Fe) to the photic zone (predominantly derived from sea ice and shelf-derived material (Aguilar-Islas et al., 2008) ). During times of low Fe input and incomplete nutrient utilisation, when increased upwelling will decrease nutrient utilisation (as the nutrient pool is larger) while productivity will show a smaller increase (from the small amount of upwelled Fe which will be scaled proportionally with macronutrient increase), causing the upwelling index to increase at a lower rate than for times of higher Fe delivery. The result of the upwelling index can be classified by three potential states which reflect interglacial, glacial and deglacial changes in biogeochemical cycling (Fig. 5 ).
section_name: Upwelling Index

section_name: Results

section_content: Over the last 850 kyr benthic foraminiferal δ 18 OU1343 varied between +2. 47 and +5. 02 ‰ (mean = +3. 58 ‰), with an amplitude of change similar to LR04 and ODP Site 849 in the eastern equatorial Pacific (range from +3. 08 to +5. 10 ‰; mean = +4. 12 ‰; Fig. 3a, 4a ). The highest Δδ 18 O849-U1343 values occur during glacial maxima of MIS 4, 6, 8, 12, 16 and 20 (Fig. 3b ). As the two δ 18 O records are used to create the age models, with rapid deglacial changes in δ 18 O as the tie points, rapid Δδ 18 O849-U1343 shifts at deglacials may be artificially amplified. However, Δδ 18 O849 -U1343 values are significantly higher in glacial periods (mean = 0. 05) than interglacial periods (mean = 0. 00; p < 0. 05), with values gradually increasing as glacial periods developed and became more intense, suggesting an increasing influence of GNPIW in the Bering Sea slope region. During deglacial periods GNPIW collapses, and the Δδ 18 O849-U1343 becomes more negative during interglacials, as δ 18 OU1343 is more similar to the deep water record from the eastern equatorial Pacific Site 849, suggesting NPDW influence at the 2 km water depth at Bering Sea slope Site U1343.
section_name: Δδ 18 O849 -U1343: GNPIW formation

section_content: Today Bering Sea green belt primary productivity is dominated by diatoms which primarily bloom in spring/summer following seasonal sea ice melt, which provides a source of dissolved iron and promotes water column stability, allowing diatoms to remain in the photic zone and utilise available nutrients (Aguilar-Islas et al., 2008). Estimates of annual primary productivity in the green belt are >170 g C m -2 yr -1 (Sambrotto et al., 2008) with blooms sustained by tidal mixing and lateral supply of iron rich waters from the shelf (Aguilar-Islas et al., 2008). 

Conditions, however, are reversed during winter months (November to April) when sea ice expansion and associated reductions in light availability reduces primary productivity, leading to incomplete nutrient utilisation and the emergence of high-nutrient low chlorophyll (HNLC) conditions similar to the lower subarctic Pacific and central Bering Sea (Galbraith et al., 2008; Knudson and Ravelo, 2015b; Sambrotto et al., 2008). With diatoms as the primary contributors to biogenic material within the IODP Site U1343 sediment record (Takahashi et al., 2011), the opal mass accumulation (MAR) record, recalculated using sedimentation rates from our new age model, is thought to approximate changes in primary productivity (Kim et al., 2014), where productivity was significantly higher in interglacials (mean = 3. 30 g cm -2 kyr -1 ) than in glacials (mean = 2. 66 g cm -2 kyr -1 ; p < 0. 001; Fig. 3e ) (Kim et al., 2014). Although the phytoplankton community structure will likely have varied in response to changing seasonality, diatoms have remained the dominant contributors to export production in the Bering Sea (Takahashi et al., 2011). As changes in the composition of the phytoplankton community will be minor in comparison to glacial-interglacial variation in export production, we assume these changes are unlikely to have impacted our geochemical signals on millennial timescales. Furthermore, silica supply is not considered a limiting nutrient in the Bering Sea (Tsunogai et al., 1979), and hence the opal MAR likely reflects first order changes in productivity (Kim et al., 2014). 

Changes in δ 15 NU1343 are assumed to primarily reflect temporal variations in nitrate utilisation which is influenced by the amount of productivity and nutrient availability, which is in turn dependent on the rate of upwelling of macronutrient rich NPDW, driven by eddies in the BSC along the shelf break. Terrestrial input with low δ 15 N is not considered to have driven the record due to the lack of covariation of δ 15 N and sedimentary C/N ratio (Fig. S2 ), and the broad similarity with ODP Site 882 and MD2416 in the North Pacific which are far from land (Fig. 3c ) (Galbraith et al., 2008; Jaccard et al., 2010). Similarly, variable inorganic contribution is not considered to bias the bulk δ 15 N record, as the intercepts on %TOC vs. %TN do not vary significantly between glacial and interglacials (Fig. S4 ). In addition, diatom-bound δ 15 N records from the Bering Sea and subarctic Pacific show similar trends to their bulk δ 15 N counterparts (Brunelle et al., 2007; Galbraith et al., 2008; Studer et al., 2012), indicating bulk sediment nitrogen is a good proxy for nutrient utilisation in surface water. δ 15 NU1343 values range from +2. 72 to +8. 89 ‰ (mean = +5. 59 ‰) with peak values occurring in glacials, particularly during MIS 2 at ~14 kyr and MIS 8 at ~245 kyr (Fig. 2c ). 

Prior to ~ MIS 13, there appears to be less distinctive glacial-interglacial variability in δ 15 N-U1343, which could reflect sedimentary preservation, however diagenetic alteration of nitrogen is considered less of an influence in this high sedimentation rate setting (Brunelle et al., 2007; Robinson et al., 2012). Alternatively, the greater glacial-interglacial variability in nutrient utilisation over the last ~600 kyr may be expected, given that 100 kyr cycles in the late Quaternary are characterised by colder and more severe glacial periods, and warmer interglacial periods, than the preceding 41 kyr cycles (McClymont et al., 2013). 

Low δ 15 NU1343 excursions (< 4. 5‰) may be a result of inorganic or terrestrial input, as the modern deep ocean average δ 15 N is ~4. 5‰ (Sigman et al., 2009). However, the δ 15 N record is not a direct recorder of deep ocean δ 15 N and more likely reflects surface water δ 15 N, which will be a combination of the deep water signal, as well as surface water processes such as productivity. Confidence in the fidelity of low δ 15 NU1343 values (<5. 0‰) is provided by other subarctic Pacific records from open ocean sites great distances from the continent, such as from Site MD2416 (Galbraith et al., 2008), Site 882 (Jaccard et al., 2010) (Fig. 3c ) and Site 887 (Studer et al., 2012). Furthermore, diatom-bound δ 15 N results also show values less than 5‰ in the subarctic (Brunelle et al., 2007; Studer et al., 2012). Although nitrate uptake can only account for increases in δ 15 NU1343, terrestrial/inorganic nitrogen input is unlikely to impose a first-order control on the glacial-interglacial variability at U1343.
section_name: δ 15 N and Opal MAR

section_content: The upwelling index is highly correlated with the LR04 deep ocean δ 18 O record (r = -0. 72, p < 0. 001)a combination of global ice volume and deep ocean temperaturewith a significantly lower upwelling index in glacials (mean = -0. 40) than in interglacials (mean = -0. 10; p < 0. 001) (Fig. 4 ). Minimal upwelling index values during glacial maxima (due to increased rates of nitrate utilisation (δ 15 NU1343) despite lowered primary productivity (opal MAR)), suggest there was a smaller pool of nutrients in the photic zone, concomitant with expanded GNPIW (Fig. 3b-c, e ). We propose that presence of GNPIW was a pivotal mechanism in reducing nutrient supply from deep water upwelling (see Discussion). A rapid deglaciation phase occurred at all major glacial terminations (dashed lines in Fig. 4 ) when nutrient upwelling peaked, indicating ample nutrient supply from resumed NPDW upwelling and the establishment of a high productivity/low nitrate utilisation state. 

Although upwelling is also correlated with relative sea level (Elderfield et al., 2012 ) (r = -0. 55, p < 0. 001) and Antarctic air temperature (Jouzel et al., 2007 ) (r = 0. 56, p <0. 001) (Fig. 4 ), the high correlation with global surface ocean temperature (Snyder, 2016 ) (r = 0. 60, p < 0. 001) and ice core atmospheric CO2 records over the last 800 kyr are of particular interest (r = 0. 60, p < 0. 001) (Lüthi et al., 2008). These tight correlations are despite the limitations to our productivity and nutrient utilisation proxies and possible inter-site age model discrepancies, therefore suggesting a common link between upwelling, atmospheric temperature and CO2. 

Although correlation does not imply causation, the consistency in the relationship between the upwelling index and the ice-core derived CO2 provides support for a shared underlying mechanism. Exceptions to this correlation between the upwelling index and CO2 occurs during MIS 13, when the upwelling index was high but pCO2 was lower than subsequent interglacials, and MIS 11, when the upwelling index was lower than might be expected given the higher pCO2 for this interglacial (Fig. 4 ). While the former discrepancy may be a result of the MAR calculations, where sedimentation is relatively high during MIS 13, this also correspond with the Mid-Brunhes Event (MBE) where interglacial conditions became warmer (Jouzel et al., 2007). A productivity peak at MIS 13 in response to increased sea ice melt may explain the high upwelling index result at this time. Alternatively, notably warm conditions during "super interglacial" MIS 11 (Melles et al., 2012) likely caused earlier and rapid sea ice melt, which would have reduced iron availability to the green belt, reducing the duration of the spring bloom and limiting summer productivity (Aguilar-Islas et al., 2008). After post-MBE warming, from MIS 11 onward, the relationship between upwelling index and Antarctic CO2 is particularly strong (r = 0. 70, p < 0. 001).
section_name: Glacial-interglacial variability in the upwelling index

section_content: Changes in nutrient upwelling are likely to have been controlled by a combination of sea level change, sea ice presence and GNPIW extent (Kender et al., 2018). Geochemical and biological proxy evidence indicate that an expansion of sea ice occurred during glacials throughout the last 850 kyr (Detlef et al., 2018; Stroynowski et al., 2017; Teraishi et al., 2016), causing lower rates of primary productivity (Kim et al., 2014) and actively restricting air-sea gas exchange. 

Coincident with glacials and maximal sea ice extent, our Δδ 18 O849-U1343 results provide evidence for intermediate water presence to at least 2 km water depth, supporting data from IODP Site U1342 which indicates GNPIW presence to at least 800 m water depth (Knudson and Ravelo, 2015a ) (Fig. 3a-b ). Brine rejection on the Bering Sea shelf during enhanced sea ice formation in glacials is the proposed mechanism for formation of a more saline water body, which would sink to intermediate depths at the shelf slope and propagate southward as GNPIW (Knudson and Ravelo, 2015a), as observed for NPIW in the modern Okhotsk Sea. 

Although direct sea ice reconstructions from the Bering Sea are of lower resolution, they broadly support GNPIW expansion, as glacial periods during this time interval were characterised by extended thick pack ice which peaked during the mid-to-late glacials, when we see the largest offsets in Δδ 18 O849-U1343 (Detlef et al., 2018; Stroynowski et al., 2017). 

Changes in sea-ice would have been in combination with sea level falls of ~50-150 m during glacials (Elderfield et al., 2012) (Fig. 4d ), causing closure of the Bering Strait (modern depth of ~50 m) and retention of both cold/fresh surface waters and denser modified shelf-water (following brine rejection) within the Bering Sea, as the northward flow into the Arctic Ocean was prevented (Kender et al., 2018; Knudson and Ravelo, 2015a). The resultant expanded GNPIW at the shelf edge likely prevented eddies in the BSC from reaching NPDW, isolating nutrients in the photic zone without significant resupply (low upwelling index) and causing the observed high rates of nitrate utilisation but low productivity, whilst also trapping CO2 in deep waters during glacials over the last 850 kyr (Fig. 5b ). 

In addition to sea ice and GNPIW expansion during glacials triggering a reduction in productivity and nutrient upwelling, lowered sea level would have hindered Alaskan Stream water from flowing into the basin by shoaling the passes along the eastern Aleutian arc (Fig. 1b ). Although we cannot quantify this effect, it would have caused some degree of decrease in surface water flow around the Bering Sea basin. This is important because the strength of BSC flow controls eddy field fluctuation, where increased flow maintains higher primary productivity by promoting vertical nutrient supply, and vice versa (Mizobata and Saitoh, 2004). 

There is evidence for glacial reductions in the inflow of the Alaskan Stream from reduced proportions of the diatom species Neodenticula seminae near the inflow sites at the Aleutian arc in the southern Bering Sea, as well as at the slope (Stroynowski et al., 2017; Teraishi et al., 2016). The resultant decrease in BSC flow in glacials would have therefore reduced large eddy formation, further lowering NPDW upwelling of nutrients, productivity and photic zone CO2 (Mizobata and Saitoh, 2004). Simultaneously, sea ice expansion over a shallowed Beringian shelf during glacials may have increased the mobilisation of sediments and transport of shelfderived iron to the slope (Aguilar-Islas et al., 2008). Reduced vertical macronutrient supply but sustained iron supply from sea ice, which is a key source of bioavailable iron in the green belt today (Aguilar-Islas et al., 2008), could help to explain the co-occurrence of low productivity and high nutrient utilisation in glacials (Fig. 2, 5b ). Overall, our results support previous proposals that enhanced sea ice and GNPIW formation during periods of low glacial sea level were a key mechanism for controlling deep water upwelling and green belt nutrient supply over the last 850 kyr (Kender et al., 2018). Reduced upwelling of nutrient-rich NPDW due to GNPIW presence and lowered BSC strength and physical overturning in glacials could then have caused the green belt region to act as a larger net carbon sink. 

During deglaciation, conditions would have reversed with the break up and retreat of sea ice (Kim et al., 2014; Stroynowski et al., 2017) and sea level rise, causing re-strengthening of the BSC and both local and regional collapse of GNPIW (Fig. 5b ). Coupled with high shelf-derived iron from deglacial flooding of the Beringian shelf and potentially increased deglacial Yukon River flow (VanLaningham et al., 2009), this rapid restoration of regional upwelling nutrients explains the observed productivity peaks during deglaciations (Kim et al., 2014) from at least MIS 20-21 onwards. Modelling has shown that late in the last deglacial, wind strength across the Bering Sea and subarctic remained greatly increased due to the remanence of large North America ice sheets (Gray et al., 2018). This would have aided eddy formation in the Bering Sea and upwelling of NPDW across the subarctic region through promoting vertical mixing. 

Once ice sheets had declined and sea ice diminished, seasonally variable upwelling and vertical mixing across the subarctic Pacific would have become re-established in interglacials, similar to modern conditions (Fig. 5c ). 

We have demonstrated for the first time that the glacial reductions in Bering Sea slope primary productivity were driven at least in part by reduced upwelling of nutrients (Fig. 4b ), and propose that this was largely controlled by sea ice and expanded low-nutrient GNPIW (Fig. 3b ). Whilst we do not suggest that this mechanism solely explains the observed proxy evidence from across the North West Pacific Ocean during the Pleistocene, similar productivity and nutrient utilisation δ 15 N records are seen at ODP Site 882 and MD2416 (Fig. 3c, d, f ) (Galbraith et al., 2008; Jaccard et al., 2010), whilst an upwelling index derived for ODP Site 882, using a Δδ 15 N882-1012 record for regional nutrient utilisation and opal (%) as a productivity record (Jaccard et al., 2010) (Fig. 3d, f ), shows remarkably similar trends to the Bering Sea upwelling index. Although the upwelling index for Site 882 only covers a short time period, the high correlation between the upwelling index for Sites U1343 and 882 (r = 0. 77, p < 0. 001) provides support for a common controlling mechanism of GNPIW between the Bering Sea and wider subarctic Pacific region. 

We propose that reduced upwelling of nutrients in the subarctic Pacific Ocean, as recorded at ODP Site 882 and MD2416, was the cause of reduced glacial export production, forced by the expansion of low nutrient GNPIW, similar to the Bering Sea slope (at U1343). The formation of dense GNPIW in the Bering Sea likely propagated south through the deep western Kamchatka Strait (Fig. 1 ), although this does not discount additional contributions of NPIW from the Okhotsk Sea (Max et al., 2014). This mechanism was proposed for the last glacial cycle in the subarctic Pacific Ocean, where model and proxy evidence from the Bering and Okhotsk Seas indicated that enhanced GNPIW during Heinrich Stadial 1 increased intermediate-to-deep ocean stratification, isolating abyssal carbon in NPDW (Gong et al., 2019). Furthermore, the disappearance of GNPIW during the latter part of last deglaciation allowed renewed upwelling of deeper water rich in nutrients and is considered responsible for increases in surface water proxies for CO2 (Gray et al., 2018). Our new data support GNPIW production as the underlying mechanism driving the strength of nutrient upwelling both along the Bering Sea slope and across the open subarctic Pacific Ocean over the past 850 kyr. 

Changes in the Southern Ocean are considered to be a dominant component in driving glacialinterglacial changes in pCO2 (Gottschalk et al., 2016; Jaccard et al., 2013). For example, glacial reduction of atmospheric pCO2 has been linked to an increased rate of biological sequestration in the Southern Ocean and/or reduced supply of nutrients and CO2 upwelled from the deep ocean as a result of intermediate water formation and expanded sea ice cover (Sigman et al., 2010). Although net CO2 flux in the Southern Ocean is regarded as the largest single control on glacial-interglacial variability in pCO2 (Sigman et al., 2010), we demonstrate here that similar processes occurred in the subarctic Pacific region. With modelling studies suggesting that ~30 ppm atmospheric CO2 was released into the atmosphere due to increased overturning in the subarctic region at the last deglaciation (Rae et al., 2014), we surmise that similar magnitude ocean-atmospheric CO2 fluxes may have occurred over the last 850 kyr in response to changes in deep water nutrient upwelling across the Bering Sea and subarctic Pacific Ocean. This aids ongoing efforts to fully explaining observed trends in atmospheric pCO2 over this interval (Gray et al., 2018; Rae et al., 2014).
section_name: Discussion

section_content: Overall, we constructed a new high resolution benthic oxygen isotope stratigraphy for IODP Site U1343, to aid the development of a new nutrient upwelling index which uses δ 15 Nbulk and opal MAR to assess glacial-interglacial variability in the supply of nutrients from NPDW upwelled at the Bering shelf edge. We find a tight correlation between the upwelling index results and global temperature, sea level and atmospheric pCO2, as well as the upwelling index calculated for the short term record at ODP Site 882 in the western subarctic Pacific Ocean. 

We posit that increased sea ice volumes in the Bering Sea and the resultant expansion of GNPIW caused suppression of both nutrient and CO2 upwelling during periods of Bering Strait closure, both in the Bering Sea and across the subarctic Pacific Ocean region. 

Results here suggests that the regional significance of sea ice and deep water ventilation in the subarctic Pacific Ocean has been generally underestimated when modelling Quaternary climate change, where GNPIW growth and propagation from the Bering Sea influences both local and, importantly, regional upwelling of NPDW and export production. The strong correlation between high latitude upwelling and glacial cycles in atmospheric CO2 (Fig. 4b-c ), and the bipolar nature of Quaternary glaciations, also supports suggestions that sea ice was a primary control on the past global climate state (Sigman et al., 2010). Overall, improved understanding of this mechanistic relationship between the Bering Sea and the subarctic Pacific will aid future efforts to constrain the role of glacial-interglacial oceanographic variability in global climate, in complement to Southern Ocean-driven dynamics. 

Geosciences Facilities Steering Committee grant IP-1674 IP- -1116 (to S. K. ) (to S. K. ). SW performed the sample preparation, statistical analyses and led the writing of the manuscript. All authors assisted in writing and contributed to interpretations on the manuscript. SK conceived the overall project and isotope analyses were overseen by MJL and CR. We are also grateful to Dyke Andreasen and Colin Carney from the Ravelo laboratory, as well as Christopher Kendrick and Jack Lacey from the NERC Isotope Geosciences Facility, for support with instrumentation and sample preparation. Further thanks is extended to Sunghan Kim and colleagues, for providing their nitrogen and carbon data from Site U1343. Finally, we would like to thank the reviewers for their helpful suggestions and comments throughout the submission process. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. 882 (50. 35, 167. 58; 3244 m) and MD2416 (51. 27 N, 167. 73 E; 2317m) are also shown, together with surface (red arrows) and deep water (blue arrows) circulation. The high productivity 'green belt' region (Springer et al., 1996), stimulated by shelf-adjacent Bering Slope Current (BSC) is also marked (green patterned band). b) Cross sections of the Aleutian passes and volumes of transport through which Alaskan Stream water exchanges with the Bering basin (Stabeno et al., 1999; Takahashi et al., 2011). (Kender et al., 2018). b) Raw δ 18 O data (red dots) are plotted against δ 18 OU1343 (black line) which has been tuned to E. batialis. c) Comparison of δ 18 OU1343 from this study (black line) against previous data from IODP Site U1343 (Asahi et al., 2016 ) (blue dots), which were combined in d) to create the updated age model for IODP Site U1343 (orange dots), which has been smoothed using linear interpolation to a 2 kyr resolution (black line). from this study and previously published lower resolution data (Kim et al., 2017), smoothed using a linearly interpolation to create a 2 kyr resolution dataset (black line). b) C/N ratios, which are low and show a lack of covariation to δ 15 N, suggesting that terrestrial input of nitrogen was low and did not drive the δ 15 N record. c) New sedimentation rates calculated for IODP Site U1343 using the updated age model presented in this study, which were then applied to d) the biogenic opal record (Kim et al. 2014) to calculate e) the opal mass accumulation rate (MAR) data. (Kender et al., 2018). The total number of samples for each group is also noted.
section_name: Conclusion

section_name: Acknowledgements

section_content: We would like to thank the International Ocean Drilling Program for providing samples as well as the Expedition 323 staff and crew of the JOIDES Resolution, and the curators at the Kochi Core Centre. This research was funded by Natural Environment Research Council (NERC) Envision DTP, CASE funding from the British Geological Survey and NERC Isotope
section_name: Acknowledgements