Mid-Holocene sea-ice dynamics and climate in the northeastern Weddell Sea inferred from an Antarctic snow petrel stomach oil deposit

Understanding past variability in Antarctic sea ice is of critical importance to determine how it regulates global climate processes, biogeochemistry, and Southern Ocean marine ecosystems. Records of changes in Holocene sea-ice conditions in the Weddell Sea are limited to a few marine sediment cores and inferences from continental ice cores. Here we present a novel record of sea ice and climate from ∼ 6390 to 1830 cal. yr BP based on accumulation rates and multi-proxy geochemical analyses of a snow petrel stomach-oil deposit from the Heimefrontfjella Range, Dronning Maud Land, East Antarctica. Three different sea-ice configurations are interpreted from the record. In the first interval, from 6390 to 5960 cal. yr BP, we see evidence of high productivity and inputs of krill, which suggests foraging both at the continental shelf edge in the marginal ice zone (MIZ) and extending offshore over pelagic waters. We infer that the MIZ lay within the foraging range of Heimefrontfjella. In the second interval, from 5960 to 4320 cal. yr BP, productivity remained high, but there was a reduced influence of krill and likely more fish in the diet. This is consistent with foraging both over the continental shelf edge and offshore, supportive of the summer sea ice retreat reaching the shelf edge more frequently. Finally, in the final interval, between 4320 to 1830 cal. yr BP we infer very low productivity where increased sea-ice conditions restricted access to foraging grounds and open waters, with a less accessible MIZ resulting in a more dense sea-ice pack. Our results highlight how specific Holocene sea-ice configurations can be interpreted from the geochemical composition of snow petrel stomach-oil deposits. We also show, for the first time, the utility of phytol and cholesterol analysis for tracking past avian diet.

1 Introduction

Antarctic sea-ice conditions are highly variable and closely coupled with continental, oceanic and atmospheric processes which both interact with and influence global climate (Brandon et al., 2010). Mechanistically, the formation of sea ice results in brine rejection which can directly contribute to the formation of Antarctic Bottom Water (AABW) (Crosta et al., 2022), helping drive ocean circulation, including the deep overturning cells (Ferrari et al., 2014) and supporting important circulation systems such as the Weddell Gyre (Vernet et al., 2019). Within the sea ice, upwelling can result in the formation of open ocean polynyas (areas of open water), whereas near the coast or at ice shelf fronts polynyas can form as a result of katabatic winds (Comiso and Gordon, 1987). These polynyas can support high levels of primary productivity in the ocean (Smith et al., 2010; Sarmiento et al., 2004; Arrigo et al., 2003).

Antarctic sea-ice records reveal pronounced declines in extent since 2016, associated with recent warming (Eayrs et al., 2021), with historic lows in 2023 and 2024 (Ionita, 2024; Purich and Doddridge, 2023; Gilbert and Holmes, 2024; National Snow and Ice Data Center (NSIDC), 2024; Wang et al., 2024) and projections of up to 67 % decline by 2100 (Collins et al., 2013). Improved palaeoenvironmental reconstructions of Antarctic sea ice are vital to put these instrumental observations (< 50 years) into a longer-term context. Such reconstructions provide a historical basis to understand the interactions between climate and sea-ice conditions, and interactions between sea ice and the extent of floating ice shelves and grounded ice.

Existing knowledge of Holocene Antarctic sea-ice evolution suggests there were three distinct phases, but these can be out of phase geographically due to regional responses and uncertainty associated with dating (Crosta et al., 2022). For coastal records, a cooler early Holocene between 11.5 to 8 ka BP (e.g. Barbara et al., 2010; Denis et al., 2010; Peck et al., 2015; Mezgec et al., 2017; Nichols et al., 2019), was followed by a warmer mid-Holocene (∼ 7 to ∼ 4–3 ka BP) with higher sea surface temperatures and longer ice-free summers (Crosta et al., 2022) and then a cooler late Holocene or “neoglacial” phase ∼ 5–3 to 1–0 ka BP marked by increased sea-ice extent as surface water temperatures reduced (Barbara et al., 2016; Taylor et al., 2001). However, records retrieved from beyond the coastal regions or those influenced by CDW (Circumpolar Deep Water) show an opposite pattern, with a warmer late Holocene (Bianchi and Gersonde, 2004; Nielsen et al., 2004; Etourneau et al., 2013; Vorrath et al., 2023), potentially in relation to the latitudinal and insolation and thermal gradients specifically associated with wind stress and upwelling (Denis et al., 2010). The integration of multiple records from the Atlantic sector of the Southern Ocean suggest there was some late Holocene cooling, driven by enhanced cold-water export from the Weddell Gyre as a cavity developed under the Ronne Filchner Ice Shelf, combined with a northward shift of the Southern Hemisphere westerly wind belt (Xiao et al., 2016). However, the changes within the Weddell Sea remain poorly understood and there is a lack of data to reconstruct past sea-ice evolution (Verleyen et al., 2011), particularly over the continental shelf.

To address this, we analysed a snow petrel (Pagodroma nivea) stomach-oil deposit, accumulated at a nesting site in the Heimefrontfjella region of Dronning Maud Land, northeastern Weddell Sea. Such deposits have been previously demonstrated to record palaeoenvironmental information (McClymont et al., 2022; Berg et al., 2019; Ainley et al., 2006) by tracking the biochemical signature of changes in snow petrel diet, which is in turn related to environmental characteristics of their feeding grounds including sea-ice extent and ocean productivity. Snow petrels have a close affinity with sea-ice during the breeding season (Delord et al., 2016; Ainley et al., 2006), as they feed in areas of intermediate sea-ice cover and switch between neritic (close to shore) and pelagic (offshore) feeding grounds (Ainley et al., 1998; Ainley et al., 1984).

Modern end-member studies have shown that specific prey species (e.g. krill, fish and squid) can be separated based on their fatty acid composition and can be used to infer different sea-ice conditions in the foraging area. For example, Antarctic (Euphausia superba) krill are high in C_14:0, C_16:0 and C_18:1 fatty acids (e.g. Cripps et al., 1999), whereas squid are dominated by C_16:0 and fish by C_18:0 fatty acids (Lewis, 1966); reviewed in McClymont et al. (2022). However, some fish can also be high in C_18:1, especially ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (e.g. the myctophid Electrona carlsbergi) (Connan et al., 2008). This broad framework allows us to distinguish between an Antarctic krill-rich diet, reflecting an offshore (pelagic) habitat, and a fish-rich diet, reflecting near-shore and continental shelf environments, based on the modern spatial distributions of fish (Ran et al., 2022; Liu et al., 2024; Freer et al., 2019), Antarctic krill (Mcbride et al., 2021) and cephalopod (mainly squid) (Xavier et al., 2016) species. Stable isotopes in stomach oils are also potentially excellent indicators of ecological baseline nutrient availability, trophic status (δ¹⁵N) and the productivity and location of past foraging habitats (δ¹³C) (McClymont et al., 2022; Ainley et al., 2006).

Using a radiocarbon-based age-depth model, we analyse a range of organic geochemical biomarkers (including the first analyses of phytol and cholesterol in stomach oil deposits), bulk elemental chemistry and stable isotopes in the stomach oil deposit from Heimefrontfjella. These record three distinctive periods of Holocene climate and sea-ice cover, which we compare with existing records of environmental changes in Antarctica and the Southern Ocean.

