Introduction

Record lows of Arctic sea ice have been repeatedly renewed in recent decades, and an ice-free summer is projected to occur within the coming four to five decades or even earlier1,2. Thermal control of the Arctic sea-ice decline has been linked to the enhanced heat transport from lower latitudes into the Arctic and its marginal seas, in association with background warming within the atmosphere-ocean coupled system3,4,5. Similar to the modern situation, paleoceanographic evidence has indicated a remarkable Arctic sea-ice retreat in the summer6,7 and extended melting seasons8,9 during the Early-Holocene episode. Such retreat of the paleo-Arctic sea ice is likely attributed to higher solar insolation10,11, increasing riverine heat flux into the Arctic Ocean12, temporal atmospheric circulation pattern resembling the modern negative-phase Arctic Oscillation13,14, and enhanced Atlantic/Pacific Water inflows6,15. Particularly, previous studies have respectively suggested that the Pacific Water influences millennial sea-ice variation in the Chukchi Sea6,7,16, and that Holocene poleward oceanic heat transport is regulated by the North Pacific ocean circulation system from the low-latitude towards the subarctic regions17, which is potentially coupled with the changes in the atmospheric circulations7,18,19, e.g., the Aleutian Low20. However, the role of the poleward oceanic heat transport via Bering Strait Inflow (BSI), known as a key thermal regulator in the modern Arctic Ocean21,22, in triggering the Holocene sea-ice retreat and its spatial heterogeneity over the Arctic marginal seas remains poorly understood.

Under modern climatic conditions, the BSI serves as a major ocean throughflow between the North Pacific Ocean and the Arctic Ocean22, while the North Atlantic Ocean and the Arctic Ocean are mainly connected via the Barents Sea23 and Fram Strait24. Within the Arctic Ocean and its marginal seas, the upper-layer waters sourced from the BSI and North Atlantic Ocean horizontally confront each other in the East Siberian Arctic Shelf (ESAS)25,26 (Fig. 1a), which thus is one of the most characteristic regions of seasonal sea-ice coverage and the crossroad of oceanological and climatological processes connecting with North Pacific Ocean and North Atlantic Ocean. Vertically, the BSI water stays in shallower depths and more directly affects sea-ice retreat in the ESAS27, because of relatively lower densities compared to the North Atlantic Water22,28,29,30. Geographically, the ESAS is a typical study area that clearly reflects the important impact of the BSI on the Arctic sea-ice cover. For example, the higher sea surface temperature in the Pacific sector and sea-ice retreat in the ESAS in September 2007 is much more severe than that in September 2008 (Fig. 1c, d).

Fig. 1: Oceanographic setting of the modern Arctic Ocean and Northern Hemisphere area.
figure 1

a, b Summer sea surface temperature (SST)80 in the Northern Hemisphere (NH) and ESAS, respectively. Position of the surface currents are shown as white arrows, including Kuroshio Extension (KE), North Pacific Current (NPC), Alaska Current (AC) and California Current (CC) in North Pacific Ocean17, as well as BSI, North Atlantic Water (NAW), Transpolar Drift (TPD), Beaufort Gyre (BG) and Siberian Coastal Current (SCC) in the circum-Arctic regions25. Yellow star marks the location of core LV77-36-1 (this study), black rhombs show locations of other cores with published paleoclimate and paleo-sea ice records used and discussed in this study (cores numbered as 1 to 22, Supplementary Table 1). c September 2007 and d September 2008 SST in the circum-Arctic regions from the ERA5 satellite data. Mean September sea-ice extent is taken from the National Centers for Environmental Prediction, representing the second-strongest observed sea-ice minimum in 200781. b represents the study area range (black rectangles in a, c, d).

Full size image

In this study, we reconstruct the Holocene sea-ice conditions in the ESAS, using a combination of sea ice biomarker namely ‘IP25’ (ice proxy with 25 carbons)31 and phytoplankton biomarkers (dinosterol and brassicasterol) as semi-quantitative sea ice proxy ‘PIP25’ (see “Methods” for procedure) for the sea-ice coverage in spring and summer32. Based on a published age model12, our sea ice record of core LV77-36-1 (155.66°E, 74.10°N, 36 m water depth, Fig. 1b) covers the past 8.3 ka BP with a high time-resolution of approximately 30 years. Together with a compilation of existing paleoclimatic and paleoceanographic evidence (Supplementary Table 1 and Supplementary Table 2), modern observation data and simulation results from multiple ocean general circulation models (Supplementary Table 3), we propose that an intensification of BSI with associated enhancement of poleward oceanic heat transport from the Pacific into the Arctic Ocean (named as Pacificization effect33) played a dominant role in driving the pronounced sea-ice retreat in the Pacific-side Arctic Ocean during the Late Holocene.

