Benthic foraminifera indicate Glacial North Pacific Intermediate Water and reduced primary productivity over Bowers Ridge, Bering Sea, since the Mid-Brunhes Transition

The Mid-Brunhes Transition (MBT) saw an increase in the amplitude of glacial cycles expressed in ice core and deep ocean records from about 400 ka, but its influence on high-latitude climates is not fully understood. The Arctic Ocean is thought to have warmed and exhibited reduced sea ice, but little is known of sea ice marginal locations such as the Bering Sea. The Bering Sea is the link between the Arctic and Pacific Ocean and is an area of high productivity and CO₂ ventilation; it hosts a pronounced oxygen minimum zone (OMZ) and is thought to be the location of Glacial North Pacific Intermediate Water (GNPIW) formation in the Pleistocene. To understand palaeoceanographic change in the region, we analysed benthic foraminiferal faunas from Bowers Ridge (Site U1342, 800 m of water depth) over the past 600 kyr, as they are uniquely well preserved and sensitive to changes in deep and surface ocean conditions. We identified and imaged 71 taxa and provide a full taxonomy. Foraminiferal preservation is markedly higher during glacials, indicating the presence of less corrosive GNPIW. The most abundant species are Bulimina exilis, Takayanagia delicata, Alabaminella weddellensis, Gyroidina sp. 2, Cassidulina laevigata, Islandiella norcrossi, and Uvigerina bifurcata, consistent with broadly high net primary production throughout the last 600 kyr. Correspondence analysis shows that the most significant Assemblage 1 comprises B. exilis, T. delicata, Bolivina spissa, and Brizalina, which occur sporadically within intervals of laminated, biogenic-rich sediment, mostly during glacials and also some deglacials, and are interpreted as indicating very high productivity. Other assemblages contain the phytodetritivore species A. weddellensis, I. norcrossi, and C. laevigata, indicative of seasonal phytoplankton blooms. Before the MBT, more numerous intervals of the very high-productivity Assemblage 1 and A. weddellensis occur, which we suggest reflect a time of more sea-ice-related seasonal stratification and ice edge blooms. Our inference of a decrease in sea ice meltwater stratification influence in the central Bering Sea after the MBT is consistent with records showing that the Arctic and Pacific Ocean warmed during glacials and suggests that high-latitude productivity and sea ice changes were an important feature of this climate event.

1 Introduction

The Mid-Brunhes Transition (MBT) was a period of climate change that saw the emergence of greater glacial–interglacial variability in temperature and CO₂ from ∼400 ka (Jouzel et al., 2007; Lüthi et al., 2008). Glacials became slightly colder over Antarctica, interglacials warmer, and global deep ocean δ¹⁸O records (Lisiecki and Raymo, 2005) show greater variability after the MBT, indicating changes to deep ocean temperature and global ice volume. The cause of the MBT is under debate but may have been related to changes in high-latitude Southern Ocean sea ice, windiness, upwelling, and deep ocean ventilation (Kemp et al., 2010; Yin, 2013). However, there is a need for high-latitude palaeoceanographic records over the MBT to test the global expression of this climate event (Barth et al., 2018) and determine if the amplification of climate change occurred in the higher latitudes (Cronin et al., 2017). Intermediate water temperature and sea ice records from the Arctic Ocean have shown that the MBT was associated with a pronounced warming and reduction in sea ice influence (Cronin et al., 2017). The Bering Sea is a critical region, as there is a pronounced OMZ, deep ocean upwelling, CO₂ ventilation to the atmosphere (Stabeno et al., 1999), sea ice, and an oceanographic link between the Arctic and Pacific Ocean.

