The North Pacific is characterized by two large-scale wind-driven gyre systems: the cyclonic Subarctic Gyre in the north and the anticyclonic Subtropical Gyre in the south (Huang, 2015). The Subarctic Gyre has a strong western boundary current. The Oyashio Current originates from the East Kamchatka Current and is further fed by cold, low saline nutrient-rich water from the Okhotsk Sea. Close to the coast of northern Japan the Oyashio Current has an annual mean temperature of ~5°C and a salinity of ~33 (Locarnini et al., 2018; Zweng et al., 2019). At ~42°N the Oyashio Current splits and one path continues as the northeastward flowing Subarctic Current (SAC) (Qiu, 2019). The western limb of the North Pacific Subtropical Gyre is also formed by a strong western boundary current. The Kuroshio Current originates from the warm, saline North Equatorial Current (NEC) and transports oceanic heat from the WPWP to the northern regions (Figure 1). The location where the North Equatorial Current (NEC) bifurcates into the Kuroshio Current and the Mindanao Current (Metzger and Hurlburt, 1996; Qiu and Lukas, 1996; Ujiié et al., 2003) is located at ~15.5°N for the annual average, but varies with time and depth. The meridional migration of the NEC bifurcation is strongly influenced by ENSO (Kim et al., 2004), which is the main driver of upper ocean circulation in the tropical Pacific, (Hu et al., 2015; Joh and Di Lorenzo, 2019). In El Niño years the NEC’s bifurcation point migrates to the north and the Kuroshio velocity and volume transport decrease. During La Niña years, in contrast, the NEC’s bifurcation migrates to the south and the Kuroshio velocity and volume transport increase (Qiu and Lukas, 1996; Yuan et al., 2001; Kim et al., 2004). Close to the coast of southern Japan, the Kuroshio Current has an annual mean temperature of ~19°C and a salinity of ~35 (Locarnini et al., 2018; Zweng et al., 2019). After separating from the Japanese coast at ~35°N, the Kuroshio Current enters the open basin of the North Pacific and continues as Kuroshio Extension (Qiu, 2002). At Shatsky Rise (~159°E), the Kuroshio Extension bifurcates, where the main body of the Kuroshio Extension continues eastward. A secondary branch, the Kuroshio Bifurcation Front (KBF), extends northwards (Qiu, 2002).

Along the northern boundary of the Kuroshio Extension/Kuroshio Bifurcation Front, numerous meanders and mesoscale eddies occur (Figure 1) due to the dia- and isopycnal mixing between Kuroshio and Oyashio water with vertically different temperature, salinity and velocity structures (Yasuda et al., 1996; Shimizu et al., 2001; Yasuda et al., 2002; Isoguchi et al., 2006; Mitnik et al., 2020). This mixing zone of water originating from the Kuroshio and Oyashio currents is termed the Kuroshio-Oyashio transition zone or mixed water region and is characterized by a pronounced latitudinal sea surface temperature, salinity, and nutrient gradient (Yasuda, 2003; Garcia et al., 2018; Locarnini et al., 2018; Zweng et al., 2019) (Figure 1). The northern end of the transitional zone is marked by the Subarctic Front (SAF) or Oyashio Front, which is defined as the 4°C isotherm at 100 m depth (Favorite, 1976). Subarctic surface water may drift southward across the SAF, and may cover the surface in the transition zone between the SAF and the Subarctic Boundary (SAB). The SAB is a near surface salinity front south of the SAF (Yasuda, 2003). East of 180°W/E the distinction between the Kuroshio Extension and the Subarctic Current is no longer clear and together they form the broad eastward-moving North Pacific Current (Qiu, 2019). It has been suggested, that there is a linkage between the latitudinal displacement of the Kuroshio-Oyashio transition zone and the climatic conditions of the tropical Pacific e.g. long term ENSO-like variability (Yamamoto et al., 2005).

Intermediate water masses (~400-900 m) are dominated by North Pacific Intermediate Water (NPIW) in the Subtropical Gyre and the mixed water region (Talley, 1993), and by Pacific Subpolar Intermediate Water (PSIW) in the subarctic Gyre (Emery, 2001; Fuhr et al., 2021).