Results and discussion

Pronounced Pacific-side Arctic sea-ice melting during the Late Holocene

The IP25, a biomarker produced by sea ice algae during the ice microalgae growing seasons and deposited into sediments following the ice melting seasons (mainly from spring to summer)31, shows concentrations in core LV77-36-1 ranging from 0.22 to 4.63 μg g−1 OC−1 with an average value of 0.80 μg g−1 OC−1 (Fig. 2a). This varying magnitude aligns with Holocene IP25 records in other Arctic shelf areas6,34,35. Relatively stable IP25 concentrations with a lower average of 0.5 μg g−1 OC−1 were observed between 8.3 and 2.2 ka BP with cyclic variations (Fig. 2b), while the Neo-Glaciation interval (i.e., the last 2.2 ka BP) is characterized by about 4-times higher values. Phytoplankton biomarkers such as dinosterol and brassicasterol are predominantly produced during the spring blooms under reduced ice and ice-free conditions32. In our study, the record of dinosterol concentrations is similar to brassicasterol concentrations, with average values of 23.26 μg g−1 OC−1 and 16.43 μg g−1 OC−1, respectively (Fig. 2d, e). As both ice-free and permanent ice cover can result in lower IP25 concentrations (Supplementary Fig. 1), it’s necessary to combine IP25 with open-water phytoplankton biomarkers and calculate PIP25 values (see “Methods” for procedure) to semi-quantitatively reconstruct Arctic sea-ice conditions. Studies based on the surface sediments have proven an evidently positive correlation between sedimentary PIP25 values and modern satellite-derived spring and summer sea-ice concentrations32,36,37 (see “Supplementary Note 1”). In the present study, both brassicasterol-based PBIP25 and dinosterol-based PDIP25 exhibit a general increasing trend throughout the Holocene (Fig. 2f). Notably, our high-resolution PIP25 record highlights a temporal sea-ice retreat in the ESAS between 4.5 and 2.2 ka BP (PIP25 < 0.5), which contrasts with the long-term trend of sea-ice expansion throughout the Holocene (Fig. 2f).

Fig. 2: Holocene (8.3 to 0 ka BP) evolution of sea ice and phytoplankton biomarkers in core LV77-36-1.
figure 2

a IP25 concentration record. b Amplified IP25 concentration during the interval between 8.3 and 2.2 ka BP. c–e Records of total organic carbon (TOC) contents, brassicasterol and dinosterol concentrations, respectively. f Records of brassicasterol-based PBIP25 and dinosterol-based PDIP25 values. Red dashed lines show the boundaries between 4.5 and 2.2 ka BP. Green, cyan, blue, and gray shadows indicate ice-free, variable/less ice, seasonal ice and extended/permanent ice cover, respectively32,36.

Full size image

Here, the Holocene-scale sea-ice expansion was reflected both in our core and other PIP25 records from circum-Arctic marginal seas, including the Chukchi Sea6, East Siberian Sea6, Laptev Sea34, Kara Sea35, the seas around Fram Strait15 and the Beaufort Sea38. Previous studies have linked the long-term, circum-Arctic sea-ice expansion to the contemporaneous decrease in solar radiation6,10,39 (Fig. 3b) through the Holocene. In contrast, the intercalated sea-ice retreat between 4.5 and 2.2 ka BP appeared to present exclusively in the Pacific-side Arctic marginal seas, e.g., the Bering Sea40 (Fig. 3l), Chukchi Sea6,7 and East Siberian Sea (Fig. 3k, this study). While this occurrence was absent in the Atlantic-side Arctic, e.g., the Laptev Sea34 (Fig. 3d) and the Fram Strait area15 (Fig. 3c), where sea-ice dynamics may be dominantly forced by Atlantic climate15. Notably, Stein et al.6 documented a similarly dramatic sea-ice retreat in the core PS72/350 and core ABA2B/1B from the ESAS, which likely occurred earlier than record in our study6. However, these mismatches might be explained by the differences in age models and core locations, respectively. Core PS72/350 has a still very tentative age model with no AMS14C chronological fix points in the mid-late Holocene6 (Fig. 3n). Core ABA2B/1B, on the other hand, is located in the central Chukchi Sea much closer to the BSI pathway than core LV77-36-1, which possibly contributes to the more sensitive and earlier sea-ice loss triggered by the intensified BSI (Fig. 3m). By integrating existing published PIP25 records (Fig. 3), our results suggest a geographic uniqueness of the pronounced sea-ice retreat between 4.5 and 2.2 ka BP within the Pacific-side Arctic Ocean and its marginal seas, named as ‘Pacific-side Arctic sea-ice melting (PASIM)’ event.