2 Materials and methods

2.1 Regional context

The Antarctic seasonal sea ice cycle is characterised by growth to maximum extent in September, followed by retreat to a minimum in February. Between the minimum and maximum sea-ice extents lies the Marginal Ice Zone (MIZ), characterised by sea-ice concentrations between 15 % and 80 %. In the northeastern Weddell Sea, the median sea-ice edge (15 % sea-ice concentration) retreats rapidly with spring melt. In October/November, when the snow petrels start to return to their nesting sites, the sea-ice edge is located north of 60° S. Rapid sea-ice retreat through December brings the sea-ice edge to ∼ 70° S by early January, spanning the snow petrel incubation phase, then retreats further westward until the minimum in February (Fig. 1). The Weddell Sea gyre advects sea ice clockwise across the basin where it outflows in the northwest (Gupta et al., 2025). Sea ice tends to accumulate to the greatest extent in the south of the Weddell Sea due to this gyre (Hutchings et al., 2012).

Figure 1 Location of Holocene snow petrel stomach oil record 3012MUM2, from Heimefrontfjella (74°34.14^′ S, 11°15.02^′ W; yellow circle) in the context of the Weddell Sea, and Antarctic Peninsula (left panel) and Antarctica and the Southern Ocean (right panel). Dashed black lines indicate sectors within 600, 900 and 1100 km snow petrel foraging range from the stomach oil deposit. Also shown are median modern sea ice extents in February and September (Fetterer et al., 2017), RAISED grounded ice sheet extents at 5 and 10 ka (Bentley et al., 2014) and inferred modern Antarctic grounding lines (Rignot et al., 2022). Other locations and core sites mentioned in the text include: ice cores from EPICA DML (Masson-Delmotte et al., 2011); Dome Fuji (Masson-Delmotte et al., 2011); EPICA Dome C (Masson-Delmotte et al., 2011); Vostok (Masson-Delmotte et al., 2011); Taylor Dome (Masson-Delmotte et al., 2011) and James Ross Island (JRI) (Mulvaney et al., 2012); marine sediment cores U1357 (Ashley et al., 2021); TTN057-13PC4 (Hodell et al., 2001); TN057-13PC4 (Divine et al., 2010; Nielsen et al., 2007); TN057-17PC1 (Divine et al., 2010; Nielsen et al., 2007); ANT28/D1–7 (Nie et al., 2022); NBP0003-JPC38 (Barbara et al., 2013); NBP0502-Site 2, Herbert Sound (Totten et al., 2015); NBP0602A and NPP0703-JPC02, Firth of Tay (Michalchuk et al., 2009; Majewski and Anderson, 2009); SS01 (Bak et al., 2007); and JKC36 and JPC36, Perseverance Drift (Kyrmanidou et al., 2018). Map developed using Quantarctica (https://www.npolar.no/quantarctica, last access: 15 October 2025).

The preferred foraging habitat of snow petrels has been observed to track the MIZ, where the presence of open waters within the sea ice facilitates access to prey (Wakefield et al., 2025). This foraging behaviour can see snow petrels foraging over and beyond the continental shelf, reaching ∼ 700–1400 km from the nest site depending on the time of the breeding cycle and sea-ice extent (Honan et al., 2025). Based on regional tracking studies to the east, in Dronning Maud Land (Honan et al., 2025; Wakefield et al., 2025) three foraging options are available for snow petrels at Heimefrontfjella: (1) they head north to the MIZ lying to the south of the outer ice edge; (2) they head north east to the Maud Rise, where an intermittent open-ocean polynya and earlier sea ice melt provides foraging habitat; (3) they forage in the MIZ closer to the continental shelf, which can be associated with coastal polynyas. Tracking studies show that the most likely scenario is that early in the breeding season snow petrels forage in association with the northern edge of the MIZ as the sea ice retreats (Wakefield et al., 2025). As the sea-ice edge retreats so does the foraging range, so that later in the season snow petrels are more likely to be found in coastal waters (Honan et al., 2025; Wakefield et al., 2025).

Table 1 Radiocarbon dates including raw ¹⁴C, F¹⁴C and Bayesian calibrated ages for the 3012MUM2 sequence. Samples analysed by BETA Analytic (Florida) ¹⁴C-AMS were prepared without acidification. Repeat analyses at 5 and 10.5 cm were prepared at SUERC Radiocarbon Laboratory using acidification and run for ¹⁴C-AMS at UC Irvine Keck Carbon Cycle laboratory. A reservoir offset based on closest Holocene ΔR of 670 ± 50 years (Björck et al., 1991) was used, following previous studies (McClymont et al., 2022) and converted to calendar ages using the MARINE20 calibration (Heaton et al., 2020). All dates (¹⁴C and calibrated) rounded to nearest 10 calendar years for presentation, original measured dates are available at Pangaea (https://doi.org/10.1594/PANGAEA.980519, Stevenson et al., 2025). The Bayesian model output is shown in Fig. 2.

^∗ Bayesian P-Sequence model applies single age to paired date due to R_Combine function in OxCal.

2.2 Heimefrontfjella stomach oil deposit 3012MUM2

Stomach-oil deposit 3012MUM2 was collected in season 2014–2015 from the Boysennuten nunatak in the Heimefrontfjella Range (74°34.14^′ S, 11°15.02^′ W) (Fig. 1). The ∼ 32 cm × 23 cm deposit with a maximum thickness of 19 cm had an irregular, mamillated outer surface (McClymont et al., 2022). It was located immediately beneath a sheltered rock crevice, a typical habitat for nesting snow petrels, at an elevation of 1336 m above sea level (Fig. S1a and b in the Supplement). It was kept in the dark and frozen at −20 °C throughout the transportation processes to Durham University where sampling was carried out from 2021. The deposit was sliced using a circular saw while still frozen to preserve the internal millimetre-scale laminae which, when oriented, spanned a depth of ∼ 19 cm along the cutting axis (Fig. S1c). Sub sampling was carried out at 2.5 mm resolution with 3.0 mm biopsy punches for organic geochemistry, isotopes and radiocarbon analyses. We selected a sampling line which represented the maximum accumulation rate and a continuous sequence through the stratigraphy; noting some heterogeneity either side of this line (Figs. 2c and S1c). A 19 cm long slice from the opposite face of the cut was mounted in plastic trunking for high-resolution X-ray fluorescence (XRF) analysis of elemental composition.

2.3 Radiocarbon analyses and age-depth modelling

An age-depth model for 3012MUM2 was constructed from twelve ¹⁴C ages (Table 1). The top and bottom ages were sampled at 0 and 18.5 cm, immediately above and below the first geochemical samples (0.5 and 18.25 cm, respectively). Most of the radiocarbon ages were carried out by Beta Analytic (Miami, USA) using ¹⁴C-AMS via graphitization on untreated samples, which were oxidised to CO₂ by combustion in O₂ and converted to graphite with Co powder as a catalyst. To assess for the effects of acid removal on ¹⁴C ages samples at 5.0 and 10.5 cm were repeat-sampled and analysed at SUERC (Scottish Universities Environmental Research Centre) Environmental Radiocarbon Laboratory by digestion in 1 M HCL (hydrochloric acid) (at 80 °C, 2 h), washed free of mineral acid with deionised water, dried and homogenised. Carbon was recovered from the residue as CO₂ by heating with CuO in a sealed tube and converted to graphite by Fe/Zn reduction. Samples were then analysed by AMS at the Keck Carbon Cycle AMS Facility, University of California, Irvine, USA. For calibration of radiocarbon ages to calendar ages (Table 1) the MARINE20 radiocarbon age calibration (Heaton et al., 2020) was used, taking into account the regional marine reservoir of ΔR = 670 ± 50 years, measured at Hope Bay in the western Weddell Sea (Björck et al., 1991). This approach has previously been applied to snow petrel stomach-oil deposits (McClymont et al., 2022). A Bayesian age-depth model was then built in OxCal v4.4 (Bronk Ramsey, 2009) using the default settings (applied with a general outlier model, except for paired dates where we used the SSimple outlier model). Our choice of a Bayesian approach to age-depth modelling means that the age uncertainties of all dates are considered for the entire age model.

2.4 Bulk stable isotope analysis and organic matter elemental composition

Carbon and nitrogen stable isotope analysis was performed using a Costech ECS400 elemental analyser coupled to a Thermo Scientific Delta V Advantage isotope ratio mass spectrometer in the Stable Isotope Biogeochemistry Laboratory (SIBL) at Durham University. The method used is described in McClymont et al. (2022). Results are reported in standard delta (δ) notation in per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB) and atmospheric nitrogen (AIR), respectively. The linear range for δ¹³C was between −46 ‰ and +3 ‰ and for δ¹⁵N between −4.5 ‰ and +20.4 ‰, based on daily analysis of international (e.g. IAEA-600, IAEA-CH-3, IAEA-CH-6, IAEA-N-1, IAEA-N-2, NBS 19, USGS24, USGS40) and in-house standards, enabling a 2-standard-deviation analytical uncertainty of ± 0.1 ‰ for international standards (replicated) and < ± 0.2 ‰ on replicated samples. Total carbon (wt % C) and nitrogen (wt % N) were obtained simultaneously using an internal standard of glutamic acid (40.8 wt % C; 9.5 wt % N).

2.5 Biomarker analyses

The biomarker sub-samples (0.02–0.05 g) were extracted in 4 mL dichloromethane (DCM) : hexane (3:1) after addition of internal standards (nonadecane, heptadecanoic acid, 5α-androstane, 5α-androstanol) and then sonicated for 15 min. Extracts were decanted and the procedure was repeated three more times. Extracts were combined and taken to dryness using rotary evaporation and N₂. The entire sample was then saponified using 1 mL KOH (8 %) in methanol (95 %), heated for 1 h at 70 °C and left overnight. The neutral fraction was extracted with 3 × 3 mL hexane. The remaining sample was acidified to pH < 3 using drops of 2 M HCL, followed by extraction of fatty acids with hexane and evaporation to dryness with N₂. Fatty acid methyl esters (FAMEs) were generated by methylating the fatty acid fraction using 3 mL methanol:HCL (95:5) heated for 12 h and left to cool. Samples were rinsed with 4 M DCM-cleaned H₂O and then FAMEs were extracted with at least three DCM:hexane rinses (4:1), before evaporation to dryness under N₂. Neutral fractions were separated in up to four fractions using 4 cm deactivated silica (heated 140 °C for 16 h) columns in glass pipettes plugged with extracted cotton wool (silica pore size 60 Å, 220–240 mesh particle size; 35–75 µm particle size; Sigma-Aldrich 60738-1KG). Hexane was used to condition the columns (x3 column volumes) followed by injection of the sample dissolved in 500 µL DCM directly into the column. Elution order into separate fractions used three column volumes each of hexane, DCM, DCM:methanol (1:1) and methanol. All fractions were decanted and transferred into GC vials and evaporated to dryness using N₂. Fraction F3 (DCM:methanol (1:1), which recovered n-alcohols, phytol and sterols, was then further derivatized to trimethylsilyl esters prior to analysis using 50 µL DCM and 50 µL BSTFA (N,O-Bis(trimethylsilyl)trifluoroacetamide) (with 1 % TMCS (chlorotrimethylsilane)) heated for 1 h at 70 °C and left overnight prior to analysis. Samples were evaporated to dryness and dissolved in hexane prior to analysis.

All extracts were analysed using a Thermo Trace 1310 gas chromatograph coupled to an ISQ LT single quadrupole mass spectrometer (GC-MS). FAMEs extracts were separated using a Restek FAMEWAX (crossbond polyethylene glycol) column (30 m × 0.25 mm × 0.25 µm), similar to McClymont et al. (2022) but with some modifications. Briefly, samples were injected (0.8 µL) into a programmable temperature vapouring (PTV) injector in CT Splitless mode with inlet temperature at 250 °C, carrier gas in constant flow and with helium carrier gas set to 1.5 mL min⁻¹ (split flow 15.0 mL min⁻¹; splitless time 1.5 min; purge flow 5.0 mL min⁻¹). GC oven temperature was set to 100 °C for 3 min followed by 2 °C min⁻¹ to 230 °C; hold of 12 min. Prep-run timeout was 10 min and equilibration time 0.5 min. MS settings included: transfer line, 230 °C; ion source temperature, 230 °C, mass range 38 to 600 $m/z$ (every 0.5 s). Compounds were identified from their respective mass spectra and retention times, with quantities calculated relative to the peak area of the internal standard heptadecanoic acid and an assumption of a 1:1 response (validated by comparison with Supelco 37 component FAME mix; CRM47885, Merck).

Fraction F3 (previously eluted in DCM:methanol (1:1)) was separated using a Restek Rxi-5ms (crossbond 5 % diphenyl$/$95 % dimethyl polysiloxane) column (60 m × 0.25 mm × 0.25 µm). Similarly, samples were injected (0.8 µL) into a programmable temperature vapouring (PTV) injector in CT split-less mode but at 280 °C inlet temperature, with helium carrier gas set to 2.3 mL min⁻¹ set in constant flow mode (split flow set to 23 mL min⁻¹; splitless time 1 min; septum purge flow 5.0 mL min⁻¹). GC oven temperature was set to 50 °C hold for 2 min followed by 10 °C min⁻¹ to 200 °C; followed by a slower ramp of 3 °C min⁻¹ to 300 °C and a hold of 20 min (prep run timeout was 10 min, equilibration time 0.5 min). MS settings included: transfer line, 310 °C; ion volume, 300 °C, mass range 50 to 550 $m/z$ (every 0.5 s). Compounds were identified from their respective mass spectra and retention times, with quantities of trimethylsilyl esters (TMS) calculated relative to the peak area of the internal standard 5α-androstanol and an assumption of a 1:1 response (validated by identically derivatized standard mix which included cholesterol-TMS).

Figure 2 (a) Bayesian ¹⁴C-AMS age-depth model for stomach oil deposit 3012MUM2, generated in OxCal (Bronk Ramsey, 2009) applying a reservoir offset (ΔR^Holocene 670 ± 50 years (Björck et al., 1991) calibrated to calendar ages using MARINE20 (Heaton et al., 2020). Dark blue ellipse indicates 68.3 % range, light blue ellipse indicates 95.4 % range. (b) Accumulation rate between age control points calculated from median ages in the Bayesian model in (a). (c) Photograph of stomach oil deposit 3012MUM2 after sectioning, complete with approximate location of ¹⁴C ages in blue dots. Constrained hierarchical cluster analysis was used to determine three significant clusters (Org A–C) using organic parameters in rioja (Juggins, 2020) and the broken-stick model (Bennett, 1996).

2.6 XRF analyses

XRF (X-ray fluorescence) analysis was carried out at Durham University, Department of Geography using a GEOTEK XRF Core Workstation (MSCL-XYZ) equipped with a rhodium source X-ray generator with a 10 mm cross-core slit width and a 1 mm downcore window. During operation the XRF scanner was set to analyse four different beam conditions, with a counting time of 10 s per beam. Beam conditions applied to the generator included: (1) no filter, 10 kV; (2) 25 µm silver filter, 20 kV; (3) 125 µm silver filter, 30 kV; (4) 625 µm copper filter, 50 kV. Detector measurement of photons ranged 2–35 keV. To pre-screen complex data, including removal of missing values and selecting elemental compositions from the most appropriate beam, the R 3.6.0 package “tidyverse” (Wickham et al., 2019) was used to produce a master XRF dataset.

2.7 ICP-OES analysis

To determine local bedrock chemistry, a sample of gneissic granite rock (consistent with local geology; Juckes, 1972) attached to deposit 3012MUM2 was soaked in ethanol to remove stomach-oil residue (repeated 3–4 times). It was then crushed using a fly press, freeze dried for 48 h, and ground to a fine, homogenous powder using a ball mill. Organic matter was removed by adding 4 mL of 30 % hydrogen peroxide to a ∼ 0.5 g aliquot of rock. The sample was then digested for 4 h in 16 mL of Aqua Regia using a DigiPREP digestion block, and subsequently diluted to 50 mL with deionised water and filtered at 0.45 µm. Elemental composition was determined using an Agilent Technologies 5100 Inductively Coupled Plasma Optical Emissions Spectrometer (ICP-OES).

2.8 Statistical and numerical analyses

Cluster analysis was carried out to highlight changes in geochemistry between neighbouring units in both inorganic (XRF) and organic (stable isotopes, TOC, $\mathrm{C}/\mathrm{N}$ ratio and biomarkers) parameters with depth. For each dataset (organic (Org A–C) and inorganic (XRF A–C), separately) a constrained hierarchical cluster analysis based on sample order was performed using the rioja package in R (Juggins, 2020). A broken stick analysis (Bennett, 1996) was used to identify the maximum number of statistically significant clusters.

Principal components analysis (PCA) was carried out in Canoco V.5.51 (Ter Braak and Smilauer, 2002) on log ₁₀-transformed and centred data: the inorganic (XRF) and organic (bulk organic geochemistry and biomarkers) parameters were treated separately. Input data for the organic PCA included bulk organics and biomarker concentrations, rather than fluxes (which included TOC normalized data and ratios). Input data for the inorganic PCA included XRF units in counts per second, rather than fluxes. As most XRF parameters had samples with counts < 500 we chose to retain all parameters that had passed the pre-screening process (see Sect. 2.5 XRF analyses). The lists of elements and compounds used in inorganic and organic PCA are shown in Tables S1 and S2 in the Supplement.

3 Results

3.1 Stomach-oil deposit 3012MUM2 age model

The stomach-oil deposit 3012MUM2 spans 1260 (1060–1420) cal. yr BP at 0 cm to 6410 (6250–6610) cal. yr BP at 18.5 cm (Table 1, Fig. 2a). When taking biomarker and isotope samples we avoided the margins of the deposit, which were easily deformable. As a result, the oldest and youngest biomarker and isotope samples lie at 0.5 cm (1830 (1170–2530) cal. yr BP) and 18.25 cm (6390 (6210–6590) cal. yr BP), respectively in the Bayesian model, or ∼ 1800 to 6400 cal. yr BP when rounded to closest 100 years. The accumulation rate based on the median Bayesian modelled age between radiocarbon dating depths varied between 8.6 and 125.0 mm kyr⁻¹ (Fig. 2b).

Figure 3 Organic parameters as fluxes and ratios measured in stomach oil deposit 3012MUM2, plotted against the age-depth model (Fig. 2). Smooth lines in bold are 3-point moving averages. Constrained hierarchical cluster analysis was used to determine three significant clusters (Org A-C) using organic parameters in rioja (Juggins, 2020), compared with the broken-stick model (Bennett, 1996). (a) Fatty acid concentrations (C_14:0, C_18:0, ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$); (b) nitrogen stable isotopes (δ¹⁵N); (c) Other biomarkers (phytol and cholesterol); (d) Fatty acid ratios (${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$, ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$, ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$); e) other measures of bulk organic matter composition (C:N_atomic ratio, δ¹³C, % C and % N).

3.2 Compositional changes in bulk organic matter and stable isotopes

The 3012MUM2 samples were high in organic C (36.6 %–68.3 %; mean 49.8 %) and total N (2.9 %–16.8 % (mean 9.7 %) (Fig. 3). The C:N_atomic ratio varied between 3.4 and 20.8 (mean 7.5). Bulk δ¹³C had a narrow range from −31.0 ‰ to −29.5 ‰, with a mean of −30.3 ‰ and a very small standard deviation (SD) of 0.3 ‰ (Fig. 3). In contrast bulk δ¹⁵N had a wide range between 9.3 ‰ and 19.1 ‰ (mean 12.2 ‰, SD 1.6 ‰) (Fig. 3).

Figure 4 Key inorganic elements measured in stomach-oil deposit 3012MUM2, plotted as qualitative counts per second (cps) normalised to accumulation rate, to approximate element fluxes against the age-depth model (Fig. 2). Elements include S, Cl, Si, Ba, Al, Fe, Ca, Br P and Cu. Plots are coloured according to groupings on the PCA. Black lines indicate 7-point moving averages. Cluster boundaries (XRF-A, B and C) and XRF PCA are uniquely calculated from entire XRF dataset (Fig. S3).

3.3 Compositional changes in biomarkers

Here and in the discussion, we present the majority of biomarkers as fluxes (Fig. 3). Given similar patterns with biomarker concentrations (Figs. S2–S4) we consider the trends we have identified to be robust. The 3012MUM2 deposit samples likely comprised wax esters consistent with existing stomach oil studies (Imber, 1976; Lewis, 1966, 1969; Warham et al., 1976; Watts and Warham, 1976), and contributions from free lipids. Once saponified and derivatised, the extracts were dominated by fatty acids (FA) and alcohols (FAlc) (Fig. S3). In 3012MUM2, C_16:0 was the most abundant fatty acid (mean 38.9 % ± 8.8 % SD), followed by ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) (mean 23.6 % ± 10.9 % SD), C_14:0 (FA) (mean 22.1 % ± 5.2 % SD), C_18:0 (FA) (mean 8.4 % ± 2.1 % SD) and C_16:1 (FA) (mean 7.1 % ± 3.5 % SD).

By flux, FA were most abundant (C_14:0, C_16:0, C_16:1, C_18:0, ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$) (total mean 182.25 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, SD 116.69 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) followed by n-alkanols (C_14:0, C_16:0, C_18:0, C_20:0 and C_22:0) (FAlc) (total mean 47.95 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, SD 44.33 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$). The other compounds were less abundant, with cholesterol the highest (mean 9.95 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, SD 7.47 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) and phytol in lower concentrations (mean 0.31 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$; 0.22 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$). Phytol is formed from the ester-linked side chain of chlorophyll a and can therefore primarily be considered a biomarker for phytoplankton (Rontani and Volkman, 2003). Cholesterol is a ubiquitous marker but in this context could be considered a krill marker (both Antarctic and ice krill), since it can account for more than 76 % of total sterols in krill (Ju and Harvey, 2004) and is typically lower in concentration in fish (e.g. cholesterol in Dissostichus mawsoni ranges 4.7 %–14.3 % of total lipids; Clarke, 1984). We also assessed the distributions of sterols (e.g. 22-dehydrocholesterol), stanols (e.g. coprostanol) and cholesterol derivatives (e.g. cholesterol α/β-epoxide); cholesterol was the most abundant contributor (Fig. S4).

During saponification of the relevant wax esters, n-alkanols (FAlc) were also formed with saturated even chain lengths; C_14:0 (FAlc) to C_22:0 (FAlc) n-alcohols were the most abundant (Fig. S3). Amongst the n-alcohols C_16:0 (FAlc) was the most abundant (mean 63 % ± 2.9 % SD), followed by C_14:0 (FAlc) (mean 22.4 % ± 1.8 % SD) and C_18:0 (FAlc) (mean 13.0 % ± 2.1 % SD). C_20:0 (FAlc) (mean 1.3 % ± 0.4 % SD) and C_22:0 (FAlc) (mean 0.3 % ± 0.1 % SD) n-alcohols were relatively minor contributors.

Broadly, PC axis 1 (30 % variance explained) for the organic geochemistry indicators had high positive loadings in C:N ratio and ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA), while PC axis 2 (24 % variance explained) had strong positive loadings in C_16:0 (FA), C_18:0 (FA), δ¹⁵N and C:N ratio (Table S2 and Fig. S9).

3.4 Changes in biomarkers through time

Cluster analysis identified three statistically significant organic zones (Org-A, B and C) which are used as a framework to discuss the changes in key biomarkers through time. Although clusters were identified from organic analyses, there were 8 zones of colour visible in 3012MUM2 (Fig. S1c). In terms of lithological units, zone Org-C closely matched with zone 8 (darkest layer, black/brown), with zones 4–7 matching most closely with zone Org-B (lighter, yellowy/orange) and zones 1–3 matching with Org-A (medium-dark, grey brown to black brown).

3.4.1 Organic Zone – C (∼ 6390–5960 cal. yr BP)

The base of the deposit is marked by relatively high fatty acid (FA) fluxes: C_14:0 (FA), C_18:0 (FA) and ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) increased to maxima around ∼ 6300 cal. yr BP and declined by the top of the zone (5960 cal. yr BP) (Fig. 3). This is broadly coincident with a pulse of high δ¹⁵N (reaching ∼ 20 ‰, ∼ 6100 cal. yr BP), which follows the trend in C_14:0 (FA) but is slightly offset (Fig. 3). Ratios of fatty acids ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ and ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ remained low, with a pulse in ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ (FA) reaching a maximum at ∼ 6300 cal. yr BP (Fig. 3). A peak of 60 % C ∼ 6100 cal. yr BP is reflected in the maximum C:N_atomic ratio value of ∼ 21 (Fig. 3). δ¹³C fluctuated between −29.9 ‰ and −30.8 ‰. Phytol and cholesterol were relatively high with phytol fluctuating between 0.06 and 0.6 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ and cholesterol between 0.3 and 1.6 $\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 3). All n-alkanols (FAlc) increased to higher concentrations around ∼ 6200 cal. yr BP (Fig. S3). Most other biomarkers of sterols, stanols and cholesterol derivatives were generally high in this zone (Fig. S3).