Study Site U1342 (Fig. 1) is located on Bowers Ridge, an extinct arc system extending 300 km north from the Aleutian Island arc, and was the shallowest site cored during the Integrated Ocean Drilling Program (IODP) Expedition 323 (818 m of water depth; Expedition 323 Scientists, 2011). Deep water connections are restricted to the Aleutian Island passes (<4000 m of water depth; Hood, 1983; Stabeno et al., 1999) such that the majority of Bering Sea deep water is sourced from aged low-oxygen North Pacific Deep Water (NPDW). Nutrient-rich intermediate nutricline water upwells over Bowers Ridge, particularly along the northern slope, inducing high primary productivity and associated organic carbon flux to the sea floor (Stabeno et al., 1999; Takahashi, 2005). There is also significant slope-shelf water exchange that brings nutrients up from the slope (Stabeno et al., 1999) and a very minor component of deep water (∼4000 m) formed within the Bering Sea, probably from sea ice brine rejection (Warner and Roden, 1995). Today, seasonally sea-ice-covered areas of the Bering Sea are found mainly over the northern continental shelf (Fig. 1). Bowers Ridge and the southern Bering Sea are affected by the relatively warm and high-salinity Alaskan Stream surface water current (Fig. 1). The onset and termination of summer stratification, from the solar warming of surface water and also sea ice melt, causes spring and autumn phytoplankton blooms (Niebauer et al., 1995; Eslinger and Iverson, 2001; Kuroyanagi et al., 2002) and associated high seasonal phytodetrital flux to the sea floor. High productivity, combined with the low-oxygen NPDW-sourced water bathing Bowers Ridge (Stabeno et al., 1999), causes the OMZ to be pronounced (∼0.6 mL L⁻¹ at the depth of Site U1342; Conkright et al., 2002).

Figure 1 (a) Map of the Bering Sea with modern surface ocean flow (white arrows), showing the position of Site U1342 (this study) and other sites referred to in the text. Bathymetry above 100 m of depth is marked as grey, approximating glacial coastlines. Last Glacial Maximum (∼30 kyr ago) sea ice conditions (blue lines and text) are from Katsuki and Takahashi (2005), interpreted from diatom assemblages. Map drawn in Ocean Data View (Schlitzer, 2019). (b) Bering Sea north–south section showing dissolved oxygen concentrations (Schlitzer, 2019) and the position of core sites discussed in this study.

Studies have shown that sea ice was present in the Bering Sea through at least the last 2.5 Myr, became more prevalent during the mid-Pleistocene (Stroynowski et al., 2015, 2017; Detlef et al., 2018), and likely resulted in GNPIW expansion (Horikawa et al., 2010; Knudson and Ravelo, 2015; Kender et al., 2018; Worne et al., 2019) and a reduced OMZ. However, high-resolution Bering Sea records of organic carbon flux and OMZ presence have not yet been produced over the MBT. Existing low-resolution records at the Bering Sea slope do suggest that sea ice was reduced at ∼500 ka (Stroynowski et al., 2017), but there are no existing records with the resolution required to test the sea ice influence on productivity at a marginal location in the central Bering Sea. Diatom assemblage records from several piston cores in the Bering Sea and North Pacific (Katsuki and Takahashi, 2005) have been used to reconstruct Quaternary sea ice variation (see Leventer et al., 2007). Katsuki and Takahashi (2005) showed that during the Last Glacial Maximum (LGM), when sea level was over 100 m lower than today and much of the continental shelf was exposed (Fig. 1), seasonal sea ice covered the northwestern slope, drift ice covered the majority of the central Bering Sea and much of Bowers Ridge, and open water covered the south–central Bering Sea and parts of northern Bowers Ridge, including the location of Site U1342 (Fig. 1). However, there is currently little information on the evolution of intermediate water properties over Bowers Ridge during this time period or of sea ice, organic carbon flux, and OMZ and GNPIW presence over the MBT.

Deep sea benthic foraminifera respond to changes in organic carbon flux and oxygen availability in modern settings (e.g. Kaiho, 1994; Thomas and Gooday, 1996; Jorissen et al., 2007; Gooday and Jorissen, 2012; Kaminski, 2012) and have been shown to respond to environmental variability in the Bering Sea during the last glacial cycle (Gorbarenko et al., 2005; Okazaki et al., 2005). For instance, the phytodetritivore species Alabaminella weddellensis dominated faunas during the LGM over the Umnak Plateau (Site UMK-3A; Fig. 1), and the high-productivity–low-oxygen genera Rutherfordoides and Bulimina dominated faunas during the deglacial when surface water productivity peaked (Okazaki et al., 2005). Coretop “mudline” samples from U1342 of sub-recent age reveal an assemblage dominated by the high-productivity–low-oxygen taxa Bulimina, Brizalina, Globobulimina, Stainforthia, Rhumblerella, and Martinottiella (Expedition 323 Scientists, 2011; Kender and Kaminski, 2017), which is consistent with the position of U1342 within the core of the modern OMZ. To better understand surface and deep water evolution over glacial–interglacial cycles of the last 600 kyr at ∼800 m of water depth (within the modern OMZ), we analysed benthic foraminiferal assemblages in high resolution at Site U1342 to constrain palaeoceanographic changes.