Fig. 3: Proxy records for the Arctic sea-ice cover, North Pacific and North Atlantic climate changes during the Holocene.
figure 3

a Neodymium isotopic composition (εNd) in ODP Site 1063 (site 1 in Fig. 1) from the North Atlantic, as index of AMOC variation56. b Summer insolation at 75°N39. c, d PIP25 in core MSM712-2 (site 6 in Fig. 1) and core PS51/159 (site 9 in Fig. 1) representing Holocene sea-ice conditions in the Fram Strait15 and Laptev Sea34, respectively. e, f Records of diatom abundance in ice-related species (T. hyperborea) and warm-water species (P. sulcata) in the study core LV77-36-1 (site 10 in Fig. 1)41, respectively. g, h Records of chlorite plus muscovite in core ARA2B/1B (site 14 in Fig. 1)6 and ostracoda abundance in core 2-PC1 (site 12 in Fig. 1)43, respectively, both as indicators for BSI transport into the Chukchi sea. i, j Uk’37-SST records in core GC33 (site 19 in Fig. 1) in the Eastern Bering Sea48 and in core 85JC (site 20 in Fig. 1) along the Alaska Current path49, respectively. k–n PIP25 records in the Western East Siberian Sea (core LV77-36-1, this study), Bering Sea (core BR07, site 18 in Fig. 1)40, Chukchi Sea (core ARA2B/1B, site 14 in Fig. 1) and Eastern East Siberian Sea (core PS72/350-2, site 11 in Fig. 1)6, respectively. The paleo-records (a–d) were dominated by the North Atlantic climate trend, while the paleo-records (e–j) represent the North Pacific circulation changes. The orange and dark triangles represent the chronological fix points of study core LV77-36-112 and core ARA2B-1A6, while the blue triangles may act as specific age points of core PS72-350 based on the consistent PIP25 record with our PIP25 record in core LV77-36-1 (see text for further details). The dark and light pink shadows show the sea-ice loss during PASIM period (4.5–2.2 ka BP) and other periods.

Full size image

Enhanced Pacific oceanic heat triggers PASIM event via BSI

In the present study, we attribute the PASIM occurrence to a stronger BSI, which transported larger amount of warmer North Pacific Water into the Arctic Ocean. The associated oceanic heat poleward transported through the ocean circulation system and caused the pronounced sea-ice melting in the Pacific-side Arctic. This interpretation aligns with paleoceanographic evidence for intensified BSI during the PASIM period6,41,42,43 (Fig. 3e–j and Supplementary Fig. 3, see ‘Supplementary Note 2’). Thalassiosira hyperborea is a typical ice-related diatom species that abounds with sea-ice cover, whereas Pyxidicula sulcata is a eurytherm species preferring nutrient-rich and warm waters41. Therefore, the anti-phase variation in core LV77-36-1, marked by higher abundances of warm-water species and lower abundances of ice-related species between 4.5 and 2.2 ka BP41, is interpreted as enhanced Pacific-derived oceanic heat input triggering the subsequent sea-ice melting (Fig. 3e, f). A contemporaneously increased BSI is also reflected by lower abundance of cold-water ostracode species Acetabulastoma dunelmensis in the Chukchi Sea between around 5 and 2 ka BP43 (Fig. 3h), as well as by elevated chlorite contents in the Chukchi-Alaskan margin sediments6,42,44,45 (Fig. 3g and Supplementary Fig. 3b, c). The neodymium isotope composition (ɛNd) in the Chukchi Sea further supports this trend, indicating a progressive increase in Pacific Water fluxes during the Holocene46,47 (Supplementary Fig. 3d). All these proxy records varying during the period between 4.5 and 2.2 ka BP seem to collectively support a stronger BSI at that time and our PASIM hypothesis.

