Summer sea-ice variability on the Antarctic margin during the last glacial period reconstructed from snow petrel (Pagodroma nivea) stomach-oil deposits

Antarctic sea ice is a critical component of the climate system affecting a range of physical and biogeochemical feedbacks and supporting unique ecosystems. During the last glacial stage, Antarctic sea ice was more extensive than today, but uncertainties in geological (marine sediments), glaciological (ice core), and climate model reconstructions of past sea-ice extent continue to limit our understanding of its role in the Earth system. Here, we present a novel archive of past sea-ice environments from regurgitated stomach oils of snow petrels (Pagodroma nivea) preserved at nesting sites in Dronning Maud Land, Antarctica. We show that by combining information from fatty acid distributions and their stable carbon isotope ratios with measurements of bulk carbon and nitrogen stable isotopes and trace metal data, it is possible to reconstruct changing snow petrel diet within Marine Isotope Stage 2 (ca. 24.3–30.3 cal kyr BP). We show that, as today, a mixed diet of krill and fish characterizes much of the record. However, between 27.4 and 28.7 cal kyr BP signals of krill almost disappear. By linking dietary signals in the stomach-oil deposits to modern feeding habits and foraging ranges, we infer the use by snow petrels of open-water habitats (“polynyas”) in the sea ice during our interval of study. The periods when consumption of krill was reduced are interpreted to correspond to the opening of polynyas over the continental shelf, which became the preferred foraging habitat. Our results show that extensive, thick, and multiyear sea ice was not always present close to the continent during the last glacial stage and highlight the potential of stomach-oil deposits as a palaeoenvironmental archive of Southern Ocean conditions.

1 Introduction

Antarctic sea ice is globally important: large seasonal fluctuations in its spatial extent influence planetary albedo, ocean–atmosphere exchanges of heat and climatically active gases including CO₂, and the formation of intermediate water masses and deep water masses which create the world's largest sink of heat and carbon (e.g. Ackley et al., 2015; Delille et al., 2014; Arrigo et al., 2008; Frölicher et al., 2011). Antarctic sea-ice extent is projected to decline by 16 %–67 % by 2100 depending on greenhouse gas emission scenarios, but with low confidence due to a wide range of model responses (Collins et al., 2013).

The geological record offers an opportunity to set the relatively short (∼ 50 years) instrumental observations of sea ice into a longer-term context. During the Last Glacial Maximum (LGM, ∼ 19–23 ka) (Tierney et al., 2020), the development of thick, multiyear, and more extensive sea ice is inferred from the disappearance of, or diagnostic changes in, assemblages of marine microfossils from many Southern Ocean sediment cores (e.g. Hillenbrand and Cortese, 2006; Grobe and Mackensen, 1993; Bonn et al., 1998; Lucchi et al., 2002; Gersonde et al., 2005; Benz et al., 2016; Collins et al., 2012; Allen et al., 2011). An expanded “sea-ice cap” is proposed to have increased deep-ocean storage of CO₂ by limiting the air–sea gas exchange and by enhancing CO₂ export through increased Antarctic Bottom Water production (e.g. Stephens and Keeling, 2000; Ferrari et al., 2014).

Figure 1 Sea-ice records of the last glacial stage for the Atlantic sector of the Southern Ocean. The stomach-oil deposit presented here (WMM7) is shown by the yellow star (Untersee Oasis; 71^∘22^′ S, 13^∘19^′ E). The likely snow petrel foraging range is indicated by the shaded ellipse: although snow petrels can forage 2600 ± 1000 km from their nests (pale shading) they are found in their highest densities within 700 km of the nest site during the breeding season (darker shading) (Delord et al., 2016). Boundaries of the maximum summer sea-ice extent (22–30 ka; Allen et al., 2011; Collins et al., 2012) and LGM summer sea-ice extent (19–23 ka; Gersonde et al., 2005; Xiao et al., 2016) are marked. Sites where polynya activity has been inferred for the last glacial stage are shown in red (Bonn et al., 1998; Smith et al., 2010; Sprenk et al., 2014; Weber et al., 1994), noting that the KC sites (Smith et al., 2010) and south-eastern Weddell Sea sites (PS1599, 1789, 1790, and 1974; Weber et al., 1994; and PS1795; Sprenk et al., 2014) infer polynyas opening “upstream” of the core locations. Note that the continental shelf edge is effectively also delineated by the solid red line of LGM maximum ice-sheet extent.

Antarctic sea-ice environments today are dynamic and complex at a range of temporal and spatial scales (Parkinson, 2019; Turner et al., 2020), and they include open waters within the sea-ice pack (polynyas), which span a large size range (1000–∼ 400 000 km²; Arrigo and van Dijken, 2003), are poorly represented in models (Mohrmann et al., 2021), and yet impact ocean circulation, formation of sea ice, and air–sea gas exchange (Mohrmann et al., 2021; Morales Maqueda et al., 2004). A similarly complex picture has emerged for the last glacial period: in the Atlantic sector of the Southern Ocean, spanning the Weddell, Lazarev, and Scotia seas (Fig. 1) there is fragmentary empirical evidence for polynyas during and before the LGM, as detailed by intervals of high productivity or bottom-current activity in marine sediment cores, and occupation of nesting sites by seabirds that require open water within their foraging range (e.g. Smith et al., 2010; Sprenk et al., 2014; Mackensen et al., 1989; Thatje et al., 2008; Berg et al., 2019; Weber et al., 1994). Millennial-scale variability in the extent of the seasonal sea-ice zone has also been described (Collins et al., 2012; Gersonde et al., 2003; Rae et al., 2018), but large (10 %–20 %) uncertainties remain in both models and geological datasets (Collins et al., 2012; Roche et al., 2012; Gersonde et al., 2003). Maximum summer sea-ice extent was likely reached at 30–22 ka (Collins et al., 2012; Gersonde et al., 2003; Allen et al., 2011; Xiao et al., 2016), raising the possibility that Southern Ocean sea ice–climate feedbacks were more important before the global LGM (Allen et al., 2011; Xiao et al., 2016).

Here, we investigate changes in the sea-ice environment during the time of proposed maximum sea-ice extent (30–22 ka) by analysing a sequence of preserved stomach oils of snow petrels (Pagodroma nivea). Although sometimes referred to as “Antarctic mumiyo” (e.g. Berg et al., 2019; Thor and Low, 2011), this is a misnomer because its biological origin is so different from the original use of “mumiyo”, which was used to describe an organic, tar-like substance of unknown origin found in high-altitude rocks and caves, especially in Asia (Hiller et al., 1988, and references therein; Aiello et al., 2011). We therefore refer to Antarctic “stomach-oil deposits” here.

Progressive accumulation of regurgitated stomach oils at the entrance to snow petrel nests in rock crevices and scree slopes (e.g. Johansson and Thor, 2004) leads to stratified deposits which may span several millennia (Hiller et al., 1988; Wand and Hermichen, 2005; Hiller et al., 1995; Thor and Low, 2011; Berg et al., 2019). The deposits include a mixture of stomach oils, guano, aeolian dust, and other mineral matter (Berg et al., 2019). Snow petrels forage in open waters at the margins of, or amongst, the sea ice within a restricted range (several 100 km) from their nesting sites (Delord et al., 2016; Marchant and Higgins, 1990). Snow petrels feed on fish, krill, and squid (Ainley et al., 1984; Delord et al., 2016), and the stomach oils preserve a biochemical fingerprint of these prey. For example, krill tend to have high abundances of C_14 : 0, C_16 : 0, and C_18 : 1 fatty acids and a variety of carotenoid pigments (e.g. Cripps et al., 1999; Färber Lorda and Ceccaldi, 2020), as well as Cu, Zn, and Fe (Rainbow, 1989; Palmer Locarnini and Presley, 1995; Liu et al., 2013), whereas squid and fish are dominated by C_16 : 0 and C_18 : 0 fatty acids, respectively (Lewis, 1966). The presence of Cu in crustaceans (including krill) is linked to hemocyanins in invertebrates (Bridges et al., 1983), which provide an alternative oxygen carrier to the hemoglobins of vertebrate blood. Stable isotope ratios provide useful proxies for trophic level (δ¹⁵N) and foraging habitat (δ¹³C), enabling diet comparisons between snow petrels and other seabirds (Delord et al., 2016; Cherel et al., 2010; Steele, 2005; Phillips et al., 2011).