2 Methods and materials

The sediments at Site U1342 are a mixture of biogenic material (largely diatom frustules and foraminifera, with minor nannofossils, silicoflagellates, sponge spicules, and radiolarians), intermittent laminated intervals (predominantly parallel but occasionally cross-bedded), volcaniclastic material (fine to coarse ash) occurring in both discrete layers and scattered in the sediment, and rare clay to pebble-sized siliciclastics (Expedition 323 Scientists, 2011). This site is chosen as it is on a ridge at ∼800 m of water depth, and the lysocline is shallow in the North Pacific such that calcite dissolution is strong below ∼1000 m. As the Bowers Ridge region is far from continental areas and was submerged during at least the Pleistocene, the coarser siliciclastic component was likely derived from sea ice and/or icebergs (Aiello and Ravelo, 2012). The laminated sediments have a distinct biogenic composition dominated by diatom and foraminiferal tests. Conversely, the homogenous intervals contain a variety of sedimentary components that include both biogenic and terrigenous particles (Aiello and Ravelo, 2012). Expedition 323 Scientists (2011) report that the majority of laminated intervals are bedding-parallel, while cross-bedding and scouring indicative of bottom water currents are less common.

The age–depth tie points used to plot our data against age are derived from detailed benthic foraminiferal oxygen isotope (δ¹⁸O) analyses carried out on Uvigerina peregrina in the same samples by Knudson and Ravelo (2015). The isotope stratigraphy of U1342 has a correlation coefficient of 0.64 with the global reference stack LR04 (Lisiecki and Raymo, 2005) after a three-point running mean was applied. The accuracy of the age model is probably better than ±10 ka, which is taken into account when evaluating the data with respect to global glacial–interglacial cycles.

A total of 160 samples were collected from the top ∼20 m CCSF (composite depth scale; Expedition 323 Scientists, 2011) of Site U1342 at ∼14 cm spacing, corresponding to the last ∼600 kyr, and an average time resolution of ∼3.7 kyr between each sample. Sample processing was carried out at the University of California, Santa Cruz. Samples were freeze-dried, gently washed over a 63 µm mesh screen, and stored in glass vials before picking. Slides are housed in collections at the University of Exeter, UK. Due to the small size of the samples (ranging between 10 and 20 cc), we aimed to pick a minimum of 100 specimens of benthic foraminifera per sample into cardboard reference slides. Samples that contained many more than 300 specimens were subdivided using a micro-splitter. Benthic foraminiferal fragmentation percentage was determined by calculating the percentage of broken specimens (fragments) to the sum of whole (unbroken) and fragmented specimens in each sample.

Correspondence analysis (CA) (using the software PAST; Hammer et al., 2005) was carried out on samples containing >50 specimens and with species with fewer than 50 specimens in the whole dataset removed. Using a reciprocal averaging algorithm (Greenacre, 1984; Hammer and Harper, 2006) to compare species occurrences across samples, similar CA scores for samples indicate comparable faunal characteristics (e.g. Kuhnt et al., 2002; Hammer and Harper, 2006; Kender et al., 2008). CA attempts to show multidimensional data in fewer dimensions, with axis 1 (describing the largest variance in the dataset) being shown to represent the strongest ecological parameter in many studies (Hammer and Harper, 2006).

3 Results

We identified and counted a total of 71 taxa (Table S1 in the Supplement) and completed detailed taxonomic analysis and imaging of each one (see Taxonomy section in the Supplement). In total, 106 of the 160 samples contained >100 specimens (Fig. 2c). Diversity averaged 12 taxa per sample but ranged up to over 20 in several samples (Fig. 2c). The number of taxa typically increased in samples with greater abundance (Fig. S1 in the Supplement), with samples having fewer than 50 specimens (n=43) showing the most reduced diversity (Fig. 2c). In these samples, the foraminifera showed a higher degree of fragmentation (Figs. 2b, S1 in the Supplement), which has been used as a qualitative proxy for dissolution (Metzler et al., 1982; Thunell, 1976; Berger, 1970). We therefore removed these 43 samples from assemblage analysis. Intervals of less dissolution and good preservation, as defined by abundance, diversity, and fragmentation, appear to be restricted to glacial maxima as defined by benthic foraminiferal δ¹⁸O (Fig. 2a, grey bars). A tight correlation with abundance, diversity, and glacial cycles can be identified, with the highest abundances in the glacial periods before 400 ka (Fig. 2). The species that make up the majority of the assemblages are Bulimina exilis, Takayanagia delicata, Alabaminella weddellensis, Gyroidina sp. 2, Cassidulina laevigata, Islandiella norcrossi, and Uvigerina bifurcata, which change in their dominance multiple times throughout the core (Fig. S2 in the Supplement).