At a basin-wide scale of the North Pacific, a similar intensification during the PASIM period was observed not only in the BSI, but also in the Alaska Current48,49 (Fig. 3i, j), the Kuroshio Extension and the North Pacific Subtropical Gyre17,50,51 (Supplementary Fig. 3f–h), as indicated by alkenone-based sea surface temperature (SST) records. This synchronous enhancement of the upper-ocean circulation system may reveal a stronger oceanic heat transport from the lower-latitude Pacific to the Arctic Ocean. Particularly, the enhanced North Pacific Subtropical Gyre, indicated by the Uk’37-SST gradient between west Pacific warm pool and Kuroshio Current path17, may contribute to increased poleward heat transport that ultimately reinforced the BSI. In addition, the Aleutian Low pressure cell over the subarctic Pacific may have exhibited an eastward shift and/or enhancement during the PASIM period, as suggested by higher z-score values of diatom δ18O in the Arctic-Alaska region20 (Supplementary Fig. 3e), which could increase sea level pressure over the Aleutian Basin, raise dynamic sea surface height and thereby result in the BSI intensification42,52. This variation suggests stronger lower-latitude Pacific winds53,54 and warmer water inflows51 poleward transport through the intensified BSI, within the atmosphere-ocean coupled system20. Once entering the Arctic Ocean, the increased Pacific oceanic heat via the intensified BSI and related ocean circulations contributes to more sea-ice melting in the Pacific-side Arctic marginal seas (Fig. 3 and Supplementary Fig. 3).

In contrast to the enhanced Pacific-derived oceanic heat facilitating the PASIM event, the reduction in Atlantic Meridional Overturning Circulation (AMOC) and North Atlantic Water promoted sea-ice expansion in the Atlantic-side Arctic and circum-Arctic marginal seas15,34,38,55. Sedimentary εNd record in the North Atlantic suggests the persistent weakening of AMOC, with outstanding weakening during the PASIM period56 (Fig. 3a). The cooling scenario induced by the reduction in AMOC is also reflected by multiple paleoclimatic records, including ice-rafted debris in the North Atlantic sediments10, δ18O in ice cores (Supplementary Fig. 4b)57, diatom δ18O record in Kola Peninsula (Supplementary Fig. 4c)58, and pollen records near Lake Baikal (Supplementary Fig. 4d)59. In addition, proxy-based SST reconstructions along the North Atlantic Water pathway indicate the reduced Atlantic oceanic heat poleward transport between 4.5 and 2.2 ka BP (Supplementary Fig. 4e–g)60,61,62, which served as part of the Holocene cooling climate (Fig. 3a, b) but had limited effect on the PASIM event. At this point, our results suggest that the PASIM event was resulted from the warming effect due to temporally intensified Pacific oceanic heat and exceeded the effect of long-term cooling climate, which explains the geographic uniqueness of PASIM event in the Pacific-side Arctic Ocean.

Dipole variation in Holocene Arctic sea ice and its implication for the warming future

At the circum-Arctic scale, multiple independent paleo-sea ice reconstructions synchronously revealed a sea-ice loss in Pacific-side Arctic during the PASIM period, based on other paleoceanographic proxies including sea ice-related diatom species41, benthic foraminifera species typically associated with seasonal sea-ice cover43, ice-rafted debris12,16, and dinocysts27,63 (Supplementary Fig. 5b, see “Supplementary Note 3”). Moreover, according to the compilation of Holocene sea-ice reconstructions using biomarker records (PIP25 and IP25, Fig. 4a) and paleoceanographic evidence mentioned above (Supplementary Fig. 5 and Supplementary Table 2), we identify the spatial seesaw variation in Arctic sea-ice conditions during the PASIM period, characterized by sea-ice loss in the Pacific side (e.g., from the Bering sea to the East Siberian Sea) contrasting with gradual sea-ice expansion in the Atlantic side (e.g., from the Laptev Sea to the Icelandic Sea). Such pattern conceptually resembles the modern ‘Arctic Dipole’ mode of sea-ice thickness change, based on the sea level pressure anomalies between Atlantic and Pacific sectors64,65,66. Our findings suggest that a similarly competition of Pacific-Atlantic influences may trigger the dipole regime of circum-Arctic sea-ice variation over the past geological time-scale67.

