Intermediate- and Deep-Water Oxygenation History in the Subarctic North Pacific During the Last Deglacial Period

Deglacial dissolved oxygen concentrations were semiquantitatively estimated for intermediate and deep waters in the western Bering Sea using the benthic foraminiferal-based transfer function developed by Tetard et al. (2017), Tetard et al. (2021a). Benthic foraminiferal assemblages were analyzed from two sediment cores, SO201-2-85KL (963 m below sea level (mbsl), the intermediate-water core) and SO201-2-77KL (2,163 mbsl, the deep-water core), collected from the Shirshov Ridge in the western Bering Sea. Intermediate waters were characterized by an oxygen content of ∼2.0 ml L−1 or more during the Last Glacial Maximum (LGM)–Heinrich 1 (H1), around 0.15 ml L−1 during the middle Bølling/Allerød (B/A)–Early Holocene (EH), and a slight increase in [O2] (∼0.20 ml L−1) at the beginning of the Younger Dryas (YD) mbsl. Deep-water oxygen concentrations ranged from 0.9 to 2.5 ml L−1 during the LGM–H1, hovered around 0.08 ml L−1 at the onset of B/A, and were within the 0.30–0.85 ml L−1 range from the middle B/A to the first half of YD and the 1.0–1.7 ml L−1 range from the middle to late Holocene. The [O2] variations remind the δ18O NGRIP record thereby providing evidence for a link between the Bering Sea oxygenation at intermediate depths and the deglacial North Atlantic climate. Changes in the deep-water oxygen concentrations mostly resemble the deglacial dynamics of the Southern Ocean upwelling intensity which is supposed to be closely coupled with the Antarctic climate variability. This coherence suggests that deglacial deep-water [O2] variations were primarily controlled by changes in the circulation of southern-sourced waters. Nevertheless, the signal from the south at the deeper site might be amplified by the Northern Hemisphere climate warming via an increase in sea-surface bioproductivity during the B/A and EH. A semi-enclosed position of the Bering Sea and sea-level oscillations might significantly contribute to the magnitude of oxygenation changes in the study area during the last deglaciation. Interregional correlation of different proxy data from a wide range of water depths indicates that deglacial oxygenation changes were more pronounced in the Bering and Okhotsk marginal seas than along the open-ocean continental margin and abyssal settings of the North Pacific.