Figure 2 (a) U1342 benthic foraminiferal δ¹⁸O record used to construct the age model (Knudson and Ravelo, 2015). (b) The degree of foraminiferal fragmentation at Site U1342 as a proxy for bottom water dissolution. (c) Abundance (red) and diversity (blue) of benthic foraminifera at Site U1342 as a proxy for preservation. Grey bars indicate glacial maxima as defined by the positive δ¹⁸O.

CA axis 1 describes 19% of the variance within the dataset, with positive values signifying times of abundant B. exilis, T. delicata, Epistominella exigua, Bolivina spissa, Brizalina alata, and Brizalina earlandi. These species make up Assemblage 1, which correlates very well with CA axis 1 (Fig. 3b). CA axis 2 (14 % of the variance) shows a good correlation with A. weddellensis, Cassidulinoidies parkerianus, Globobulimina auriculata and Stainforthia fusiformis (Assemblage 2) for positive values, and C. laevigata (Assemblage 3) for negative values (Fig. 3c). Other axes do not define further meaningful assemblages, although I. norcrossi and U. bifurcata are important constituents of the fauna outside these three assemblages (Fig. S2 in the Supplement) and plot near the CA axes 1 and 2 origin (Fig. 3a), signifying no particular abundance preference for the three defined assemblages. Both CA axes 1 and 2 show considerable variability with time when plotted against age (Fig. 4), with numerous intervals of dominant Assemblage 1 co-occurring with laminated sediment but no clear glacial–interglacial variability. Of the 15 samples analysed from laminated sediments, 13 (87 %) have CA axis 1 values above 0. Of the 98 samples analysed from non-laminated sediments, 78 (80 %) have CA axis 1 values below 0. In contrast, CA axis 2 shows a weak preference for unlaminated, homogenous samples.

Figure 3 (a) Correspondence analysis (CA) of the dataset for samples with >50 specimens and for species with >50 individuals. Samples are marked with grey dots, and grey circles are the approximate position of the interpreted assemblages. (b) Percentage counts for Assemblage 1 (B. exilis, T. delicata, E. exigua, B. spissa, B. alata, and B. earlandi) against CA axis 1 scores for each sample (with >50 specimens). (c) Percentage counts for Assemblage 2 (A. weddellensis, C. parkerianus, G. auriculata, and S. fusiformis) and Assemblage 3 (C. laevigata) against CA axis 2 scores for each sample (with >50 specimens). Linear regression for Assemblage 2 is only for samples with positive CA scores, and Assemblage 3 is only for negative CA scores.

Figure 4 Benthic foraminiferal dataset CA axis 1 and 2 scores compared with various other records plotted on age. Marine Isotope Stages (MISs) are shown at the top. (a) The Site U1342 benthic foraminiferal δ¹⁸O record (blue) used to construct the age model (Knudson and Ravelo, 2015) against the global benthic δ¹⁸O composite record (Lisiecki and Raymo, 2005). (b) Site U1342 magnetic susceptibility record (Expedition 323 Scientists, 2011). (c) Benthic foraminifera CA axis 1 scores, with high values indicating Assemblage 1 (Fig. 3). (d) Benthic foraminifera CA axis 2 scores, with high values indicating Assemblage 2 and low values Assemblage 3 (Fig. 3). (e) Bottom water temperature (BWT) record from Mg∕Ca of ostracods in various Arctic Ocean intermediate water depth sites (Cronin et al., 2017). (f) Sea surface temperature (SST) record from alkenones at North Pacific Site 882 (Martínez-Garcia et al., 2010). Green bars indicate laminated sediment at Site U1342 (Knudson and Ravelo, 2015).

4 Discussion and conclusions