Here, we undertake elemental scans, organic geochemistry, and stable isotope analysis to investigate snow petrel diet during the last glacial stage from a stomach-oil deposit collected at the Untersee Oasis in central Dronning Maud Land (DML). Since snow petrels have a restricted foraging range during the breeding season (Delord et al., 2016), our analyses enable us to reconstruct changes in their diet and foraging habitat relatively close to the Antarctic margin (Fig. 1). In turn, we hypothesize that the biochemistry of the stomach-oil deposits provides diagnostic signatures of snow petrels foraging in sea ice over or beyond the continental shelf and in polynyas, offering novel insights into the evolution of sea-ice environments during the last glacial stage.

2 Materials and methods

2.1 Untersee Oasis sequence WMM7 and its regional context

The stomach-oil deposit WMM7 was collected during the GeoMaud expedition (1995–1996) from the Untersee Oasis (71^∘21.6^′ S, 13^∘18.96^′ E) in DML (Fig. 1) (Wand and Hermichen, 2005). The deposit was retrieved from under a boulder on a steep slope west of the Untersee, ca. 320 m above the lake surface, at an elevation of 880 m above sea level. WMM7 was kept in dark, cold storage (4 ^∘C, 60 % humidity) at the Alfred Wegener Institute, Germany, until analysis. The sample had a waxy consistency, so to preserve its internal structure it was frozen immediately before slicing and subsampling at Durham University. A central slab of WMM7 was sectioned for non-destructive, high-resolution X-ray fluorescence (XRF) analysis (40 mm thick, 50 mm wide, and 150 mm long; marked in Fig. 2a). Samples for stable isotope and lipid analysis were taken adjacent to the XRF slab (Fig. 2a).

Figure 2 Stomach-oil deposit WMM7 stratigraphy and the two age–depth models explored in Sect. 2.2 and Table 1. (a) Cross-section of the section used for subsampling, scaled with panels (b) and (c). Note the clear laminae, which slope to the right away from the main slab. A zone of darker brown and less well-defined laminae is found in the central section (Unit II in panel b). A 5 cm wide slab was cut for XRF analysis. Data from the centreline scan are shown in Fig. 3 and used for statistical analysis. The results of additional scans to the left and right of centre and from two new cuts (3 cm left, 3 cm right; see arrows) are shown in Fig. A1. Biomarker and bulk stable isotope samples were taken immediately to the left of the XRF slab. (b) Three-unit stratigraphy and two age–depth models constrained by six bulk radiocarbon dates, calibrated using either a Holocene “no sea ice” ΔR or a glacial-stage “maximum sea ice” ΔR (see Table 1 for ¹⁴C calibration). Figures 3–6 apply the Holocene “no ice” ΔR age model; (c) accumulation rate between age control points for the two age–model scenarios. Linear interpolation was applied to generate an age–depth model for all sampling points.

Before subsampling, WMM7 weighed ca. 1 kg and was 155–194 mm thick, 144 mm wide, and 120 mm deep (Fig. 2a). It had an irregular, mammillated outer surface but was characterized internally by millimetre-scale laminae that were traced through the deposit. The laminae visibly slope away from the centre of the deposit, showing that the deposit progressively draped over the underlying rock. Three distinct units were observed in the sequence during sampling (Fig. 2b): the oldest part of the sequence (160–110 mm depth, Unit III) was characterized by laminated, yellow–brown deposits overlain by a zone of relatively dark, brown–black deposits with more poorly defined laminae (Unit II). The Unit II–I transition was gradual (80–70 mm), and Unit I contained yellow–brown deposits with sub-millimetre-scale laminae (Unit I). No hiatuses were visible in the stratigraphy; these would be indicated by breaks in the structure or visible sediment, which would have been deposited and concentrated during an interval when stomach oils were not being deposited. The linear age model (Fig. 2) also indicated continuous accumulation.

2.2 Radiocarbon analysis

Previous ¹⁴C analysis of the top of WMM7 provided an age constraint of 21 551 ± 110 years (Berg et al., 2019). To establish an age–depth model, six additional bulk samples were dated (Table 1), and linear interpolation was applied between dating points. Each sample was digested in 2 M HCl (80 ^∘C, 8 h), washed free from mineral acid with deionized water, and then dried. Samples were graphitized using an automated graphitization system (Rethermeyer et al., 2019) and analysed for ¹⁴C by accelerator mass spectrometry (AMS) at CologneAMS, Germany. Radiocarbon ages were converted to calendar ages (Table 1) using the MARINE20 radiocarbon age calibration (Heaton et al., 2020) by adopting two approaches to explore the potential influence of changing sea-ice extent on ΔR (Heaton et al., 2021). First, radiocarbon ages were corrected for the Southern Ocean marine reservoir effect by applying a ΔR of 670 ± 50 years from pre-bomb ¹⁴C ages at Hope Bay in the western Weddell Sea (Björck et al., 1991). As this ΔR reflects Holocene conditions, wherein sea ice is low or absent from the region, this calibration is described as a “no sea ice” scenario (Table 1; Heaton, 2021). As variable sea-ice extent introduces additional uncertainty to radiocarbon age calibration in polar regions (Heaton et al., 2020), a second scenario was explored by applying a modelled glacial-stage maximum sea ice ΔR of 2480 ± 50 years (Table 1). As snow petrels forage in surface waters within sea ice or close to the sea-ice edge, we suggest that the glacial-stage ΔR is likely an overestimate. For the results presented here, the datasets are described according to the Holocene “no sea ice” age model (Table 1), recognizing that this provides a maximum age for the WMM7 record. This approach also ensures that our results can be incorporated into studies using stomach-oil deposits for constraining ice-sheet thinning histories, for which they are used to indicate the onset of ice-free conditions (e.g. Hiller et al., 1988; Hillenbrand et al., 2014, and references therein). However, we acknowledge that the application of a constant, modern ΔR to a region affected by greater sea-ice extent during MIS 2 leads to larger uncertainty and potential bias in our calendar ages than indicated by the results in Table 1 (Heaton et al., 2020, 2021).

Table 1 Radiocarbon dates and calibrated ages. All ¹⁴C analyses were performed on bulk samples at CologneAMS, Germany. COL3022 was previously published (Berg et al., 2019). All calibrations to calendar ages used MARINE20 (Heaton et al., 2020). To explore the likely range of impacts of sea ice on our ¹⁴C calibrations, we first apply the nearest Holocene ΔR of 670 ± 50 years (Björck et al., 1991), which assumes low or absent sea ice sampled by WMM7 (ΔR^Holocene). A second calibration (ΔR^{glacial stage}) assumes enhanced sea-ice cover, as suggested for the last glacial stage, generated by adding 1800 years of additional ocean ¹⁴C depletion as suggested by Heaton et al. (2020, 2021).

2.3 XRF analysis

XRF analysis was performed using the ITRAX core scanner at the National Oceanography Centre, Southampton, UK (Croudace et al., 2006). A molybdenum X-ray source was used (45 kV, 40 mA) at a step size of 200 µm and exposure time of 400 s per increment. To test for potential contamination during subsampling and to explore internal consistency of the recovered signals, a further five profiles were measured on the central slab (Figs. 2a, A1). Multiple XRF scans confirmed that the major patterns outlined below are consistently recorded in WMM7, with differences accounted for by changes in the orientation and spatial continuity of the laminae (Fig. A1).

2.4 Bulk organic matter elemental composition and stable isotope analysis

Directly adjacent to the left margin of the XRF slab, 15 contiguous 10 mm subsamples were taken by scalpel for bulk stable isotope analysis (Fig. 2a). Carbon and nitrogen stable isotope analyses were performed using a Costech elemental analyser (ECS 4010) connected to a Thermo Scientific Delta V Advantage isotope ratio mass spectrometer. Carbon isotope ratios are corrected for ¹⁷O and reported in standard delta (δ) notation in per mil (‰) relative to Vienna PeeDee Belemnite (VPDB). Nitrogen isotope ratios are reported against atmospheric nitrogen (AIR). Isotopic accuracy was monitored through multiple 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, which provided a linear range for δ¹³C between −46 ‰ and +3 ‰ and for δ¹⁵N between −4.5 ‰ and +20.4 ‰. Analytical uncertainty was typically ±0.1 ‰ (2 SD) for replicate analyses of the international standards and < 0.2 ‰ (2 SD) on replicate sample analysis. Total organic carbon (TOC) and nitrogen data were obtained as part of the isotopic analysis using the internal standard: glutamic acid (40.82 wt % C, 9.52 wt % N).

2.5 Biomarker distributions and stable carbon isotope analysis