To determine the sedimentary elemental composition, the archive half of each core segment was scanned using an Avaatech X-ray Fluorescence (XRF) core scanner at the Alfred‐Wegener‐Institute (AWI) Bremerhaven. For preparation, the uppermost sediment layer was removed to provide a fresh and undisturbed surface, which was then covered with SPEXCerti Prep Ultralene foil. Each core segment was triple scanned at a one-centimeter resolution with a current of 0.15, 0.175 and 1 mA, tube voltages of 10, 30, and 50 kV, and acquisition times of 10, 15, and 20 sec, respectively. Element intensities were determined with the proprietory Aavatech WinAxil (Batch) software covering elements from aluminum to barium and are reported as counts per second (cps). As proxy for changes in paleoproductivity we use the barium/titanium ratio (Ba/Ti) (biogenic barium reflecting diatom productivity; e.g. Nürnberg et al., 1997). Based on the assumption that the covariance between aluminum (Al) and Ti content of terrigenous material remained constant in space and time (Galbraith et al., 2007), the XRF-derived barium counts normalized against terrigenous background elements like Al or Ti are well-established paleoproductivity proxies in the North Pacific area (Nürnberg et al., 2004; Jaccard et al., 2010; Korff et al., 2016). Al and Ti are both exclusively of detrital origin and correlate well at our study site (Supplementary Figure 6). Yet the detection of light elements like Al during XRF scanning shows a higher risk to be influenced by water trapped between the sediment and the scanning foil than heavier elements (Tjallingii et al., 2007). Thus, we use Ti for normalizing our Ba record.

Fe has been widely used to document variations in terrigenous sediment delivery and often relates to aeolian dust flux and thus wind strength (Jahn et al., 2003; Nürnberg et al., 2004; Jahn et al., 2005; Mohtadi et al., 2007; Helmke et al., 2008; Kaiser et al., 2008; Lamy et al., 2014; Abell et al., 2021). Therefore, we use the iron (Fe) record derived via XRF scanning as proxy for relative changes in terrigenous Fe input. Because the use of raw intensities or peak integrals of single element counts cannot be interpreted in terms of sediment composition we report our Fe-record as log-values (Weltje et al., 2008). The very good correlation of our Ti and Fe records (R²=0.98) (Supplementary Figure 6) implies that a potential diagenetic imprint on Fe is negligible (Croudace and Rothwell, 2015).

For geochemical studies, we use the most abundant planktic foraminiferal species foraminifera Globigerina bulloides (G. bulloides) and Neogloboquadrina pachyderma (N. pachyderma) (sinistral-coiling, Darling et al., 2006). In the Northwest Pacific, G. bulloides is a near-surface dweller and mainly found above the thermocline (Iwasaki et al., 2017; Schiebel et al., 2017; Taylor et al., 2018), which occurs at ~40-50 m at the study site (Locarnini et al., 2018) (Figure 1 and Supplementary Figure 1). In the study region, G. bulloides calcifies throughout the year (Kuroyanagi et al., 2008; Sagawa et al., 2013; Taylor et al., 2018). In the North Pacific, N. pachyderma is mainly found at thermocline depth (Sarnthein et al., 2004; Riethdorf, 2013; Taylor et al., 2018), thus, we assume a depth habitat of 50-130 m for our study site (Locarnini et al., 2018). N. pachyderma is a subpolar species and prefers water temperatures below 7°C (Reynolds and Thunell, 1986), thus in modern settings it mainly occurs from autumn to spring at our study site, which is shown by sediment trap studies from the North Pacific (Kuroyanagi et al., 2008; Sagawa et al., 2013; Taylor et al., 2018). For a detailed discussion about seasonality and habitat depth, see Supplement 1.

For converting Mg/Ca of foraminiferal calcite into water temperature we used well-established species-specific equations from Elderfield and Ganssen (2000) for G. bulloides (Mg/Ca=0.56e^(0.10*T)) and from Kozdon et al. (2009) for N. pachyderma (Mg/Ca=0.13*T+0.35) (Supplement 2). The temperature derived from G. bulloides is referred to as annual sea surface temperature (SST_Mg/Ca: ~20-60 m) and the temperature from N. pachyderma as subsurface temperature (subSST_Mg/Ca: ~50-130 m). Late Holocene temperatures of both species match the modern annual mean temperatures in the defined depth habitats (Locarnini et al., 2018; Supplementary Figure 1). For our interpretations, we assume that the foraminiferal species did not significantly change their habitat depth over time.