3.4.2 Organic Zone – B (∼ 5960–4320 cal. yr BP)

The transition to zone B (∼ 5960 cal. yr BP) is marked by a rapid decrease in C_14:0 (FA), C_16:0 (FA) and C_18:0 (FA) fatty acid fluxes which thereafter remained broadly stable throughout the zone (Figs. 3, S2). Zone Org-B has markedly lower δ¹⁵N than zone Org-C, fluctuating between ∼ 9 ‰ and 13 ‰. In contrast, the C_16:1 (FA) and ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) fatty acid concentrations decreased from around 5400 cal. yr BP to the top of zone Org-B, which is also observed in phytol and cholesterol from ∼ 5100 cal. yr BP (Figs. S2 and S3). The major n-alkanols were markedly lower in zone Org-B compared to zone Org-C (Fig. S3), with a similar trend to lower values in C:N_atomic ratio (Fig. 3). Notably, at ∼ 4900 cal. yr. BP there was a short-lived unusual interval of high δ¹⁵N, high ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ (FA), ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ (FA) and a small increase in C:N ratio (Fig. 3). Organic carbon composition was variable in zone Org-B, fluctuation between −30.7 ‰ to −29.5 ‰ for δ¹³C, and between ∼ 36.6 % to 66 % for %C (Fig. 3). Other biomarkers of sterols, stanols and cholesterol derivatives were generally high at the bottom of zone Org-B (∼ 5960 cal. yr BP), typically decreasing by the top (∼ 4320 cal. yr BP) (Fig. S4).