A 2 mm stainless-steel dermal punch was used to extract samples for biomarker analysis from the same line as the bulk stable isotope samples (Fig. 2a), including triplicates at 1 cm depth and 8 cm depth. Lipids were extracted from 0.03–0.3 g of each sample using repeated ultra-sonication (3 × 10 min) in 10 mL dichloromethane (DCM)–methanol (3 : 1), following addition of two internal standards of known concentration (androstanol, hexatriacontane). Ultrasonic extraction yielded between 16 and 43 mg of total lipid extract.

An aliquot of the total lipid extract was dissolved in acetone and analysed by UV–Vis spectrophotometry using a Dionex HPLC Quaternary pump and photo-diode array detector (McClymont et al., 2007). A grass standard and blank acetone were injected regularly to monitor potential instrument drift. The relative absorbance at wavelengths characteristic of chlorophyll derivatives (“chlorins”) and potential carotenoid derivatives, at 410, 435, and 665 nm (Jeffrey et al., 1997), was calculated for each sample:

$$\begin{array}{}\text{(1)}& P\mathit{\lambda }={\displaystyle \frac{(A\mathit{\lambda }\times \mathrm{DF})}{M}},\end{array}$$

where Pλ is the relative magnitude of absorbance for a given wavelength (abs. g⁻¹​​​​​​​), Aλ is the integrated area for that wavelength averaged over three repeat measurements, DF is the dilution factor (i.e. the aliquot of the total sample injected), and M is the mass of material which was extracted. To take into account the potential influence of variable organic matter deposition or preservation, we also normalized all photosynthetic pigment absorbance to TOC content (generated from Sect. 2.4). There may be multiple sources of our target pigments (Jeffrey et al., 1997), so we refer to them as P410, P435, and P665 to identify the trends in absorbance for each specific wavelength.

The remaining extracts were saponified using 2 mL KOH (8 %) in methanol (95 %) and heated for 2 h at 70 ^∘C. Neutral lipids were extracted using hexane; the remaining extracts were acidified using 2 M HCl, and the fatty acids were extracted using hexane. An internal standard of known concentration (heptadecanoic acid, 0.4 mg mL⁻¹) was added to the fatty acid fractions before generating fatty acid methyl esters (FAMEs) by methylating with 3 mL methanol–HCl (95 : 5) for 12 h at 70 ^∘C, then allowing it to cool to room to temperature. After adding 4 mL of DCM-rinsed H₂O to each sample, FAMEs were recovered sequentially using hexane followed by hexane–dichloromethane (4 : 1) and pooled. FAME fractions were then taken to dryness in a stream of N₂. The isotopic value of the methanol was determined through methylation of a phthalic acid with known isotopic value (Lee et al., 2017).

The FAMEs were identified and quantified with a Thermo Trace 1310 gas chromatograph linked to an ISQ LT single quadrupole mass spectrometer (GC-MS). Chromatographic separation was performed with a Restek Famewax 30 m × 0.25 mm × 0.25 µm column. Sample extracts were injected (0.8 µL) in CT split mode (80 : 1 ratio) into a programmable temperature vaporizing (PTV) injector at 230 ^∘C with a constant helium carrier gas flow rate of 1.5 mL min⁻¹. The oven temperature was initially set to 100 ^∘C for 3.0 min, then ramped at 2 ^∘C min⁻¹ to 230 ^∘C, and held at 230 ^∘C for 10 min. The transfer line temperature was set to 230 ^∘C and the ion volume temperature to 220 ^∘C. A mass range of 38 to 600 $m/z$ was scanned every 0.5 s, giving at least 20 data points per compound peak. Samples were identified and quantified by comparison to a Supelco 37 component ME mix (CRM47885, Sigma-Aldrich) with peak area ratios calculated with reference to the peak area of the internal standard. Spectral confirmation was performed using a NIST EI reference library.

The carbon isotopic compositions (δ¹³C) of individual saturated FAMEs were determined on eight samples (22–132 mm depth) using a Thermo GC-C-IRMS system at the Department of Geography, Durham University. All samples were run in duplicate. The system was composed of a Trace 1310 GC coupled to a Thermo Delta V Plus via a GC IsoLink II and a Restek Famewax 30 m × 0.25 mm × 0.25 µm column. Samples were injected (2 µL) in splitless mode into an S/SL injector set to 240 ^∘C with a helium carrier gas flow rate of 1.5 mL min⁻¹. The oven temperature was initially set to 50 ^∘C for 1 min, then ramped to 100 ^∘C at 10 ^∘C min⁻¹, then ramped to 240^∘ at 3 ^∘C min⁻¹, and finally held for 10 min at 240 ^∘C. The alumina (with CuO, NiO, and Pt wires) combustion reactor was operated at 1000 ^∘C and conditioned with oxygen each day immediately before use. CO₂ reference gas pulses were introduced at the start and end of each chromatogram to provide an isotope ratio reference point and to check the system stability during the run. All the quantified FAME peaks in the selected carbon range were baseline-resolved apart from the cis- and trans-isomers of C_18 : 1. A FAME standard (CRM47885 – Supelco 37 component mix, Sigma-Aldrich) was run with each batch of isotope reference standards to confirm the retention time of the FAME peaks. Individual FAME isotope ratio values were corrected using a multipoint linear normalization of a C₁₄–C₂₀ FAME reference material (F8-3 standard provided by Arndt Schimmelmann​​​​​​​, Indiana University, Bloomington). Reference standard FAMEs from C₁₄–C₂₀ were used to generate the normalization curve, covering δ¹³C values from −23.24 ‰ to −30.92 ‰. Reference material F8-3 was also used to determine the usable amplifier signal range, which minimized the residuals, and gave an r² value of at least 0.995 for the normalization plot. The concentration of the F8-3 FAME standard used for the linear normalization was adjusted to obtain reference amplifier intensities within this range (1000 to 6000 mV). Each sample was then diluted and pre-run to determine the optimum solvent volume required to fit within the amplifier signal range of the reference standards. The long-term pooled standard deviations of the ME F8-3 reference compounds were all < 0.35 ‰ (ranging from 0.21 ‰ to 0.34 ‰, n = 21). Fatty acid δ¹³C was calculated through mass balance corrections of the measured FAME δ¹³C and known δ¹³C of the added methyl group (Lee et al., 2017). It is reported relative to Vienna PeeDee Belemnite (VPDB).

2.6 Statistical analysis

Principal component analysis (PCA) was performed to investigate the potential for co-evolving relationships in the stomach-oil geochemistry using the PAST3 software (Hammer et al., 2001). We focussed on the XRF data for PCA since it was analysed at the highest resolution and with the highest number of variables, and we used the variance–covariance matrix since all variables were measured in the same units (i.e. counts). Elements which consistently recorded counts below 500 were excluded from the analysis.

Cluster analysis was performed to identify units of similar geochemical composition using the rioja package in R 3.6.0 (Juggins, 2020), whereby a hierarchical clustering is performed, constrained by sample order, and compared to a broken-stick model of a random distribution of zones within a sequence (Bennett, 1996). Due to the different sampling resolutions of the geochemical methods employed here, we resampled the XRF data by averaging element counts across the same depth window sampled by each of the discrete geochemical measurements (10 mm diameter). This resampling approach does not alter the main signals of the principal components, although additional clusters are identified using the original XRF data (Fig. A2). We present the resampled XRF PCA here as an independent measure of geochemical change that is comparable to the resolution of the lipid distribution and stable isotope signatures.

3 Results

3.1 WMM7 age model

Stomach-oil deposit WMM7 spans 23 987 (23 668–24 366) cal yr BP through 30 307 (29 949–30 685) cal yr BP using the Holocene “no sea ice” ΔR for calibration (Table 1). The unit transitions are well constrained by ¹⁴C dates (Fig. 2b). The accumulation rate is remarkably stable at 26.6 ± 3.7 mm cal kyr⁻¹ (Fig. 2c).

3.2 Elemental composition (XRF)