We are aware that changes in pH and salinity of the surrounding seawater, calcite dissolution and diagenetic effects can influence foraminiferal Mg/Ca, and hence the reconstructed absolute temperature values. Yet, we assume that these effects are either neglectable for our interpretation or lie within the assigned error range. A detailed discussion can be found in the Supplement 2–4.

Via the combined δ¹⁸O and Mg/Ca measurement of the foraminiferal calcite, we calculate the δ¹⁸O of seawater (δ¹⁸O_sw) by removing the temperature effect from the initial foraminiferal δ¹⁸O signal. For G. bulloides we use the equation from Shackleton (1974) (T=16.9-4.38*(δ¹⁸O_calcite-δ¹⁸O_sw)+0.1*(δ¹⁸O_calcite-δ¹⁸O_sw)²) and for N. pachyderma the species-specific equation from Mulitza et al. (2003) (T=3.55*(δ¹⁸O_calcite-δ¹⁸O_sw)+12.69). The δ¹⁸O_sw data is further corrected for changes in global ice volume (δ¹⁸O_sw-ivc reported in ‰ versus SMOW; ivc = ice-volume corrected) by using a data set from De Boer et al. (2014). This was calculated based on the global LR04 stack (Lisiecki and Raymo, 2005) and simulations of continental ice sheets (De Boer et al., 2014). Modern δ¹⁸O_sw is positively correlated to salinity (e.g. δ¹⁸O_sw = 0.44*S-15.13; LeGrande and Schmidt, 2006), yet previous studies have shown that this relationship can vary over time e.g. through changes in regional freshwater budgets, ocean circulation, and sea ice regimes (Caley and Roche, 2015; Holloway et al., 2016). Therefore, we refrain from converting the calculated δ¹⁸O_sw-ivc into absolute salinity units. However, we interpret high values in terms of relative high saline conditions and low values in terms of relative freshening. We present δ¹⁸O_sw-ivc as relative changes deviating from the calculated modern value.

Accelerated mass spectrometer radiocarbon (AMS¹⁴C) dating provides an age estimation for the sediment surface (0–1 cm) from station SO264-45-1 (MUC). We chose the undisturbed surface sample of the MUC for this measurement, because the surface of gravity cores is often disturbed during core recovery. Approximately 90 µg C of the planktic species G. bulloides and the benthic species Uvigerina peregrina were measured at the MICADAS (Mini Carbon Dating System) facility at the AWI in Bremerhaven, Germany (AWI sample nr. 3976.1.1 and 3976.2.1). The raw ¹⁴C age was converted into calendar age, using the Calib8.20 software (Stuiver et al., 2020) along with the marine calibration dataset MARINE20 (Heaton et al., 2020) and a regional reservoir age correction of ΔR = 273 +- 70 years. In this case, ΔR is the mean reservoir age of six sites from the North Pacific published by Kuzmin et al. (2001; 2007) and Yoneda et al. (2007). The generated AMS¹⁴C age range (2 sigma) for U. peregrina is 263-745 years BP and 1871-1301 years BP for G. bulloides. This age difference lies within the age range for one sample (~1500 years/sample) as discussed below.

As benthic foraminiferal δ¹⁸O values from the MUC SO264-45-1 fit rather well to those of the uppermost centimeters of the gravity core SO264-45-2 (Figure 2), it is likely that both have similar surface ages. Thus, based on the ¹⁴C ages from the MUC’s surface we assume that the gravity core top sediments are of late Holocene age.

The stacked benthic δ¹⁸O isotope record of SO264-45-2 was visually correlated to the global benthic isotope stack LR04 (Lisiecki and Raymo, 2005) by using QAnalySeries v1.5.0 (Kotov et al., 2018). Our benthic record shows smaller amplitudes than LR04, which is most likely caused by low sedimentation rates, which mute extreme values. By applying 21 tie lines and linear interpolation between the tie lines, a correlation coefficient of r = 0.82 was achieved. Further support for the age model is provided by the spectral analysis of the benthic δ¹⁸O stack that reveals dominant cyclicities with a frequency of 0.01 and 0.025 as a response to the cyclic fluctuations in the Earth’s orbital parameters eccentricity and obliquity (Figure 2). The spectral analysis reveals that all obliquity cycles of the past 650 ka are reflected in our benthic δ¹⁸O stack. The calculated sedimentation rates are rather low with an average of 0.66 cm/ka for the discussed interval, thus one sample of 1 cm width reflects on average ~1500 years. Yet due to the observed bioturbation of the sediment (Nürnberg, 2018), the absolute age range of foraminifera in one sample might even be higher.