3.4.3 Organic Zone – A (∼ 4320–1830 cal. yr BP)

For most biomarkers zone Org-A is a relatively stable period, with fluxes and concentrations of C_14:0 (FA), C_18:0 (FA) and ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) remaining low, broadly similar to the previous zone Org-B. The fatty acid ratios ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ and ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ are stable throughout (Figs. 3, S2 and S3). Between 4000 and 4320 cal. yr BP, δ¹⁵N was slightly higher (reaching 13.7 ‰) before decreasing towards the top of the deposit (Fig. 3). C_14:0 (FAlc), C_16:0 (FAlc) and C_18:0 (FAlc) n-alkanols also reached maxima at ∼ 3800 cal. yr BP (Fig. S3). C_16:1 (FA) and ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) concentrations were markedly lower than zone Org-B, and decreased further to the top of zone Org-A (Figs. S2 and S3), with a similar trend in ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) fluxes (Fig. 3). Of note, after ∼ 3300 cal. yr BP ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) ratios increased, and reached a maximum of ∼ 4 at 2500 cal. yr BP. Similarly, C:N_atomic gradually increased to a maximum of ∼ 17 by ∼ 2900 cal. yr BP, before decreasing to the top of the sequence (∼ 11.7), while δ¹³C displayed the reverse trend (increasing from ∼ −31 ‰ at ∼ 2700 cal. yr BP to −29.9 ‰ at the top of the sequence) (Fig. 3). Phytol and cholesterol were relatively low in zone Org-A, and both displayed trends of decreasing concentrations to the top of the sequence (Fig. 3). Other biomarkers of sterols, stanols and cholesterol derivatives were similarly also low throughout zone Org-A (Fig. S4).