Each glacial period of the past 600 kyr experienced a remarkable increase in the preservation of calcareous benthic foraminifera at Site U1342 (fragmentation, abundance, and diversity; Fig. 2). The most likely reason is the presence of a different, less corrosive glacial water mass over Bowers Ridge, which has been identified as GNPIW from offsets in oxygen and carbon isotopes at U1342 (Knudson and Ravelo, 2015) and other sites (Cook et al., 2016; Kender et al., 2018; Worne et al., 2019), and neodymium isotope changes within the Bering Sea (Horikawa et al., 2010). Positive benthic δ¹³C and depleted benthic δ¹⁸O during glacials at U1342 were interpreted as the export of surface waters to depth via sea ice brine rejection and intermediate water formation (Knudson and Ravelo, 2015). Surface waters were likely depleted in ¹²C, CO₂, and DIC compared with aged deeper NPDW due to primary productivity, which could have led to an increased carbonate ion concentration and better preservation of CaCO₃. Modern NPIW is formed partially in the Sea of Okhotsk, which has typical δ¹³C values of −0.2 ‰ to −0.4 ‰ at 800 m of water depth compared with −0.6 ‰ to −0.7 ‰ in the Bering Sea (Keigwin, 1998; Cook et al., 2016) where intermediate water does not currently reach. Alternative possible causes of increased glacial preservation are decreased primary productivity, which is not supported by the distribution of laminated sediment and benthic foraminiferal assemblages (see below), and lower sea level, which is unlikely as the modern δ¹³C chemocline is at ∼250 m of water depth in the Bering Sea (Keigwin, 1998).

To assess the palaeoecological changes at Site U1342, we removed samples heavily affected by poor preservation and examined the underlying benthic foraminiferal assemblage changes. Surprisingly, although preservation shows strong glacial cyclicity, foraminiferal assemblages do not show a glacial trend but a preference for laminated sediment that also does not follow glacial cyclicity (Fig. 4). Benthic foraminiferal faunas are controlled by a range of factors that include primary productivity and organic carbon supply, sediment heterogeneity, bottom water currents, organic carbon type, quantity and seasonality, bottom water and sediment pore water oxygen levels (Sun et al., 2006; Jorissen et al., 2007) and the ability of organic carbon to reach the sea floor after passing through different water masses with varying degradative constituents (Arndt et al., 2013). In the deeper Bering Sea, organic carbon flux and oxygen levels are the major two variables, controlled by a combination of export production and intermediate water ventilation (Kender and Kaminski, 2017). Of our most dominant species, C. laevigata, Uvigerina peregrina, and species from the genera Bulimina, Bolivina, Stainforthia and Globobulimina are known denitrifiers able to survive in anoxic water in the modern ocean by respiring nitrate (Piña-Ochoa et al., 2010). We interpret the assemblages described here as generally tolerant of low oxygen, even though the majority of our species have not yet been analysed for nitrate reduction. Whilst E. exigua, Cibicidoides, and other minor constituents found here are unable to respire nitrate (Piña-Ochoa et al., 2010), E. exigua is abundant in the OMZ of the Peru–Chile Trench (Sen Gupta and Machain-Castillo, 1993; Erdem and Schönfeld, 2017), and if the majority of species at U1342 are able to respire nitrate it appears as though oxygen limitation alone is unlikely the predominant control on benthic faunas. Some species may still have a competitive advantage over others in low-oxygen settings, and indeed various OMZs today contain faunas that have characteristic assemblages (Hermelin and Shimmield, 1990; Sen Gupta and Machain-Castillo, 1993; Kaminski et al., 1995; Kaiho, 1999; Gooday et al., 2000; Schumacher et al., 2007). In the core of the OMZ of the Arabian Sea (∼0.1 mL L⁻¹), C. laevigata, B. exilis, Uvigerina semiornata (a morphologically similar species to U. bifurcata), U. peregrina, Globobulimina, and Bolivina species all show high abundance at different depths (Schumacher et al., 2007). Bulimina spp., T. delicata, Bolivina spissa, and U. peregrina were abundant in OMZs of offshore California (Sen Gupta and Machain-Castillo, 1993), and I. norcrossi was abundant in the Okhotsk Sea where oxygen levels were 0.3 mL L⁻¹ (Bubenshchikova et al., 2008). The only abundant species in our dataset that does not typically occur in high abundance in OMZs at lower latitudes is A. weddellensis, although it is associated with low-oxygen-tolerant species in this study, and it is present in six of our laminated samples with an abundance greater than 50. Thus, relatively high organic carbon flux conditions were likely present throughout much of the past 600 kyr, as the defined assemblages and the residual fauna, consisting largely of U. bifurcata, are all indicative of modern high-productivity OMZs, and typical abundant deep water benthic foraminifera from oligotrophic well-oxygenated settings, such as Planulina wuellerstorfi and Cibicidoides mundulus, are all but absent.