Deposit WMM7 was dominated by Fe, Ca, Cu, and Ti. Elements likely to be indicative of minerogenic inputs (e.g. Ti, Al, K, Ca), for example from windblown particulates, had low signals overall (Fig. A3). There were no horizons of elevated minerogenic inputs, which might be expected during a hiatus in stomach-oil accumulation. Elements with both a biogenic and minerogenic (e.g. Cu, Fe) or principally biogenic source (e.g. As, Zn) (Huang et al., 2009; Liu et al., 2013) varied both at the millimetre-scale (likely reflecting individual laminae) and with depth in the deposit (Fig. A3). To account for the influence of minerogenic contributions to the elemental composition of WMM7, all data were normalized to Ti (Figs. 3, A4). There were no clear down-core trends in Fe $/$ Ti and Si $/$ Ti (Fig. 3a, b) indicating a dominant minerogenic source for Fe and Si. Cu $/$ Ti was generally elevated at 30.3–28.6 ka with minima at ca. 29.3 ka (Unit III) and very low values between 28.5 and 27.4 ka (Unit II). The highest Cu $/$ Ti was recorded between 27.4 and 26.9 ka, then remained elevated until 25.7 ka (Unit I), until values decreased to the top of the deposit (Fig. 3c). In contrast, Br $/$ Ti and S $/$ Ti were low from 30.3 to 28.7 ka (Unit III), then elevated between 28.7 and 27.4 ka (Unit II; Fig. 3d, e). A gradual increase in Br $/$ Ti and S $/$ Ti was then recorded until a second maximum at 25.5 ka (Unit I); both Br and S then decreased to 24.5 ka. A dominant biogenic source for the trends in Si, Cu, Br, and S is inferred for WMM7, noting that the patterns in Br $/$ Ti and S $/$ Ti were different to those of Cu $/$ Ti and may include a sea-salt origin given some similarities to Cl $/$ Ti (not shown; r²=0.36).

Figure 3 XRF characterization of stomach-oil deposit WMM7. All XRF data are normalized to Ti. The age model uses Holocene “no sea ice” ΔR (Table 1). The ∼ 100-year smoothing is shown in coloured lines. Cluster boundaries are determined by constrained hierarchical clustering analysis in the rioja package of R (Juggins, 2020) compared to a broken-stick analysis (Bennett, 1996). Three significant clusters were identified (h), which broadly coincide with Units I–III, noting the gradual Unit II–I boundary.

The PCA of the XRF data confirms that there were two main patterns in the elemental composition of WMM7 (Fig. 3f, g; Table 2). PC1 accounted for 93 % of the variance and is driven by positive loading from Fe, Cu, and Ca, which could have both minerogenic and biogenic origins. PC2 accounted for 6 % of the variance, driven by Cu (positive loading), Ca, and Fe (negative loading), but PC2 was not a statistically significant factor using the broken-stick threshold in PAST3. Peaks in PC1 occurred before 28.5 ka (Unit III), between 27.1 and 26.3 ka (lower Unit I), and at 24.7 (uppermost Unit I). PC1 minima occurred in Unit II (28.5–27.6 ka) when high Br $/$ Ti and S $/$ Ti confirmed enhanced organic matter inputs. In contrast, PC2 broadly followed the Cu $/$ Ti trend, recording declining values from 29.9–27.6 ka (Units III and II) and peak values at 27.1 ka (lower Unit I, Fig. 3). Three clusters were identified in the resampled XRF data (Fig. 3i), broadly aligned to Units I–III.

Table 2 Principal component analysis. The elements contributing to each principal component are shown; bold font highlights Fe, Cu, and Ca, which dominate the loadings. PC1 explains 93 % of the variance (statistically significant), whereas PC2 explains only 6 % of the variance and is not statistically significant. PC1 and PC2 through time are shown in Fig. 3.

3.3 Organic matter composition

All WMM7 samples were highly organic-rich (27 %–50 % C, 11 %–13 % N, Fig. A5), with C $/$ N ratios between 6 and 12 (Fig. 4a). The C $/$ N ratios were consistent with high lipid and low N contributions as observed in modern Weddell Sea fish, bird, and seal tissues (Rau et al., 1992). The extracted lipids were dominated by fatty acids and alcohols (not presented here), consistent with the presence of triglycerides and wax esters observed today in stomach oils of procellariform seabirds, i.e. albatrosses and petrels (e.g. Horgan and Barrett, 1985). The saturated fatty acids C_16 : 0 (46 ± 4 %) and C_14 : 0 (24 ± 1 %) dominated all samples, with additional minor contributions from C_18 : 1 (16 ± 3 %), C_18 : 0 (8 ± 1 %), and C_16 : 1 (5 ± 1 %) (Fig. 4b, c).

Figure 4 Organic matter composition of stomach-oil deposit WMM7. The age model uses Holocene “no sea ice” ΔR (Table 1). Unit II is highlighted with brown shading, noting the gradual Unit II–I boundary. Cluster boundaries are determined by constrained hierarchical clustering analysis in the rioja package of R (Juggins, 2020) compared to a broken-stick analysis (Bennett, 1996). Three significant clusters were identified (j), which broadly coincide with Units I–III, noting the gradual Unit II–I boundary. The Cu $/$ Ti data from Fig. 3 are shown in (g), noting the reverse x axis to align with the interpretation fatty acid indicators of krill in (e) and (f).

Table 3 Examples of snow petrel prey biochemistry and oceanographic conditions related to snow petrel diet and prey distributions. Snow petrels forage in close association with the sea ice, either at the sea-ice edge or in leads or polynyas within the sea-ice pack. The prey biochemistry information is used as a framework to interpret the chemical signatures recorded in stomach-oil deposit WMM7, including our use of fatty acid ratios to assess relative contributions of prey (Figs. 4 and 5).

The dominant fatty acids in WMM7 are consistent with a snow petrel diet of krill (C_14 : 0, C_16 : 0), squid (C_16 : 0), and fish (C_16 : 0, C_18 : x) (Cripps et al., 1999; Lewis, 1966) (Table 3). To explore whether prey contributions changed through time, the relative abundance of the major fatty acids were compared as ratios (Fig. 4d–f), assuming a krill (C_14 : 0), fish (C_18 : x), or mixed source (C_16 : 0) (Table 3). Minor variations to C_18 : 0 $/$ C_16 : 0 occurred (0.18 ± 0.03, not shown), whereas C_18 : 0 $/$ C_18 : 1 broadly decreased from 30.3 to 25.9 ka, then increased towards the top of the deposit (Fig. 4d). C_16 : 0 $/$ C_14 : 0 values were generally low (∼ 1.75) throughout the deposit but peaked at 29.1–28.8, 27.6, and 25.5 ka. C_18 : 0 $/$ C_14 : 0, broadly increased from 30.3–27.4 ka (Units III and II), decreased between 27.4 and 25.9 ka (lower Unit I), and then increased to reach a maximum at the top of the deposit (24.3 ka).

The major pigments in WMM7 were P410 and P435 (Fig. 4h), consistent with chlorophyll derivatives (chlorins; Harris and Maxwell, 1995) and some carotenoid sources (Jeffrey et al., 1997). P665​​​​​​​, recording exclusively chlorins, was very low (3 orders of magnitude below P410, Fig. A5). P410 was elevated at 29.8–28.8 ka (Unit III) and showed a sustained maximum between 28.5 and 27.2 ka (extending across the Unit II–I boundary, Fig. 4h). P435 remained low from the base of the sequence until a small peak at 27.2 ka. A rapid decrease in P410 and P435 resulted in pigment minima at 26.3 ka. From 25.9–24.3 ka, P410 oscillated at intermediate values, and there is a slight increase in P435 (Fig. 4h).

3.4 Stable isotope composition

Bulk δ¹³C ranged from −27.0 ‰ to −29.9 ‰. Bulk δ¹⁵N ranged between 11.1 ‰ and 12.6 ‰. There was a positive correlation between δ¹⁵N_bulk and δ¹³C_bulk between 30.3 and 25.9 ka, but they were negatively correlated after 25.9 ka (Fig. 5a, b). The two major fatty acids were offset from δ¹³C_bulk by up to −5.8 ‰ (−3.3 ± 1.4 ‰, C_14 : 0) and up to −3.9 ‰ (−1.3 ± 1.6 ‰, C_16 : 0; Fig. 5c). The most enriched values were recorded by δ¹³C_18 : 1 (−26.0 ± 1.2 ‰), which is offset to δ¹³C_18 : 0 by +1.27 ± 1.3 ‰ and to δ¹³C_bulk by +2.1 ± 1.6 ‰. Although lipid δ¹³C is usually depleted relative to δ¹³C_bulk as a result of fractionation during lipid synthesis (e.g. Rau et al., 1992), the dominance of lipids in WMM7 resulted in δ¹³C_bulk reflecting the combined influences of the main lipid components.

