Silicic Acid Cycling in the Bering Sea During the Mid‐Pleistocene Transition

The rate of deep‐ocean carbon burial is considered important for modulating glacial‐interglacial atmospheric CO2 concentrations and global climate during the Quaternary. It has been suggested that glacial iron fertilization and increased efficiency of the biological pump in the Southern Ocean since the Mid‐Pleistocene Transition (MPT) was key in lowering atmospheric pCO2 and facilitating rapid land ice accumulation. There is growing evidence that a similar mechanism may have existed in the subarctic Pacific Ocean, although this has not yet been assessed. Here, the silicon isotope composition of diatoms (δ30Sidiatom) from the Bering Sea upwelling region is used to assess the role of nutrient cycling on the subarctic Pacific biological pump during the MPT. Results show that during and after the “900 kyr event,” the high productivity green belt zone was characterized by low silicic acid utilization but high supply, coincident with the dominance of diatom resting spores. We posit that as nutrient upwelling was suppressed following pack ice growth and expansion of glacial North Pacific Intermediate Water (GNPIW), primary productivity became nitrate‐limited and enhanced opal remineralization caused a relative increase in silicic acid supply. However, preferential preservation and higher cellular carbon content of diatom resting spores, as well as increased supply of iron from expanded sea ice, likely sustained the net efficiency of the Bering Sea biological pump through the MPT. Remnant iron and silicic acid may also have propagated into the lower subarctic Pacific Ocean through GNPIW, aiding a regionally efficient biological pump at 900 kyr and during post‐MPT glacials.

Although the biological pump in the Southern Ocean, where ∼30%-50% of silicon export production occurs (Gnanadesikan, 1999;Reynolds, 2009), is considered the most significant region for high-latitude ocean-atmosphere CO 2 exchange, there is growing evidence that the subarctic Pacific Ocean could also play an important role in Quaternary climate dynamics (Gray et al., 2018;Haug et al., 1995Haug et al., , 2005Jaccard & Galbraith, 2018;Rae et al., 2014Rae et al., , 2020Sigman et al., 2010). The Bering Sea is considered a key region for understanding how water column dynamics impact the subarctic biological pump in the North Pacific Ocean (Takahashi, 1998). At the Bering shelf slope today, nutrient-rich North Pacific Deep Water (NPDW) upwells due to large eddies in the shelf-adjacent Bering Slope Current. This fuels a region of high productivity known as the "green belt" (annual primary production of 175-275 gC m −2 yr −1 ; Springer et al., 1996), dominated by spring and summer diatom blooms and supported by the release of iron-rich, ice-entrained sediments from the shelf during spring sea ice melt (Figure 1; Aguilar-Islas et al., 2008). Sea ice also modulates the biological pump through its control on light availability and water column stability, as melt waters stratify the water column and allow diatoms to remain in the photic zone for longer (Aguilar-Islas et al., 2008). However, productivity becomes iron-limited in the summer when a high nutrient low chlorophyll (HNLC) environment forms, reducing the efficiency of the biological pump and transforming the region into a seasonal source of CO 2 to the atmosphere (Aguilar-Islas et al., 2007;Takeda, 2011).
Beyond the modern day, the upwelling of nutrient-rich NPDW is also proposed to have been influenced by expanded glacial North Pacific Intermediate Water (GNPIW) event, as a result of brine rejection following sea ice growth (Kender et al., 2018;Knudson & Ravelo, 2015). This GNPIW is proposed to have propagated across the subarctic Pacific region, isolating deep waters, contributing to lower atmospheric pCO 2 and increasing the duration of post-MPT glacials. However, the influence of these oceanographic changes on Bering Sea green belt export production and the regional biological pump have not been constrained. Here, diatom silicon isotope (δ 30 Si diatom ) data from IODP Site U1343 in the green belt region near the Bering slope ( Figure 1) are presented to document changes in rates of silicic acid utilization and the efficiency of the biological pump between 0.6 and 1.2 Ma. In combination with diatom assemblages, diatom productivity (opal mass accumulation rate (MAR)), nitrate utilization (δ 15 N bulk , the nitrogen isotope composition of bulk organic matter) and nutrient upwelling data from the same site are used to investigate the cause and effect of MPT oceanographic changes on the green belt biological pump. Through this and quantification of relative changes in silicic acid utilization and supply through the MPT, we show that sea ice expansion, coupled with the preferential preservation and higher cellular carbon content of diatom resting spores during the 900 kyr event, increased regional carbon export and contributed to atmospheric CO 2 removal.