3.5 Inorganic composition (XRF)

Cluster analysis identified three major XRF clusters: XRF-C (6390–4570 cal. yr BP); XRF-B (4570–4180 cal. yr BP); XRF-A (4180–1830 cal. yr BP). Zone XRF-C broadly coincides with organic zone Org-C and B, with XRF-B aligning (although slightly offset) with the uppermost part of zone Org-B. Zones XRF-A and zone Org-A are broadly aligned.

Based on XRF mean counts per second (cps) the largest elemental contributors to stomach-oil deposit 3012MUM2 were Cl, Ca, Fe, S, K, Br and P. Key XRF-derived inorganic compositions normalised to accumulation rate are presented in Fig. 4. Fluxes are characterised by higher variability in XRF-C with peaks in multiple elements at ∼ 6300, ∼ 5300–4900 cal. yr BP and at the top of the zone (∼ 4570 cal. yr BP). In XRF-B values reach a peak at ∼ 4300 cal. yr BP. Elemental fluxes are lower and relatively stable in XRF-A, with a small peak in some elements such as Fe ∼ 2400 cal. yr BP.

Based on the PCA of counts per second data, PC1 reflected 28 % of the variance (Table S1; Fig. S10) and included Si, Ba Al and Fe, which tend to peak around ∼ 6300 cal. yr BP in zone XRF-C (Fig. S5). PC2 reflected 22 % of variance (Table S1; Fig. S10), and included S, Cl, Br, P and Cu, which shift towards lower values and lower variability in XRF-A (Fig. 4).