Figure 5 Stable isotope ratios in organic matter from stomach-oil deposit WMM7. The age model uses Holocene “no sea ice” ΔR (Table 1). (a) Bulk δ¹⁵N; (b) bulk δ¹³C; (c) fatty acid δ¹³C and bulk δ¹³C. The standard deviations of replicates are shown: where no error bar is visible the standard variation is smaller than the width of the symbol. Grey shaded data points are those for which the standard variation of the repeat measurements indicated poor reproducibility (> 1 ‰). (d) C_18 : 0 $/$ C_14 : 0 fatty acid ratio as an indication of relative krill contribution, from Fig. 4; (e) cluster boundaries from Fig. 4.

Fatty acid δ¹³C showed several short-term oscillations (Fig. 5c). Both δ¹³C_16 : 0 and δ¹³C_14 : 0 decreased from 29.4–28.4 ka (Unit III, lower Unit II), whereas δ¹³C_18 : 0 increased and δ¹³C_18 : 1 showed no clear trend. All fatty acid δ¹³C values were low at 28.4 ka (Unit II). The δ¹³C_14 : 0 remained low from 28.4–27.4 ka (Unit II), whereas δ¹³C_16 : 0 and δ¹³C_18 : 0 increased until 26.7 and 27.1 ka, respectively (Unit I). Unit I showed highly variable fatty acid δ¹³C. Minima in δ¹³C_14 : 0 and δ¹³C_18 : 0 were recorded at 26.7 ka and subsequently in δ¹³C_16 : 0 and δ¹³C_18 : 1 at 26.2 ka. All fatty acid δ¹³C increased to 25.4 ka and showed a slight decrease at the top of the sample (25.0 ka). The relationships between the four fatty acid δ¹³C changed through time: from 29.4–28.4 ka there was a close relationship between δ¹³C_14 : 0 and δ¹³C_16 : 0. In contrast, from 28.4 ka the δ¹³C_14 : 0 correlated with δ¹³C_18 : 0, and δ¹³C_16 : 0 correlated with δ¹³C_18 : 1.

Cluster analysis on the lipids, pigments, and stable isotope data identified three zones (Figs. 4j, A3). The lower zones aligned with Units III and II, but the upper zone boundary occurred at ∼ 27.0 ka, slightly above the Unit II–I and XRF zone boundaries and likely driven by the sustained maxima in P410 (Fig. 4h).

4 Discussion

Stomach-oil deposit WMM7 is characterized by biochemical variability, which can be attributed to changing biogenic composition through time. In this section, we first evaluate the likely sources of this variation, infer changes in snow petrel diet through time, and investigate how the sea-ice environment may have varied. By doing this we aim to provide greater insight into both the biochemistry of the stomach-oil deposit and critically evaluate its use as a palaeoenvironmental archive.

4.1 Contributions to stomach-oil biochemistry

The stomach oil regurgitated by snow petrels at nesting sites provides an integrated archive of diet during the most recent foraging trip in the summer breeding season. These trips may be for several days (Barbraud et al., 1999). During this time, the snow petrel concentrates the lipid components of their prey into stomach oil; this provides an energy-rich food source for their chicks (Warham, 1977; Watts and Warham, 1976) and can be spat in defence against predators or against other snow petrels in disputes over access to suitable crevices for nesting. Unlike muscle, feather, or adipose tissues, which involve biosynthesis (e.g. Rau et al., 1992), the basic biochemical composition of the prey is not altered during stomach-oil formation (Watts and Warham, 1976). Thus, the stomach-oil deposits provide a window into the biochemistry of the prey consumed by snow petrels in open waters at the margins of, or within, the sea-ice zone.

We have considered non-diet influences on the stomach-oil deposit geochemistry. A contribution of snow petrel guano to stomach-oil deposits has been indicated previously (Berg et al., 2019; Hiller et al., 1988). However, WMM7 does not show the negative relationship between C $/$ N and other elements (Cl, P, S) which would have been expected if there was a strong residual signature of guano contributions to the deposits (Berg et al., 2019). The C $/$ N ratios are consistent with particulate organic carbon (∼ 10) in modern sea ice (Henley et al., 2012) and fish collected in the Weddell and Lazarev seas (∼ 3–11) (Rau et al., 1992; Friedrich and Hagen, 1994). We suggest that any contribution of guano to WMM7 is lower than for other published sequences derived from DML stomach-oil deposits (Berg et al., 2019), perhaps due to local topography (including aspect and microclimate), the distance from the nest that the sample was taken, or the history of occupation.

We did not observe formation of K- and C-bearing phosphates by weathering (Berg et al., 2019); rather, the highest P counts were recorded in Unit II where Ca counts were lower and there was no change in K (Fig. A1). Progressive oxidation of organic matter was not evident: no systematic trends in labile pigment abundances (Fig. 4h) or fatty acid δ¹³C were observed (Fig. 5). Furthermore, bulk δ¹³C values are negatively correlated with C $/$ N, as observed today in Antarctic fish (and other predators), reflecting the variable contributions of ¹³C-depleted and low-N lipids to the bulk signal (Rau et al., 1992; Friedrich and Hagen, 1994; Cherel et al., 2011). Again, this difference in weathering signal between deposits may reflect local factors at the nesting sites (discussed above).

As observed in other glacial-stage stomach-oil deposits (Berg et al., 2019), WMM7 differs from some modern snow petrel stomach oils (e.g. Warham et al., 1976) by substantially lower contributions of C_18 : 1 relative to C_16 : 0 and C_14 : 0. The lower contributions of polyunsaturated (< 1 %) and monounsaturated fatty acids (e.g. C_18 : 1 9 ± 2 %, C_16 : 1 2.6 ± 0.3 %) than in fresh stomach oil of other procellariiform seabirds (Connan et al., 2007) including snow petrels (Warham et al., 1976; Watts and Warham, 1976) could indicate post-depositional oxidation (Berg et al., 2019). However, fatty acid contributions similar to WMM7 have been recorded in a late Holocene DML stomach-oil deposit (Aiello et al., 2011), in prey of snow petrels (Cripps et al., 1999), and in stomach oils from other Procellariiformes (Wang et al., 2007). As previously noted (Berg et al., 2019), the samples of fresh stomach oils (Warham et al., 1976; Watts and Warham, 1976) were from snow petrels foraging in the Ross Sea, where prey availability, and hence stomach-oil biochemistry, may also be different. We therefore suggest that the fatty acid signatures in WMM7 primarily signal a dietary intake rather than variable preservation.

Analyses of stomach contents and observations at sea indicate that snow petrel diet is dominated by fish (12 %–95 %) and krill (35 %–71 %) (e.g. Ridoux and Offredo, 1989; Ainley et al., 1984; Falla, 1937; Fijn et al., 2012), including the Antarctic silverfish Pleuragramma antarctica (family: Nototheniidae) and Antarctic lanternfish Electrona antarctica (family: Myctophidae) (Ainley et al., 2006, 1984). Variable contributions from squid have also been recorded (Ridoux and Offredo, 1989; Ainley et al., 1984; Bierman and Voous, 1950). The dominant fatty acids in WMM7 are consistent with the main snow petrel prey (Table 3): Antarctic krill (Euphausia superba, high abundances of C_14 : 0, C_16 : 0, sometimes C_18 : 1; Cripps et al., 1999; Raclot et al., 1998), squid (dominated by C_16 : 0, plus longer-chain fatty acids C_20 : 5 and C_22 : 6; Raclot et al., 1998), and both notothenoid and myctophid fish (high concentrations of C${}_{\mathrm{18}:\mathrm{1}(n-\mathrm{9})}$, C_16 : 0, and several monounsaturated and polyunsaturated C₂₀ and C₂₂ fatty acids; Imber, 1976; Raclot et al., 1998; Mayzaud et al., 2011). Although only recorded as a minor (∼ 2 %) contributor to snow petrel diet (Ridoux and Offredo, 1989), in coastal waters ice krill (Euphausia crystallorophias) is found in very high densities (La et al., 2015), characterized by high abundances of the C_{18 : 1ω9} and C_16 : 0 fatty acids as well as C_14 : 0 and C_16 : 0 alcohols (Ju and Harvey, 2004; Bottino, 1975). The C_16 : 0 fatty acid thus has a mixed origin from krill, fish, and squid in contrast to C_14 : 0 (Antarctic krill) and C_18 : x (fish, ice krill). The offsets in δ¹³C for C_14 : 0, C_16 : 0 and C_18 : x (and their different trends through time) confirm that they have different sources in WMM7 (Fig. 5). The offsets between individual fatty acids exceeds the ∼ 1 ‰ enrichment which might be expected through an increase in trophic level (Post, 2002), suggesting an alternative control through temporal or spatial variations in foraging habitat of the snow petrels and/or their prey (Sect. 4.3).