INTRODUCTION
Over the past few decades, modern climate change has led to a steady decrease in oxygen concentrations in the World Ocean due to temperature increases and the slowdown of ocean-atmosphere gas exchange associated with global warming (IPCC, 2014). Declines in oceanic oxygen content result in the expansion of oxygen minimum zones (OMZs), which are defined as [O 2 ]deficient layers in the water column (Paulmier and Ruiz-Pino, 2009) and characterized by oxygen concentrations lower than 0.50 ml L −1 (Gilly et al., 2013). As the expansion of OMZs significantly affects global marine species distribution, their habitats, and fisheries' productivity (IPCC, 2014), a deep understanding of the potential climate-associated drivers of the OMZs development is crucial. Here, we present the first timeseries of semi-quantitative reconstructions of the intermediateand deep-water oxygen concentrations in the Bering Sea and, in combination with previously published data, discuss the causes of abrupt changes of oxygen content in the water column in the North Pacific since the Last Glacial Maximum (LGM).
The triggers of OMZ variability are still under debate. Many studies propose that the OMZ temporal and spatial variations are associated with the interplay of three main processes: 1) the strength in upper ocean productivity and associated downward fluxes of carbon, 2) subsurface oxygen consumption by biota during biological remineralization of organic matter and 3) changes in water mass circulation renewing the oxygen content (e.g., Wyrtki, 1962;Mix et al., 1999;McKay et al., 2005;Paulmier and Ruiz-Pino, 2009;Moffitt et al., 2014;Moffit et al., 2015;Praetorius et al., 2015). Biological oxygen utilization is related to bacterial activity during the organic matter degradation process both in the water column and on the seafloor after its deposition. Variability in ocean circulation is determined by intermediate-or deep-water formation through the sinking/ upwelling of water masses from the sea surface/deep ocean or its lateral advection from the other areas. Besides, abrupt ocean warming is proposed to be the driver of the OMZ expansion via the decrease in oxygen solubility and the increase in sea surface bioproductivity in the Gulf of Alaska during the last deglaciation (Praetorius et al., 2015).
The intensity of North Pacific Intermediate Water (NPIW) formation has been considered a major player in explaining ventilation changes in the Bering and Okhotsk marginal seas. Multiproxy data suggest the Bering Sea was a major source of the NPIW precursor during severe cold intervals of the Last Glacial Maximum [LGM: 23-19 ka BP (Hughes et al., 2013)], Early Deglacial [ED: 19-17.7 ka BP (Andrews and Voelker, 2018;Praetorius et al., 2020)] and Heinrich I [H1: 17.7-14.6 (Andrews and Voelker, 2018]; Tanaka and Takahashi, 2005;Rella et al., 2012), whereas the contribution of intermediate waters from the Okhotsk Sea to NPIW played a larger role during warmer times, like today (Shcherbina et al., 2003). This switch between the source areas has reportedly coupled with changes in prevailing wind direction that created favorable conditions for the sea-ice formation and associated brine rejection since at least the last glacial times (Rella et al., 2012).
Previous studies provide evidence for non-synchronous intermediate-and deep-water oxygenation changes in the northeastern Pacific (Mix et al., 1999;Jaccard and Galbraith, 2012;Belanger et al., 2020). In the northwestern Pacific, variations in radiocarbon-derived ventilation ages together with the epibenthic δ 13 C Cibicidoides wuellerstorfi signal demonstrated the antiphase ventilation history of intermediate (700-1,750 m) and deep (>2,100 m) waters during the last deglaciation (Keigwin, 1998;Matsumoto et al., 2002;Max et al., 2014).
Here, we study BF assemblages in comparison with previously obtained results from two sediment cores retrieved in 2009 during the SO201-KALMAR Leg 2 cruise of the R/V "Sonne" (Max et al., 2012;Ovsepyan et al., 2013;Max et al., 2014;Ovsepyan et al., 2017;Riethdorf et al., 2013a, Riethdorf et al., 2013b. Based on our results, available published data, and transfer functions (Tetard et al., 2017;Tetard et al., 2021a), we estimate millennial-scale variations in oxygen concentrations at present-day water depths (w.d.) of 963 and 2,163 m on the Shirshov Ridge in the western Bering Sea, from 22 to 9 ka BP and from 21 to 1 ka BP, respectively. In combination with several independent proxies, our data suggest a link between OMZ dynamics and sea surface bioproductivity and ventilation changes in the western Bering Sea. In addition, we analyze spatiotemporal oxygen variations in the North Pacific during the last deglaciation.

OCEANOGRAPHIC SETTING
The Bering Sea is the northernmost marginal sea of the Pacific Ocean located between the Far East of Russia and Alaska and bounded by the Commander and the Aleutian Islands in the FIGURE 1 | (A) Map of the Pacific Ocean with modern (arrows) water mass circulation (after (Tomczak and Godfrey, 1994;Kawabe and Fujio, 2010;Rella and Uchida, 2014)), modern austral winter (orange shaded area), and modern austral summer (brown dashed line) position of the Southern Hemisphere Westerlies Wind Belt (after (Rella and Uchida, 2014)) and the key core locations: 85KL, 77KL (this study) and MD07-3088 (Siani et al., 2013). AABW, Antarctic Bottom Water; AAIW, Antarctic Intermediate Water; NPIW, North Pacific Intermediate Water; NPDW, North Pacific Deep Water. Pink arrow marks the flow direction of the NPIW precursor. (B) Core locations, surface circulation (red arrows; Stabeno et al., 1999), and modern position of sea ice margin (dashed lines) in the Bering Sea  and the Sea of Okhotsk (Rostov et al., 2001). AA', BB', and CC' are the profiles of oxygen concentrations across the Bering Sea and northwestern Pacific, Sea of Okhotsk, and northeastern Pacific, respectively. AS, Alaskan Stream; AC, Alaska Current; BSC, Bering Slope Current; EKC, East Kamchatka Current; OG, Okhotsk Gyre; ESC, East Sakhalin Current; WSAG, Western Subarctic Gyre. (C) Oxygen concentrations in the water column in the Bering Sea and northwestern Pacific along the profile AA'. (D) Oxygen content in the water column in the Sea of Okhotsk along the profile BB'. (E) Oxygen concentrations in the water column in the northeastern Pacific along the profile CC'. The map was created using Ocean Data View software (Schlitzer, 2020). Bathymetry according to GEBCO Compiling Group (2020). south. Approximately half of the sea belongs to the shelf area, with depths of 0-200 m, whereas the other part is represented by a deep basin, the Aleutian Basin, with a maximal water depth of 3,900 m. Two submarine rises, Shirshov Ridge and Bowers Ridge, are located in the western and southern parts of the Bering Sea, respectively, ( Figures 1A,B). Numerous passages between the islands enable warm surface waters from the North Pacific to reach the Bering Sea. These surface waters are involved in the counterclockwise flowing gyre stirring the water masses along the continental margin around the deep part of the sea. The shallow Bering Strait (∼50 m) is the gateway for exporting water and heat to the Arctic. The modern transport of the Bering Sea surface waters via the Bering Strait is about 1 Sv (Woodgate, 2018).
The interaction of the Aleutian Low and Siberian High atmospheric pressure cells results in a large contrast between summer (8-10°C) and winter (0-2°C) sea surface temperatures (Arseniev, 1967;Miura et al., 2002;Locarnini et al., 2006). As a result, a mixed layer with seasonal thermocline depths of 25-50 m forms in the Bering Sea (Arseniev, 1967;Miura et al., 2002). The occurrence of an underlying cold dichothermal layer is related to brine release during the winter sea ice formation. Its lower modern boundary marks a permanent thermocline and halocline at the water depth of 150-200 m along 57°N (Miura et al., 2002). Due to the low salinity of surface and subsurface waters, no deep-water formation occurs in the North Pacific today (Warren, 1983).
Just below the permanent halocline, the water mass with the potential density (σ θ ) of 26.8 kg m −3 occurs at 250-300 m below sea level (mbsl). The potential density mentioned here is a value calculated with in situ salinity, potential temperature, and pressure equal to zero, minus 1,000 kg m −3 . The potential density value of 26.8 kg m −3 is typical of the low saline NPIW occupying 300-800 m of the water column in the North Pacific today (Talley, 1993;Yasuda et al., 1997). Below 250-300 mbsl, the Bering Sea interior is filled by a homogeneous nutrient-rich mesothermal layer, represented by North Pacific Deep Water (NPDW, Kawabe and Fujio, 2010), which is originated from a diapycnal mixing of Antarctic Bottom Water (AABW, Talley, 2013) with overlying waters in the northeastern North Pacific.
Diatoms are the most important primary producers in spring (April-May) and coccolithophores-in early autumn (September) (e.g., Springer et al., 1996;Harada et al., 2012). Maximal values of primary production (175-275 g C m 2 y −1 ) are concentrated along the shelf-slope front within the so-called "Green Belt" (Springer et al., 1996). The sea surface above the deep western Bering Sea is a typical high-nutrient lowchlorophyll-a area due to iron limitations; therefore, primary production usually does not exceed 50-100 g C m 2 y −1 .
In November, sea ice appears within the coastal polynyas near the islands on the northern and eastern shelves (Niebauer et al., 1999), while the deep-sea area remains ice-free during the entire year. Sea ice completely disappears between July and September (Niebauer, 1980). The modern OMZ is developed between 500 and 1,700 mbsl, with the core depth at 900-1,000 mbsl being more pronounced in winter (Figures 1C-E, Paulmier and Ruiz-Pino, 2009;Garcia et al., 2014). Today, the average oxygen concentrations are about 0.34 ml L −1 at 960 mbsl and about 1.6-1.7 ml L −1 at 2,100 mbsl (Garcia et al., 2014).

Core Locations and Stratigraphy
We used two piston cores, 170°24.79' E,w. d. 968 m,170°41.97' E,w. d. 2,163 m,hereafter 77KL), recovered from the Shirshov Ridge in the western Bering Sea during the R/V Sonne cruise SO201 in 2009 (Figure 1; Dullo et al., 2009). Cores 85KL and 77KL were collected from the intermediate-and deep-water levels, respectively. In Core 85KL, the uppermost 16 cm-interval was disturbed during the core processing. We excluded this part from further studies.
The age models of the deglacial part of the sections are provided by Max et al. (2012). In general, both records are part of a stratigraphic framework of sediment records from the northwest Pacific realm, which is based on radiocarbon dating and regional correlation of high-resolution XRF (Ca intensities) and spectrophotometric (color b*) data. Seven AMS 14 C dates have supported this correlation for Core 85KL and five AMS 14 C dates for Core 77KL (Max et al., 2012). The authors transferred the AMS 14 C ages into calendar ages by assuming the surface reservoir age of 700 years and using the Intcal09 atmospheric calibration curve (Reimer et al., 2009).
In this study, we recalculated radiocarbon ages using the IntCal20 calibration curve (Reimer et al., 2020), assuming the variable reservoir corrections derived from oceanic radiocarbon simulations (Butzin et al., 2017; Table 1). One additional radiocarbon date from Core 77KL was measured at the MICADAS radiocarbon laboratory, Alfred Wegener Institute (AWI) ( Table 1). Sample preparation methods and operational procedures are described in Mollenhauer et al. (2021). The age of the core top sample is considered to be 700 years BP, according to the modern reservoir age (Max et al., 2012). Linear interpolation was used to calculate ages between calibrated radiocarbon dates. Two visible tephra layers were found at 116-117 and 145 cm in the upper part of Core 77KL, whereas one dark tephra layer is reported from the 47-55 cm interval in Core 85KL (Dullo et al., 2009;Ponomareva et al., 2013;Derkachev et al., 2018). It is strongly disturbed by coring and reaches down to 67 cm (Dullo et al., 2009). Although an abrupt deposition of the tephra layers interrupts normal marine sedimentation, the 1-cm thick deposits hardly affect our results obtained with the lower depth resolution from Core 77KL. The strongly disturbed tephra layer is 8 cm thick (47-55 cm) in Core 85KL. However, the thickness of the same ash layer in the nearby cores does not exceed 4 cm (Derkachev et al., 2018). This thickness is also lower than the depth resolution of the proxy data used in this study. Thus, we determine ages for proxy time series from the 47-55 cm interval by linear interpolation.
Age models for previously published cores from the North Pacific were revised by the recalibration of available radiocarbon dates implemented in PaleoDataView (Lagner and Mulitza, 2019) using IntCal20 (Reimer et al., 2020). For most cores, we used  Table S7). The plots of proxy records against ages were created using the revised age-depth models.

Microfossil Assemblages
BF assemblages have been previously studied in >63 µm size fraction from Core 85KL with a 5 cm-step (time resolution of 200-600 years; Ovsepyan et al., 2013). For each sample, 250-300 specimens were counted in the whole sample or suitable aliquots if necessary (Murray, 2006). Planktic foraminifers (PF) were determined in the >100 µm size fraction from the same core depths. At least one hundred tests (shells) were counted for each sample. We follow the same approach with Core 77KL to make the results comparable for both cores. In Core 77KL, the sampling interval of 5 cm provides a time resolution of 100-600 years. The BF absolute abundances were estimated as tests per gram of dry bulk sediment. The accumulation rates of benthic tests were calculated using the dry bulk density values sediment and sedimentation rates. Absolute abundances and accumulation rates of BF were used as proxies of organic matter flux supplied from the photic layer (Gooday, 2003;Jorissen et al., 2007). BF-to-PF ratio was applied as a dissolution indicator (Dittert et al., 1999).