δ 30 Si diatom
δ 30 Si diatom has been used previously in paleoceanographic studies from the high latitude Southern Ocean (De La Rocha et al., 2011;Egan et al., 2012;Hendry & Brzezinski, 2014) and subarctic Pacific (Bailey et al., 2011;Maier et al., 2015;Swann, 2010) to the eastern equatorial Pacific (Pichevin et al., 2009(Pichevin et al., , 2012 and tropical Atlantic Oceans (Griffiths et al., 2013). As diatoms preferentially utilize 28 Si relative to 30 Si, the ratio of these isotopes can be used to assess variation in the rate of silicic acid utilization and/or supply to the photic zone, with diatoms from periods of higher utilization/lower supply having a higher proportion of 30 Si preserved in their frustules (i.e., higher δ 30 Si diatom ; De La Rocha et al., 1997).
Using the high resolution age model for IODP Site U1343, based on the orbital tuning of benthic foraminiferal δ 18 O (Asahi et al., 2016;Kender et al., 2018;Worne et al., 2020), 30 raw sediment samples were selected between 0.6 and 1.2 Ma, with an average resolution of ∼23 kyr. Diatoms were extracted from these samples with 30% H 2 O 2 and 5% HCl to remove organic material and carbonates, respectively (Swann et al., 2013). Samples were then repeatedly mixed with sodium polytungstate (SPT) at a series of specific gravities from 2.25 to 2.10 g/ml 10.1029/2021PA004284 3 of 16 to isolate a pure diatom sample before being sieved at two size fractions (>38 μm and 20-38 μm), to remove all non-diatom contaminants. Sample purity was confirmed through a combination of X-ray fluorescence (XRF), light microscopy (using a Zeiss Axiovert 40 C inverted microscope) and Scanning Light Microscopy (SEM), with contaminated samples subjected to another round of SPT cleaning or disregarded for isotope analysis. Once cleaned, 28 of the original samples contained sufficient clean material to be analyzed for δ 30 Si diatom from the >38 μm fraction (δ 30 Si large ), and 29 samples from the 20-38 μm fraction (δ 30 Si small ). Samples were prepared for analysis using a step-wise fluorination procedure and δ 30 Si diatom was measured using a Thermo Finnigan MAT 253 at the NERC National Environmental Isotope Facility at the British Geological Survey (Leng & Sloane, 2008). Results are present in δ-notation in per mille (‰) deviations relative to NBS28, using within-run aliquots of NBS28 (1σ = 0.02‰), with replicate analysis of sample material from Site U1343, across the study interval, indicating an analytical reproducibility (1σ) of 0.01‰.

Diatom Assemblage Counts
To ensure that any differences between δ 30 Si large and δ 30 Si small are not an artifact of sieving or species-specific vital effects, diatom assemblage counts were completed on the >38 μm and 20-38 μm samples, with a minimum of 300 valves identified to species level. An exception was made for Chaetoceros RVS, which was grouped at genus level. Using the mean dimensions for each species (Leblanc et al., 2012), relative diatom species biovolumes were calculated for each sample in order to assess the biological origin of individual isotope values. Here, we assume that the biovolume of T. antarctica var borealis resting and vegetative spores (RVS) are the same as those listed for T. antarctica vegetative cells in Leblanc et al. (2012).

Silicic Acid Utilization
Silicon cycling at IODP Site U1343 in the green belt is characterized by a continuous supply of silicic acid (Si(OH 4 )) to the photic zone, from the upwelling contribution of NPDW as well as from the remineralization of Figure 1. Oceanographic map of the Bering Sea. IODP Site U1343 (this study), situated adjacent to the Bering slope is marked by a yellow star, and is situated in the high productivity green belt region (green patterned shape). Productivity is fueled by turbulent eddies in the surface waters (red arrows) of the shelf-adjacent Bering Slope Current, originating from the Alaskan Stream which enters in the southern Bering basin. ODP Site 882, a site of significant previous paleoceanographic study, is marked by a yellow circle. The modern winter sea ice edge is also marked here by a blacked dashed line. diatoms as they sink through the water column (Coachman et al., 1999;Egan et al., 2012). Following this, δ 30 Si diatom is determined by the original δ 30 Si of dissolved silicic acid upwelled to the photic zone (δ 30 Si(OH) 4 ), the fraction of silicic acid which is unutilized and remains in the photic zone (f) and the fractionation factor between diatoms and the dissolved silicic acid (ε).
Previous studies from the subarctic Pacific Ocean and Bering Sea have used a δ 30 Si(OH) 4 value of +1.23‰ Reynolds et al., 2006;Swann et al., 2016), which is derived from the nearest surface value in the subarctic Pacific Ocean (50°N, 167°E; Reynolds et al., 2006). However, a modeling study which considers global opal export fluxes to reconstruct silicon cycling indicates that δ 30 Si(OH) 4 in the subarctic Pacific Ocean can be as high as +1.63‰ (Reynolds, 2009). In this study, both these values are initially used for δ 30 Si(OH) 4 to explore how they impact calculated rates of silicic acid utilization. However, given that deep waters in the Bering Sea are known to be particularly silica rich (Tsunogai et al., 1979), we ultimately assume that the δ 30 Si(OH) 4 value of upwelling deep water at the Bering Slope is more likely represented by the higher modeled value (+1.63‰) than values measured in the lower latitude subarctic Pacific Ocean (+1.23‰; see Section 3.4).
An ε of ∼−1.1 ‰ is often used in paleoceanographic studies Reynolds et al., 2006;Swann et al., 2016), based on experimental work which demonstrates a constant fractionation factor that is independent of species and temperature (De La Rocha et al., 1997). However, culture studies indicate that ε can vary between different polar and sub-polar marine species (Sutton et al., 2013), although there is no supporting evidence from similar studies in natural settings (Cardinal et al., 2007;Fripiat et al., 2012). During the analyzed interval, the diatom assemblages at IODP Site U1343 are dominated by two species: Thalassiosira antarctica var. borealis and Chaetoceros RVS (Teraishi et al., 2016;Worne et al., 2021). Although ε for these species have not been investigated, the Southern Ocean "bipolar twin" of Thalassiosira antarctica var. borealis (known as Thalassiosira antarctica var. antarctica), which thrives in sea ice associated waters and blooms in autumn (Pike et al., 2009), was found in culture studies to have an ε of ∼−0.74 ‰ (Sutton et al., 2013). Some species within the Chaetoceros genus such as Chaetoceros brevis have been shown in culture studies to have a significantly lower ε of −2.09‰ (Sutton et al., 2013), although the impact of this on the IODP Site U1343 δ 30 Si record is likely to be low given that these frustules are generally <20 μm and so mostly removed during the sieving process. To account for potential inter-species variations in ε, a combination of diatom assemblage data, biovolume estimates and mass-balance equations were used to determine the weighted average value of ε used in Equation 2 for each sample and size fraction. Equation 3 demonstrates this, using the aforementioned culture study estimates of ε for T. antarctica var. borealis and Chaetoceros RVS (Sutton et al., 2013), as well as a generic value of −1.1‰ for other taxa (de la Rocha et al., 1997). In doing this, we assume that the ε of −2.09 ‰ for C. brevis (Sutton et al., 2013) is representative of the whole genus.