Our interpretation of C_14 : 0 as an indicator of Antarctic krill inputs to snow petrel diet is broadly supported by elevated Cu (measured as Cu $/$ Ti) (Rainbow, 1989; Palmer Locarnini and Presley, 1995; Liu et al., 2013) when C_14 : 0 is also relatively high (by C_16 : 0 $/$ C_14 : 0 and C_18 : 0 $/$ C_14 : 0) (Fig. 4). Thus, we identify two independent measures of the relative contributions of krill which can be applied to WMM7 to investigate broad-scale changes in snow petrel diet through time.

4.2 Snow petrel diet 30.3–24.3 ka

The variations in biogenic contributions to stomach-oil deposit WMM7 reveal changing snow petrel diet at decadal (XRF data), centennial, and millennial timescales for the early part of Marine Isotope Stage (MIS) 2, which includes the start of the LGM (30.3–24.3 ka). We note that although each stomach-oil regurgitation provides a snapshot of snow petrel diet during the summer breeding season (November to February), prey biochemistry (and thus stomach-oil composition) may also provide information over varying timescales, reflecting life cycle and tissue turnover rates of the fish, krill and squid consumed by the snow petrel.

Our multi-proxy analysis confirmed that three zones could be identified. Between 30.3 and 28.7 ka (∼ Unit III) elevated Cu $/$ Ti and C_14 : 0 contributions (low C_16 : 0 $/$ C_14 : 0 and C_18 : 0 $/$ C_14 : 0) identified Antarctic krill as an important component of snow petrel diet, but likely decreasing through time. Although offset by ∼ 2.0 ‰, supporting different sources, both δ¹³C_14 : 0 and δ¹³C_16 : 0 decreased across this time interval, suggesting either a common environmental (prey habitat) signal or that Antarctic krill contributed to both C_14 : 0 and C_16 : 0. In contrast, the much higher (and increasing) δ¹³C_18 : 0 and δ¹³C_18 : 1 confirmed a different prey source and/or habitat signal contributing to the C₁₈ fatty acids. Thus, the prey other than Antarctic krill, which were most likely to be fish, occupied a different foraging habitat (see Sect. 4.3).

At 28.7 ka (Unit II) there was a shift in the dominant contributions to the snow petrel stomach oils. Low Cu $/$ Ti indicate a prolonged (∼ 1100 years) interval wherein Antarctic krill was not a major component of snow petrel diet (Figs. 3, 4), supported by increasing C_16 : 0 $/$ C_14 : 0 and C_18 : 0 $/$ C_14 : 0. Reduced minerogenic inputs (Fig. 3f) and increased marine organic matter (Fig. 4h, i) characterize Unit II: low PC1, elevated S, Br (e.g. Leri et al., 2010), chlorins, and C_18 : 1 (Jónasdóttir, 2019; Cripps et al., 1999). Although this biochemical signature suggests a dominant phytoplankton fingerprint in snow petrel diet between 28.5 and 27.2 ka, phytoplankton are too small to be consumed directly by seabirds. The elevated chlorins (Fig. 4h) seem more likely to indicate that the snow petrels consumed prey with a phytoplankton-dominated diet. Decreasing δ¹⁵N_bulk in Unit II may also reflect a shorter food chain, although the shift is small (< 1 ‰, Fig. 5). Preservation of intact phytoplankton in stomach oils is also feasible through secondary ingestion (i.e. in the stomachs of snow petrel prey), since undigested algal cells have been reported in penguin guano (Mychra and Tatur, 1991) and diatoms have been isolated in P. antarctica larvae (Vallet et al., 2011) and other stomach-oil deposits (Berg et al., 2019).

The similarity in the trends between δ¹³C_16 : 0 and δ¹³C_18 : 1 (Fig. 5c) confirms an increased importance of prey with C₁₈ fatty acids in their tissues between 28.7 and 27.2 ka, consistent with incorporation of a phytoplankton signal in predator tissues, through their consumption of copepods, squid, or fish, which can occur with minimal alteration (e.g. Lee et al., 1971). Elevated (∼ 20 %) C_18 : 1 through Units I and II supports an increased contribution of fish to the snow petrel diet between 28.7 and 27.4 ka, but identifying the particular fish species is more challenging. Both P. antarctica (Mayzaud et al., 2011) and several species of myctophid fish including E. antarctica (Imber, 1976; Raclot et al., 1998) contain C_18 : 1, and their δ¹³C and δ¹⁵N values suggest that P. antarctica and E. antarctica have overlapping trophic niches (Rau et al., 1992). High abundances of C_16 : 0 in myctophids (Mayzaud et al., 2011) may account for the similarity in δ¹³C_16 : 0 and δ¹³C_18 : 1 in Unit II. Although the main myctophid prey today are euphausiids (including krill; Saunders et al., 2019), this diet is not consistent with the low Cu $/$ Ti and low C_14 : 0 contributions to WMM7 in Unit II, nor with the > 4 ‰ offset between δ¹³C_14 : 0 and δ¹³C_18 : 1. Alternatively, an increased contribution from nototheniid fish (e.g. P. antarctica), the diet of which includes copepods as well as krill, might explain the apparently elevated phytoplankton signature and the reduced contribution of Antarctic krill to snow petrel diet throughout Unit II.

From 27.4–25.9 ka (lower Unit I), Antarctic krill returned as an important component (high Cu $/$ Ti, higher C_14 : 0 contributions), although the positive relationship between δ¹³C_16 : 0 and δ¹³C_18 : 1 suggests that fish was also important in the diet. Reduced input of marine organic matter was apparent in Br, S, and chlorins (Figs. 3 and 4). Fatty acid δ¹³C showed higher variability than in Units II and III. From 25.2 ka to the top of the deposit (24.3 ka; Upper Unit I), increasing minerogenic contributions may explain the simultaneous decrease in Cu $/$ Ti, Br $/$ Ti, and S $/$ Ti, which contrasts with the anti-phase behaviour between Cu and Br or S discussed previously, and the Cu $/$ Ti minima, whilst C_14 : 0 inputs increase slightly. Fatty acid δ¹³C values also show concomitant peaks and lows from 25.5–25.0 ka, suggesting there may be a strong baseline (i.e. primary producer) control over the stable isotope values.

We note that the changes in relative contributions of Antarctic krill and fish in the diet that we inferred for WMM7 are not reflected in the range of δ¹⁵N_bulk (∼ 1.2 ‰), which lies below modern trophic-level offsets of > 2 ‰ (e.g. Seyboth et al., 2018) including some krill–fish offsets of ∼ 5 ‰–7 ‰ in the Southern Ocean (e.g. Polito et al., 2011), and suggest minimal changes to snow petrel trophic level between 30.3 and 24.3 ka. However, δ¹⁵N_bulk reflects δ¹⁵N in prey tissues, which can themselves vary over space and time (Quillfeldt and Masello, 2020) so that trophic values can only be calculated if baseline (primary producer) δ¹⁵N is known (Post, 2002). Since baseline δ¹⁵N is poorly constrained for our interval of study and since snow petrel prey integrate δ¹⁵N across different temporal and spatial scales, it seems likely that the small range of δ¹⁵N_bulk in WMM7 reflects either baseline δ¹⁵N_bulk or mixing of multiple signals over time, which we are unable to disentangle with our analyses.

In summary, the stomach-oil deposit WMM7 indicates several changes to snow petrel diet between 30.3 and 24.3 ka. Antarctic krill and fish provided important contributions between 30.3–28.7 ka (Unit III) and 27.4–24.3 ka (Unit I). Between 28.7 and 27.4 ka (Unit II) a prolonged interval of low Antarctic krill contributions suggests a shift towards a more fish-dominated diet. The strong phytoplankton fingerprint during this time also suggests that the fish themselves had a reduced krill contribution to their diet, highlighting a substantial change to the food chain in the foraging habitat of the snow petrels. In the following section, we consider the likely palaeoenvironmental drivers and implications of these signals.

4.3 Evolution of sea-ice environments (30.3–24.3 ka)

Seabird diets vary over time and space, depending on the availability of prey species, in turn reflecting environmental conditions within their foraging range (Mills et al., 2020, 2021). For example, it has been proposed that myctophid fish become more important to seabird diet when krill availability is low, during winter when sea ice is more extensive, and/or in response to temporal and spatial changes in summer production (Nicol, 2006; Watanuki and Thiebot, 2018). Prey biochemistry may also reveal environmental information: δ¹³C decreases with increasing latitude (Francois et al., 1993; Trull and Armand, 2001) or between coastal and offshore waters (Trull and Armand, 2001), whereas δ¹³C increases within the sea ice during spring melt (Dunbar and Leventer, 1992). These differences are reflected in the δ¹³C of crustaceans, fish, and squid in these habitats and in turn in their predators including seabirds (Cherel et al., 2011; Jaeger and Cherel, 2011; Delord et al., 2016; Phillips et al., 2009).