4 Discussion

We found pronounced variations in geochemical proxies in stomach-oil deposit 3012MUM2 showing snow petrel diet was responding to changes in sea-ice conditions and ocean productivity related to changes in ice sheet extent and climate. In this discussion we first outline the rationale and caveats in our proxy interpretations, then explain the links between diet changes and sea ice environments through time.

4.1 Proxy interpretations

High levels of organic C and N in 3012MUM2 were similar to previously measured snow petrel stomach oil deposits (e.g. McClymont et al., 2022; Berg et al., 2023; Hiller et al., 1988), while bulk δ¹³C values were within the range of previous measurements of a Holocene deposit (Ainley et al., 2006). The fatty acid (FA) distributions found in the stomach oils are commensurate with previous stomach-oil and source end-member studies which suggest a diet of krill (mainly C_14:0 (FA), C_16:0 (FA)), squid (C_16:0 (FA)) and fish (C_16:0 (FA), ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA)) (Cripps et al., 1999; Lewis, 1966) as summarised in McClymont et al. (2022) and further explored in Berg et al. (2023). While some FAs are labile and susceptible to degradation, particularly unsaturated FAs (Stefanova and Disnar, 2000), the similar profiles shown by C_14:0 and C_18:0 FA with ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ FA suggests that changes in diet, not preservation, is the main driver for the latter FA (Fig. S3).

In the PCA biplots for the XRF analysis (Table S1 and Fig. S10) negative loadings from elements which are commonly associated with seabirds (e.g. P, Zn, Sr, Ni), including Cu (Shatova et al., 2016; Shatova et al., 2017; Castro et al., 2021; Sparaventi et al., 2021) support our interpretation of dietary sources to the deposit (Duda et al., 2021). In contrast, the positive loading of Fe, Al, Si and Ca on PC1 aligns with the main elements measured in the rock specimen taken from 3012MUM2 (Table S3), which reflects the local gneissic granite bedrock (Juckes, 1972). These elements are interpreted to reflect locally-derived, probably bedrock, erosional contribution. Cu was also present in 3012MUM2, which is known to be a key Antarctic krill marker as a Cu backbone structure is found in hemocyanin (Bridges, 1983) and has been observed in previous snow petrel deposits to be associated with Antarctic krill (McClymont et al., 2022). Broadly elevated Cu in zone Org-B with elevated C_14:0 FA and cholesterol (Fig. 3) also supports an Antarctic krill source for these components.

A wide range of δ¹⁵N values is recorded, which could reflect environmental change (e.g. changing circulation, degradation) or dietary change (e.g. baseline values in phytoplankton, trophic level). Past variations in the δ¹⁵N of circumpolar deep water (CDW) which upwells to the sea surface can contribute to high δ¹⁵N (Kemeny et al., 2018), but it is so far unknown if this occurred during the timescale of our deposit between 6390 and 1830 cal. yr BP (as opposed to glacial-interglacial timescales). Fragmented coral δ¹⁵N records of upper CDW indicate mid- to late Holocene values spanning ∼ 7 ‰–11 ‰ (Chen et al., 2023) but without the resolution to compare shorter-term changes with the variations we observed in 3012MUM2. Diatom-bound records of baseline δ¹⁵N show a long-term decline through the Holocene as nutrient availability has increased (e.g. Ai et al., 2020; Studer et al., 2018), which aligns in part with our overall trend. Previous work has cautioned that guano could contribute to stomach-oil deposits (Berg et al., 2019), which would also elevate δ¹⁵N (e.g. Wainright et al., 1998; Bokhorst et al., 2007). Finally, weathering would also elevate δ¹⁵N values due to microbial biosynthetic pathways causing ¹⁴N release (Macko and Estep, 1984). However, we only identified minor contributions from microbial fatty acids in 3012MUM2 to support this interpretation. In our discussion below, we focus on the environmental and dietary information provided by δ¹⁵N in combination with other proxy indicators.

We focus mainly on organic compounds for interpretation of palaeoclimate because dietary lipids and organic isotopes reflect the main snow petrel dietary sources which can, in turn, provide secondary information on sea ice distribution (e.g. McClymont et al., 2022). In particular, we focus on C_14:0 (FA) as an indicator of past Antarctic krill contributions to the diet, C_18:0 (FA) for fish, ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) for mainly fish, bulk δ¹⁵N as primarily an indicator of trophic status and productivity, and phytol as an indicator of past productivity. In terms of past ice configurations, we interpret a diet with more Antarctic krill and high productivity as recording snow petrels feeding in the open ocean (within the MIZ or the open ocean) (e.g. zone Org-C). Zones of higher productivity with reduced Antarctic krill contributions are inferred to represent enhanced snow petrel feeding closer to shore, in neritic zones, either due to a more proximal MIZ or in coastal polynyas along the Antarctic ice sheet margin (e.g. zone Org-B). While zones with lower productivity and accumulation rates reflect sea ice expansion, a less accessible MIZ and a more dense ice pack (e.g. zone Org-A).

4.2 Interpreting changes in sea ice in mid-Holocene stomach-oil deposits