Silicic Acid Supply
To evaluate the cause of changes in silicic acid utilization through the MPT, records of biogenic silica productivity (opal MAR) at Site U1343 (Kim et al., 2014) can be used to constrain the supply of silicic acid to the photic zone, relative to a value of 0 for the youngest sample at 0.605 Ma (Equation 4). Opal MAR values were linearly interpolated to match the sample ages of the δ 30 Si samples. Through this, we assume that the rate of opal dissolution is constant over time and does not impose a first order control on variability in the opal MAR record.

Silicic Acid Versus Nitrate Utilization
The respective utilization of nitrate in comparison to silicic acid is influenced by iron availability/limitation, as diatoms utilize a larger proportion of silicic acid relative to nitrate when there is an iron deficiency, causing high diatom Si:N ratios (Hutchins & Bruland, 1998). δ 15 N bulk records have been previously used in paleoceanographic studies in the Bering Sea and wider subarctic Pacific Ocean to assess variability in nitrate utilization rates, where phytoplankton will preferentially assimilate the lighter 14 N over 15 N. At a given rate of nitrate supply, the ratio of these isotopes (δ 15 N) is therefore a reflection on the rate and efficiency of nitrate utilization, where higher δ 15 N values reflect periods of enhanced nitrate utilization. Similarly, at a given rate of productivity, δ 15 N will vary depending on the amount of nitrate available, where higher δ 15 N values represent periods of lower nitrate supply.
On the other hand, the δ 15 N bulk record can also be influenced by several other factors, for example low δ 15 N bulk excursions (<+4.5‰) may result from an influx of inorganic or terrestrial material. This is particularly relevant for Site U1343 given its proximity to the paleo-coastline, where terrigenous nitrogen may have been transported from riverine sources and/or coastal shelf erosion as a result of sea-level rise during deglaciation (Khim et al., 2010;Tesdal et al., 2013). However, diatom fossil assemblages do not exhibit deglacial peaks in species associated with shallow coastal waters (benthic), riverine input (freshwater) or those dwelling on the eastern shelf species, hence we expect terrigenous input from these sources was minimal (Teraishi et al., 2016;Worne et al., 2021). Furthermore, previous studies which present C/N and organic carbon isotope (δ 13 C org ) data from IODP Site U1343 confirm that any contribution from inorganic nitrogen does not drive glacial-interglacial variability in δ 15 N bulk (Kender et al., 2018;Worne et al., 2019). δ 15 N bulk has also been proven to show similar trends to diatom-bound δ 15 N (δ 15 N db ), where δ 15 N bulk often exhibits the same temporal trends but with lower amplitude of change (Galbraith et al., 2008;Studer et al., 2012), although this is not always the case, for example offsets between foraminiferal δ 15 N and bulk records were observed in the Southern Ocean (Martínez-García et al., 2014;Robinson et al., 2005).
When using either δ 15 N bulk or compound specific δ 15 N proxies, it is typically assumed that the δ 15 N of the nitrate in the source waters remains constant on geological timescales, however, it is thought that denitrification signals observed in the eastern tropical North Pacific (ETNP) may propagate northward influencing the δ 15 N of nitrate in NPDW over time (Galbraith et al., 2008;Liu et al., 2005). To account for this, the δ 15 N record from ODP Site 1012 (δ 15 N 101 ) in the ETNP is often subtracted from δ 15 N bulk records from the subarctic Pacific (such as δ 15 N U1343 ) to control for changes in the δ 15 N of nitrate in sources waters (Galbraith et al., 2008;Knudson & Ravelo, 2015). Therefore, Figure 3 presents the δ 15 N U1343-1012 record, to facilitate a broad comparison between silicic acid and nitrate utilization. Averages and statistical test results presented below have been calculated using the raw δ 15 N U1343 data ( Figure S1 in Supporting Information S1).

