The sea-ice proxy IP₂₅ is a highly branched isoprenoid (HBI) biomarker produced by certain species of marine Arctic sea-ice diatoms during the spring algal bloom^32,36,37. Its occurrence in marine sediments is now widely used as a reliable proxy for the presence of seasonal (spring) sea ice in paleo records from polar regions^32,38. We identified IP₂₅ in all samples studied (except for one sample at ca. 17.5 ka), consistent with a near-continuous presence of seasonal sea ice, albeit with variable duration or extent.

During late MIS 4 and early MIS 3 (64–45 ka) IP₂₅ values remain generally low, with IP₂₅ peaks occurring mostly during the stadials, except for the end of GI 14 (52–49 ka), one of the longest interstadials, when IP₂₅ increases. This is preceded by a peak in concentration of HBI III, a molecular biomarker produced by some pelagic diatoms in the open waters of the retreating sea-ice edge or in the marginal ice zone (MIZ)³⁹. This may be indicative of a progressive increase in seasonal sea-ice cover after the peak warmth at ca. 52 ka. Highest IP₂₅ concentrations occur during the mid-late MIS 3 (45–29 ka), during the Last Glacial Maximum (LGM; 24–19 ka), and the late glacial period (19–17.5 ka) (Fig. 2B).

The high IP₂₅ values during the mid-late MIS 3 may have resulted from a gradual retreat of an otherwise extensive sea-ice cover, thus favoring the growth of IP₂₅-producing diatoms, or from seasonal sea ice within the dynamic MIZ in the study area. For the same period, the concentrations of HBI III remain generally low (Fig. 2C). The combination of elevated IP_25, yet relatively low concentrations of HBI III, suggests the occurrence of a seasonal, but persistent, sea-ice layer in the study area. This interpretation is supported further by estimates of high SpSIC (and supported by the biomarker-based CT model) and low HBI T₂₅ indicating only intermittent spring phytoplankton blooms (Fig. 2D,E). The previously published results on the benthic foraminiferal faunas show dominance of the benthic species Stainforthia spp. and Nonionella spp., which feed on short-lived seasonal pulses of phytoplankton deposited on the sea floor also pointing to presence of an extensive seasonal sea-ice cover²⁴ (Fig. S1).

Similarly, the LGM and the late glacial are characterized by substantial shifts between high and low IP₂₅ concentrations, accompanied by low/absent HBI III (Fig. 2B,C), likely as a result of rapid fluctuations in the degree of seasonal sea-ice cover. In contrast to the mid-late MIS 3 period, the previously published benthic foraminiferal concentration data for the LGM show maximum concentrations of benthic foraminifera, and assemblages dominated by Cassidulina neoteretis²⁴, a species that thrives under the influence of cooled Atlantic water and is found in stratified waters with dense sea ice or short-lasting seasonally ice-free conditions (see Cage et al.⁴⁰ and references therein) (Fig. S1).

The results from core HH15-1252PC show that seasonal sea-ice (i.e., SpSIC) cover was more extensive than at present during the last glacial period. These long-term severe sea-ice conditions appears to have masked the expected millennial-scale signals in the record. However, a clear millennial-scale pattern can be seen in the distribution of IP₂₅ (seasonal sea ice) and HBI III (retreating sea-ice edge/MIZ) for some GSs and GIs (Figs. 2 and S1). The most clearly defined stadials are characterized by generally low IP₂₅ (with the exception of events between 30 and 40 ka; see below) and a very low, quasi-absent HBI III, indicating overall extensive sea-ice conditions during most of these periods (Figs. 2B,C and S1). This pattern of low IP₂₅ is particularly pronounced for HSs, and especially for HS 6, HS 5, HS 2 and HS 1, when sea-ice cover was likely so significantly advanced and dense that the production of both IP₂₅ and HBI III was very restricted. The concentration of IP₂₅ increases slightly towards the end of these HSs, as does HBI III, which reaches a maximum at the beginning of the subsequent GI. In addition, SpSIC declines towards the end of each of these HSs. Spring phytoplankton blooms (HBI T₂₅ > 1) occurred mainly during GIs, although not exclusively (Fig. 2D). This likely indicates a rapid increase in phytoplankton production due to the close vicinity of the ice margin/MIZ to the core site, and/or longer-lasting seasons of open water concomitant with the initial warming of GIs (Fig. 2C). However, the progressive opening of the ocean surface that started within the HSs (HS 6, HS 5, HS 2 and HS 1) and culminated at the beginning of GIs, suggests that the transition from (1) dense sea-ice cover to (2) a retreating sea-ice edge to (3) open-ocean conditions may not be caused exclusively by the abrupt atmospheric warming, recorded in Greenland ice cores at the beginning of the GIs^5,41 (Figs. 2 and 3).Figure 3