By the start of the 3012MUM2 record at ∼ 6400 cal. yr BP, the Filchner and Ronne Ice Shelf had reached its modern position (Nichols et al., 2019; Johnson et al., 2019; Hillenbrand et al., 2017; Grieman et al., 2024). In other parts of Antarctica, cooling of the surface ocean and freshening of shelf surface waters since the middle-Holocene (Ashley et al., 2021) led to increased sea-ice concentrations and sea-ice duration over the continental shelf (Crosta et al., 2008; Mezgec et al., 2017; Johnson et al., 2021). Further offshore, sea ice became less extensive since ∼ 6500 cal. yr BP (Nielsen et al., 2004; Xiao et al., 2016), reflecting shifts in the latitudinal insolation and thermal gradients (Denis et al., 2010) as well as the multi-centennial expression of climate modes, such as ENSO (El Niño-Southern Oscillation) (Crosta et al., 2021). For the wider Southern Ocean, exemplified in East Antarctica, after 4500 cal. yr BP although sea ice and turbulence persisted in locations proximal to the ice shelf, it did not increase further offshore due to northward transport of subpolar surface waters as a response to southern Westerlies reinforcement (Denis et al., 2010).

The evolution of Holocene sea ice has been separated into three distinct phases around Antarctica (Early Holocene ∼ 11.5 to ∼ 8 ka BP; mid-Holocene ∼ 7 to ∼ 4–3 ka BP; the late Holocene ∼ 5–3 to 1–0 ka BP), with the phasing of these periods differing depending on regional response to long-term forcing (Crosta et al., 2022). Despite this wider knowledge there is an absence of any paleoenvironmental reconstructions of sea ice for the Holocene within the Weddell Sea region. Here we explore how each of the three statistically significant geochemical zones in the stomach oil deposit align with the timing or direction of the environmental shifts recorded by Crosta et al. (2022) and similar studies elsewhere in Antarctica.

4.2.1 Zone Org-C (6390–5960 cal. yr BP)

Between 6390–5960 cal. yr BP we infer that the snow petrels at Heimefrontfjella had a mixed diet which reflected access to areas of high productivity similar to today, with a MIZ situated in pelagic waters but with the potential for foraging over the continental shelf when the MIZ was situated closer to the coast. This interpretation is based on geochemical evidence of a diet sporadically high in either Antarctic krill (Euphausia superba) (evidenced by FA C_14:0, Cu: pelagic waters) or fish (evidenced by C_18:0 (FA): continental shelf waters). While we are unable to constrain sea-ice extent precisely, we interpret the high but variable snow petrel stomach oil accumulation rates to reflect repeated nest occupation from an accessible foraging habitat. The MIZ must therefore have been located within the modern snow petrel foraging range i.e. within ∼ 1200 km from Heimefrontfjella (Fig. 1; Wakefield et al., 2025).

We acknowledge that there are intervals in zone Org-C when the C_14:0 FA record does not align with Cu and cholesterol, which are instead positively correlated with indicators of guano: S and N (Cheng et al., 2016; Tatur et al., 1997; Roberts et al., 2017; Sun et al., 2000; Liu et al., 2005). The negative relationship between C:N and other elements (Cl, P, S) (Figs. S7 and S8) also supports a contribution from guano (Berg et al., 2019; Hiller et al., 1988). Although the presence of guano is unlikely to be an issue for interpreting dietary changes, as different biomarkers are involved (e.g. FA C_14:0, FA C_16:0, FA ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ for diet), it may have contributed to the elevated bulk δ¹⁵N values in zone Org-C. The δ¹⁵N values are exceptionally high in zone Org-C (∼ 20 ‰), exceeding those observed in modern Southern Ocean top predators (∼ 12 ‰–14 ‰) (Hückstädt et al., 2012; Valenzuela et al., 2018; Reisinger et al., 2016; Van Den Berg et al., 2021). The increase in δ¹⁵N is unlikely to reflect higher trophic status (Hodum and Hobson, 2000) or coastal foraging (St John Glew et al., 2021), because the fatty acid distributions suggest high Antarctic krill contributions to the diet, which should introduce lower δ¹⁵N than fish (Rau et al., 1992) and are associated with pelagic foraging (Brault et al., 2018). In zone Org-C, there could be potential for microbial degradation or enhanced guano inputs causing elevated δ¹⁵N. However, we also note that high δ¹⁵N are coincident with peak C_14:0 FA fluxes ∼ 6400 cal. yr BP, suggesting foraging could have occurred in pelagic waters with high nitrate utilisation (cf. Studer et al., 2018), consistent with nutrient limitation in highly productive polynyas. The presence of Antarctic krill markers suggests foraging in an open-ocean polynya, such as the intermittent Maud Rise polynya (Fig. 1) (Jena and Pillai, 2020; Turner et al., 2020; Holland, 2001). Feeding at the Maud Rise polynya is possible as snow petrels are able to forage at long ranges > 700 km (Honan et al., 2025), although they preferentially feed in broken ice in proximity to their nesting site during chick-rearing (Wakefield et al., 2025). Further work is required to identify whether other records with elevated δ¹⁵N can be identified, or to rule out an impact from guano or degradation on the unusual δ¹⁵N values in zone Org-C, especially as low %N could indicate degradation. At present, we interpret the data from zone Org-C to reflect foraging in a MIZ situated in pelagic waters but with the potential for foraging over the continental shelf during intervals when the MIZ was situated closer to the coast.

Figure 5 Summary plot comparing the marine and ice-core climate proxies (a–e) with 3012MUM2 geochemical proxy data (f–j). (a) James Ross Island (JRI) ice core temperature anomaly relative to 1961–1990 mean based on 100 year averages in Mulvaney et al. (2012); (b) Total yearly insolation at 65° S (Laskar et al., 2004); (c) Common δ¹⁸O signal in 5 Antarctic ice cores, interpreted by PCA axis 1 scores from Hodgson and Bentley (2013) based on common and residual signals in five Antarctic δ¹⁸O records; Vostok, Taylor Dome, EPICA Dome C, EPICA Dronning Maud Land and Dome Fuji (Masson-Delmotte et al., 2011); (d) Antarctic Peninsula sea ice from the ratio of highly-branched isoprenoid (HBI) diene to triene and HBI diene individually (core JPC38) (Barbara et al., 2016); (e) Adélie region sea ice from the ratio of HBI diene to triene and HBI diene individually (core U1257) (Ashley et al., 2021); (f) C_14:0 (FA) flux as a krill maker; (g) nitrogen stable isotopes (δ¹⁵N); (h) ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$ (FA) flux as an ice krill and fish marker; (i) phytol flux as a productivity marker; (j) cholesterol flux; and (k) Holocene time ranges for palaeoclimate periods from Crosta et al. (2022). Smooth lines are calculated from moving averages. Cluster boundaries (vertical black lines) are based on zones Org-A–C in Fig. 3.

4.2.2 Zone Org-B (5960–4320 cal. yr BP)

Between 5960–4320 cal. yr BP, we infer snow petrel foraging which is in a similar or slightly more productive MIZ than in zone Org-C due to elevated contributions from phytol. A mixed diet is inferred from intermediate contributions from C_14:0 (FA) and ${\mathrm{C}}_{\mathrm{18}:\mathrm{1}(\mathrm{n}-\mathrm{9})}$, from Antarctic krill and fish, respectively, with no change in diet compared to zone Org-C based on the fatty acid ratios ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ and ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ (Fig. 3). This mixed diet continues to indicate feeding both at the continental shelf edge and offshore. The main difference between zones Org-C and Org-B appear to be in the higher %N and sustained inputs of phytol, suggesting a more stable interval of higher productivity (and potentially better preservation) compared to the more variable zone Org-C. High accumulation rates in zone Org-B (Fig. 2) suggest frequent nest occupation during this zone, supported by elevated and sustained accumulation rates for the productivity marker (phytol) (Figs. 3 and 5) and inputs of potential guano-related elements (Fig. 4). Phytol is formed from chlorophyll a (Rontani and Volkman, 2003) and is brought into snow petrel diet through krill gut contents (Sargent and Falk-Petersen, 1981). As for zone Org-C, foraging in productive waters relatively proximal to the nest site could be explained either by access to foraging in coastal polynyas if the MIZ was situated further offshore, or by a MIZ relatively close to the continent e.g. if winter or spring sea-ice extents are relatively low, enabled in part by ice shelf fronts which were relatively stable during this period (Grieman et al., 2024). The variable but continued presence of Antarctic krill markers implies that there were several intervals where foraging occurred in a MIZ which was situated in pelagic waters between 5960 and 4320 cal. yr BP.