The ages of five tephra layers mentioned in the core description (Nürnberg, 2018) are shown in Figure 2. We assume that T4 might be the prominent Pauzhetka ash found in several records from the Northwest Pacific and Okhotsk Sea and dated to 421.2 ± 6.6 ka BP (Ponomareva et al., 2018). Yet according to our age model, the depth of T4 refers to an age of ~390 ka BP, which is likely within the range of dating uncertainty.

Over the last ~650 ka BP, δ¹⁸O values of G. bulloides fluctuate between 1.5 and 4 ‰, and those of N. pachyderma between 2.4 and 3.8 ‰ (Figure 3). With only few exceptions in MIS 6-8 the δ¹⁸O of N. pachyderma are consistently heavier than those of G. bulloides. The δ¹⁸O fluctuations of both species follow the glacial-interglacial pattern of the LR04 reference record (Lisiecki and Raymo, 2005) with heavy δ¹⁸O values during glacials and light δ¹⁸O values during interglacials (Figure 3). In MIS 11, however, G. bulloides shows lighter δ¹⁸O values than in the adjacent interglacials. In MIS 7, the δ¹⁸O are not light throughout the entire interglacial, yet decrease steadily. δ¹⁸O of G. bulloides is lightest in MIS 15 and 13 compared to the other interglacials in our record, while the δ¹⁸O of LR04 is heavier in MIS 15 and 13 compared to the other interglacials.

Mg/Ca ratios from G. bulloides vary between 0.82 and 2 mmol/mol, which leads to SST_Mg/Ca estimates between 3.8 and 12.7°C. The SST_Mg/Ca development of G. bulloides can be divided into two phases (Figures 3, 4). Between ~600 and 280 ka BP (phase B), the surface ocean temporally reached SST_Mg/Ca of up to ~12°C (with a mean of 8.7°C), thereby being up to ~6.5°C warmer than modern annual mean SST at 30 m (Figure 3). Within this time interval, the SST_Mg/Ca broadly follow glacial-interglacial cycles except for MIS 10, during which the SST_Mg/Ca is high (~10°C) throughout the glacial and MIS 11, where temperatures are ~2°C lower than in the adjacent glacials. In contrast to the cold MIS 16 (min ~ 7°C), MIS 14 remains rather warm with minimum temperatures of ~9.5°C. At 480 ka BP, the SST_Mg/Ca abruptly decreases by ~4°C (calculated as the temperature difference between three data points prior and after the decrease; the decrease itself occurs between two data points) coming close to the modern-day SST value (Figure 3). This decrease is followed by an overall increasing SST_Mg/Ca trend for the next 200 kyrs. Phase B ends at ~ 280 ka BP after a second abrupt SST_Mg/Ca decrease of ~3.5°C (calculated as the temperature difference between three data points prior and after the decrease; the decrease itself occurs between two data points). From here on (phase A) the SST_Mg/Ca remain rather low with an average of ~ 6.1°C and show only minor variations (max ±2°C) around the modern annual mean temperature at 30 m of ~5.5°C (Locarnini et al., 2018). The warmest SST_Mg/Ca occur at the beginning and end of MIS 5 and in MIS 2, while in the middle of MIS 7 and 6 the SST_Mg/Ca are coldest.

N. pachyderma shows consistently lower Mg/Ca ratios than G. bulloides, which vary between 0.41 and 1.14 mmol/mol with subSST_Mg/Ca estimates between 0.5 and 6.1°C (Figure 3). Throughout the entire record, the subSST_Mg/Ca show amplitude variations of ±3°C around an average of ~2.9°C. This average matches the modern annual temperature at thermocline depth of ~3°C (Locarnini et al., 2018). The highest subSST_Mg/Ca of up to ~6°C is reached at the transition to the Holocene and the lowest with ~0.5°C during MIS 12. In contrast to the SST_Mg/Ca the subSST_Mg/Ca do not show a shift towards warmer temperatures prior to 280 ka BP, yet rather cold subsurface conditions between MIS 13 and MIS 10. From MIS 5 to MIS 1 subSST_Mg/Ca and SST_Mg/Ca follow a similar low-amplitude pattern.