Stacks of Heinrich Stadial (HS) showing normalized bottom water temperature (BWT; red dashed line) and Spring sea-ice concentration (SpSIC; green line) in time. Stack including HS 6 to 1 (A), stack including HS 6, HS 5, HS 2, HS 1 (B), and the interpreted water column structure for three time slices represented in the stacks (pre-, syn- and post-transition from high to low BWT) (C). The thickness of the main water masses is not scaled and is exaggerated. The stacking was done setting time 0 at the mid-point of an abrupt BWT drop during HS (see Material and Methods). AW Atlantic Water; DW Deep Water; PwSI Polar water with Sea Ice; SI Sea Ice.

Between 30 and 40 ka, HS 4 and HS 3 exhibits higher production of IP₂₅ and absence of HBI III, indicating presence of a lasting sea-ice cover (high IP₂₅) with only infrequent or transient seasonal open-water conditions at the core site (very low/absence of HBI III). Indeed, conditions may have been similar to the LGM, as indicated by high (but variable) % of C. neoteretis (Fig. S1).

Bottom water temperature estimates based on Mg/Ca measurements in benthic foraminifera in core HH15-1252PC²⁴ show generally high BWTs (up to 5 °C) during HSs, consistent with studies from the southern Nordic Seas^23,42. The BWT and sea ice show coherent patterns during most HSs, with the exception of HS 4 and HS 3. We believe this to be due to the greater variability in the BWT data (Figs. 2 and S1) and/or uncertainties in the age model²⁴. To investigate the relationship between BWT and sea ice during the HS intervals further, we normalized and stacked our data for all HSs (Figs. 3a; S3A,C; see Material and Methods) and then focused our attention and discussion on HS 6, HS 5, HS 2 and HS 1, while excluding HS 4 and HS 3 (Fig. 3B).

For each of HS 6, HS 5, HS 2, and HS 1, decreases in SpSIC and a transition from extensive to intermediate/marginal sea-ice conditions (Fig. 2) occur alongside reductions in BWT (Figs. 2 and 3), demonstrating concomitant sea-ice retreat and cooling of the intermediate-depth ocean. It is worth noting that non-HS stadials do not show the clear coupling between BWT and SpSIC observed during HSs and have thus not been included in the stacking of data. During HSs, we interpret higher temperatures at deep and intermediate depths in the Fram Strait to be the result of the presence of a strong halocline causing surface stratification, cessation/reduction of deep water formation linked to major AMOC reductions²² and a larger southward extent of the sea ice cover in the Nordic Seas^16,18. This would have reduced the ocean–atmosphere heat exchange and, in turn, allowed warm Atlantic water to flow at greater depths^22–24. Consistent with these hypotheses, we observe extensive spring sea-ice cover (SpSIC > 50%) at the beginning of HSs during periods with high BWT (Figs. 2E,F and 3). Such conditions would have been unfavorable for both open-ocean and sea-ice related diatom flora, as reflected by low HBI III and IP₂₅, respectively (Fig. 2B,C). Although extensive sea-ice conditions dominated during the entire last glacial period, a decline in sea ice from extensive (SpSIC > 50%) to intermediate extent (SpSIC 50–10%) occurred concomitant with reductions in BWT during late HSs (Figs. 2E,F and 3). These results therefore support the presence of a sea-ice lid, which would have limited heat loss from the subsurface to the atmosphere, thus contributing to higher BWT at high latitudes during the early part of HSs^22–24.