Statistical Analyses and Semiquantitative Oxygen Reconstructions
To assess the past dissolved oxygen concentrations, we used the transfer function technique developed by Tetard and co-authors (2017), which is based on the BF data from core MD02-2508 retrieved off Baja California (Tetard et al., 2017) and later recalibrated using 45 worldwide core-tops of which the modern [O 2 ] is known (Tetard et al., 2021a). This approach can be divided into three main steps. First, a principal component analysis (PCA) was conducted on matrices of relative abundances of BF species to obtain the dominant taxa in the assemblages. As oxygen is the main factor controlling the distribution of BF in the OMZ areas (e.g., Jorissen et al., 1995), the dominant species obtained by the PCA were then assigned to dysoxic (0.1-0.5 ml L −1 ), suboxic (0.5-1.4 ml L −1 ) and oxic (>1.4 ml L −1 ) groups (oxygen ranges after (Moffitt et al., 2014;Moffitt et al., 2015)) based on their oxygen affinities found in the literature. Next, a correspondence analysis (CA) was carried out to check if the dominant species are oriented along Axis 1, assuming that oxygen content is the major factor controlling the species distributions at our sites. For every sample, the reconstructed oxygen value was determined by a transfer function that depends on the relative abundance of the dysoxic and oxic assemblages (Tetard et al., 2017;Tedard et al., 2021a). All the statistical treatments were carried out using the software package PAST (PAleontological STatistics; Hammer et al., 2001). Based on the approach described above, deglacial oxygen concentrations were also calculated for the Sea of Okhotsk using published BF data obtained from core MD01-2415 (Bubenshchikova et al., 2015) to compare the deglacial evolution in oxygenation histories between the Bering Sea and the Sea of Okhotsk, the two major sources of NPIW precursors.  (Max et al., 2012) and new AMS 14 C dates from the sediment cores SO201-2-85KL and SO201-2-77KL calibrated using IntCal20 (Reimer et al., 2020) with applied variable reservoir corrections (Butzin et al., 2017). The calendar dates are given as calibrated weighted mean age within ±2σ range and calculated using PaleoDataView software (Langner and Mulitza, 2019). Previous calibrated ages from Max et al. (2012) are also shown. Other Proxies for Bioproductivity,

Bottom-Water Oxygenation, and Ventilation
The deglacial changes in element composition along the sediment records were measured using the Avaatech X-ray Fluorescence (XRF) core scanner at the Alfred Wegener Institute for Polar and Marine Research, Bremerhaven. The cores were triple-scanned at 1 mA and tube voltages of 10 kV (Al, Si, S, K, Ca, Ti, Fe), 30 kV (Cu, Zn, Br, Rb, Zr, Sr, Mo) and 50 kV (Ag, Cd, Sn, Te, Ba), using a sampling resolution of 1 cm and 30 s count time (Max et al., 2012). The Mn/Ti ratio derived from XRF measurements was applied as a proxy of changes in bottom-water redox conditions (e.g., Rothwell and Croudace, 2015).
To support foraminiferal-based reconstructions, total organic carbon (TOC) content was used as an independent indicator for assessing the intensity of organic matter flux to the seafloor (e.g., Riethdorf et al., 2013a). Variations in δ 13 C were measured in tests of BF species C. wuellerstorfi in Core 85KL. At Core 77KL, carbon isotopes were measured on mixed Uvigerina akitaensis which is taxonomically similar to U. peregrina (Bubenshchikova et al., 2008)) and U. auberiana. The δ 13 C records were used as a proxy for changes in bottom-water ventilation (Riethdorf et al., 2013b;Max et al., 2014).

Revised Deglacial Chronology of Cores 85KL and 77KL
The modified age model of Core 85KL differs from that previously published by Max et al. (2012), especially within the B/A and H1 intervals where the new calibrated weighted mean ages become ∼700-1,000 years younger than previously published ones (Max et al., 2012;Supplementary Figure S1). During the EH, new ages, on the contrary, appear to be 200-300 years older than the ages from the initial stratigraphy. Revised chronology of Core 77KL demonstrates younger ages (by 300-400 years) within B/A, YD, and EH and older ones (by 400-500 years) during the H1 (Supplementary Figure S1). An increase in age model differences up to 1 kyr occurs in the uppermost and lowermost part of the section, probably due to age extrapolation. Thus, the revised age model of Core 77KL has not been changed substantially for the key millennial-scale intervals (H1, B/A, YD, and EH) of the current study.
Sedimentation rates are calculated between the calendar ages by linear interpolation. Sedimentation rates range within 3-19 cm kyr −1 and 8-55 cm kyr −1 within the deglacial interval of Cores 85KL and 77KL, respectively. According to modified age models, our studied intervals of 147-16 cm in 85KL and 265-0 cm in 77KL correspond to 22-9 ka BP and 21-1 ka BP, respectively.

Benthic Foraminiferal Assemblages and Downcore Trends
The application of PCA reveals that, in Core 85KL, twelve main species with high positive and negative loadings are distributed along the first axis representing principal component (PC) 1 and explaining 66.3% of the variance (Supplementary Table S8 Table S8). Ten dominant species are considered the main taxa, including E. tenuata, Rutherfordoides erectus, T. delicata, Bolivina pacifica, Nonionellina labradorica, C. norvangi, Epistominella arctica, I. norcrossi, A. weddellensis and S. fusiformis. The assignment of these species to oxygen-related groups is shown in Table 2, and the downcore distributions are displayed in Supplementary Figure S2, S3.
Correspondence analysis was applied to each matrix to plot the main species on two axes representing environmental parameters (Tetard et al., 2017). CA Axis 1 explains 44% of the variance in core 85KL and 20% of Core 77KL, whereas CA Axis 2 explains 21% of the variance in Core 85KL and 16% of the variance in Core 77KL. As shown on the plots, all the main species are well grouped according to their oxygen tolerance along the first axes ( Figure 2).
In Core 85KL, the dysoxic group consisting of B. seminuda, E. tenuata and T. delicata displays score values between −1.03 and −1. The suboxic group containing P. naraensis, S. fusiformis, B. spissa and U. akitaensis ranges from −1.09 to 0.42, thereby overlapping the interval of the dysoxic group occurrence. The distributions of the single species throughout Core 85KL demonstrate that P. naraensis varies similarly to the species of the suboxic group; thus, assignment of this species to the suboxic assemblage seems reasonable (Supplementary Figure S2). The oxic group consisting of A. angulosa, P. takayanagii, I. norcrossi, A. weddellensis, and C. norvangi demonstrates score values within the 0.61-0.87 interval. The samples are divided into two groups. The first group contains most samples and is dominated by oxic species, whereas the remaining samples are associated with the suboxic and dysoxic taxa. In Core 77KL, BF groups are evenly distributed along the CA Axis 1 without any overlaps. The score values of the dysoxic (E. tenuata, R. erectus, T. delicata), suboxic (B. pacifica, N. labradorica, S. fusiformis), and oxic (A. weddellensis, C. norvangi, E. arctica, I. norcrossi) group varies within −1.54-−1.24, −0.85-−0.30, and 0.11-1.34, respectively. The samples are divided into three groups. The first corresponds to the samples with the prevalence of I. norcrossi, E. arctica, and C. norvangi. The second group of samples is dominated by A. weddellensis, whereas the third mainly contains the suboxic and dysoxic species. The distributions of the main species for Core 77KL are shown in Supplementary Figure S3.
In Core 85KL, CA Axis 1 demonstrates maximal positive values before the onset of the B/A, tending to decrease toward the end of the B/A and then displaying minimal values during the late B/A, YD, and EH ( Figure 3A). In Core 77KL, CA Axis 1 shows clear minima during the onset of the B/A and the YD-EH timespan ( Figure 3C). Values around zero are documented within the LGM-H1 interval, whereas positive values between 0 and 2 are found during the Holocene.
Semi quantitatively reconstructed near-bottom oxygen concentrations exhibit relatively high and stable values around 2 ml L −1 during LGM-the onset of the B/A at the core site 85KL, whereas they show high-amplitude variations within about 0.90-2.50 ml L −1 during LGM-H1 at the deeper core site 77KL ( Figures 3A,C). Very low oxygen content between 0.07 and 0.50 ml L −1 is documented for the middle B/A-EH interval, with a slight increase in oxygenation at the beginning of the YD at the core site 85KL. At the deeper location (site 77KL), the short-term minimum of dissolved oxygen concentrations (0.08 ml L −1 ) is estimated at the beginning of B/A. Then the oxygen content increases up to 0.85 ml L −1 in the middle of YD and, after the second minimum of 0.30 ml L −1 during the end of YD-EH, oxygen concentrations gradually increase to present-day values of 1.60-1.70 ml L −1 during the Holocene.
Based on BF data from Core MD01-2415 (Bubenshchikova et al., 2015), oxygen concentrations in the Sea of Okhotsk oxygen vary between 1.20 and 2.0 ml L −1 during H1, 0.50 and 1.20 ml L −1 during the B/A-EH, and 0.50 and 2.0 ml L −1 during the remainder of the Holocene. Foraminiferal Abundance and Accumulation Rates BF abundance and accumulation rates display maxima within the B/A and EH in Core 85KL and during the B/A and YD in Core 77KL ( Figure 4). The BF abundance is more pronounced during the B/A (2.5 × 10 3 tests g −1 ) than during the YD (1 × 10 3 tests g −1 ) in Core 77KL; however, it is more distinct during the EH (2.6 × 10 3 tests g −1 ) than during the B/A (1.5 × 10 3 tests g −1 ) in Core 85KL. During the LGM-H1 and mid-late Holocene intervals, BF abundance and accumulation rates are low and do not exceed 0.6 × 10 3 tests g −1 and 3.7 × 10 3 tests cm −2 kyr −1 , respectively, in both cores. BFto-PF ratio demonstrates two maxima during the B/A and EH in Core 85KL. Maximal values of this proxy are associated with the H1, YD, and mid-Holocene ( Figure 4). Variable values of BF-to-PF ratio are documented from Core 77KL.