Diatom Assemblages
Forty-three diatom species were identified in both the larger >38 μm fraction (species large ) and smaller 20-38 μm fraction (species small ). The most common species found in both size fractions was Thalassiosira antarctica var. borealis RVS (henceforth T. antarctica), with a mean of 23% relative abundance in species large and 48% relative abundance in species small (Figure 2). This species is known to be dominant during periods of thick pack ice cover (Sancetta, 1982;Stroynowski et al., 2017), most notably during the 900 kyr event (∼0.84-0.91 Ma), when maximal abundances occur in both species large and species small (up to ∼90%; Figure 2), consistent with previous diatom assemblage work at IODP Site U1343 (Teraishi et al., 2016;Worne et al., 2021).
Other prominent species in both size fractions include Neodenticula seminae and Chaetoceros RVS, where the former is used as an indicator of Alaskan Stream water inflow into the Bering Sea gyre, and shows a decreasing trend through the MPT as sea levels decline (Stroynowski et al., 2017;Figure 2), and the latter is most dominant in the spring season, thriving in turbulent waters with high nutrient availability due to upwelling of NPDW (Caissie, 2012;Stroynowski et al., 2017). Key species which show consistently higher abundances in species large than species small , include Coscinodiscus marginatus and Actinocyclus curvatulus ( Figure 2). Both taxa primarily bloom in the summer/autumn season, where the former requires high nutrient availability and warmer waters and reflects high seasonality, while the latter is more commonly associated with low nutrient and stratified waters which follow sea ice melt (Stroynowski et al., 2017;Takahashi, 1986). Increased relative abundance and biovolume of these taxa in species large therefore reflects samples derived from periods of high seasonality with a prominent summer bloom ( Figure S2 in Supporting Information S1). Overall, comparison of δ 30 Si large and δ 30 Si small can aid investigation into the seasonal variability of silicic acid utilization and supply as well as its impact on export production/biological pump through the MPT.
Additional sporadic peaks in species large include Proboscia curvirostris (maximum = 43%), Stellarima microtrias (maximum = 37%), Paralia spp. (maximum = 25%), Actinocyclus spp. (maximum = 17%) and Shionodiscus spp. (maximum = 11%; Figure 2, Figure S3 in Supporting Information S1). Peak abundance of Paralia spp. occur with increased abundance of T. antarctica in species large (e.g. 0.85 and 0.91 Ma; Figure 2); these species are considered littoral and are typically found in higher abundance on the Bering shelf where they thrive in low light coastal environments (McQuoid & Nordberg, 2003;Sancetta, 1982;Stroynowski et al., 2017). Studies from the neighboring Sea of Okhotsk which has a similar oceanographic regime to the Bering Sea, have suggested increased abundance of Paralia spp. is indicative of periods with increased sea ice transport of shelf-transported materials (Wang & Wang, 2008). Therefore, combined increases in T. antarctica and Paralia spp. represent intervals of increase sea ice volume, where IODP Site U1343 becomes more coastal following decreased sea levels and shallowing of Bering shelf waters. In general, species small showed less diversity with T. antarctica RVS consistently dominant over the analyzed interval. Other than isolated peaks of Stephanopyxis turris (maximum = 29%), increased abundances of other taxa in species small were mostly a result of fragmentation, for example at ∼0.94 Ma a high abundance of A. curvatulus was found in both species large and species small as fragments of larger valves (as well as smaller whole valves) had fallen through the sieve.