The δ¹⁸O_sw-ivc records reveal a different pattern than the SST_Mg/Ca and subSST_Mg/Ca records (Figure 3). Both, the seasurface and subsurface δ¹⁸O_sw-ivc records have rather similar amplitudes and fluctuations. Between ~650 and 350 ka BP G. bulloides shows constantly higher δ¹⁸O_sw-ivc values than N. pachyderma. Between 0 and ~350 ka BP, however, both records lie within error range of each other and alternately show higher values. Prior to 280 ka BP the δ¹⁸O_sw-ivc records of both species tend to be higher in interglacials and lower in glacials.

The XRF-based Ba/Ti record and thus the reflection of marine productivity shows a clear glacial-interglacial cyclicity. Ba/Ti values start to increase during the deglaciations and reach their maxima during the early interglacials (Figure 4). During the late interglacials and the subsequent early glacials low Ba/Ti values predominate. This pattern even holds for the period MIS 15-14, although the Ba/Ti ratios remain on an exceptionally high level during glacial MIS 14, which even exceeds typical deglacial to early interglacial values. The deglaciation towards MIS 13 is again characterized by a continuous increase in Ba/Ti. In the late interglacial MIS 13 (~500-480 ka BP), we observe a rapid and significant decrease in Ba/Ti to typical glacial values, synchronous to significant changes in Fe and SST_Mg/Ca (Figure 4).

The XRF-based normalized log Fe record, which is used as a proxy for terrigenous input commonly increases at the transitions from late interglacials to early glacials and decreases during the glacials. Prior to MIS 12, no distinct glacial-interglacial pattern is recognizable. Between ~560 and ~500 ka BP, Fe values are exceptionally low, which is followed by a sudden increase from ~500-480 ka BP (Figure 4). Ba/Ti and Fe appear to be anticorrelated with low Fe values in times of high Ba/Ti and vice versa. Thus, between MIS 13 and 5 Ba/Ti maxima occur on average ~30 ka earlier than Fe maxima.

At the study site, glacial-interglacial changes are well documented in the benthic and planktic δ¹⁸O records of all species with high δ¹⁸O values during glacials and low δ¹⁸O values during interglacials. This indicates that climate variations did have an effect on both the upper and deep ocean. Yet in our SST_Mg/Ca and subSST_Mg/Ca records changes between glacial and interglacial periods are neither particularly dominant nor continuously reflected. Thus, the observed temperature variations and especially the mode shift between phases A and B, recorded in the SST_Mg/Ca record must have been caused by different forcing mechanisms. A shift from higher SST_Mg/Ca prior to 280 ka BP to colder SST_Mg/Ca in more recent times, as reflected in our proxy data, does occur in few records from certain areas in the North Pacific (cf. 4.1.1) but it is neither a global phenomenon nor significant in the entire Pacific region (Morley and Heusser, 1997; Herbert et al., 2001; Lang and Wolff, 2011; Bordiga et al., 2013). Thus, our observed SST_Mg/Ca shift cannot be explained by a general cooling of the (North) Pacific climate after 280 ka BP yet must be linked to changes in local hydrography, most likely the complex interplay of the Kuroshio and Oyashio Current.