The presence of multiple stomach-oil deposits in DML during MIS 3 and the early LGM (Berg et al., 2019) confirms that there must have been accessible prey within range of the breeding snow petrels at the Untersee during the spring and summer. WMM7 has a very stable accumulation rate (Fig. 2), and we found no evidence for hiatuses, suggesting that the site was not affected by long periods of nest-site abandonment related to unfavourable climatic conditions (e.g. snow cover over the site, perennial sea ice preventing access to prey within foraging range) or local factors (e.g. burial by rock) (Olivier and Wotherspoon, 2006; Olivier et al., 2005). We consider it unlikely that the snow petrels were foraging at an oscillating spring–summer sea-ice margin because WMM7 deposition coincides with an expanding then maximum summer sea-ice extent in the Atlantic sector of the Southern Ocean between 29 and 22 ka (Fig. 6) (Collins et al., 2012; Gersonde et al., 2005, 2003; Allen et al., 2011; Fischer et al., 2007; Xiao et al., 2016; Goto-Azuma et al., 2019), which places the summer sea-ice margin at 55–60^∘ S, > 1500–2000 km away from WMM7 and far beyond the typical maximum foraging range for this species (Fig. 1).

Figure 6 Variations in snow petrel diet across the MIS 3–2 transition, including the interval of maximum summer sea-ice extent (SSI-max) in the Scotia Sea from Allen et al. (2011) and plotted in Fig. 1. (a) Air temperature recorded by δD in an EPICA–DML ice core (Jouzel et al., 2007); (b) atmospheric CO₂ from a West Antarctic Ice Sheet (WAIS) ice core (Bauska et al., 2021); (c) sea-salt Na flux from an EPICA–DML ice core (Fischer et al., 2007); (d) biogenic sulfate flux to the EPICA–DML ice core, generated by correcting for terrestrial contributions to non-sea-salt Na flux (Goto-Azuma et al., 2019); (e) WMM7 chlorin pigment abundance plotted using the glacial-stage (enhanced sea ice) ΔR from Table 1; (f) WMM7 Cu $/$ Ti ratio plotted using the glacial-stage (enhanced sea ice) ΔR from Table 1.; (g) WMM7 chlorin pigment abundance plotted using the Holocene (low or absent sea ice) ΔR from Table 1, as shown in Fig. 4; (h) WMM7 Cu $/$ Ti ratios plotted using the Holocene (low or absent sea ice) ΔR from Table 1, as shown in Fig. 3. As discussed in the text, we prefer to apply the Holocene ΔR calibrations of (g, h) and interpret high chlorin values (g) as indicative of an enhanced phytoplankton and/or reduced krill signal, whereas high Cu $/$ Ti ratios (h) are interpreted as evidence of enhanced krill inputs. As discussed in the text, we infer the loss of krill from the snow petrel diet at ∼ 27–28 ka to represent polynyas opening over the continental shelf. (i, j) Cu $/$ Ti signals in other DML stomach-oil deposits from analysis in Berg et al. (2019) recalibrated using the Holocene ΔR from Table 1.

Today, snow petrel foraging in high (> 80 %) sea-ice concentrations occurs in the margins of polynyas (Ainley et al., 1984), which are regions of open water within the winter or multiyear sea-ice pack critical for enabling high levels of primary productivity, carbon cycling, and air–sea gas exchange (Arrigo et al., 2008; Sherrell et al., 2015). Several marine sites in the Weddell and Lazarev seas have indicated sporadic polynya formation during MIS 2 and 3 (Fig. 1) (Sprenk et al., 2014; Smith et al., 2010; Thatje et al., 2008; Weber et al., 1994). Deposit WMM7 confirms previous proposals that polynyas were important for supporting snow petrel colonies through the last glacial (Thatje et al., 2008; Berg et al., 2019). We interpret our new records of changes in snow petrel diet in MIS 2 to have occurred as a result of changes to conditions within, and/or the locations of, polynyas offshore of DML. Unlike today, when polynyas are largely winter–spring phenomena which impact spring–summer productivity during or after sea-ice melt as “post-polynyas” (Arrigo and van Dijken, 2003), summer sea surface temperatures below freezing in the Weddell Sea (Gersonde et al., 2005) would have supported more extensive sea-ice (and polynya) formation in summer during the last glacial stage (Fischer et al., 2007).

Diet studies indicate that snow petrels consume a mixed diet of krill, squid, and myctophid fish when feeding beyond the continental shelf, whereas fish are most important in shelf waters, at least in the Ross Sea (e.g. Ainley et al., 1992; Fijn et al., 2012;​​​​​​​ Ridoux and Offredo, 1989; Ainley et al., 1984; Falla, 1937). A similar pattern in Antarctic petrels (Thalassoica antarctica) has been attributed to a lack of Antarctic krill inshore during the breeding season (Nicol, 1993). Although adult krill and larvae over-winter within the sea ice, post-larval Antarctic krill are mostly oceanic (Atkinson et al., 2008). During the austral spring and summer, corresponding to the breeding season of the snow petrels, adult Antarctic krill move to deeper waters for egg development (Nicol, 2006), which may account for the observed increase in krill in snow petrel diets when they forage beyond the continental shelf. Given the importance of deep waters for the Antarctic krill life cycle, we infer that the mixed diet of krill and fish in Units I and III represents snow petrels foraging in polynyas located beyond the continental shelf. In contrast, the observed shift in fatty acid and element profiles in Unit II suggests that fish became more important to snow petrel diet, suggesting that polynyas had opened up over the continental shelf between 28.7 and 27.4 ka. We hypothesize that these shifts in foraging habitat reflect changes in sea-ice conditions by influencing either prey distributions or access to surface waters for feeding.

We discount a latitudinal control over δ¹³C in WMM7 because the seasonal and spatial ranges of δ¹³C within sea ice can exceed the Southern Ocean latitudinal δ¹³C gradient (Dunbar and Leventer, 1992; Kennedy et al., 2002) due to changing CO_2(aq) supply and productivity (Henley et al., 2012; Kennedy et al., 2002). Organic matter in the sea-ice zone tends to have elevated δ¹³C compared to open waters beyond the ice edge (Henley et al., 2012; Kennedy et al., 2002), but δ¹³C_14 : 0 production in open waters for WMM7 is not consistent with the expanded sea-ice extent of MIS 3–2 outlined above. Rather, low δ¹³C_14 : 0 is consistent with (krill) production in the winter–spring sea ice (Dunbar and Leventer, 1992; Henley et al., 2012) and/or more oceanic (rather than coastal) habitats (Cherel et al., 2011), with ¹³C-enriched C₁₈ fatty acids in WMM7 reflecting primary production and fish tissues from summer and/or more coastal habitats. No long-term trends are recorded in δ¹³C_18 : 0 or δ¹³C_18 : 1, and all four of the major fatty acids have highly variable δ¹³C (ranges > 2 ‰) after 27.0 ka (Unit I). In contrast, a long-term decrease in δ¹³C_14 : 0 between 29.4 and 28.4 ka (Units III and II) occurs alongside declining Antarctic krill contributions (Cu $/$ Ti and C_18 : 0 $/$ C_14 : 0). The cause of the long-term decline in δ¹³C_14 : 0 is unclear but could reflect prolonged winter (low δ¹³C) sea ice, which leads to ¹³C depletion (Dunbar and Leventer, 1992) and also to restricted access to krill (Watanabe et al., 2020), even though increased winter sea ice is important for Antarctic krill recruitment (Atkinson et al., 2004; Jaeger and Cherel, 2011). Alternatively, declining δ¹³C_14 : 0 could reflect enhanced upwelling of CO_2(aq)-rich but low-δ¹³C circumpolar deep water (Bostock et al., 2004; Menviel et al., 2015) during the season of Antarctic krill production. The decreasing δ¹³C_14 : 0 and krill contributions from 30.3 ka culminate in the shift to a diet rich in fish by 28.7 ka (Unit II), which we infer to represent coastal polynya development, reducing krill recruitment in shallower waters and/or reducing snow petrel access to krill with more widespread spring sea ice. Elevated phytoplankton contributions between 28.5 and 27.2 ka in WMM7 are also consistent with the highest primary productivity and/or low grazing observed today in shallow (continental shelf) waters, including coastal polynyas (Delmont et al., 2014; Karnovsky et al., 2007; Arrigo and van Dijken, 2003).