Silicic Acid Utilization
The potential impact of inter-species variations in ε on rates of silicic acid utilization was checked using mass balance equations and relative diatom biovolumes for two taxa (T. antarctica and Chaetoceros RVS). These species are consistently present and dominate the bulk assemblages during the 900 kyr event (Teraishi et al., 2016;Worne et al., 2021) and have ε values in culture experiments that diverge from the typical value of −1.1 found in marine diatom studies (Sutton et al., 2013;Equation 3). Utilization rates (which were found to be normally distributed) calculated using variable ε values were statistically compared using a Welch's t-test in R (R Core Team, 2017). Results show that varying ε to consider the relative abundance and biovolume of T. antarctica does not significantly alter the rate of utilization large (p > 0.5 and p > 0.5, respectively) or utilization small (p > 0.5 and p > 0.5, respectively), nor does varying the ε value to consider the combined relative abundance and %biovolume of T. antarctica and Chaetoceros RVS for utilization large (p > 0.5 and p > 0.05, respectively) or utilization small (p > 0.5 and, p > 0.5, respectively; Figure 3, Figures S2 and S3 in Supporting Information S1). Therefore, the standard marine ε value of −1.1‰ is a suitable choice for assessing MPT changes in utilization at IODP Site U1343 and is used through the rest of the discussion.
Through selecting a single North Pacific δ 30 Si(OH) 4 value, it is assumed that all nutrients are supplied from NPDW. However, additional inputs may derive from the Alaskan Stream, originating in the north eastern Pacific Ocean, where the nearest modern δ 30 Si(OH) 4 values have been measured in the Gulf of California and peak at 1.7‰ (De La Rocha et al., 2000). The potential contribution of Alaskan Stream surface waters in this study can be inferred from changes in the abundance of N. seminae (Sancetta, 1982;Stroynowski et al., 2017). As N. seminae abundance never exceeds a relative biovolume of ∼2.5% in either size fraction ( Figure S2 in Supporting Information S1), it suggests that any δ 30 Si(OH) 4 signal from the Alaskan Stream would have contributed minimally to the δ 30 Si diatom record at Site U1343. Furthermore, the relative abundance of N. seminae varies similarly in both species large and species small , so this signal will not impact any variation between δ 30 Si large and δ 30 Si small . Therefore, given that using the different North Pacific Ocean δ 30 Si(OH) 4 values of +1.23‰ and +1.63‰ showed no temporal difference in the calculated rates of silicic acid utilization ( Figure S4 in Supporting Information S1), the higher estimate of 1.63‰ is used through the rest of this study as it best agrees with modeled values for Bering Sea deep waters. It also has the advantage of accounting for any potential signal from the influx of higher δ 30 Si(OH) 4 values from the Alaskan Stream (Reynolds, 2009).
Using a δ 30 Si(OH) 4 of 1.63‰, utilization large show low values during early MIS 35 (∼1.2 Ma) in the early Pleistocene, coinciding with low upwelling index values and peak abundance of summer/autumn species Actinocyclus spp. (Figures 2 and 3). Utilization rates remain broadly similar for both fractions between ∼1.1 and 0.95 Ma, followed by a subsequently low utilization large values during MIS 23 (0.9 Ma), consistent with low upwelling index values and peak abundance of T. antarctica in species small (Figures 2 and 3).
Species large also shows increased abundance of T. antarctica, as well as Paralia spp.. Similarly, low utilization large rates occur during early MIS 21 (0.84 Ma), when T. antarctica abundance is maximal in both fractions. After the 900 kyr event, utilization rates calculated from both size fractions are aligned during MIS 19-20 but appear to diverge during MIS 17/18 (∼0.70 Ma), although sampling resolution through the latter interval is low.