XRF Data
Variations in the Mn/Ti ratio exhibits two spikes corresponding to CA Axis 1 minima at the onset of B/A and during the YD in Core 77KL and in the middle of B/A and YD in Core 85KL ( Figures 3C,D). The most pronounced Mn/Ti maxima are related to the YD interval in both cores.

Oxygenation Changes on the Shirshov Ridge During Termination I
Interpretation of the millennial-scale oxygenation changes in the intermediate and deep waters strongly depends on age model uncertainties, which are found to be large, especially at the H1-B/A boundary (1,260 years) and during the YD (1,000 years) in the age model of Core 85KL (Figure 3). About 1,200 and 1,300 years uncertainties are documented within the B/A and YD, respectively, in the revised stratigraphy of Core 77KL (Figure 3). These large uncertainties are derived from the large surface reservoir ages modeled for these core sites for each millennial-scale period (Butzin et al., 2017). As a result, two main scenarios of intermediate- Red stars correspond to modern-day values. Black bars indicate the calendar age uncertainty ranges which have been recalculated from the AMS 14 C dates from Max et al. (2012). Numbers indicate the weighted means of the calendar ages; raw AMS 14 C dates are given in brackets. Color bars mark warm (light orange) and cold (light blue) intervals, according to the North Atlantic climatostratigraphy. Note that the AMS 14 C dates were converted into calendar ages by assuming variable surface reservoir ages, according to Butzin et al. (2017), and the Intcal20 atmospheric calibration curve (Reimer et al., 2020). The means and AMS 14 C dates are written in kiloyears. Red number with a star highlights age reversal. The red and light blue shades represent the 23% uncertainty of the estimates for cores SO201-2-85KL and SO201-2-77KL, respectively, (Tetard et al., 2021a). B/A, Bølling/Allerød; YD, Younger Dryas; EH, Early Holocene.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 and deep-water oxygenation history might occur in the western Bering Sea during the last deglacial interval. The first scenario includes rather synchronous deoxygenation events at both sites during the mid-B/A and EH intervals, with a slight increase in oxygen content during the YD. This scenario might be possible if the calendar age of the deoxygenation spike would be the oldest within its uncertainty range in Core 85KL, while the age of the oxygen minimum from Core 77KL should be the youngest ( Figure 3). If so, the synchronicity of oxygen-deficiency events at different depths in the Bering Sea during the B/A and EH is in line with most North Pacific reconstructions (e.g., Zheng et al., 2000;van Geen et al., 2003;Cook et al., 2005;Ikehara et al., 2006;Shibahara et al., 2007;Sagawa and Ikehara, 2008;Ishizaki et al., 2009;Itaki et al., 2009;Kim et al., 2011;Rella et al., 2012;Schlung et al., 2013;Kuehn et al., 2014;Knudson and Ravelo, 2015;Asahi et al., 2016;Tetard et al., 2017). The second scenario assumes different oxygenation histories at the various water depths with non-synchronous deoxygenation events. Considering the fact, uncertainty ranges of the radiocarbon dates measured within the CaCO 3 spikes within the B/A are significantly shifted relative to each other, the second scenario seems to be more realistic, at least for this millennialscale event. Indeed, the uncertainty intervals indicate that the dysoxic event most likely appeared within the second half of the B/A at the core site 85KL, whereas oxygen deficiency was pronounced earlier at the core site 77KL. Moreover, the shapes of the curves are dissimilar as well. We could also see the disagreement between oxygenation records when the previous age models were applied (Max et al., 2012; Supplementary Figure S1). The second scenario concurs with an antiphase intermediate-and deep-water ventilation history in the North Pacific during the last deglaciation (Keigwin, 1998;Matsumoto et al., 2002;Max et al., 2014) and seems to us more plausible.
The chronologies of our cores are based on the calibrated weighted mean ages. According to these age models, bottom water conditions were mostly oxic during the LGM-the onset of the B/A, suboxic at the mid-B/A, and dysoxic during the end of B/A-EH with a slight increase in [O 2 ] (up to 0.20 ml L −1 ) at the B/A-YD boundary and EH on the top of the Shirshov Ridge (core site 85KL) (Figure 3). High values of the Mn/Ti ratios suggest the precipitation of manganese in a solid phase during the B/A and FIGURE 4 | Oxygen estimates, productivity, ventilation, and preservation proxies from Core SO201-2-85KL (A) and Core SO201-2-77KL (B) relative to the Northern Greenland ice core record (NGRIP δ 18 O, AICC2012 timescale, Veres et al., 2013) and the Antarctic isotopic signatures (EDML δ 18 O, EPICA community members, 2006), TOC records after (Riethdorf et al., 2013a), benthic foraminiferal abundance (BF) and accumulation rates (BFAR) from Core SO201-2-85KL, according to (Ovsepyan et al., 2013(Ovsepyan et al., , 2017, BF-to-PF ratio, δ 13 C C.wuellerstorfi from Core SO201-2-85KL after (Max et al., 2014) and δ 13 C U. akitaensis/auberiana from Core SO201-2-77KL according to (Riethdorf et al., 2013b). Note that the North Atlantic climatostratigraphy has been applied to both groups of proxy records. Color bars mark warm (light orange) and cold (light blue) intervals, according to the North Atlantic climatostratigraphy (AICC2012 timescale, Veres et al., 2013). The red shades represent the 23% uncertainty of the estimates for cores SO201-2-85KL and SO201-2-77KL (Tetard et al., 2021a). Acronyms as in Figure 3.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 YD. One could interpret this manganese enrichment as evidence of an increase in oxygenation of the bottom waters (e.g., Brumsack, 2006). However, Mn could be adsorbed to biogenic carbonate or incorporated into the calcareous shells, also in the restricted basins, like the Bering Sea, with a widespread oxygen deficiency (Brumsack, 2006). Therefore, low BF-based oxygen concentrations and coeval increases of the Mn/Ti ratios might point to a strongly oxygen-depleted water mass during the B/A and most of the YD at the intermediate depths in the western Bering Sea. Like the intermediate waters, the deep-water mass of the western Bering Sea was characterized by presumably oxic conditions during the LGM-H1, from the middle B/A to the first half of YD, and during the middle-late Holocene (Figure 3). The dysoxic environmental conditions with oxygen concentrations of 0.30-0.40 ml L −1 occurred at the present-day depth of 2,163 mbsl on the Shirshov Ridge during the second half of the YD-EH interval, whereas extremely dysoxic (0.08 ml L −1 ) conditions appeared for a short time at the onset of the B/A (Figure 3). The maxima of the Mn/Ti ratios are coeval with the estimated past oxygen concentration minima, providing additional evidence of the low-oxygenated environment in the deep western Bering Sea at the onset of B/A and during the YD.
Variations in organic carbon flux to the seafloor and changes in oceanic circulation are known to be the key factors in regulating the subsurface oxygen contents (e.g., McKay et al., 2005). In addition, ocean temperatures affecting the [O 2 ] solubility have reportedly played an important role in a global ocean oxygenation state during Termination I (e.g., Jaccard and Galbraith, 2013;Praetorius et al., 2015). The colder the water, the more oxygen can be dissolved. In the marginal seas of the subarctic Pacific, intermediate-water formation is related to brine rejection during the development of winter sea-ice cover (e.g. Arseniev, 1967;Miura et al., 2002;Shcherbina et al., 2003). Since sea ice constantly forms at the freezing point, the oxygen solubility effect cannot be considered the main reason for explaining oxygenation changes in the newly formed intermediate waters from the Bering Sea. Thus, it is very likely that the dissolved oxygen content in the Bering Sea mainly results from the interplay of organic matter flux from the surface layer and the state in intermediate and deep water circulation.
High values of TOC, BF abundance, and accumulation rates point to a high export production during the B/A and EH at both sites (Figure 4), suggesting that the organic matter flux partly controlled the oxygenation state during these periods. However, these proxies cannot explain the oxygen deficiency of the intermediate-and deep-water masses during the YD. The ocean circulation dynamics may have played a leading role in forming a low-oxygenated environment during the YD, which will be discussed below. Mix et al. (1999) suggested a hypothesis about different mechanisms (Northern versus Southern Hemisphere influence) controlling the intermediate (980 m) versus deep (2,700 m) water ventilation in the northeastern Pacific since the last glacial period (Mix et al., 1999). Authors showed that enhancements of intermediate-water ventilation during the North Pacific warming events are synchronous with warm B/A and EH episodes in the North Atlantic, suggesting atmospheric transmission of climate signals across the Northern Hemisphere (Mix et al., 1999). Better ventilation of the deep layer reportedly coincided with the Antarctic Cold Reversal (ACR, 14.7-13 ka BP (Pedro et al., 2015)), providing evidence for interhemispheric teleconnection via intensively formed AABW during the Southern Ocean cooling (Mix et al., 1999). Proximal and distal sources of the intermediate-and deepwater masses are proposed to explain the discrepancy of oxygenation changes at different water depths in the northeastern Pacific (Kennett and Ingram, 1995). Finally, Max and co-authors (2014) (Figure 4). On the one hand, the 85KL oxygen record resembles the Northern Hemisphere climate variability recorded in δ 18 O ice core records from Greenland, displaying an extreme oxygen depletion during the mid-B/A and EH warming intervals. Episodes of the weakening or disappearing of the OMZ in the subarctic Pacific at intermediate depths are partly associated with cold periods in the North Atlantic. These trends are described in previous studies from the northeastern (Kim et al., 2011) and southern Bering Sea (Schlung et al., 2013). On the other hand, dysoxic conditions during the YD make the 85KL record not a typical "Greenland"-like pattern. The main drivers of the dysoxic environment formation in the intermediate waters during the YD will be discussed below.

Changes in oxygen concentrations at intermediate depths in the western Bering Sea demonstrate dysoxic conditions during the B/A, YD, and EH with extreme deoxygenation events during the B/A and EH
At the deeper site, a pronounced [O 2 ] minimum is detected at the very beginning of the B/A, whereas during the EH, oxygen depletion is noticeable but modest compared to the B/A in terms of absolute values. One could attest this record with two lowoxygenated events to the "Greenland"-type variability; however, the younger oxygen minima appeared in the middle of the YD. Moreover, the oxygenation record does not match TOC data representing the "Greenland"-like pattern. The evidence mentioned above suggests another primary driving mechanism of the deep-water [O 2 ] changes that might interact with Greenland's signal during the warm B/A and EH.
The Greenland temperature rose. On the one hand, the offset in the timing might be explained by stratigraphic uncertainties which could be even larger than those assumed in this study. On the other hand, this shift might be linked to the gradual intensification of the oxygen consumption in the Bering Sea's intermediate waters due to the abrupt warming in the Northern Hemisphere. However, sediment records from the northeastern Bering Sea (Kim et al., 2011) provided no corresponding shift in oxygenation proxies during the B/A, which might be explained by the non-synchronous reaction of the intermediate-water oxygenation to the climate warming in the western and eastern parts of the sea. However, this discrepancy demands further investigation. The TOC variations indicate low export production at the core sites 85KL and 77KL during the LGM-H1 and YD and enhanced organic matter flux to the seafloor during the B/A and EH (Figure 4; Riethdorf et al., 2013a). Previously published BF and PF data concur with the TOC record and further point to high/low sea-surface bioproductivity during the warm/cold intervals in the Bering Sea (Kohfeld and Chase, 2011;Riethdorf et al., 2013a;Ovsepyan et al., 2013, Ovsepyan et al., 2017. During the warmer intervals, high sea surface bioproductivity might be explained by high sea surface temperatures (Max et al., 2012) and the absence of sea ice cover (Max et al., 2012;Méheust et al., 2016). A breakdown of deep-ocean stratification and the subsequent upward flux of nutrients (nitrates and phosphates) within nutrient-rich deep waters (Galbraith et al., 2007;Jaccard and Galbraith, 2012;Max et al., 2014) might have further stimulated the primary production during the B/A and EH. Numerical model simulations indicate the shoaling of the global thermocline during a reinvigoration of the North Atlantic Deep Water (NADW, Broecker, 1991) formation in the North Atlantic during the B/A and EH (Schmittner, 2005) and also support the hypothesis of reduced ventilation of the North Pacific intermediate depths during the deglacial warm intervals. The enhanced freshwater supply from the continental runoff (Kim et al., 2011;Rella et al., 2012) and intensive delivery of the Alaskan Stream warm waters through numerous passages likely maintained the surface ocean stratification. Thus, these factors collectively created favorable conditions for phytoplankton blooms and high export production events to the seafloor in the western Bering Sea during the B/A and EH.
BF-to-PF ratio displays high values indicating carbonatecorrosive environment at the intermediate depths during high organic matter input events, whereas this proxy suggests good preservation and less corrosive conditions at the deeper site ( Figure 4). As we discussed above, a breakdown of the North Pacific deep-water stratification occurred during the B/A and EH (Max et al., 2014), leading to outgassing of the carbon dioxide from the deep ocean and the uplift of nutrients to the sea surface (Galbraith et al., 2007;Jaccard and Galbraith, 2012). It might be possible that CO 2 outgassing resulted in an increase in the preservation of calcareous microfossils at both sites; however, this effect was compensated by intensive biological altering of sinking organic matter at the shallower depths.
During cold episodes, nutrient and light limitations led to reduced sea surface bioproductivity and a diminished carbon flux to the seafloor. Expanded sea ice cover and a long winter season may have prevented a massive phytoplankton bloom in the euphotic layer during spring. Low export production might partly justify the high oxygenation of the intermediate depths during cold intervals (e.g., Kim et al., 2011). Epibenthic δ 13 C C.wuellerstorfi record and the radiocarbon-based ventilation ages suggest the occurrence of the nutrient-poor, wellventilated intermediate water mass in the northwestern Pacific during the cold intervals of H1 and YD (Max et al., 2014). Matsumoto et al. (2002) reported expanding the glacial mode of NPIW in the North Pacific down to 2 kmbsl during glacial times. This expansion further explains relatively high BF-based oxygen concentrations during the H1 (Figure 4). The appearance of newly formed intermediate waters, as a precursor of glacial NPIW, was proposed to be a result of brine rejection during winter sea-ice formation around the polynyas on the northeastern Bering Sea shelf forming under the influence of northerly winds during the cold climates (Rella et al., 2012;Knudson and Ravelo, 2015). An enhancement of intermediate-water formation might also be favored by increasing sea surface salinity due to a decrease in precipitation over the basin under severe conditions (Rella et al., 2012).
However, neither low TOC content nor high δ 13 C C. wuellerstorfi cannot fully explain low [O 2 ] values at the core site 85KL during most of the YD (Figure 4). Indeed, TOC values were lower during the YD than during the LGM-H1, but oxygen deficiency was an order of magnitude greater during the YD (about 0.20 ml L −1 ) than during the onset of deglaciation (about 2.0 ml L −1 ) (Figure 3). At the same time, δ 13 C C.wuellerstorfi values during the YD are lower than those during the H1, possibly suggesting the less intensive formation of intermediate waters. In addition, the higher alkenone-based sea surface temperatures during the YD than at the H1-B/A boundary (Max et al., 2012), as well as relatively light NGRIP oxygen isotope signatures, might serve as evidence for the shortened winter season and consequent restriction of the intermediate-water formation in the western Bering Sea during the YD compared to the LGM-H1. Less intensive production of the intermediate waters might lead to an intrusion of old oxygen-poor NPDW at least to the present-day depth of 963 msbl leading to reduced oxygen concentrations at the core-site.
Thus, changes in ventilation might contribute to intermediatewater oxygenation in the western Bering Sea. A rough covariation of the NGRIP record, the Bering Sea oxygenation, export productivity, and ventilation data (Figure 4) might confirm the previously discussed hypothesis about a close coupling between the climate in the North Atlantic and the intermediate-depth of the North Pacific (e.g., Mikolajewicz et al., 1997;Mix et al., 1999;Max et al., 2012;Schlung et al., 2013;Ovsepyan et al., 2017;Tetard et al., 2017).

Evidence for Southern Hemisphere-Sourced Signal
Downcore TOC records and BF abundance and accumulation rates display two maxima indicating a two-step increase in export production during the B/A and EH intervals at core site 77KL (Riethdorf et al., 2013a; Figure 4). The timing and the amplitudes of changes in BF-based oxygen concentrations partly match the TOC variations, as seen in Figure 4. In particular, the abrupt Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 deep-sea deoxygenation is reconstructed for the onset of B/A, whereas the TOC maximum is documented later (i.e., at the middle of this interval). Moreover, the amplitude of the TOC spike during the EH is equal to that during the B/A in the western Bering Sea; however, the deoxygenation of the deep waters was stronger during the B/A compared to the EH. These findings suggest that organic matter flux variations representing the signal from Greenland contributed to the dysoxic environment formation but seemed to play a minor role in oxygenation state changes at the deep-sea site in the western Bering Sea Termination I. BF-based oxygen concentrations and Mn/Ti records at both sites as well as carbon isotope data from Core 77KL provide evidence for an oxygen-deficient environment in the deep waters during the YD (Figures 3, 4). These low-oxygen conditions cannot be explained by the export production changes and the Northern Hemisphere climate oscillations. Thus, we address the hypothesis about the contribution of deep-ocean circulation to the oxygenation state in the western Bering Sea.
The δ 13 C U. akitaensis/auberiana records from Core 77KL roughly covary with the past oxygen concentrations based on BF assemblages (Figure 4). The periods of stronger ventilation (high δ 13 C U. akitaensis/auberiana values) are associated with higher oxygenation, while the contrary (low δ 13 C U. akitaensis/auberiana values) is observed in periods of weaker ventilation. High δ 13 C U. akitaensis/auberiana values indicate the presence of relatively well-ventilated deep water conditions during the LGM-H1 and middle-late Holocene and coincide with relatively high oxygen concentrations in the bottom waters at the present-day depth of 2,163 mbsl, during the same intervals ( Figure 4). Within the B/A-EH interval, low δ 13 C U. akitaensis/auberiana values are associated with nutrient-rich and poorly oxygenated deep waters; however, they do not match marked [O 2 ] minima, shown by BF-based estimates at the onset of B/A and during the late YD-EH. This discrepancy might be explained by the shallow infaunal living strategy of U. akitaensis and U. auberiana. Accordingly, the δ 13 C values possibly reflect the isotopic composition of pore waters rather than a signal related to rapid changes in bottom water [O 2 ]. Thus, short-term changes in ventilation at the sediment-water interface might be underestimated.
Nevertheless, reconstructed absolute values of the deep-sea oxygen concentrations demonstrate a higher similarity to the ventilation record than the export production variations. Based on benthic-planktic ventilation ages, Max et al. (2014) reported a breakdown of the deep northwestern Pacific stratification and upwelling of poorly oxygenated nutrient-rich waters at least to the intermediate water level during the B/A and EH. Thus, old oxygen-depleted waters most likely occur at the deeper site during these time intervals. Low values of δ 13 C U. akitaensis/auberiana indicate that oxygen-poor NPDW occupied the present-day water depth of 2,163 m during the YD as well. At the intermediate depths, δ 13 C C.wuellerstorfi records demonstrate higher values during the YD (Max et al., 2014), but they are not as high as during the H1, likely pointing to the less intensive intermediate-water formation during the YD and admixing of old poor-ventilated NPDW from below.
It has been suggested that the Southern Ocean is the primary source of 70% of the deep waters in the Pacific (Gebbie and Huybers, 2010); thus, variations in the southern-sourced water formation might be a potential trigger of ventilation and oxygenation changes in the deep subarctic North Pacific in the past, as was previously proposed by many authors (Mix et al., 1999;Galbraith et al., 2007;Jaccard and Galbraith, 2012;Rella and Uchida, 2014;Du et al., 2018;Gorbarenko et al., 2020). A  Veres et al., 2013), benthic foraminiferal oxygen estimates versus δ 13 C C.wuellerstorfi from Core SO201-2-85KL after (Max et al., 2014) from Core SO201-2-85KL (B) and Core SO201-2-77KL (C), ventilation ages ( 14 C B-P ) from Core MD07-3088 from the southeastern Pacific as a proxy of Southern Ocean upwelling intensity (D), (Siani et al., 2013), opal flux from Core TN057-13PC as a proxy of Southern Ocean diatom productivity (E), (Anderson et al., 2009), and EDML δ 18 O as a proxy of Antarctic air temperature changes (F), (EPICA community members, 2006). Light blue bars indicate cold intervals in the Northern and Southern Hemisphere (AICC2012 timescale, Veres et al., 2013;Pedro et al., 2015). Light orange bars highlight warm intervals, according to the North Atlantic climatostratigraphy (AICC2012 timescale, Veres et al., 2013). Yellow bars show episodes of enhanced upwelling in the Southern Ocean, according to modified stratigraphy from Siani et al. (2013). The age models for Cores MD07-3088 and TN057-13PC are based on recalculated AMS 14 C ages using SHCal20 (Hogg et al., 2020) and consistently estimated reservoir ages (For details, see Supplementary Information). The red shades represent the 23% uncertainty of the estimates for cores SO201-2-85KL and SO201-2-77KL (Tetard et al., 2021a). Acronyms as in Figure 3. H1, Heinrich 1; ACR, Antarctic Cold Reversal.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 comparison of the [O 2 ] concentration record at core site 77KL with the Antarctic oxygen isotope record (EPICA community members, 2006) demonstrates that episodes of deep-water deoxygenation coincide with the warming intervals in Antarctica (Figures 4, 5). A coupling between the deglacial deep Bering Sea oxygenation changes and the Antarctic climate variations has been recognized and generally supports an idea about the significant impact of the Southern Oceanrelated processes on [O 2 ] concentrations in the subarctic North Pacific ( Figure 5). If so, different timing of the events might be applied to describe and explain paleoceanographic changes in the deep subarctic Pacific ( Figure 5). We scaled our 77KL oxygenation record to the Antarctic climate chronology. In particular, we considered the two Antarctic warming events of ∼18 to 14.7 ka BP and 13 to ∼10 ka BP, subdivided by a cooling event of ACR which mostly corresponds to the B/A warming in the North Atlantic. As previously discussed, higher temperatures in Antarctica generally coincide with a weakening of deep-sea oxygenation in the subarctic Pacific ( Figure 5). As increasing air temperature decreases oxygen solubility, northward-flowing newly formed Antarctic Intermediate Water (AAIW, Talley, 2013), as a southern-sourced component of the deep Pacific interior, should contain less dissolved oxygen during ∼18-14.7 and 12.3-10 ka BP than during the LGM because the formation of this water mass is not linked to the freezing process. The AAIW is reportedly slightly depleted in oxygen during the Antarctic deglacial warming events than during the cold episodes (Muralti et al., 2010). This oxygen depletion of AAIW was explained by an intensification of upwelling of [O 2 ]-depleted waters in the Southern Ocean, which subsequently downwelled at the Subpolar Front in the Southern Hemisphere (Muralti et al., 2010). The AABW, as a source of NPDW, is also known to have been less oxygenated during the warm intervals as inferred from deglacial Southern Ocean redox-sensitive element records (Jaccard et al., 2016). Lower oxygenation in the upper Pacific Ocean (>1.5 kmbsl) during H1-B/A transition has been reported in a review paper by Jaccard and Galbraith (2012). Du et al. (2018) revealed an acceleration of abyssal Pacific circulation during the Antarctic warming events based on neodymium isotopes and model results. In particular, AABW transport into the Pacific increased from 8 Sv during the LGM to 20-25 Sv during the deglacial Antarctic warmings and then decreased to 14 Sv during the recent time (Du et al., 2018). Assuming the findings mentioned above, the decrease in the dissolved oxygen content in the southern-sourced components combined with the intensified northward transport along the abyssal Pacific might have collectively contributed to the deoxygenation of deep waters in the Bering Sea during the Southern Hemisphere's warming events.
Enhancement of the Southern Ocean's surface bioproductivity would further reduce the oxygen concentrations in the southernsourced waters via the biological respiration of sinking organic particles in the water column and after its deposition on the seafloor. Two clear episodes of higher primary productivity conditions associated with two Antarctic warming intervals were previously inferred from the opal flux records ( Figure 5E, Anderson et al., 2009). The authors hypothesized that sea surface bioproductivity maxima resulted from the twostep increase in the Southern Ocean upwelling rates when nutrients (including Si) were supplied from the deep ocean to the photic layer and produced favorable conditions for phytoplankton blooms (Anderson et al., 2009).
The strength of the Southern Ocean upwelling, which is partly linked to the intensity and position of the Southern Hemisphere westerlies (Toggweiler et al., 2006), might amplify the effect of a decrease in oxygen solubility and increase in primary productivity via changes in the intensity of AAIW and AABW export to the Pacific interior. Based on radiocarbon-derived ventilation ages from Core MD07-3088 retrieved from the southeastern Pacific, Siani and co-authors (2013) clearly showed three episodes of enhanced Southern Ocean upwelling (∼18.2, 17.7-14.5 and 12.5-11.5 ka BP, according to the revised timescale) corresponding to warming episodes in Antarctica. A direct comparison of the Bering Sea oxygenation curves and the Southern Ocean upwelling proxy records reveals significant synchronicity between the time series. Thus, reduced oxygen solubility in AAIW formation, the acceleration of abyssal Pacific circulation, enhanced Southern Ocean upwelling, and increased sea surface bioproductivity might be responsible for a supply of relatively deoxygenated deep waters to the intermediate and deep North Pacific during the Antarctic warming events.
The absence of a synchronous Bering Sea response to the earliest enhancement of the Southern Ocean upwelling (∼18.2 ka BP) might be explained by several factors. First, this Southern Ocean signal could be detected in the Bering Sea, but with a phase shift. If so, one could attribute the response to this signal to the Bering Sea deoxygenation increase at 18 ka BP (i.e., ∼0.2 kyr later than the earliest episode of the Southern Ocean upwelling). The offset might be in the stratigraphic uncertainties within the H1, which might be enlarged due to age extrapolation during age model revision. A delay in the Bering Sea reaction to oceanographic changes in the Southern Ocean was previously proposed by (Ovsepyan et al., 2018). The possible existence of the time shift is supported by the low AABW transport into the Pacific abyss between 20 and 18 ka BP compared to the following deglacial warming intervals (Du et al., 2018). Second, the earliest enhancement of upwelling seemed too weak and too short to be transmitted to the North Pacific.
Furthermore, low opal fluxes indicate a reduced bioproductivity in the Southern Ocean at 18.2 ka BP; thus, the decay of the export production could not affect the oxygenation state of northward transported waters. Third, sedimentation rates in the deep Bering Sea were probably not high enough to fix the response to this short event. Finally, more oxygen could be dissolved in the southern sourced waters during relatively cold climate conditions at this time (EPICA community members, 2006). Considering the age model uncertainties, this short deoxygenation event of the deep Bering Sea might further occur synchronously with the beginning of the early Antarctic warming interval (∼17 ka BP).
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 The appearance of a highly oxygenated environment followed the short dysoxic episode in the Bering Sea might be explained by an active intermediate-water formation in the North Pacific region during the severe H1 interval when newly formed [O 2 ]rich precursor of glacial NPIW reached the water depth of 2.5-3.0 km (Matsumoto et al., 2002;Okazaki et al., 2010). We suppose that the southern-sourced signal represented by lowoxygenated NPDW was suppressed by the Northern Hemisphere climate signal reflected by entrainment of the well-oxygenated NPIW precursor.
The second Antarctic warming and strengthening of the Southern Ocean upwelling from 12.5 to 11.5 ka BP corresponded to most of the YD and the onset of the EH is coeval with a decrease in oxygenation at both sites in the western Bering Sea. This coincidence supports the active transmission of southern-sourced water from the Southern Hemisphere to the North Pacific and its upwelling there to the intermediate depths.
As it was mentioned above, the upwelling of southern-sourced waters might compensate for a weak intermediate water formation in the Bering Sea during the YD compared to the H1. Along with the Southern Ocean, the Eastern Equatorial Pacific might also contribute to the deoxygenation of northward-propagating deep waters by enhancing sea surface bioproductivity during deglaciation (e.g., Bradtmiller et al., 2006;Ivanova et al., 2012;Hoogakker et al., 2018). In particular, the lowest oxygen concentrations were reconstructed for NPDW at 1.36 and 2.65 kmbsl in the Eastern Equatorial Pacific, suggesting an expansion of the eastern tropical OMZ during 17-15 ka BP (Hoogakker et al., 2018).
Thus, we argue that deep-water Pacific circulation might significantly control the oxygenation state of the deep western Bering Sea. Southern Hemisphere signal could also be detected at the intermediate depths, particularly during the YD-EH interval.

Deglacial Oxygenation Changes in the Open North Pacific and its Marginal Seas
A previously published global review of the available proxy records from the different parts of the World Ocean demonstrates a decrease in oxygenation of the ocean interior above 1.5 kmbsl and an increase in oxygen content below this depth between the LGM and EH (Jaccard and Galbraith, 2012). A decrease in oxygen solubility reportedly explains the global upper ocean deoxygenation as the sea-surface temperature rose. In contrast, growing deglacial oxygenation of the deep ocean occurred due to the respired carbon release from the ocean to the atmosphere even in case of ocean warming (Jaccard and Galbraith, 2012).
To compare deglacial oxygenation changes in the North Pacific area, including its marginal seas, we analyze our records in combination with previously published data from the study area within a wide depth range (Table 3, Figure 6). We selected those based on BF assemblages and uranium trace elements known to be sensitive indicators to oxygenation changes (Jaccard et al., 2014). We also consider bulk δ 15 N record from Core U1340 from the southern Bering Sea as a proxy of core site denitrification that occurs at low-oxygen environments (Schlung et al., 2013). All records, except for ODP882 (Jaccard et al., 2009) and ODP887 (Galbraith et al., 2007) collected from water depth below 2.3 kmbsl, demonstrate an oxygenated bottom-water environment during the LGM-H1 interval. High [O 2 ] content seems to be related to the sea-ice formation and associated brine rejection during this period rather than to changes in temperature-driven oxygen solubility in the Pacific interior. This statement is consistent with the hypothesis about the expansion of glacial NPIW in the North Pacific during the last glacial time (e.g., Matsumoto et al., 2002).
At the beginning of the B/A, all records, except for 85KL, display an intense deoxygenation event that was previously described and interpreted in many published works (e.g., Zheng et al., 2000;Crusius et al., 2004;Cook et al., 2005;Galbraith et al., 2007;Shibahara et al., 2007;Kim et al., 2011Kim et al., , 2017Jaccard et al., 2009;Jaccard and Galbraith, 2012;Schlung et al., 2013;Tetard et al., 2017). Several reasons might result in this abrupt change in oxygenation at this time. First, high organic matter flux to the sea floor (e.g., Crusius et al., 2004;Cook et al., 2005;Kim et al., 2011;Belanger et al., 2020) might impact the oxygenation state via an increase in biological respiration. Second, Northern Hemisphere climate warming most likely contributed to decreased oxygen solubility in the remote icefree areas (Jaccard and Galbraith, 2012;Jaccard et al., 2014). Third, shoaling of the global thermocline reportedly resulted in nutrient flux from the deep ocean (possibly by upwelling of NPDW) to the sea surface during the reinvigoration of Atlantic meridional overturning circulation and thereby maintained the high sea-surface bioproductivity (Schmittner, 2005). Finally, the closed Bering Strait prevented leakage of meltwaters from the North Pacific; this led to an enhancement of sea surface stratification and favored an increase in sea surface productivity and organic matter drawdown during the B/A (Davies-Walczak et al., 2014;Kim et al., 2017).
During the rest of the B/A, some shallower records collected from the depths of ∼0.7-1.0 kmbsl) display evidence for low oxygen bottom-water conditions probably due to high sea-surface productivity. On the contrary, the time-series obtained from the deeper cores taken within 1.0-2.3 kmbsl indicate an increase in oxygenation probably due to the shallowing of the central organic matter respiration depths. For the sites deeper than 2.3 kmbsl, highly oxygenated environments seemed to be related to a depletion of the global nutrient inventory (Jaccard and Galbraith, 2012) and supply of relatively young well-ventilated bottom waters (Galbraith et al., 2007;Du et al., 2018;Gorbarenko et al., 2020). The records obtained from core sites between 1 and 2.3 kmbsl demonstrate that the oxygen content was slightly higher during the major part of B/A than during its onset but lower than during the LGM-H1 interval at this depth range ( Figure 6). It seems that even relatively well-ventilated southern-sourced waters supplied from the Southern Ocean during the deglaciation were less oxygenated than intermediate waters formed in the northern latitudes during the cold intervals. Furthermore, the colder climate in Antarctica during the ACR (warm B/A in the North Atlantic time scale) favored more oxygen to be dissolved in the southern-sourced waters.
At the beginning of YD, a modest short-term increase in oxygenation was detected at the core sites with depths ∼0.7-2.3 kmbsl due to an intensification of intermediate water formation during cold periods. During the rest of the interval, oxygenation decreased again, probably due to an increase in southern-sourced water supply via the vigor Pacific abyssal circulation (Du et al., 2018) at times of reinforcement of the Southern Ocean upwelling (Siani et al., 2013). These southernsourced waters likely penetrated the shallow depth (∼0.7-1.0 kmbsl), explaining the oxygen minimum at the corresponding core sites during the YD.
During the rest of deglaciation, the deepest sites experienced a further increase in oxygenation due to continued reduction of oxygen solubility (Jaccard and Galbraith, 2012;Jaccard et al., 2014), whereas the upper ocean underwent oxygen deficiency due to high sea-surface bioproductivity conditions over the EH (Kim et al., 2011;Ovsepyan et al., 2013Ovsepyan et al., , 2017Bubenshchikova et al., 2015) Enhancement of deoxygenation due to highly productive sea surface conditions during the EH became less pronounced with greater depth, demonstrating a very modest maximum between 1.0 and 2.3 kmbsl. The latter might be partly explained by the opening of the Bering Strait at 11 ka BP (Jakobsson et al., 2017). This link between the Pacific and the Arctic Ocean allowed a meltwater leakage to the high latitudes, thereby weakening the sea surface stratification. Decreasing stratification, in turn, might reduce the sensitivity of the North Pacific environment to producitivity-related water column dysoxia (Davies-Walczak et al., 2014;Kim et al., 2017).
According to the analyzed proxy distributions, the most pronounced deoxygenation event occurred in the North Pacific at depths greater than 1.0 kmbsl during the onset of the B/A. This depth range of the oxygen-deficient environment is slightly wider than in the previously published review, which reported the lower boundary of the oxygen-deficient layer at ∼2.9 kmbsl during the B/A (Moffitt et al., 2015). Above 1.0 kmbsl in the marginal seas, the bottom waters deoxygenation during the EH was stronger than during the B/A at some sites (Bubenshchikova et al., 2015; this study) probably due to higher productivity conditions during the EH.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 638069 Ingram (1995) emphasized that any restriction in water mass circulation amplifies even small changes in oxygenation in such basins compared to the open water conditions. The sensitivity of enclosed and semi-enclosed basins to environmental changes is related to its physiography that increases the water residence time in the interior and prevents the dilution of nutrients and landoriginated freshwater, responsible for the development of high productivity conditions (Rabalais et al., 2010). Thus, slight variations in oxygenation at the given depth in the open North Pacific might be significantly amplified in its marginal seas.