During the major GIs that follow HS 6, HS 5, HS 2 and HS 1, the reactivation of deep water formation could have resulted in generally low BWT, similar to modern day values (ca.  − 0.5–0 °C)^23,24. In these periods, prominent peaks in HBI III generally occur at the beginning of each GI that follows a HS, are likely the result of a rapid retreat of the sea-ice edge coincident with initial warming (Fig. 2C). Open ocean and reduced sea-ice conditions (relatively higher HBI III) occurred during GIs, although extensive spring sea-ice persisted during the GIs of the cold mid-late MIS 3 and MIS 2 (Fig. 2E).

Core HH15-1252PC (79.04°N; 6.89°E; 1257 m water depth; 9.35 m core recovery) was retrieved north of Vestnesa Ridge, on the western Svalbard continental slope on a cruise with the R/V Helmer Hanssen in summer 2015. After retrieval, magnetic susceptibility was measured with a Bartington MS2 loop sensor. The core was subsequently split longitudinally in two halves. The archive halves were X-rayed with a GEOTEK 7.9 Multi Sensor Logger and color imaged with a Jai L-107CC 3 CCD RGB line scan camera installed on an Avaatech XRF at the laboratory of the Department of Geosciences at UiT The Arctic University of Norway in Tromsø. The work halves were subsampled for benthic foraminiferal Mg/Ca analyses to investigate the evolution of bottom water temperature (BWT), as well as benthic and planktic foraminiferal stable isotopes (δ¹⁸O and δ¹³C), faunal distribution of benthic foraminifera and ice-rafted debris counts (see²⁴ for details).

The age-depth model of core HH15-1252PC was constructed using planktic foraminiferal δ¹⁸O, supported by the magnetic susceptibility and the down-core distribution of benthic foraminifera Cassidulina neoteretis²⁴. Planktic foraminiferal δ¹⁸O from our core were correlated to the North Greenland Ice Project (NGRIP) ice core δ¹⁸O with the GIC05modelext timescale b2k (i.e. before 2000 AD)^5,60,61, assuming that massive meltwater events occurred during Greenland Stadials⁷. Thereafter, the age-depth model was built by linear interpolation of the identified tie-points and validated using seven AMS-¹⁴C dates (see²⁴ for details on the construction of the age model).

We used the previously identified tie-points²⁴ to build an age-depth model with quantified uncertainties. To account for uncertainties resulting from the spaced sampling intervals (average spacing between adjacent planktic foraminiferal δ¹⁸O samples of 3.32 cm) and the visual identification of tie-points, we used a fixed uncertainty around the tie-points of 150 years following the approach by Hughen and Heaton⁶². The uncertainties for each tie-point were then used to build the age-depth model using the Bayesian-based package BACON⁶³ in R Software.

Based on the new age model, core HH15-1252PC spans from 63.8 to 12.9 ka, covering Marine Isotopic Stages (MIS) 4 to MIS 2.

For biomarker analyses, a total of 201 samples were subsampled from the archive halves of the core in 1-cm thick slices in intervals of 1 to 5 cm (one sample subsampled at a 7-cm interval), with higher resolution sampling carried out for intervals with high BWT. Samples were freeze dried, homogenized, and stored in glass vials.