4.2.3 Zone Org-A (4320–1830 cal. yr BP)

We infer that between 4320 and 1830 cal. yr BP there was a shift towards a cooler, or neoglacial, phase which is characterised by more extensive sea ice. Our interpretation is based on the lower concentrations and fluxes in most of the organic proxies, aligning with sustained low accumulation rates in zone Org-A (Figs. 2–4). Although we observe no significant changes to diet according to the fatty acid ratios ${\mathrm{C}}_{\mathrm{16}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ and ${\mathrm{C}}_{\mathrm{18}:\mathrm{0}}/{\mathrm{C}}_{\mathrm{14}:\mathrm{0}}$ (Fig. 3), the loss of phytol indicates very low primary productivity, and the loss of cholesterol and Cu suggests that contributions from Antarctic krill were also reduced. In zone Org-A we suggest that the low rates of nest occupation and reduction in Antarctic krill and productivity markers indicate limited accessibility to snow petrel prey from the Heimefrontfjella site, which is consistent with the development of a more extensive or more dense sea ice pack in the north-east Weddell Sea which displaced the MIZ further offshore.

4.3 Coherence with other records of Holocene environmental change

Our records show associations between inference of past sea ice configurations made from our snow petrel stomach oil deposit and regional records of environmental change. For example, there are similarities and differences with records in the wider Weddell Sea region, likely due to the influence of the Weddell Gyre and atmospheric processes that determine sea ice configurations. Our interpretations of high productivity and krill between 6390 and 5960 cal. kyr BP (zone Org-C, Figs. 3 and 5), correspond with a relatively productive and an open-ocean (pelagic) type sea ice configuration from 7.2 cal. kyr BP at Herbert Sound and Croft Bay, in response to a wider regional mid-Holocene warming (e.g. core NBP0502-Site 2) (Totten et al., 2015). In the Firth of Tay (Fig. 1), this mid-Holocene climatic optimum spans 7800–6000 cal. yr BP (core NBP0602A) (Michalchuk et al., 2009) and 7750 and 6000 yr BP (cores NBP0602A-8B and NBP0703-JPC02) (Majewski and Anderson, 2009), and was attributed to seasonally-open marine conditions. Scotia Sea diatom records are also in phase with our results, with seasonally open water assemblages from 8300 to 2400 yr BP (core SS01) (Bak et al., 2007).

Our interpretation of zone Org-A having more sea ice corresponds with HBI records at the Vega Drift in JPC38, when interpreted carefully (Fig. 5d) (Barbara et al., 2016). Higher HBI diene levels in our deposit zone Org-A (4320–1830 cal. yr BP) correspond to increased sea ice cover. The fact that HBI:2 $/$ HBI:3 ratios vary and are low in zone Org-A are probably due to limitations of the proxy and the issue of normalization, linked with variable pelagic productivity. At Bransfield Strait (core D1–7) persistent sea ice is observed between 5.8 and 3.8 ka BP (core ANT28/D1–7) (Nie et al., 2022).

The neoglacial, sea-ice expansion and cooling phase suggested by zone Org-A is coherent with a sediment record on the Firth of Tay (core NBP0602A-8B and NBP0703-JPC02) which shows cooling from ∼ 3500 yr BP based on foraminifera (Majewski and Anderson, 2009). At Perseverance Drift, the warm interval persists longer than in our deposit, with high abundance of the foraminifera Globocassidulina spp. between 3400 and 1800 yr BP indicating incursions of Weddell Sea Transitional Water and a period of “freshening” consistent with open-marine or seasonally open marine conditions (Kyrmanidou et al., 2018) (cores JKC36 and JPC36), indicative of cooling.

The boundary between zone Org-B and A (4320 cal. yr BP) aligns with the transition to neoglacial conditions (e.g. Crosta et al., 2022). This mid-Holocene neoglacial is consistent with the James Ross Island (JRI) ice core temperature decline (Fig. 5a) (Mulvaney et al., 2012) and decreasing δ¹⁸O signal in 5 Antarctic ice cores (Hodgson and Bentley, 2013; Masson-Delmotte et al., 2011) (Fig. 5c). We find that the onset of the neoglacial is remarkably in phase with a switch in sea ice biomarkers (both HBI diene and ($\mathrm{HBI}:\mathrm{2}/\mathrm{HBI}:\mathrm{3}$) ratio) indicative of more sea ice in the distant Adélie region in core U1357 (Ashley et al., 2021) (Fig. 5e). A marked transition is also seen in South Atlantic cores (Hodell et al., 2001) where IRD (as % lithics) increase markedly from ∼ 5 ka BP suggesting cooling waters (together with concomitant changes in diatoms SST index and δ¹⁸O on diatoms) which have been linked to the arrival of more sea ice from the Weddell Sea region (core TTN057-13-PC4). Total yearly insolation (e.g. at 65° S) features an ongoing decline during this period (Fig. 5b) but is unlikely to be the sole driver of the neoglacial (Divine et al., 2010; Renssen et al., 2005).

5 Conclusions

Analyses of accumulation rates and a range of biomarkers in a snow petrel stomach-oil deposit from the Heimefrontfjella Range have been used to infer sea ice conditions in the foraging area in the northeastern Weddell Sea over the Holocene. The record has three significant zones. In the first zone between 6390 and 5960 cal. yr BP, high stomach oil accumulation rates and high concentrations of both fish and krill fatty acids suggest easy access to productive foraging grounds. This is consistent with foraging both at the continental shelf in the MIZ and offshore. In the second zone, between 5960 and 4320 cal. yr BP productivity remained high but there was evidence of a more mixed diet suggesting foraging both at the continental shelf edge and in open water, with high productivity due to summer sea ice retreat reaching the shelf edge more frequently. In the third stage, between 4320 and 1830 cal. yr BP low stomach oil accumulation rates and reductions productivity markers indicate that increased sea ice conditions consistent with the transition to neoglacial conditions seen in a range of records from the northwestern Weddell Sea (Barbara et al., 2016; Nie et al., 2022), from ice records on the Antarctic continent (Hodgson and Bentley, 2013), and changes elsewhere in the wider Antarctic region including Adélie Land (Ashley et al., 2021) and the South Atlantic (Hodell et al., 2001; Divine et al., 2010). This study has also shown, for the first time, the utility of phytol and cholesterol for tracking past snow petrel diet and thus interpretation of environmental conditions in their foraging areas.