CONCLUSION
Based on benthic foraminiferal assemblages and previously established transfer function (Tetard et al., 2017;Tetard et al., 2021a) The deglacial deep water [O 2 ] concentrations mostly follow the Antarctic air temperature oscillations and changes in the intensity of the Southern Ocean upwelling. The deglacial acceleration of Pacific abyssal circulation might be responsible for the active southern-sourced water supply to the North Pacific. The deep Pacific waters are suggested to be highly oxygenated compared to those during the glacials (Jaccard and Galbraith, 2012;Gorbarenko et al., 2020). However, in the western Bering Sea at the present-day water depths of 963 and 2,163 mbsl, they remained more depleted in oxygen than glacial waters. In addition, the southern-originated signal might be amplified by the Northern Hemisphere climate warming via an increase in sea-surface bioproductivity in the Bering Sea during the B/A and EH.
A comparison of previously published benthic foraminiferal data, redox-sensitive elemental records, and geochemical time-series from the North Pacific cores reveals three depth-related patterns of oxygenation histories reflecting the dominance of the North Atlantic climate signal at ∼0.7-1.0 kmbsl, the mixed Northern and Southern Hemisphere signals at 1.0-2.3 kmbsl and presumably the Southern Hemisphere signal at >2.3 kmbsl.
The majority of the available deglacial record confirmed previously published conclusions about the most intensive OMZ conditions at the beginning of the B/A. The strongest OMZ at the onset of the B/A might be explained by both influence of the North Atlantic and Antarctic climate signals, which induced the maximum of sea surface bioproductivity; availability of nutrients due to the shallowing of the thermocline, deoxygenation of the upper ocean during the air temperature rise, and strong supply of southern-sourced deep waters to the deep North Pacific. The closed Bering Strait might have further amplified the deoxygenation of the water column during the B/A.
The interregional correlation of oxygenation proxy records indicates that oxygen deficiency was stronger in the marginal seas (Bering and Okhotsk Seas) than the open North Pacific in the same intervals, probably due to increased water residence time and a higher concentration of nutrients delivered by rivers and meltwater to the marginal seas.
Our knowledge and understanding of past variations in ocean's oxygen concentrations and its relationship to natural climate variability are essential to infer and model the ocean's response to modern climate change. Millennial-scale quantitative oxygen estimates extent the existing observational datasets further in the past and, thus, might have the potential for filling gaps in predicting the future ocean OMZ behavior.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article is available in the Supplementary information.

AUTHOR CONTRIBUTIONS
EO designed the research and wrote the paper, EO, EI, and LM performed the analyses, EO and MT analyzed data, RT provided samples. All authors discussed the results and helped to improve the article at all stages.

ACKNOWLEDGMENTS
EO appreciates Christian März for fruitful discussions and valuable recommendations, Eric Galbraith for providing a published age model for Core ODP882 and Consuelo Martínez Fontaine for giving the revised age model for Core MD07-3088. We thank Lester Lembke-Jene for his support in AMS 14 C sample handling. The English editing was provided by N. Day (ServiceScape Inc.). We thank the reviewers NS and YO and the topic editor KN for thorough revision and constructive comments that helped to improve this paper. We also greatly appreciate YO for English editing of the revised version of the article.