Higher marine productivity during deglaciations and interglacials and lower productivity during glacials as observed in our Ba/Ti record is common in the North Pacific and its marginal seas (Narita et al., 2002; Gorbarenko et al., 2004; Kienast et al., 2004; Nürnberg and Tiedemann, 2004; Jaccard et al., 2005; Brunelle et al., 2007; Shigemitsu et al., 2007; Galbraith et al., 2008; Gebhardt et al., 2008; Jaccard et al., 2010; Riethdorf et al., 2013). In our record we further observe, that the productivity already starts to increase during the deglaciations and decreases during the late interglacials. This is similar to observations at ODP Site 882 (Jaccard et al., 2010) located north of our site (Figure 1) which also shows productivity increases during deglaciations and its highest values in the early interglacials (Figure 4). A similar feature is observed for site SO202-39-3 (Korff et al., 2016) located south of our study site in the area of the Kuroshio Extension (Figures 1, 4). The forcing mechanisms for these productivity patterns, however, are still a matter of debate (Jaccard et al., 2010; Knudson and Ravelo, 2015; Korff et al., 2016). For the North Pacific, most authors suggest a change of nutrient supply, and in the subarctic marginal seas light limitation through sea ice cover and changes in stratification as main drivers for productivity changes (Narita et al., 2002; Gorbarenko et al., 2004; Kienast et al., 2004; Jaccard et al., 2005; Brunelle et al., 2007; Shigemitsu et al., 2007; Galbraith et al., 2008; Gebhardt et al., 2008; Jaccard et al., 2010; Riethdorf et al., 2013; Davis et al., 2020). As our site is located south of the area where sea ice would have a direct influence on productivity and the average winter SST exceeds 3°C, we exclude it as a forcing mechanism (Supplementary Figure 1). Therefore we conclude that the availability of nutrients or their utilization causes the observed glacial-interglacial productivity pattern.

Iron fertilization has been invoked as an important driver for nutrient sequestration efficiency and thus, an enhanced supply of the micronutrient Fe could increase productivity (Boyd et al., 2004; Harrison et al., 2004; Tsuda et al., 2003). One of the sources of Fe is aeolian transport (Hovan et al., 1991; Kawahata et al., 2000; Boyd et al., 2007; Shigemitsu et al., 2007). To test whether the Fe input via dust caused our observed productivity pattern, we compared our XRF-Fe record of core SO264-45, used as a proxy for terrigenous input via wind to our Ba/Ti record (Figure 4). The proxy records are apparently anticorrelated suggesting that Fe fertilization via dust is not the (principal) driver of marine productivity at our study site, although, we cannot say to what extent it still contributed to the observed pattern.

The amount of nutrient transport via water masses could have contributed to changes in productivity (Jaccard et al., 2005; Galbraith et al., 2007; Nishioka et al., 2011; Costa et al., 2018; Gray et al., 2018). Today, nutrient-rich water from the Oyashio/Subarctic Current dominates the upper ocean at our study site. The transported nutrients reach the mixed layer through vertical wintertime mixing (Nishioka et al., 2011). Our SST_Mg/Ca record, which is linked to changes in the Kuroshio-Oyashio transition zone does not show the same fluctuations as our productivity record. Thus, although it could have contributed to changes in productivity we assume that the interplay of the Kuroshio and Oyashio Current is not the main driver for the observed productivity pattern at our study site. We do, however, consider that the amount of transported nutrients might have changed on glacial-interglacial timescales. Lembke-Jene et al. (2017) show that the export of nutrients from the Okhotsk Sea increased at the beginning of the deglacial warm phases Allerød and Preboreal, because of an enhanced input of iron and nutrient-rich terrestrial material from the Siberian hinterland via the Amur River caused by melting processes. They argue that such an increased export of nutrients from the Okhotsk Sea could have caused temporary nutrient-replete conditions in the Subarctic North Pacific. Riethdorf et al. (2013) also observe an enhanced input of terrestrial-derived organic matter from flooded shelf areas during early deglacial phases. As water from the Bering and Okhotsk Sea feeds the Oyashio/Subarctic Current, these nutrients could have reached the study site and caused an increase in productivity. Assuming that an enhanced terrestrial sourced nutrient supply from marginal seas like the Okhotsk Sea and the Bering Sea via subsurface water did not only occur on millennial- but also on glacial-interglacial timescales, this could contribute to our observed productivity pattern. Yet to further clarify on this additional proxy records on longer time scales and in the open North Pacific would be needed.

Further it has been proposed for subarctic and Antarctic sites that strong stratification during glacials limits nutrient availability at the surface while weak stratification in interglacials enables upwelling of nutrients-rich deep water thereby enhancing productivity (Jaccard et al., 2005; Brunelle et al., 2007; Galbraith et al., 2008). Yet this topic alongside the efficiency of nutrient uptake is still a matter of debate (e.g. Knudson and Ravelo, 2015) for which further studies with detailed nutrient records and very reliable control would be needed.