The cause of the millennial-scale switching between foraging in coastal and open-ocean polynyas is not clear but could reflect ocean temperatures, wind speeds and directions (Morales Maqueda et al., 2004), and availability of continental shelf waters. Colder Southern Ocean waters during glacial stages (Gersonde et al., 2005) are expected to limit formation of the open-ocean polynyas which rely on ocean-driven sea-ice melt (Ferrari et al., 2014), and reduced vertical mixing may have resulted from sea-ice-driven stratification in the upper water column (Crosta and Shemesh, 2002). However, upwelling driven by wind and/or terraces on the continental slope may explain LGM polynya formation at site PS1506 (Thatje et al., 2008) and upstream of site PS1795 (Sprenk et al., 2014), within the foraging range of snow petrels breeding at the Untersee (Fig. 1). Enhanced wind speeds during glacial stages are also expected to have encouraged sea-ice break-up and coastal (ice-sheet- or ice-shelf-proximal) polynya formation (Morales Maqueda and Rahmstorf, 2002; Sprenk et al., 2014; Smith et al., 2010; Weber et al., 1994), as observed today (Morales Maqueda et al., 2004).

Unlike other parts of the Antarctic ice sheet (Bentley et al., 2014), access of snow petrels to continental shelf waters during MIS 2 may also have been facilitated by relatively minor changes to vertical and lateral ice-sheet extent in DML over the last 100 000 years (Hillenbrand et al., 2014; Mackintosh et al., 2014). The exact position of the LGM ice-sheet limit on the DML continental shelf is not well constrained (Fig. 1) (Hillenbrand et al., 2014), but our data from WMM7 suggest that at least part of the continental shelf was free of glacial ice, enabling the development of coastal polynyas within MIS 2. Times of ice-sheet advance may then have reduced the area of available continental shelf and increased wind-driven polynya formation over the continental slope, as suggested by Weber et al. (1994) and Sprenk et al. (2014), which would account for the increased contribution of Antarctic krill to snow petrel diet. The relatively restricted extent of the continental shelf at DML (< 100 km) may make coastal polynya formation in this region particularly sensitive to ice-sheet advance.

The sustained presence of polynyas through MIS 3–2, as recorded in WMM7 and in other DML sequences (Fig. 6) (Berg et al., 2019), may also offer an explanation for the small amplitude of millennial-scale oscillations in sea-salt Na to the EPICA ice core and the weak relationship with Antarctic temperatures (Fischer et al., 2007). Although we show that polynya properties changed through time (Fig. 6), the open waters within the sea-ice pack would have provided a continuous supply of sea-salt Na to the ice core site. Polynyas may have also affected the strength of the sea ice–climate feedbacks during MIS 2: introducing only 2 %–8 % open waters into the LGM sea-ice pack (compared to 10 %–20 % for winter today) reduces the Southern Ocean contribution to the LGM CO₂ drawdown from ∼ 80 % to 15 %–50 % via enhanced ocean–atmosphere CO₂ transfer (Morales Maqueda and Rahmstorf, 2002). In contrast, increasing brine formation transfers dense water and carbon to the deep ocean (Bouttes et al., 2011) and could have been enhanced by polynya formation (Paillard and Parennin, 2004). Brine formation over the continental shelves or at the ice-sheet margin has been proposed as conducive to formation of dense glacial Antarctic Bottom Water (AABW) and the associated deep-ocean storage of CO₂ (Paillard and Parrenin, 2004; Adkins, 2013; Adkins et al., 2002). It is currently difficult to evaluate the relationship between the proposed variability in polynya positions and the millennial-scale oscillations in atmospheric CO₂ (Fig. 6), in part because it is unclear whether the variations in surface ocean productivity observed in the stomach-oil deposits are related to changes in the efficiency of the biological pump and CO₂ drawdown (e.g. for Unit II with high chlorin inputs). Furthermore, the relative impact of polynyas compared to other Southern Ocean carbon cycle processes is unclear, given that a combination of brine formation related to sea-ice growth, changes in deep-ocean stratification, and iron fertilization of subantarctic waters is invoked to account for the observed CO₂ drawdown (e.g. Bouttes et al., 2011; Marzocchi and Jansen, 2019; Sigman et al., 2021). The relative impact of the polynyas between winter and summer seasons during the last glacial stage is also uncertain, since the large changes in sea-ice extent (Fig. 1) will likely also have affected the air–sea gas exchange on seasonal timescales. Our age model uncertainties also limit confident correlation between WMM7 and the ice core CO₂ record, so further testing is required to explore whether polynya development along the DML coastline contributed to observed changes in atmospheric CO₂.

5 Conclusions

Here, we present a multi-proxy analysis of stomach-oil deposits of snow petrels from Dronning Maud Land, Antarctica. We show that variation in trace metals and fatty acid distributions strongly suggests changes in snow petrel diet between 30.3 and 24.3 ka. We show that, as today, a mixed diet of Antarctic krill and fish characterizes much of the record. However, between 28.7 and 27.4 cal kyr BP signals of Antarctic krill in the diet almost disappear. By linking dietary signals in the deposits to modern feeding habits and foraging ranges, we highlight the presence of open water (polynyas) within more extensive summer sea-ice cover during MIS 2. The reduced contribution of Antarctic krill to the snow petrel diet between 28.7 and 27.4 ka suggests restriction of polynyas to the continental shelf, limiting krill recruitment or access to waters where krill were present. Our results show that extensive, thick, and multiyear sea ice was not always present close to the continent during MIS 2. These results add to a growing body of evidence which shows that seasonal changes in sea-ice extent and the presence of polynyas were likely important drivers of sea ice–climate feedbacks during glacial stages, whilst also highlighting the potential of stomach-oil deposits as a palaeoenvironmental archive of Southern Ocean conditions.

Appendix A

Figure A1 Consistency of XRF scanning signals in deposit WMM7. The XRF data presented in Fig. 3 were recovered from a central slab of WMM7, where we identified the longest sequence of horizontal laminae (Fig. 3). Here we show additional scans to investigate the internal consistency of the signals. Panels (a)–(e) correspond to lines a–e in panel (f). Scans were undertaken to the left (b) and right (d) of the centreline (c), and then two further cuts were used to investigate signals perpendicular to the cut surface (a, e). Away from the central line the signals are less strong and sometimes intermittent, but the laminae also dip sharply away from the centre (f), suggesting more heterogenous accumulation, perhaps due to distance from the nest and/or nest morphology.

Figure A2 Comparison of cluster analysis results for original XRF data, resampled XRF data (shown in Fig. 3), and organic indicator data (shown in Fig. 4). XRF data between 0–10 and 155–160 mm were removed before analysis. The left panel shows the original Cu $/$ Ti data (grey line) with ∼ 100-year smoothing (orange line, as in Fig. 3), alongside resampled Cu $/$ Ti data (red squares) used for cluster analysis. The right panel shows the P410 pigment absorbance signal, which may account for the offset in the O1 $/$ O2 cluster boundary compared to $\mathrm{1}/\mathrm{2}$ boundary in other analyses.

Figure A3 Element contributions to stomach-oil deposit WMM7, recorded by ITRAX XRF scanning (Croudace et al., 2006). Note that the amplitude of the signal decreases from left to right. Potential minerogenic (Ti, K, Al; brown), biogenic (Br, Zn, S, As; green), and mixed (Fe, Ca, Cu, Si; black) sources are noted. Cluster analysis was performed after removing data from 0–10 and 155–160 mm where several elements record values of zero.

Figure A4 Effect of normalization for minerogenic inputs on the Fe and Cu signals in WMM7. Due to low or absent contributions of Al (Fig. A3), data were normalized to Ti. Fe $/$ Ti and Cu $/$ Ti are shown here for comparison as the original ratios and on a log scale. Fig. 3 presents all normalized XRF data on a linear scale, excluding 0–10 and 155–160 mm, where very low counts of one or more element lead to anomalous values.

Figure A5 Additional indicators of organic matter composition of WMM7. The %N and %TOC were used to generate the C $/$ N ratios presented in Fig. 4. The dominant pigments (P410, P435) are shown in comparison to the much lower contribution of P665, which has an exclusive chlorin origin but which is also 3 orders of magnitude lower in abundance than P410.