Lipid analysis was carried out according to Belt et al.⁶⁴, but with a slight modification to the extraction method. Thus, freeze-dried samples (ca. 5 g) were saponified in a methanolic Potassium Hydroxide (KOH) solution (ca. 10 mL Methanol:MilliQ water (9:1); 5% KOH) for 60 min (70 °C). Hexane (3 × 2 mL) was added to the cooled (room temperature) saponified content, with non-saponifiable lipids (NSLs) transferred to clean vials and dried over Nitrogen (N₂, 25 °C). NSLs were then further fractionated using silica (SiO₂; 0.5 g) column chromatography with non-polar fractions containing HBIs eluted with hexane (6 mL). Prior to extraction, samples were spiked with an internal standard (9-octylheptadec-8-ene, 9-OHD, 100 ng) to permit quantification.

Analysis of purified fractions containing HBIs was carried out using an Agilent 7890A GC coupled to a 5975 series mass selective detector (MSD) and operating conditions specified in Belt et al.⁶⁴. The identification of individual HBIs was based on their characteristic retention indices and mass spectra (Belt, 2018), while quantification was achieved by comparison of mass spectral responses of selected ions (m/z 350 (IP₂₅); 348 (HBI II, also referred to as IPSO₂₅³²); 346 (HBI III and IV)) with those of internal standard (9-OHD, m/z 350), and normalized according to their respected instrumental response factors and the mass of sediment extracted⁶⁴.

Semi-quantitative measures of spring sea-ice concentrations (SpSIC) and HBI T₂₅ index (a measure of phytoplankton blooms at the productive marginal ice zone) were derived according the Eqs. (1) and (2) respectively^33,65.1SpSIC%=IP25IP25+HBIIII×0.63-0.0692/0.0107\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$SpSIC\left( \% \right) = \left( {\left( {\frac{{IP_{25} }}{{\left( {\left( {IP_{25} + HBI III} \right) \times 0.63} \right)}}} \right) - 0.0692} \right)/0.0107$$\end{document}2T25=HBIIIIHBIIII+HBIIV0.62\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{25} = \frac{{\left( {\frac{HBI III}{{\left( {HBI III + HBI IV} \right)}}} \right)}}{0.62}$$\end{document}

Classification tree (CT) methods⁶⁶ were used to further categories spring sea-ice conditions into marginal (0–10%), intermediate (10–50%), and extensive (> 50%).

To measure the total organic carbon (TOC) ca. 3 g of sediment was subsampled every 5 cm. Dried and powdered samples were treated with 10% hydrochloric acid (HCl) for the removal of carbonates and subsequently measured in a Leco CS-200 at Department of Geosciences at UiT the Arctic University of Norway in Tromsø.

Due to differing sampling intervals in the biomarker analysis and the organic carbon analyses, 33 samples do not have paired carbon content measurements. However, the downcore trends and relative changes of lipid biomarker data normalized with the bulk sediment weight and with the total organic content do not show any deviation from each other (Fig. 2B) and therefore it was deemed unnecessary to carry out such analyses in the remaining samples.

To construct the stacks presented in Figs. 3 and S3 we used an approach that consist of the following steps. First, we estimated the timing of significant changes in BWT throughout the record by applying a Bayesian change point analysis method using the ‘bcp’ package in R⁶⁷. We detected changes in the mean values of the BWT timeseries that crossed a credibility level of 95%. Using visual inspection of the results, we then narrowed down the number of change points to those corresponding to Heinrich Stadials, which for our data set sums up to six transitions (i.e. HS 6 to HS 1; Fig. S3A–C). Next, we defined a time vector t that spans from -1500 to + 1500 years at 200-year steps. We set each of the six transitions defined in the BWT data to t = 0, and linearly interpolated the different records (i.e. BWT, HBI III, IP₂₅, SpSIC) onto the time vector t. The interpolation step was chosen based on the mean temporal resolution across all the proxy records considered for stacking. Finally, with all the individual events resampled to identical time spacing, we estimated the correlation (Figs. 3 and S3) and then scaled the data between  − 1 and 1, and averaged them to obtain stacked records. The resulting stacks, excluding HS 4 and HS 3 (Fig. S3D–F; see main text), were used to examine the phasing and co-variability between BWT and sea-ice proxies relative to the transitions observed in the BWT data.