It seems most likely that the assemblages in our study are controlled predominantly by changes in the supply of organic carbon to the sea floor, as most of the species are tolerant of low oxygen. The high correlation of Assemblage 1 with laminated sediments (in >80 % of the samples – green vertical bars in Fig. 4) is evidence that those species are indicative of the highest organic carbon flux to low-oxygen regime to the sea floor at this site. That is because (1) the associated low oxygen (the absence of deep bioturbation, allowing laminations to be preserved) was likely caused by high benthic respiration rates from high organic carbon flux, and (2) the laminations are enriched in biogenic material as indicated by diatom-rich layers (Expedition 323 Scientists, 2011) and consistent low magnetic susceptibility (MS) (Fig. 4) signifying non-magnetic biogenic material (Expedition 323 Scientists, 2011). Alabaminella weddellensis (from Assemblage 2) is regarded as an opportunistic species, blooming at times of high seasonal phytodetrital flux to the sea floor (Smart et al., 1994; Thomas and Gooday, 1996). On the Arctic shelf, high proportions of I. norcrossi coincide with seasonal sea ice and probable summer ice edge productivity (Steinsund, 1994; Polyak et al., 2002; Ivanova et al., 2008), and the species is considered to respond and reproduce with highly seasonal phytodetritus deposition (Wollenburg et al., 2004). A laboratory experiment by Alve (2010) showed that living C. laevigata (Assemblage 3), in contrast to many other species, disappeared from the assemblage after a halt in the supply of fresh phytodetritus, indicating that this species requires an input of fresh phytodetritus to survive. Outside Assemblage 1 (highest organic carbon flux), it appears as though I. norcrossi, A. weddellensis, and possibly C. laevigata were competing as seasonal phytodetritus species in a relatively (compared to Assemblage 1) lower organic carbon flux environment, responding to changes in the type, quantity, and quality of organic carbon supplied to the sea floor.

Diatom assemblages from nearby Site BOW-8A (Fig. 1) can be used to indicate surface ocean conditions proximal to U1342 over the past 270 kyr (Katsuki and Takahashi, 2005). BOW-8A is dominated by Neodenticula seminae, indicating that the Alaskan Stream (Sancetta, 1982) influenced the palaeoceanography at the western Bowers Ridge over this time, and low proportions of the sea-ice-living diatoms Thalassiosira gravida and Fragilariopsis cylindricus indicate that this location was similarly unaffected by persistent sea ice over much of this time period (Katsuki and Takahashi, 2005). Jin et al. (2006) modelled Bering Sea reductions in spring thermal stratification, wind mixing and tidal mixing, and found that – even though they had different effects on the timing of blooms and zooplankton–phytoplankton concentrations – all caused reduced net primary production, thus providing possible mechanisms for past changes in organic carbon flux to the sea floor. We propose three broad scenarios for the annual cycling of surface water masses at Bowers Ridge over the Quaternary in the absence of significant sea ice (as indicated by diatom records; Katsuki and Takahashi, 2005). The first scenario is the modern setting (interglacials) in which water mixing brings intermediate nutricline water, influenced by nutrient-rich NPDW, into the photic zone, and summer warming causes spring–summer blooms and high net primary production (Eslinger and Iverson, 2001; Katsuki and Takahashi, 2005). The second scenario is weaker glacial spring thermal stratification due to colder glacial conditions, coupled with enhanced glacial wind mixing as modelled due to the presence of larger North American ice sheets (Gray et al., 2018). Although mixing would have brought deeper water into the photic zone, reduced thermal stratification would have acted to reduce net spring–summer bloom productivity (Katsuki and Takahashi, 2005; Jin et al., 2006). In addition, during glacial times, wind mixing would have brought a smaller proportion of nutrients into the photic zone due to enhanced nutrient-poor GNPIW presence (Kender et al., 2018; Knudson and Ravelo, 2015). The third scenario is glacial spring–summer stratification from southward-flowing cold low-salinity surface water sourced from melting sea ice, causing a spring–summer bloom and elevated seasonal phytodetrital flux (Eslinger and Iverson, 2001; Katsuki and Takahashi, 2005; Aguilar-Islas et al., 2007). Site U1342 was probably affected by seasonal stratification over the majority of glacials (Fig. 4), as suggested by the high proportions of phytodetritivores I. norcrossi, A. weddellensis, and possibly C. laevigata, perhaps similar to scenario three. Although sea ice melt stratification would likely have been more extensive during strong glacials, we may not expect to see an increase in net primary production as the greater presence of GNPIW would have separated high-nutrient deeper water from the photic zone. These competing mechanisms at Bowers Ridge may explain why there is no clear glacial–interglacial cyclicity to the foraminiferal assemblages and laminated sediments. The lack of abundant sea ice diatoms over the last ∼270 kyr indicates that the bloom production, suggested by abundant phytodetritivores, may have been caused by stratification from low-salinity spring water from melting sea ice in the north and east (Katsuki and Takahashi, 2005). The intervals low in phytodetritivore abundance may indicate episodes of little spring stratification and associated seasonal blooms.