Silicic Acid Supply
Applying the utilization rates calculated using δ 30 Si(OH) 4 of 1.63‰ and ε of −1.1, in combination with extrapolated opal MAR data from Site U1343 (Kim et al., 2014; Figure 3), mean rates of Si(OH) 4 supply to the photic zone (henceforth "supply") were calculated. Supply rates are typically 25% higher through the study interval than  (Kim et al., 2014) and ODP Site 882 (Haug et al., 1995), and biogenic silica content (BSi%) from ODP Site 1090 (Diekmann & Kuhn, 2002), reflecting changes in siliceous productivity. (d) δ 15 N U1343-1012 data representing nitrate utilization rates in the Bering Sea green belt Liu et al., 2005;Worne et al., 2020). Blue bars represent glacial periods, as in Figure 2. (e) Upwelling index data set of Worne et al. (2019Worne et al. ( , 2020 (red) which reflects macronutrient supply to the green belt. (f) Atmospheric CO 2 concentrations from ice core records (Lüthi et al., 2008), foraminiferal boron isotopes (Hönisch et al., 2009) and modeling studies (Chalk et al., 2017). those at 0.61 Ma in the >38 μm fraction (supply large ), with maximum supply occurring at 0.904 Ma and increased supply large during the 900 kyr event (∼0.84-0.95 Ma; Figure 3). Alternatively, supply rates calculated from the 20-38 μm fraction (supply small ) are a mean of 12% higher than supply rates at 0.61 Ma, with the maximal supply small occurring at 1.16 Ma. There is a slight increase in supply small around ∼0.85 Ma, which coincides with the end of the 900 kyr event. Generally, the range of supply rate values are much lower in supply small (−60%-120%) than supply large (−56%-415%; Figure 3).

Seasonal Investigation of δ 30 Si diatom Using Assemblage Results
Given that species large generally has a variable but higher abundance of summer/autumn species than species small , the diatom assemblage data can be used to constrain the seasonal changes in silicic acid utilization and supply. When species large is dominated by key summer species, including C. marginatus and A. curvatulus, variability in δ 30 Si large likely reflects silicic acid cycling during the summer/autumn bloom. Alternatively, when species large becomes dominated by T. antarctica and/or Paralia spp., δ 30 Si large will reflect silicic acid cycling during periods of reduced seasonality when the sea ice melt-associated bloom contributes more significantly than the summer bloom to annual production. As the smaller fraction is more constantly dominated by T. antarctica, δ 30 Si small can be considered a better reflection of conditions during the early spring/sea ice season. Variation between the utilization and supply rates calculated from the δ 30 Si diatom of each size fraction can therefore provide insight into the seasonal differences in silicic acid cycling.
There is a larger range of supply large values (compared to supply small ), indicating a greater variability in the seasonality and supply/availability of silicic acid for diatom productivity in the summer months than the winter/sea ice season (Figure 3). This follows modern observations of nutrient availability/biogeochemical cycling in the Bering Sea, where the rate of summer productivity is more variable as it is largely dependent on the length/duration of the spring melt season and the success of the preceding spring bloom (Fukai et al., 2020). For example, a year with a particularly large/long spring bloom will deplete nutrient availability for the summer/autumn period causing lower supply large . Alternatively, a year with a particularly short/early sea ice melt (occurring further north on the shelf) will have a higher supply large , although this could also cause a reduction in iron delivery to the green belt, which may in turn cause enhanced HNLC conditions and limit summer productivity (Aguilar-Islas et al., 2007.

Pre-900 kyr Event
MIS 35 is considered to be a "premature" 100 kyr cycle, where MIS 36 had a prolonged and gradual termination, followed by a long interglacial and gradual return to glacial conditions in MIS 34 (McClymont et al., 2013).
Results here indicate that the MIS 36 termination and early MIS 35 are characterized by lower rates of summer silicic acid utilization (utilization large ) with little change in silicic acid supply at the green belt, although there are a low number of datapoints through this interval. In these samples, there was an increased biovolume of Actinocyclus spp. in species large ( Figure S2 in Supporting Information S1), which thrives in lower nutrient and stratified water columns which are relatively cold and follow prolonged sea ice melting (Stroynowski et al., 2017). Coincident with higher rates of utilization small , this may indicate that the sea ice-associated spring bloom was larger during this interval, causing nutrient limitation for productivity later in the summer months. This is consistent with lipid biomarker sea ice reconstructions from IODP Site U1343 which indicate that there was a pronounced seasonal advance and retreat of the sea-ice margin during this interval (Detlef et al., 2018).
From MIS 34 to 24, low resolution of the data set prevents orbital evaluation of nutrient utilization trends, however silicic acid supply and utilization appear to be broadly consistent through time (Figure 3). Slightly lower δ 30 Si large occurs during MIS 31, which has been globally recognized as a notably warm interglacial (e.g., de Wet et al., 2016;Teitler et al., 2015). Previous records from the Bering Sea indicate that this interglacial is characterized by high productivity fueled by ample nutrient supply (from deep water upwelling) and predominantly open ocean conditions (Detlef et al., 2018;Kim et al., 2014). However, calculations here do not indicate any significant increase in silicic acid supply. Overall, short temporal fluctuations in nutrient utilization through this interval are consistent with suggestions that pre-MPT conditions demonstrated high seasonality and an inconsistent sea ice margin (Detlef et al., 2018).

During and After the 900 kyr Event
The 900 kyr event at IODP Site U1343 exhibits low δ 30 Si diatom and opal MAR, as well as decreased silicic acid utilization (particularly δ 30 Si large and utilization large ) during MIS 23 (Figures 2 and 3). However, nitrate utilization rates mirror the average for preceding interglacials (δ 15 N U1343 = 5.56‰ at 0.90 Ma; Figure 3, Figure S1 in Supporting Information S1). Although high δ 15 N bulk records may reflect diagenetic alteration (Robinson et al., 2005), the high resolution of the δ 15 N U1343 record provides confidence that nitrate utilization rates did not demonstrate a similar rate of decline, as observed with silicic acid (utilization large ; Figure 3). Therefore, if diatoms utilized a larger proportion of the nitrate pool, relative to silicic acid, this may be indicative of abundant iron availability during the 900 kyr event (Hutchins & Bruland, 1998). This is supported by a shift in species large from summer-blooming species C. marginatus to T. antarctica and Paralia spp. (Figure 2). T. antarctica resting spore (RS) are believed to form under periods of nutrient stress and/or low light intensities (Peters & Thomas, 1996), hence they are dominant at the green belt during periods of thick pack ice and suppressed nutrient upwelling, when other species cannot survive (Caissie, 2012). A simultaneous increase in Paralia spp. is consistent with low light availability and increased storm frequency. Overall, this indicates that there was an increased supply of shelf-derived iron from sediments entrained in an expanded sea ice cover, which was blown southward to the green belt (Wang & Wang, 2008). Continued growth of RS could have resulted from continued iron supply from shelf-derived materials, with export production becoming limited by nitrate and low light availability rather than iron.
Although δ 30 Si large (and utilization large ) are low during the 900 kyr event, the supply large rates show maximal values at MIS 23 ( Figure 3). Coincident with low upwelling index values, this suggests that there must have been an increase in supply of silicic acid an alternative source than upwelling of nutrient rich NPDW. In the modern Bering Sea, regeneration of silica is 4-5 times higher than in the adjacent deep waters of the North Pacific Ocean (Coachman et al., 1999;Tsunogai et al., 1979). Therefore, we posit that as silicic acid supply from deep water upwelling was suppressed, and photic zone silicic acid concentrations decreased, there was an increase in remineralization of vegetative diatoms in the water column during the 900 kyr event, supported by low opal MAR ( Figure 3). Although dissolution rates were increased, as resting spores are much more heavily silicified than their vegetative cell counterparts (Doucette & Fryxell, 1983) they sank through the water column faster and were preferentially preserved, leading to maximal abundance of T. antarctica RS in species large ( Figure S5 in Supporting Information S1). This is further supported the dominance of T. antarctica and Chaetoceros resting spores in the bulk sediment during the 900 kyr event (Worne et al., 2021).
Similar trends with notably decreased δ 30 Si and utilization (particularly in the large fraction) occur at 0.85 Ma (MIS 22/21 deglaciation) and ∼0.7 Ma (MIS 17; Figure 2), both also coincident with increases in supply large and also high relative abundance of T. antarctica RS relative abundance (Figure 2), reflecting increased volumes of pack ice in the 100-kyr post-MPT world (Detlef et al., 2018;Stroynowski et al., 2017). Although higher abundance of C. marginatus (and slight increases in A. curvatulus) reoccur in species large during MIS 17, sustained prevalence of T. antarctica in species small and bulk assemblage data (Worne et al., 2021) is consistent with slightly increased supply large , indicating that remineralization of vegetative diatoms likely continued to be an important source of silicic acid after the 900 kyr event. This is in line with previous suggestions that despite partial recovery of summer productivity/seasonality, post-MPT interglacials exhibited a more sea ice-dominant environment (Detlef et al., 2018;Stroynowski et al., 2017;Worne et al., 2021).

Wider Implications -The Biological Pump
The impact of diatom remineralization and sea ice iron fertilization in the Bering Sea green belt during the 900 kyr event has significance for understanding ocean-atmosphere CO 2 exchange and the nature of glacial-interglacial periodicity through the MPT. Iron fertilization is shown to have occurred in the subantarctic zone of the Southern Ocean during the MPT (Figure 3) and is suggested to have been, at least in part, responsible for increasing CO 2 sequestration via the biological pump and leading to low atmospheric CO 2 concentrations (Figure 3; Chalk et al., 2017). Similarly, greater iron availability in the Bering Sea would have increased the efficiency of the biological pump and caused the green belt to sequester a larger proportion of atmospheric CO 2 . However, the net impact of this potential iron fertilization on carbon burial is unclear, given the evidence presented here for simultaneously increased diatom dissolution.
During the 900 kyr event, and subsequent glacial periods, T. antarctica RS dominated the preserved diatom assemblage (Teraishi et al., 2016;Worne et al., 2021). Resting spores are known to have higher ratios of cellular carbon to nitrogen than their vegetative counterparts (McQuoid & Hobson, 1996), supported by evidence from the Southern Ocean which suggests that the selective preservation of T. antarctica and Chaetoceros resting spores (given their heavy silicification) forms the key component of seasonal and annual carbon export (Rembauville et al., 2015(Rembauville et al., , 2016. Although it is beyond the scope of this paper to quantify the impact increased resting spore production and green belt nutrient utilization would have had on carbon cycling, we speculate that dominance of T. antarctica RS and Chaetoceros RVS during the 900 kyr event, alongside sustained iron availability from sea ice-derived sources, may have maintained carbon burial rates, despite increased opal dissolution. At the same time, GNPIW formation and expansion is suggested to have occurred in the Bering Sea during and after the 900 kyr event, which proliferated through the subarctic Pacific Ocean causing regional suppression of NPDW upwelling CO 2 ventilation (Kender et al., 2018;Worne et al., 2020; Figure 3). If the green belt region was iron and silicic-acid replete, it is highly plausible that excess of these nutrients would have propagated in the GNPIW as an important nutrient source to the wider subarctic Pacific (März et al., 2013;Matsumoto & Sarmiento, 2008). The subsequent transfer of this pool of silicic acid and iron to lower latitudes would have promoted the growth of diatoms in those regions, increasing the net global export of carbon into ocean interior and lowering atmospheric pCO 2 . This proposition provides a subarctic analog to the "Silicic Acid Leakage Hypothesis," where glacial increases in iron input to the Southern Ocean lowered diatom Si:N, thereby increasing the pool of unutilized silicic acid in the photic zone (Matsumoto et al., 2002;Matsumoto & Sarmiento, 2008). Furthermore, surface water samples from the western subarctic Pacific Ocean provide evidence for a similar propagation of iron in modern NPIW, which originates from the Sea of Okhotsk (Nishioka et al., 2007).
Increased iron and silicic acid availability across the subarctic would act to increase the efficiency of the biological pump across the subarctic Pacific Ocean during the 900 kyr event, supported by relatively high biogenic opal accumulation at ODP Site 882 during MIS 23 (0.5 gopal cm 2 kyr −1 ) compared with the Pleistocene average (0.42 gopal cm 2 kyr −1 ; 0-1.2 Ma; Figure 3), despite colder sea surface temperatures across the North Pacific (Haug et al., 1995;Martinez-Garcia et al., 2010; Figure 2). This provides further corroboration that GNPIW was a likely linking mechanism between the Bering Sea and wider subarctic Pacific during the mid-late Quaternary, and that the biological pump of the subarctic North Pacific is an important consideration for global ocean-atmosphere CO 2 dynamics (Jaccard & Galbraith, 2018).
Overall, although the Bering Sea green belt alone is unlikely to have placed a principal control on increasing the severity of post-MPT glacials via carbon cycle feedbacks, increased iron and silicic acid delivery to the larger subarctic Pacific region during the MPT and subsequent glacials would have fueled high deglacial productivity once sea levels began to rise and upwelling/vertical mixing of NPDW resumed (Figures 2 and 3). The resultant increase in export production would have driven higher rates of CO 2 sequestration and hence slower net release of CO 2 from the NPDW across the subarctic region, which may have contributed to increased glacial duration globally.

Conclusion
The 900 kyr event is characterized by reduced silicic acid utilization in the Bering Sea green belt, despite an increase in silicic acid supply to the photic zone during this interval. Species assemblage results also indicate that there was low seasonality and a very short/no summer bloom period, likely due to reduced light and nutrient availability following the development of thick pack ice and GNPIW expansion. We propose that expansion of GNPIW and suppression of macronutrient upwelling after the closure of the Bering Strait caused the region to become nitrate-rather than iron-limited, while photic zone silicic acid supply was maintained by diatom remineralization. Low nitrate availability also triggered an increase in the production of diatom resting spores, which are dominant in the MPT sediments at IODP Site U1343. Increased resting spore production was likely fueled by sea ice-derived iron and may have caused enhanced CO 2 export during the 900 kyr event, given their preferential preservation and higher cellular carbon content, although this may be offset by simultaneously increased opal dissolution. Further quantitative investigation into utilization rates and carbon cycling is required to test the net effect of these opposing mechanisms following suppressed NPDW upwelling.
As a result of nitrate limitation, iron and silicic acid that were left unutilized at the Bering slope could have propagated through the wider subarctic Pacific Ocean in GNPIW providing an important regional glacial and deglacial nutrient source. Low resolution evidence after the 900 kyr event indicates that more severe sea ice seasons persisted in both glacials and interglacials, with sea ice-derived iron and/or remineralized silicic acid continuing to provide an important supply for both green belt and subarctic Pacific Ocean productivity. However, future high resolution investigations are required to evaluate the relative variability in different nutrient cycles on orbital and suborbital timescales and their impact on the regional biological pump. Overall, this study highlights that understanding high latitude upwelling systems in relation to sea ice dynamics and the biological pump is of critical importance for assessing controlling mechanisms for ocean-atmosphere CO 2 exchange and global climate.

Data Availability Statement
The data presented in this paper are stored in the Pangaea data repository (https://doi.pangaea.de/10.1594/ PANGAEA.933139). BGS (GA/15S/003). δ 30 Si analysis was completed through funding provided by the NERC Isotope Geosciences Facilities Steering Committee grants IP-1740-0517 (to GEAS) and IP-1413-1113 (to SK). S. Worne performed the sample preparation, statistical analyses and led the writing of the manuscript. J. H. Lacey performed laboratory analyses of prepared sample material. All authors assisted in writing and contributed to interpretations on the manuscript. The authors would also like to thank the reviewers for their constructive comments which greatly advanced the interpretations of this manuscript.