Fig. 4: Maps of compilation of paleo-sea ice records, modern observation and modeling results.
figure 4

a The paleo-sea ice records both in the Atlantic side and Pacific side Arctic during the PASIM period between 4.5 and 2.2 ka BP, based on IP25 and PIP25 index values (see Supplementary Table 2 and Supplementary Fig. 5 for further details). Records of similarly reduced sea-ice and increased sea-ice cover are shown in red and blue symbols, respectively. b The correlation coefficient between BSI velocity and Arctic sea-ice concentration from 2008 to 2020. Red shadows represent negative correlation (mainly in the Pacific side), whereas blue shadows represent positive correlation (mainly in the Atlantic side). Black dots represent regional grid points that pass the significance test (p < 0.05). c Map of mean BSI Δvelocity between May and July under the strong La Niña conditions for the MRI-ESM2 model run. The black arrows indicate the inflows direction. The star marks the location of the study core LV77-36-1.

Full size image

To further assess this hypothesis, we compared the paleo-climate records with modern observation and numerical modeling dataset (Fig. 4). The mooring observation data indicates that intensification of BSI ranges from 0.7 to 1.2 Sv in recent decades22,68, which has a significantly negative correlation with satellite-derived sea-ice concentration in the Pacific-side Arctic Ocean (including the Chukchi Sea and the East Siberian Sea) during 2008–2020 (see ‘Modern observations’ in Methods, Fig. 4b). Similarly, numerical modeling resluts from five ocean general circulation models in the Ocean Model Intercomparison Project Phase 2 (Supplementary Table 3) show that intensified BSI transports warm-water into the Chukchi Sea and the East Siberian Sea, which has a vital effect on the sea-ice retreat between May and July under modern climate conditions of North Pacific Ocean (see “Reanalyses data of model running" in Methods, Fig. 4c and Supplementary Fig. 6). Therefore, modern satellite, physical oceanographic observation data and ocean modeling results consistently indicate the feasibility of our hypothesized physics for the PASIM event and the dipole pattern of sea-ice cover across the circum-Arctic Ocean.

Our high-resolution records of Arctic sea-ice conditions, combined with paleoclimate evidence from the Pacific Ocean and Atlantic Ocean, stress the PASIM event of anomalous sea-ice loss in the Pacific-side Arctic between 4.5 and 2.2 ka BP. It is triggered by intensified advection of oceanic heat from the lower-latitude Pacific Ocean, poleward transport via BSI and associated basin-scale upper-ocean circulation system. Our findings demonstrate that the temporary enhancement of Pacific-derived oceanic heat overcame the long-term Holocene cooling trend within the circum-Arctic climate system (Fig. 5). The PASIM event is further supported by a compilation of paleo-sea ice records and the hypothesized physics feasibility based on modern observations and multi-model results. These results provide crucial paleoceanographic insights into how North Pacific oceanic heat via the BSI drives pronounced Arctic sea-ice retreat, as partly also shown in earlier studies6,7. More importantly, our study highlights the latitudinal interaction between the North Pacific Ocean and the Arctic Ocean during the Holocene.

Fig. 5: Schematic illustration of Pacificization effect during the PASIM period.
figure 5

a During the PASIM event between 4.5 and 2.2 ka BP, the enhanced upper-ocean currents in the North Pacific Ocean resulted in increased oceanic heat transport from lower latitudes into the ESAS via the parallel stronger BSI. As such, this enhanced Pacificization effect with poleward oceanic heat contributed to the proposed PASIM event, whereas AMOC persistently weakening with less heat inflow into the Arctic Ocean (outstanding reduction during the PASIM event). b During the normal Holocene periods, the long-term trend of sea-ice expansion in the Arctic Ocean was mainly driven by decreased solar insolation in combination with the persistently deceleration in AMOC with reduced oceanic heat inflow from North Atlantic (relatively stronger than PASIM period), especially when North Pacific Subpolar Gyre circulations and BSI were also weakened, resulting in decreasing Pacificization effect and increasing sea-ice cover in the Pacific-side Arctic. Here, red arrows and orange arrows represent the North Pacific and North Atlantic inflows, respectively, and red dashed circles mark the areas under Pacificization effect.

Full size image

Conceptually, our illustration about the physics of intensified BSI in triggering the PASIM event is broadly in line with the known Pacificization process33, which has been defined as the dominated effect of Pacific Water on the Arctic rapid changes. Meanwhile, paleoceanographic evidence60,61,62 (Supplementary Fig. 4e–g) suggests the reduced Atlantic Water poleward transport, in combination with decreased solar radiation, leading to the cooling process and sea-ice expansion over the circum-Arctic region throughout the Holocene. This Holocene-scale pattern resembles the conception of Atlantification regime69,70, which is associated with anomalous advection of Atlantic Water and serves as the counterpart of Pacificization. Overall, our finding about the PASIM occurrence represents the millennial event when the warming impact via Pacificization overcame the cooling effect via Atlantification as noted above, suggesting the competition effect of “Pacificization-Atlantification” in shaping spatial heterogeneity of the Arctic sea-ice cover during the Holocene.

In particular, coupled modeling results mentioned above have projected enhanced BSI duing the strong La Niña years compared to El Niño mode in the modern North Pacific Ocean (Fig. 4c and Supplementary Fig. 6). Furthermore, studies have suggested more frequent occurrence of La Niña events71,72, accompanied with a stronger Kuroshio system due to global warming73,74. In the future, projected warming of the surface Pacific Ocean in the coming decades75 is expected to a larger poleward oceanic heat transport via BSI29. Here, our paleoclimatic findings have important implication for an intensified sea-ice melting in the Pacific-side Arctic Ocean and attributed to the enhanced Pacificization effect, on top of a basin-scale sea-ice decline due to global warming.

Methods

Sediment sampling and preparation

During cruise LV77 of R/V “Akademik Lavrentiev” in August and September 2016, the gravity core LV77-36-1 (74.10°N, 155.66°E) was collected in the ESAS, at a water depth of 36 m. This core site was located near the mean edge of sea ice between 1981 and 2010 (Fig. 1b). The length of the recovered sediment core is 376 cm. The chronology of this core, representing the last about 8.3 cal ka BP, was based on AMS14C ages obtained from shells and has been published by Dong et al.12. Downcore measurements of grain size12 and major and trace elements41 have previously confirmed that the sedimentary sequence is continuous without any disturbance and thus suited for paleoceanographic and paleoclimatic investigations. For this study, the upper 40 cm section was subsampled continuously at 1-cm intervals, while the lower section was subsampled at 2-cm intervals, resulting in a total of 209 samples for organic geochemical analysis with a high-resolution of ~30 years per sample. Sediments were stored at −20 °C until freeze-dried and homogenized for further treatment.

Analyses of IP25 and phytoplankton sterol biomarkers

Analysis of highly branched isoprenoids (HBIs, e.g., IP25 and HBI III) and phytoplankton sterol biomarkers (e.g., brassicasterol and dinosterol) were carried out at the Marine Organic Geochemistry Laboratory, Ocean University of China. The extraction and purification procedures followed previous described procedure76,77. Briefly, freeze-dried sediment samples (~5 g) were extracted 4 times using CH2Cl2/MeOH (2:1, v/v). Prior to extraction, the internal standards 7-hex-ylnonadecane and C19 n-alkanol were added to the sediments for quantification. The hydrocarbon fraction (containing HBIs) and alcohols (containing sterols) were separated via column chromatography (SiO2) using 8 ml n-hexane and 12 ml dichloromethane: methanol (95:5 v/v), respectively. Alcohols were saponified with 6 % KOH in methanol for 12 h, and further silylated with 40 μL bis-tri-methylsilyl-trifluoroacet-amide (BSTFA) and 40 μL CH2Cl2 at 70 °C for 1 h, befor sterol analysis.

HBIs was analyzed by gas chromatography (Agilent 7890 A GC; 30 m HP-5MS column, 0.25 mm i.d., 0.25 μm film thickness) coupled to an Agilent 5975 C VL mass selective detector (MSD, 70 eV constant ionization potential, ion source temperature 230 °C). The oven was kept initially at 60 °C for 3 min and then programmed to 150 °C at 15 °C per minute, followed by 10 °C per minute to 320 °C holding for 15 min. Sterols were analyzed by gas chromatography (Agilent 6890 N GC; 30 m HP-1 capillary column, 0.32 mm i.d., 0.17 μm film thickness). The oven was kept initially at 80 °C for 1 min and then programmed to 200 °C at 25 °C per minute, followed by 4 °C per minute to 250 °C, 1.7 °C per minute to 300 °C for 10 minutes, and finally 310 °C holding for 8 min. For the quantification of HBIs, the molecular ion of IP25 (m/z 350) and HBI III (m/z 346) in relation to the ion of internal standard 7-hexylnonadecane (m/z 266) was used. The content of sterols were calculated from the ratio of their GC peak integrations to that of the C19 n-alkanol IS. The biomarker concentrations are reported as μg g−1 of bulk dry weight sediment, and presented in the normalization by total organic carbon contents of the samples (i.e., μg g−1 OC−1), to compensate for different sedimentation rates in the study area. Finally, we gained 209 biomarker data of HBIs, dinosterol and brassicasterol concentrations.

The phytoplankton-IP25 (PIP25) index was defined by Müller et al.32,

$${{{{\rm{PIP}}}}}_{25}=\frac{{{{{\rm{IP}}}}}_{25}}{({{{{\rm{IP}}}}}_{25}+{{{\rm{phytoplankton}}}}* c)}$$

where c is the mean IP25 concentrations divided by the mean phytoplankton biomarker concentrations. Brassicasterol, dinosterol and HBI III were used as phytoplankton biomarkers to derive PBIP25, PDIP25 and PIIIIP25 indices, respectively (Supplementary Fig. 2). The higher PIP25 means more severe sea-ice cover. In this study, given that water salinity and turbidity changes caused by riverine input have an inhibitory effect on the growth of brassicasterol producers, we primarily used dionsterol concentrations to calculate this index (i.e, PDIP25)37,78 (see ‘Supplementary Note 1’ for details).

Analyses of total organic carbon

The analysis of total organic carbon contents were carried out at the Key Laboratory of Marine Geology and Metallogeny, First Institute of Oceanography, Ministry of Natural Resources, Qingdao, China. The detailed analytical methods of bulk organic matter proxies can be found in Hu et al.79. Briefly, portions of the freeze-dried sediment samples (~1 g) were decalcified using 2 M hydrochloric acid, and the carbonate-free samples were then analyzed for total organic carbon and total nitrogen in duplicate using a Vario EL-III Elemental Analyzer. In total, 209 samples were analyzed for total organic carbon content.

Modern observations

The summer Sea Surface Temperature (SST) dataset used in Fig. 1a, b was obtained from the World Ocean Atlas 2023 (WOA2023) in the National Centers for Environmental Information80 (NCEI), while September SST in 2007 and 2008 in Fig. 1c, d was obtained from the reanalyses data from ERA5 monthly averaged data on single levels produced by the European Centre for Medium-Range Weather Forecasts81 (ECMWF). We used the satellite observation data of sea-ice concentrations (SIC) in the circum-Arctic region from OSI SAF Global SIC climate data record (v3.0, 2022)82, at a global 1° × 1° spatial resolution. The BSI variation was observed in A3 mooring site (168°57′N, 66°19.6′W) of 2 m surface water83, with continuously monthly velocity data from only A.D 2008 to 2020. Therefore, in this study, we used the monthly data from January 2008 to December 2020 (i.e., 12 data per year) both in BSI velocity and SIC variability to ensure sufficient data points for significance test. Then we can deseasonally calculate the correlation coefficient between BSI and SIC and focus on their close relationship in the whole year.

Reanalyses data of model running

Here, our reanalyses of the modeling results were based on the Ocean Model Intercomparison Project Phase 284, in the Coupled Model Intercomparison Project Phase 6 (CMIP6). We firstly intermpolated the modeling results of MRI-ESM2, FESOM, NorESM2-LM, MIROC6, and FGOALS, whose resolutions and the amount of grid points representing Bering Strait are presented in the Supplementary Table 3. Then the strong La Niña and El Niño conditions years (Supplementary Table 4) were indentified and grouped according to Nino 3.4 index. In the end, the climatology mean of the La Niña and El Niño conditions were respectively calculated by (sqrt(u2 + v2)) of the surface inflows.