There may be different possible causes of the sporadic, relatively short-lived episodes of very high productivity and organic carbon flux at U1342, as indicated by Assemblage 1 and the laminated intervals (Fig. 4). (1) During warming (deglacials), a widespread reduction in regional subarctic North Pacific sea ice occurred, which shut down GNPIW formation (Kender et al., 2018; Worne et al., 2019) and promoted the expansion of high-nutrient–low-oxygen NPDW to intermediate and nutricline depths (Site U1342), aiding productivity. This could be comparable to the scenario proposed for the last deglacial Bølling–Allerød (e.g. Shibahara et al., 2007; Gray et al., 2018) in which laminated sediments appear to have been formed across the North Pacific and Bering Sea. In addition, melting ice sheets could have delivered an increase in surface nutrients (Katsuki and Takahashi, 2005). The deglacial collapse of the GNPIW scenario could only explain high Assemblage 1 at the termination of Marine Isotope Stages (MISs) 2, 8, 14, and 15 in our records (Fig. 4), as the remaining high-productivity episodes occurred during glacials. (2) During colder intervals (early glacials), enhanced wintertime sea ice formation in the northeast Bering Sea could have aided the production of springtime low-salinity stratification over Bowers Ridge from melting sea ice, causing an associated local increase in primary productivity (Katsuki and Takahashi, 2005). An increase in productivity large enough to produce Assemblage 1 and laminated sediment may have only been possible at times when GNPIW was not yet extensive (compared with glacial maxima) to allow for the mixing of high-nutrient intermediate water sourced from NPDW into the photic zone.

Our conceptual model indicates that both very high-productivity episodes (Assemblage 1 and laminated sediment) and longer periods of bloom productivity (phytodetritivores in Assemblages 2 and 3) required spring–summer stratification from sea ice melt elsewhere in the Bering Sea, as glacial air masses were likely a lot colder than at present (Gray et al., 2018) and less able to cause strong summer thermal stratification. We observe that Assemblage 1 is more prevalent in the earlier part of the record (from 400–600 ka), and our increased ice melt–productivity interpretation is consistent with overall higher benthic foraminifera abundance (Fig. 2c) and higher abundances of A. weddellensis (Assemblage 2) in the earlier part of the record. During the Last Glacial Maximum A. weddellensis was very abundant (60 % of the fauna) over the eastern Bowers Ridge (BOW-9A, Fig. 1) and Umnak Plateau (UMK-3A) (Okazaki et al., 2005), and both of these locations were interpreted as being influenced by sea ice (Fig. 1) (Katsuki and Takahashi, 2005). In contrast, we find very low abundances of A. weddellensis at western Bowers Ridge (U1342) during the LGM (Fig. 4), where open water conditions were likely present during glacials over the past 270 kyr (Katsuki and Takahashi, 2005). Older diatom records do not exist from Bowers Ridge, but low-resolution records from Bering slope Sites U1343 and U1344 do show a reduction in pack-ice-related diatom species from ∼500 ka (Stroynowski et al., 2017), consistent with our suggestion. This time period approximately coincides with the MBT, for which alkenone-based sea surface temperatures indicate that North Pacific (Site 882; Fig. 1) glacials warmed (Martínez-Garcia et al., 2010), and Mg∕Ca intermediate water temperatures indicate that Arctic Ocean glacials warmed and sea ice became more seasonal (Cronin et al., 2017) (see Fig. 4). Our records may therefore suggest that the Bering Sea responded to regional warming during the MBT after initial cooling and sea ice expansion in the mid-Pleistocene (Stroynowski et al., 2015, 2017).