# Sensitivity of Holocene East Antarctic productivity to subdecadal variability set by sea ice

Katelyn Johnson, Robert Mckay, Johan Etourneau, Francisco Jiménez-Espejo, Anya Albot, Christina Riesselman, Nancy Bertler, Huw Horgan, Xavier Crosta, James Bendle, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Katelyn Johnson, Robert Mckay, Johan Etourneau, Francisco Jiménez-Espejo, Anya Albot, et al. Sensitivity of Holocene East Antarctic productivity to subdecadal variability set by sea ice. Nature Geoscience, 2021, 14 (10), pp.762-768. 10.1038/s41561-021-00816-y. hal-03529697

HAL Id: hal-03529697

https://cnrs.hal.science/hal-03529697v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Sensitivity of Holocene East Antarctic productivity to subdecadal variability set by sea ice

Katelyn M. Johnson<fnref n="1" />,²✉, Robert M. McKay², Johan Etourneau³,⁴,⁵, Francisco J. Jiménez-Espejo<fnref n="1" />,³,⁶, Anya Albot², Christina R. Riesselman<fnref n="1" />,⁷,<fnref n="8" />, Nancy A. N. Bertler<fnref n="1" />,⁲, Huw J. Horgan², Xavier Crosta<fnref n="1" />,⁵, James Bendle<fnref n="1" />,⁹, Kate E. Ashley⁹, Masako Yamane¹⁰, Yusuke Yokoyama¹¹, Stephen F. Pekar¹², Carlota Escutia¹³ and Robert B. Dunbar¹³

Antarctic sea-ice extent, primary productivity and ocean circulation represent interconnected systems that form important components of the global carbon cycle. Subdecadal to centennial-scale variability can influence the characteristics and interactions of these systems, but observational records are too short to evaluate the impacts of this variability over longer timescales. Here, we use a 170-m-long sediment core collected from Integrated Ocean Drilling Program Site U1357B, offshore Adélie Land, East Antarctica to disentangle the impacts of sea ice and subdecadal climate variability on phytoplankton bloom frequency over the last ~11,400 years. We apply X-ray computed tomography, Ice Proxy for the Southern Ocean with 25 carbon atoms, diatom, physical property and geochemical analyses to the core, which contains an annually resolved, continuously laminated archive of phytoplankton bloom events. Bloom events occurred annually to biennially through most of the Holocene, but became less frequent (~2-7 years) at ~4.5 ka when coastal sea ice intensified. We propose that coastal sea-ice intensification subdued annual sea-ice break-out, causing an increased sensitivity of sea-ice dynamics to subdecadal climate modes, leading to a subdecadal frequency of bloom events. Our data suggest that projected loss of coastal sea ice will impact the influence of subdecadal variability on Antarctic margin primary productivity, altering food webs and carbon-cycling processes at seasonal timescales.

Antarctica's marine margin is a complex biological and oceanographic system in which sea-ice growth, Antarctic Bottom Water (AABW) formation and high primary productivity act as a significant CO₂ sink and ventilate the Southern Ocean<fnref n="1" />,². High primary productivity occurs where nutrients are brought to the surface, including oceanographic fronts³, polynyas⁴, upwelling near the continental shelf break⁵ and the marginal ice zone (MIZ)⁶, all of which are influenced by Antarctic wind fields. High productivity and export events around Antarctica occur with changing insolation and stratification associated with sea-ice break-up⁷. Large-scale subdecadal climate modes, specifically the El Niño–Southern Oscillation (ENSO), Southern Annular Mode (SAM) and Indian Ocean Dipole (IOD), are known to affect sea ice<fnref n="8" />,⁹ and wind fields¹⁰ around Antarctica. The teleconnection between ENSO (which varies at 2–7 year periods) and Antarctic sea-ice variability is largely driven by wind changes resulting from hemispheric-scale sea level pressure and 500 mbar height anomalies<fnref n="8" />. This teleconnection can be amplified or dampened by other subdecadal climate modes such as the IOD and SAM⁹⁻¹². Collectively, these subdecadal climate modes alter meridional and zonal wind flows⁹,¹⁰ that regulate sea-ice break-out¹¹ at 2–7 year periods, thus influencing primary productivity in Antarctica¹³⁻¹⁵. Clarifying how the annual cycle and subdecadal-scale climate modes have impacted past Antarctic coastal systems will inform models used to project future system response¹⁶.

We investigate a 170 m sediment core recovered from the Adélie Basin (Integrated Ocean Drilling Program (IODP) Site U1357B)¹⁷ along the Wilkes Land margin of East Antarctica (Fig. 1). The Adélie Basin is a region of high primary productivity near the MIZ. It also lies beneath and downstream of several large polynya systems and the westward-flowing Antarctic Coastal Current. The drill site targeted a high-sedimentation (~1.5–2 cm per year) drift deposit (Adélie Drift) dominated by pelagic biogenic deposition. It provides an ultra-high-resolution Holocene record of climate and oceanographic variability adjacent to the Mertz Polynya system, one of the largest exporters of sea ice and AABW along the East Antarctic margin². Previously collected Antarctic cores have significantly lower sedimentation rates, and alternate between massive (bioturbated) and laminated diatom ooze¹⁸,¹⁹. They cannot resolve high-frequency change at subdecadal scales. However, U1357B is continuously laminated, and high sedimentation rates afford an unprecedented opportunity to assess subdecadal and annual changes at the Antarctic oceanic margin.

## An ultra-high-resolution record of marine biogenic blooms

The east-west elongated Adélie Drift deposit formed parallel to the wind-driven Antarctic Coastal Current²,²⁰. This current influences both surface and deep waters on the continental shelf²,²⁰. Consequently, the mass accumulation rate (MAR) (Methods) in this drift is thought to reflect current strength, and only partially

<footnote n="1">¹GNS Science, Lower Hutt, New Zealand. ²Antarctic Research Centre, Victoria University of Wellington, Wellington, New Zealand. ³Andaluz Institute of Earth Sciences, CSIC-University of Granada, Granada, Spain. ⁴EPHE, PSL University, Paris, France. ⁵UMR 5805 EPOC CNRS, University of Bordeaux, Bordeaux, France. ⁶Biogeochemistry Center, JAMSTEC, Yokosuka, Japan. ⁷Department of Geology, University of Otago, Dunedin, New Zealand.</footnote>

<footnote n="8">⁸Department of Marine Science, University of Otago, Dunedin, New Zealand. ⁹School of Geography, Earth, and Environmental Sciences, University of Birmingham, Birmingham, UK. ¹⁰Institute for Space-Earth Environmental Research (ISEE), Nagoya University, Nagoya, Japan. ¹¹Analytical Center for Environmental Study, Atmosphere and Ocean Research Institute, University of Tokyo, Kashiwa, Japan. ¹²School of Earth and Environmental Sciences, City University of New York - Queens College, New York, NY, USA. ¹³School of Earth, Energy, and Environmental Sciences, Stanford University, Stanford, CA, USA. ✉e-mail: k.johnson@gns.cri.nz</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 1** | Area of study and example bloom event after sea-ice break-out. **a**, Map<sup>46-48</sup> of Adélie Land and locations of U1357B (orange circle), MD03-2601 (black circle) and MD03-2597 (pink circle) sites and Dumont d'Urville (DDU) station (black square). Primary bathymetric features<sup>48</sup>, wind directions and current locations are indicated. X-Y (black) indicates approximate location of seismic profile<sup>17</sup> in Fig. 2. **b**, Moderate Resolution Imaging Spectroradiometer (MODIS) true colour (bands 1, 4 and 3) satellite imagery capturing sea-ice break-out on 12 December 2008. **c**, MODIS true colour (bands 1, 4 and 3) satellite imagery overlaid with chlorophyll-a concentrations<sup>49,50</sup> from phytoplankton bloom event on 15 January 2009, following sea-ice break-out. Site U1357B is indicated by orange circle in **b** and **c**. Antarctic polar stereographic projection (EPSG: 3031). MODIS true-colour satellite images from NASA Worldview.

reflects changes in biological productivity (Fig. 2). Comparison of the covarying siliciclastic (detrital) and biogenic MAR²¹ (Fig. 3c) suggests that detrital and biogenic sediments are advected to the site together, under the influence of wind-driven currents and focused into the Adélie Basin by shelf bathymetry (Fig. 2). Nearby core MD03-2601 (Fig. 1) shows sedimentation rates covarying with U1357B throughout the Holocene, indicating that the sediment advection process is a regional signal (Fig. 2 and Extended Data Fig. 1).

Iceberg rafted debris (IRD) is negligible (Extended Data Fig. 2), aside from the very bottom of the core (>168 m below sea floor (mbsf)) suggesting that direct sediment supply from icebergs is limited. The geometry and location of the drift are inconsistent with deposition as part of a glaciomarine fan system. Regional bathymetric highs are characterized by poorly sorted diamicts and

muddy sands²². Grain size frequency distributions in those settings indicate the partial winnowing of the <125 µm component by bottom currents²². However, detrital siliciclastic material in the bathymetric troughs, including the Adélie Drift deposit, are consistently <125 µm with a well-defined silt-fine sand mode (Extended Data Fig. 2). This is interpreted as representing suspension settling of winnowed muds derived from diamicts on adjacent highs, suggesting that the primary control on sedimentation is from current strength and sediment advection (Fig. 2). As suspended sediment winnowed from the banks is advected towards U1357, settling of sediment will occur where current slows as it passes over deep bathymetric troughs²². Therefore, U1357 represents a sediment trap and changes in sediment grain size are a function of Antarctic Coastal Current strength, with larger grain sizes transported during stronger currents²². This is supported by the covariance of sand
<figure><img src="image_5.png" /><figcaption>b Stronger winds, less coastal sea ice and stronger current</figcaption></figure>

<figure><img src="image_6.png" /></figure>

**Fig. 2** | Simplified sediment deposition model for U1357B. **a,** Simplified deposition model of Adélie Drift during weaker winds (katabatic/zonal; blue arrows), more sea ice and subsequent weaker coastal current (yellow arrows). Biogenic and winnowed terrigenous material are selectively deposited (white arrows) into drift as water slows over basin. Light and dark laminae indicated by brown and green lines. Pink line indicates approximate location of U1357. Relative strength of winds and currents indicated by arrow size. Characteristics of this mode are reduced grain size, reduced MAR, reduced laminae thickness and increased laminae per metre. **X,** Y marks seismic profile direction as seen in Fig. 1. **b,** same as **a,** but for stronger winds (katabatic/zonal), less sea ice and stronger current.

percent with MAR curves, whereby an increase in MAR and sand percent relate to increased current speed (Fig. 3c,e and Extended Data Fig. 3).

The site also traps biogenic material produced in the Dumont d'Urville polynya above¹, as well as biogenic material advected from the Mertz Polynya to the east. It is assumed that local biogenic material dominates¹⁷. Large phytoplankton bloom events along Antarctica's coastal margins are associated with a relatively fresh and stably stratified meltwater layer originating from seasonal sea-ice melt¹,⁶. Seasonal sea-ice break-up in this region is strongly affected by changes in katabatic and zonal wind intensity²³. Diatom analysis from MD03-2597 in the Adélie Drift and nearby MD03-2601 interprets light laminae as biogenic bloom events, occurring during spring sea-ice retreat, which are rapidly exported to the seafloor²⁴,²⁵.

To determine light laminae frequency, the top and bottom of each light lamina bloom event were determined manually using X-ray computed tomography (CT) images and compared with greyscale time series extracted from raw CT data (Methods, Extended Data

Fig. 4 and Fig. 3a). Comparisons of the greyscale curve with gamma ray attenuation (GRA) bulk density, natural gamma radiation (NGR), X-ray fluorescence (XRF) silica and XRF titanium indicate that changes in the CT profile are primarily associated with changes in laminae composition (Extended Data Fig. 5). This is assessed further by independent time series analysis of the greyscale and productivity indicators in the XRF data (Fig. 4). Ice Proxy for the Southern Ocean with 25 carbon atoms (IPSO₂₅) data from U1357B (Methods) and diatom analysis from MD03-2601 (refs. ²⁶,²⁷) are used to assess the influence of sea-ice conditions on sedimentology and bloom frequency. IPSO₂₅ is interpreted as a proxy for fast ice²⁸, sea ice anchored to features along the continental margin, whereas increases in *Fragilariopsis curta* relative abundance indicate cooler temperatures and later spring sea-ice melting²⁷.

## Subdecadal drivers of coastal Adélie Land bloom events

Applying the age model to the manual laminae counts, we find that annual to biennial frequencies dominate before 4.5 ka. Subdecadal periodicities (2–7 years) dominate after 4.5 ka and are superimposed on distinct variations occurring at centennial to millennial timescales (Fig. 3d). These annual to subdecadal frequencies, as well as the lower-frequency shifts, are upheld by evolutionary spectral analysis of the raw CT-scan greyscale data and inferred XRF productivity ratio Si/Ti (Methods and Fig. 4). Comparison of these data with the MAR curve, IPSO₂₅ and diatom assemblage data indicate distinct climate states in the Holocene, noted in other Antarctic records (Fig. 3)¹⁶,²⁹. From 4.5 ka, a baseline shift occurs in coastal sea-ice proxies and sand percent, which correspond with less-frequent bloom events relative to the overall record (Fig. 3 and Extended Data Fig. 6). These laminae events occur at frequencies that are consistent with modern-day ENSO frequencies of 2–7 years (Figs. 3 and 4). However, this relationship is interrupted between 0.8 and 1.8 ka, when IPSO₂₅ is reduced and bloom events are more frequent.

Although some records argue for a baseline shift in ENSO variability at 4–5 ka (refs. ³⁰,³¹), others suggest it has been highly variable for the past 7 ka (ref. ³²). Given that a consistent pattern is not yet recognized in Holocene ENSO records³¹,³³, we remain cautious about correlating Antarctic time series with low-latitude records of Holocene ENSO variance³¹,³³. Temporally limited proxy records of other subdecadal climate modes³⁴,³⁵ (for example, IOD and SAM), which amplify or dampen the ENSO response along Adélie Land⁹⁻¹², present a similar issue. This precludes us from making a direct comparison of how variations in the intensity of these subdecadal climate modes have impacted Adélie Land. However, there is no evidence that ENSO frequencies have shifted out of the 2–7 year band into the 1–2 year band³⁰⁻³³. SAM and IOD modulate the amplitude of ENSO influence on winds¹⁰ and coastal sea ice⁹,¹¹, not the frequency. Below, we investigate how the frequency of biological bloom events has shifted through the Holocene. We interpret our data primarily in the context of local productivity drivers. We also identify whether bloom frequencies are consistent with modern subdecadal climate modes, which regulate sea-ice break-out⁸⁻¹⁰ and induce bloom events¹³,¹⁵, or the annual seasonal cycle.

## Local deglaciation-influenced bloom events (11.4 to 8 ka)

Sediments in the lowermost ~0.7 m of U1357B are poorly sorted, with IRD visible in the CT images (Extended Data Fig. 7). Upcore, IRD is largely absent from the CT images and grain size frequency distributions (Extended Data Fig. 2). Bloom events occurred at an annual frequency around 11.4 ka, before trending towards a biennial period (5–7 laminae per 10 years) between 10.8 and 9 ka (Fig. 3d). Bloom frequency was highest at ~8.2 ka, with one or multiple events occurring annually.

Before ~8.2 ka, frequent occurrences of laminae peaks are attributed to freshwater pulses from the final phase of local East Antarctic
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Holocene proxy records in Adélie Land. **a**, Raw CT greyscale data. **b**, Raw XRF linescan data of productivity ratio Si/Ti. **c**, MAR from U1357B. Green is biogenic silica MAR; brown is terrigenous²¹. **d**, Laminae frequency per 10 years smoothed in a five-point moving mean; the bold curve is an rlowess smoothing, using a 5% span of the data. **e**, Sand percentage of the light laminae, which is representative of current speed; the bold curve is an rlowess smoothing using a 5% span of the data. **f**, $IPSO_{25}$ concentration from U1357B, a proxy for fast ice conditions. **g**, Percentage of *F. curta* from MD03-2601 (ref. ²⁷), a diatom species indicative of later spring sea-ice melt. Missing data in **a** and **b** represent intervals with no core recovery.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4 |** EHA of greyscale data and XRF Si/Ti productivity proxies.</figcaption></figure>

a, Greyscale data. b, XRF Si/Ti. Both plots are overlaid with curves of laminae frequency per 10 years in white and black (same as in Fig. 3d). Normalized power is similar across both proxies, showing a distinct shift to fewer bloom events and reduced productivity at 4.5 ka. Manual laminae counts binned at 10-year intervals are consistent with the EHA. The white curve is the 10-year binned record smoothed in a five-point moving mean, while the black curve is a lowess smoothing, using a 5% span of the unsmoothed 10-year binned record. The black boxes indicate intervals with no core recovery. The 2–7 year subdecadal climate mode band is indicated by vertical black dotted lines (that is, five laminae per 10 years is a 2-year frequency).

Ice Sheet retreat<sup>36,37</sup>. Deglacial reconstructions suggest a calving bay re-entrant pattern, whereby ice retreated first in the bathymetric troughs and later from the adjacent bathymetric highs<sup>36,37</sup>. This is supported by the decline in NGR and mean grain size, and gradual increased sorting of the detrital fraction upcore, representing a declining influence of glaciomarine sediment (Extended Data Figs. 5b and 7). The low MAR during this period may indicate less lateral advection of sediments as bathymetric highs were still ice-covered, restricting sediment transport from the east.

In contrast, the sharp MAR increase at 8.2 ka likely indicates enhanced advection of material, initiated as the local bathymetric highs fully deglaciated (Fig. 3c). A high proportion of *Chaetoceros (Hyalochaete)* resting spores in diatom assemblages from MD03-2601 indicates a more stratified and stabilized water column than in later parts of the Holocene, supporting the interpretation of enhanced glacial meltwater at this time<sup>26,29</sup>. Stratified and nutrient-rich glacial meltwater<sup>38</sup> likely created conditions favourable for bloom events.

The likely dominance of a local signal on sedimentation during the deglaciation suggests that regional processes (that is, meltwater

stratification in an enclosed calving bay embayment) drove sea-ice seasonality/break-out and bloom events, not low-latitude teleconnections. A lack of fast ice, inferred from the IPSO<sub>25</sub> proxy, allowed regular bloom events to occur in most seasons, and any subdecadal climate mode influences appear to have subordinate control. In this context, evolutionary harmonic analysis (EHA) of the CT greyscale curve and XRF Si/Ti linescan data shows power throughout the 2–7 year frequency band. Laminae counts occasionally fall into this band as well, consistent with subdecadal climate mode influences (Fig. 4). Although the annual sea-ice cycle appears to regulate bloom events during this period, we propose that subdecadal climate modes were a background influence, potentially causing earlier or later break-out of seasonal sea ice.

## Annual coastal sea-ice break-out modulated blooms (8 to 4.5 ka)

By ~8 ka, regional interpretations suggest that glacial retreat was largely complete<sup>29,36</sup>, and U1357B grain size values, MARs and physical properties (for example, NGR and CT density values) stabilize, albeit with millennial-scale variations (Fig. 3 and Extended Data Fig. 5). Bloom events occur every ~1–2 years, and rarely fall into the 2–7 year subdecadal climate mode band (Fig. 4).

Sea-ice reconstructions from *F. curta* in MD03-2601 (ref. 26) suggest reduced seasonal sea-ice duration, and IPSO<sub>25</sub> data from U1357B indicate reduced fast ice cover compared with later intervals (Fig. 3f,g). Sand percent and MAR curves indicate stronger currents and high terrigenous sediment advection, implying enhanced wind stress due to reduced ice cover (Fig. 3c,e). A reduced duration of coastal ice in this region would increase the frequency of seasonal stratification from sea-ice meltwater and open water conditions. These conditions are currently observed to trigger diatom blooms in the Mertz Polynya<sup>39</sup>. Thus, during the relative warmth of the mid-Holocene<sup>29</sup>, we propose that the primary control on bloom events was the break-up and melting of seasonal sea ice. This is consistent with the observed shift towards annual to biennial laminae frequencies.

Although some studies suggest lower ENSO-related variability from Eastern Pacific equatorial records before 4.5 ka (refs. 30,31), a shift to lower variability does not explain more frequent sea-ice break-out and bloom events (Fig. 4). As with the preceding interval, spectral power in the 2–7 year band remains evident (Fig. 4), and a subordinate influence could account for break-out events' not occurring yearly. However, the annual cycle appears to be the dominant driver of coastal sea-ice break-out events throughout this interval.

## Increased coastal sea ice reduced productivity at 4.5 ka

Around 4.5 ka, a shift occurred in all records (Fig. 3 and Extended Data Fig. 6), explained as a longer period of sea-ice cover in most years. This is reconstructed by diatom assemblages and IPSO<sub>25</sub> proxies, and the decline in MARs and sand percent. Bloom events became less frequent and occurred approximately every 2–7 years. EHA analysis of the greyscale curve and Si/Ti XRF variance indicates a similar shift to the 2–7 year band (Figs. 3 and 4). Between 1.8 and 0.8 ka, there is an exception to this pattern. The IPSO<sub>25</sub> data show a drop in fast ice, sand percent increases and laminae frequency increases to near-annual to biennial events (5–8 laminae events per 10 years; Fig. 4). Albeit a qualitative measure of fast ice<sup>28</sup>, we note a consistent pattern where IPSO<sub>25</sub> values are consistently low (for example <0.2 µg g<sup>-1</sup>) and bloom events fall into the 1–2 year band (for example, 1.8 to 0.8 ka). When fast ice increases, bloom events fall into the 2–7 year band.

A reduction in primary productivity and, therefore, bloom events is expected with an overall increase in seasonal duration of sea-ice cover, due to reduced light availability and shorter growing season<sup>26</sup>. Extensive multi-year fast ice along George V Land<sup>23</sup>,
to the east of U1357, significantly influences this region today. A regional increase in multi-year ice would reduce the occurrence of bloom events. Larger-scale seasonal sea-ice break-up would occur less frequently, thereby reducing the frequency of stratification events adjacent to the Mertz and Dumont d'Urville polynyas. The mechanism for increased sea-ice duration at 4.5 ka is not a focus of this paper. It is noted around much of the Antarctic margin, being previously interpreted as a consequence of reduced local insolation forcing and enhanced ocean-ice shelf interactions²¹. Decreases in sand percent and MAR indicate an associated drop in current speed (Fig. 3c,e and Extended Data Fig. 3). This is likely due to expanded sea-ice coverage, which reduced wind stress on the ocean surface and the vigour of the coastal current. A slowdown in sedimentation rate is also observed at MD03-2601, indicating a regional slowdown in wind-driven current and sediment advection (Extended Data Fig. 1). These lines of evidence indicate that the shift in productivity and the reduced laminae frequency were due to increased presence and duration of coastal sea ice.

## Sea ice sets system sensitivity to subdecadal climate modes

Although synoptic and katabatic winds are essential for opening and maintaining the Mertz Polynya, fast ice distributions are also important. Increased fast ice extent to the west and east restrict sea-ice export and 'back-fill' the polynya, thereby limiting its size⁴⁰. Greater fast ice extent over Site U1357, which lies to the west of the Mertz Polynya (Fig. 1), would increase the probability of back-fill events, limiting bloom events. However, greater fast ice extent could also increase stratification during favourable conditions for sea-ice break-out.

Our Adélie Land record shows that the frequency of mass biogenic bloom events in coastal polynyas of Adélie Land, East Antarctica is strongly modulated by coastal sea ice. The 2-7 year variability in bloom events, consistent with subdecadal climate mode forcing, increased after 4.5 ka. This agrees with other Antarctic Holocene records which suggest increased impacts of subdecadal climate modes on westerly winds and surface temperatures in the late Holocene¹⁹,⁴¹. Changing seasonality and distribution of coastal sea ice, and shifts in zonal winds are modulated by subdecadal climate modes under modern conditions⁸⁻¹⁰,¹⁴. We propose that the increased extent of coastal ice at 4.5 ka accentuated the impact of subdecadal climate modes on sea-ice break-out. This caused biogenic blooms to shift from annual/biennial events to subdecadal-scale modulation.

This is relevant to projections of Antarctic coastal change, as Adélie Land climate anomalies associated with climate modes differ among reanalysis studies⁹,⁴²⁻⁴⁴. Climate models also struggle to capture recent sea-ice trends, due to the complexities of ocean and atmospheric feedbacks in the Antarctic⁴⁵. Thus, critical processes appear to be underrepresented in models which project the future response of Antarctic coastal systems to increased tropical and southern mid-latitude variability.

Our data highlight the importance of sea-ice dynamics in regulating the sensitivity of biological productivity to subdecadal-scale climate modes (for example, ENSO, SAM and IOD) along the Adélie Land margin. If future warming trends result in reduction or loss of coastal sea ice, as occurred during the mid-Holocene at Adélie Land, our work suggests that more frequent bloom events will result, independent of background shifts in subdecadal-scale climate modes. This has implications for future food webs in the Antarctic, and carbon cycling processes within this globally important region of AABW formation.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of

author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-021-00816-y.

Received: 14 October 2019; Accepted: 16 July 2021; Published online: 9 September 2021

## References

1. Arrigo, K. R., van Dijken, G. & Long, M. Coastal Southern Ocean: a strong anthropogenic CO₂ sink. *Geophys. Res. Lett.* **35**, L21602 (2008).

2. Rintoul, S. R. On the origin and influence of Adélie Land bottom water. *Antarct. Res.* **Ser. 75**, 151–171 (1998).

3. Moore, J. K. & Abbott, M. R. Phytoplankton chlorophyll distributions and primary production in the Southern Ocean. *J. Geophys. Res. Ocean.* **105**, 28709–28722 (2000).

4. Arrigo, K. R. & van Dijken, G. L. Phytoplankton dynamics within 37 Antarctic coastal polynya systems. *J. Geophys. Res. C.* **108**, 1–18 (2003).

5. Comiso, J. C., McClain, C. R., Sullivan, C. W., Ryan, J. P. & Leonard, C. L. Coastal zone color scanner pigment concentrations in the Southern Ocean and relationships to geophysical surface features. *J. Geophys. Res.* **98**, 2419–2451 (1993).

6. Smith, W. O. & Nelson, D. M. Importance of ice edge phytoplankton production in the Southern Ocean. *Bioscience* **36**, 251–257 (1986).

7. Rigual-Hernández, A. S., Trull, T. W., Bray, S. G., Closset, I. & Armand, L. K. Seasonal dynamics in diatom and particulate export fluxes to the deep sea in the Australian sector of the southern Antarctic Zone. *J. Mar. Syst.* **142**, 62–74 (2015).

8. Yuan, X. ENSO-related impacts on Antarctic sea ice: a synthesis of phenomenon and mechanisms. *Antarct. Sci.* **16**, 451–452 (2004).

9. Nuncio, M. & Yuan, X. The influence of the Indian Ocean dipole on Antarctic sea ice. *J. Clim.* **28**, 2682–2690 (2015).

10. L'Heureux, M. L. & Thompson, D. W. J. Observed relationships between the El Niño Southern Oscillation and the extratropical zonal-mean circulation. *J. Clim.* **19**, 276–287 (2006).

11. Stammerjohn, S. E., Martinson, D. G., Smith, R. C., Yuan, X. & Rind, D. Trends in Antarctic annual sea ice retreat and advance and their relation to El Niño–Southern Oscillation and Southern Annular Mode variability. *J. Geophys. Res.* **113**, C03S90 (2008).

12. Fogt, R. L., Bromwich, D. H. & Hines, K. M. Understanding the SAM influence on the South Pacific ENSO teleconnection. *Clim. Dyn.* **36**, 1555–1576 (2011).

13. Saba, G. K. et al. Winter and spring controls on the summer food web of the coastal West Antarctic Peninsula. *Nat. Commun.* **5**, 1–8 (2014).

14. Arrigo, K. R. & Van Dijken, G. L. Annual changes in sea-ice, chlorophyll a, and primary production in the Ross Sea, Antarctica. *Deep. Res. Part II Top. Stud. Oceanogr.* **51**, 117–138 (2004).

15. Venables, H. J., Clarke, A. & Meredith, M. P. Wintertime controls on summer stratification and productivity at the western Antarctic Peninsula. *Linnol. Oceanogr.* **58**, 1035–1047 (2013).

16. Mayewski, P. A. et al. State of the antarctic and southern ocean climate system. *Rev. Geophys.* **47**, 1–38 (2009).

17. Escutia, C., Brinkhuis, H., Klaus, A. & Expedition 318 Scientists. Site U1357. In *Proc. IODP Vol. 318 (IODP, 2011)*.

18. Shevenell, A. E., Ingalls, A. E., Domack, E. W. & Kelly, C. Holocene Southern Ocean surface temperature variability west of the Antarctic Peninsula. *Nature* **470**, 250–254 (2011).

19. Pike, J., Swann, G. E. A., Leng, M. J. & Snelling, A. M. Glacial discharge along the west Antarctic Peninsula during the Holocene. *Nat. Geosci.* **6**, 199–202 (2013).

20. Whitworth, T. et al. in *Ocean, Ice, and Atmosphere: Interactions at the Antarctic Continental Margin* (eds. Jacobs, S. S. & Weiss, R. F.) 75, 1–27 (American Geophysical Union, 1998).

21. Ashley, K. E. et al. Mid-Holocene Antarctic sea-ice increase driven by marine ice sheet retreat. *Clim. Past* **17**, 1–19 (2021).

22. Dunbar, R. B., Anderson, J. B., Domack, E. W. & Jacobs, S. S. in *Oceanology of the Antarctic Continental Shelf* (ed. Jacobs, S. S.) 309–312 (American Geophysical Union, 1985).

23. Massom, R. A. et al. Change and variability in East Antarctic sea ice seasonality, 1979/80–2009/10. *PLoS ONE* **8**, e64756 (2013).

24. Denis, D. et al. Seasonal and subseasonal climate changes recorded in laminated diatom ooze sediments, Adélie Land, East Antarctica. *Holocene* **16**, 1137–1147 (2006).

25. Maddison, E. J., Pike, J. & Dunbar, R. B. Seasonally laminated diatom-rich sediments from Dumont d'Urville Trough, East Antarctic Margin: Late-Holocene Neoglacial sea-ice conditions. *Holocene* **22**, 857–875 (2012).

26. Denis, D. et al. Holocene productivity changes off Adélie Land (East Antarctica). *Paleoceanography* **24**, PA3207 (2009).
27. Crosta, X., Debret, M., Denis, D., Courty, M. A. & Ther, O. Holocene long-and short-term climate changes off Adélie Land, East Antarctica. *Geochem. Geophys. Geosyst.* **8**, 1–15 (2007).

28. Belt, S. T. et al. Source identification and distribution reveals the potential of the geochemical Antarctic sea ice proxy IPSO₂₅. *Nat. Commun.* **7**, 12655 (2016).

29. Crosta, X., Denis, D. & Ther, O. Sea ice seasonality during the Holocene, Adélie Land, East Antarctica. *Mar. Micropaleontol.* **66**, 222–232 (2008).

30. Moy, C. M., Seltzer, G. O., Rodbell, D. T. & Anderson, D. M. Variability of El Niño/Southern Oscillation activity at millennial timescales during the Holocene epoch. *Nature* **420**, 162–165 (2002).

31. Carré, M. et al. Holocene history of ENSO variance and asymmetry in the eastern tropical Pacific. *Science* **345**, 1045–1047 (2014).

32. Cobb, K. M. et al. Highly variable El Niño–Southern Oscillation throughout the Holocene. *Science* **339**, 67–70 (2013).

33. Karamperidou, C., Di Nezio, P. N., Timmermann, A., Jin, F.-F. & Cobb, K. M. The response of ENSO flavors to mid-Holocene climate: implications for proxy interpretation. *Paleoceanography* **30**, 527–547 (2017).

34. Abram, N. J. et al. Evolution of the Southern Annular Mode during the past millennium. *Nat. Clim. Change* **4**, 564–569 (2014).

35. Abram, N. J. et al. Palaeoclimate perspectives on the Indian Ocean Dipole. *Quat. Sci. Rev.* **237**, 106302 (2020).

36. Mackintosh, A. N. et al. Retreat history of the East Antarctic Ice Sheet since the Last Glacial Maximum. *Quat. Sci. Rev.* **100**, 10–30 (2014).

37. Leventer, A. et al. Marine sediment record from the East Antarctic margin reveals dynamics of ice sheet recession. *GSA Today* **16**, 4–10 (2006).

38. Gerringa, L. J. A. et al. Iron from melting glaciers fuels the phytoplankton blooms in Amundsen Sea (Southern Ocean): iron biogeochemistry. *Deep Sea Res. II* **71–76**, 16–31 (2012).

39. Moreau, S. et al. Sea ice meltwater and circumpolar deep water drive contrasting productivity in three Antarctic polynyas. *J. Geophys. Res. Ocean.* **124**, 2943–2968 (2019).

40. Massom, R. A. et al. Effects of regional fast-ice and iceberg distributions on the behaviour of the Mertz Glacier polynya, East Antarctica. *Ann. Glaciol.* **33**, 391–398 (2001).

41. Etourneau, J. et al. Holocene climate variations in the western Antarctic Peninsula: evidence for sea ice extent predominantly controlled by changes in insolation and ENSO variability. *Clim. Past* **9**, 1431–1446 (2013).

42. Schneider, D. P. et al. Observed Antarctic interannual climate variability and tropical linkages. *J. Clim.* **25**, 4048–4066 (2012).

43. Ciasto, L. M. et al. Teleconnections between tropical Pacific SST anomalies and extratropical Southern Hemisphere climate. *J. Clim.* **28**, 56–65 (2015).

44. Marshall, G. J. & Thompson, D. W. J. The signatures of large-scale patterns of atmospheric variability in Antarctic surface temperatures. *J. Geophys. Res. Atmos.* **121**, 3276–3289 (2016).

45. IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (eds Weyer, N. M. et al.) (IPCC, 2019).

46. Mouginot, J., Scheuchl, B. & Rignot, E. MEaSURES Antarctic boundaries for *IPY 2007–2009 from Satellite Radar, Version 2*. [68°S, 65°S; 148°E, 138°E] (NASA National Snow and Ice Data Center, 2017).

47. Fretwell, P. et al. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica. *Cryosphere* **7**, 375–393 (2013).

48. Beaman, R. J., O'Brien, P. E., Post, A. L. & De Santis, L. A new high-resolution bathymetry model for the Terre Adélie and George V continental margin, East Antarctica. *Antarct. Sci.* **23**, 95–103 (2011).

49. NASA Ocean Biology Processing Group. MODIS-Aqua Level 2 Ocean Color Data Version R2018.0 NASA Ocean Biology DAAC (2017); https://doi.org/10.5067/AQUA/MODIS/L2/OC/2018

50. Hu, C., Lee, Z. & Franz, B. Chlorophyll a algorithms for oligotrophic oceans: a novel approach based on three-band reflectance difference. *J. Geophys. Res.* **117**, C01011 (2012).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2021
## Methods

**Age model.** The age model was developed from 87 ¹⁴C dates from acid-insoluble bulk sedimentary organic carbon to constrain the ages of the sediment between ~11.4 ka BP and modern day²¹ (Extended Data Fig. 1). The age model only resolves ages to 11.4 ka as ages older than this are anomalously old and assumed to incorporate reworked carbon of pre-Last Glacial Maximum age. This was also indicated by a larger terrestrial contribution observed at the lowest section of the core according to XRF data³¹. Ages younger than 11.4 ka are less likely to be affected by reworked carbon at the Adélie Drift site, as lack of IRD suggests that direct glacial influences were negligible. The consistent stratigraphic order of the 87 radiocarbon ages and sedimentation rates through the Holocene supports this interpretation. Very few Antarctic marine sediment core records presently have age models of this resolution and with this level of stratigraphic integrity. A reservoir correction age of 1,200 ± 100 was applied to the depth-to-age conversion calculated by BACON, which uses a Bayesian iteration scheme that invokes memory from dates above any given horizon and produces a weighted mean and median age-depth curve⁵². This correction is consistent with the uncalibrated age of the uppermost sample of 1,310 years.

Since C₁₆ fatty acids decompose rapidly in the water column and sediment⁵³,⁵⁴ compound-specific (CS) ages in Antarctic sediments that contain relict carbon from glacial reworked sediments often show younger ages than bulk ages⁵⁵,⁵⁶. Yamane et al. (2014) reported an age model based on CS ¹⁴C ages using C₁₆ fatty acids from core U1357A, and ages were reported with one-sigma uncertainty. In this study, the background level of the study was rigorously re-examined using the latest background evaluation method for small-scale ¹⁴C analysis developed at the Atmosphere and Ocean Research Institute, University of Tokyo⁵⁷. The modern carbon contamination (MCC) was evaluated from ¹⁴C value of IAEA-C4 (wood: Δ¹⁴C = -998.0 to -995.6‰) which was processed and measured by accelerator mass spectrometry in the same batch as other unknown samples (Supplementary Fig. 1a). The background correction was carried out differently depending on sample size using the relationship between sample size and background (Supplementary Fig. 1b). To externally evaluate the reliability of the MCC, we estimated the core-top CS ¹⁴C value using the mean sedimentation rate of lithostratigraphic unit I (0–170.25 m below seafloor). Based on the revised CS ¹⁴C ages, it is estimated that the Δ¹⁴C value of core-top sediment is about -147‰. This Δ¹⁴C value is in agreement with the pre-bomb dissolved inorganic carbon (DIC) Δ¹⁴C value of the Southern Ocean (-149.8 ± 10.4‰; ref. ⁵⁸), hence validating the CS ¹⁴C. The values are co-plotted with bulk ages with two-sigma uncertainties and show that CS and bulk organic ¹⁴C ages are consistent (Extended Data Fig. 1). This is the case for earlier values (that is, ref. ⁵¹) if we plotted values with two-sigma uncertainties, thus all ages are consistent within statistical uncertainties. Consequently, these revised compound-specific radiocarbon assessments support our inference that contamination of reworked carbon in these rapidly deposited biogenic-rich samples is minimal (Supplementary Table 1). Below, we independently assess the reliability of the bulk organic carbon age model by comparison with age models from the nearby core MD03-2601 (Extended Data Fig. 1).

The BACON methodology was applied to the ¹⁴C dates from MD03-2601 (refs. ²⁷,²⁹) to recalibrate the MD03-2601 age model²¹. The model shown in this paper is different from the one used previously³⁰, which used an inferred meteorite impact at ~15 m to determine an age of 4 ka at that depth. The old age model also removed two ¹⁴C dates at 4.4 and 5.6 ka years due to the assumption that these ages were anomalously old relative to the meteorite impact. However, the meteorite age-depth correlation cannot provide absolute age control, and the new age model presented herein indicates that the impact occurred at around 5.4 ka. Comparison between the U1357B and new MD03-2601 age model show strong covariance in sedimentation rates and suggest a regional sedimentation advection process (Extended Data Fig. 1).

**Depth scales.** Core recovery from each 9.5 m piston core run often exceeded 100% due to expansion as the core decompressed during recovery. Data derived from these initial core lengths are termed the csf-a depth scale. The standard IODP procedure to correct for expansion is to apply a linear compression algorithm to scale recovery back to 100% and create a new scale (csf-b), as it is assumed that expansion is uniform in the core. However, in U1357, expansion due to biogenic gas was particularly high and resulted in discrete sections of core being pushed apart, creating voids in the depth scale that did not represent real gaps in the stratigraphy. To account for this, the voids are numerically removed, and the depth scale adjusted, before linear compression is applied (if recovery still exceeds 100%). In this paper, we term this the csf-d scale (noting that this is not an official IODP depth scale term). Although cap expansion gaps (voids) are removed within individual core runs, the csf-d scale still contains sections with no core recovery at the base of some runs where there was less than 100% after voids within the cores were numerically removed. The sections with no core recovery were as follows: 48.82–50.0 m, 58–59.5 m, 66.19–69.0 m, 76.81–78.5 m, 86.5–88.0 m, 95.66–97.5 m, 105.16–107.0 m, 115.72–116.5 m, 124.38–126.0 m, 134.41–135.5 m, 144.23–145 m, 153.85–154.5 m and 163.94–164.0 m. Slight differences in these depths could have occurred in core storage before CT and XRF scans.

**Composite core.** Three holes (U1357A, B and C) were drilled in the Adélie Basin as part of IODP expedition 318 (ref. ¹⁷). Drilling multiple holes is standard IODP procedure for sites with palaeoceanographic focus to address core breaks and other intervals of incomplete recovery; a complete and continuous stratigraphy can usually be constructed by splicing sections from individual holes into a stratigraphic composite section. This is usually achieved in IODP cores by using physical core properties to guide placement of the least disturbed, highest recovery intervals in the spliced sections. However, cores from Site U1357 are problematic in this context as extremely high biogenic and gaseous content precluded many physical property measurements, such as magnetic susceptibility, from being registered beyond typical noise levels. This made construction of a composite core at subcentennial-scale precision extremely difficult. Given the difficulties in creating a spliced record, hole U1357B was selected as the best core for this analysis because it had less gas-related disturbances than hole A, and a more complete record than hole C, which produced a shorter core. Additionally, it also has a higher-resolution age model.

**Linear sedimentation rates.** The linear sedimentation rates were calculated for every centimetre using the age-depth model above. These were then binned every 10 cm.

**MARs.** Terrigenous and biogenic MARs were calculated using the formula below:

$$MAR = \%X \times (LSR \times BD)$$

MAR = mass accumulation rate (g cm⁻² per year)

LSR = linear sedimentation rate (cm per year)

X = percent abundance of the component of interest (that is, terrigenous or biogenic)

BD = bulk density (g cm⁻³)

Shipboard bulk density measurements were not collected on the U1357B core, which was preserved as whole-round sections until post-expedition sampling several months after collection. Moisture and density (MAD) bulk densities from the U1357A core (collected at the same site location) were used instead, with a linear fit taken though these data to derive a downhole estimate of bulk density¹⁷. The associated depths of these discrete samples were converted to age using the U1357A age-depth model. This model uses 36 bulk organic carbon dates and demonstrates the age versus depth relationship using the same Bayesian approach used in the U1357B age model. A linear fit between the age and density measurements of U1357A was interpolated to the U1357B age scale to determine the densities for U1357B. Biogenic silica and terrigenous percentage were determined using alkaline extraction spectrophotometric methods⁶⁰.

**Grain size analysis.** Grain size analysis was performed on 341 samples. Samples were treated twice with 1 M sodium hydroxide (NaOH) solution in an 80 °C water bath for 24 h to remove biogenic opal, then treated with hydrogen peroxide (H₂O₂) to remove organic material. As terrigenous material formed a minor component of the bulk sediment, post-treatment sample mass varied from ~0.035 to 0.8 g. Samples were measured on a Beckman Coulter LS 13 320 laser diffraction particle size analyser. Eighty-four sub-samples were taken after chemical treatment with NaOH and H₂O₂ to assess reproducibility. Twelve samples were split into two subsamples, and chemical treatment was performed on each subsample to test for biases relating to subsampling and chemical dissolution. Correlations calculated using least-squares regression between the original and repeat measurements were r² = 0.74 for the post chemical treatment replicates and r² = 0.75 for the pre-chemical treatment replicates.

**Computed tomography scans.** X-ray computed tomography (CT) scanners quantify the amount of X-ray energy absorbed (attenuated) by a particular object and display the resulting attenuation coefficients as a greyscale image⁶¹. Pixel values within these images are expressed as greyscale values or Hounsfield units (also known as CT number) which are calculated by comparing the sample attenuation coefficient with that of water¹⁵.

X-ray attenuation is a function of the density, porosity, chemical composition and grain size of the sample⁶². Brighter areas in the image represent higher attenuation, while darker areas represent lower attenuation. CT scans were completed on core U1357B using a Toshiba Aquilion TXL CT scanner at the Department of Petroleum Engineering, Texas A&M University. Axial scans were completed at 135 kVp and 200 mA, and coronal slices were created in open-source HOROS software⁶³. The resolution averages 1.3 pixels per millimetre, and each core was exported as its own Digital Imaging and Communications in Medicine image stack which contained 512 images. From there, the best image (accounting for cracks and other spaces in the core) from each stack was selected and exported to another CT processing software (FIJI⁶⁴) for greyscale analysis and laminae counting.

**HBI/isoprenoid/IPSO₂₅ data.** IPSO₂₅ is another name for the highly branched isoprenoid (HBI) lipid biomarker (diene II). C₂₅-HBI alkenes, in particular the di-unsaturated C₂₅-HBI with a double bond, also referred to as diene, were extracted at Laboratoire d'Océanographie et du Climat: Experimentations et
Approches Numériques, using a mixture of 9 mL CH₂Cl₂/MeOH (2:1, v:v) to which internal standards (7-hexyl nonadecane, 9-octyl heptadecane and androstanol) were added; several sonication and centrifugation steps were applied in order to properly extract the selected compounds⁴¹. After drying with N₂ at 35 °C, the total lipid extract was fractionated over a silica column into an apolar and polar fraction using 3 mL hexane and 6 mL CH₂Cl₂/MeOH (1:1, v:v), respectively. HBIs were obtained from the apolar fraction by fractionation over a silica column using hexane as eluent following the procedures reported in refs. ⁶⁵,⁶⁶. After removing the solvent with N₂ at 35 °C, elemental sulphur was removed using the tetrabutylammonium sulphite method⁶⁷,⁶⁸. The obtained hydrocarbon fraction was analysed using an Agilent 7890A gas chromatograph (GC) fitted with 30 m fused silica Agilent J&C GC column (0.25 mm i.d., 0.25 µm film thickness), coupled to an Agilent 5975 C series mass selective detector. Spectra were collected using Agilent MS-Chemstation software. Individual HBIs were identified on the basis of comparison between their GC retention times and mass spectra with those of previously authenticated HBIs (for example, from ref. ⁶⁹) using Mass Hunter software.

**Image analysis.** *Greyscale curve.* Any pixel value less than zero was converted to non-values (NaNs) by thresholding the images in FIJI⁴⁰. This eliminated noise from pervasive sub-millimetre to millimetre-scale cracks resulting from expansion due to biogenic gas in the cores (Extended Data Fig. 4).

A single greyscale curve was created by taking a line profile of the greyscale image for each core. The line profile was four pixels wide, with the pixel value of each row being the average of these four pixels. The profile was chosen to minimize core disruptions. Many CT studies choose to average all rows across the whole width of the image, but this was not possible due to the middle of this core having previously been sampled using U-Channel methods, and due to dipping of the laminae along the core liner. These image curves were then corrected for any depth offset introduced by the core liner and CT machine, and concatenated into a final dataset.

*Laminae counts.* The top and bottom of bright laminae were picked manually throughout the entire core.

Some laminae had sharp divisions between bright and dark pixels, while others had a gradual transition. In addition, some bright laminae were interspersed among a slightly lighter background, making it difficult to distinguish between multiple laminae and single events. We counted such intervals as a single lamina, and suggest that these could represent seasons when there were multiple blooms or prolonged bloom events. Visual picking of the laminae can be subjective, but was preferred over automated methods due to noise produced by gas expansion cracks, which varied from core to core. To assess this subjective aspect, laminae picks were visually overlaid on the greyscale curve to evaluate consistency throughout the length of the core (Extended Data Fig. 4). Some laminae were disrupted by cracks. We manually removed the laminae disrupted by several centimetres or more, but these accounted for less than 0.1% of laminae. Laminae were binned into 10-year intervals (Figs. 3 and 4). For bins that contained a missing interval, that is, the base of a 9.5 m core run where recovery was <100%, the binned laminae amounts were scaled to represent the actual number of years per bin. For the scaled 10-year bins, seven data points were removed because the bins contained fewer than 2 years of data. Comparison of the manually picked laminae with evolutionary spectra of the raw greyscale curve and Si/Ti values from XRF linescan data was conducted to independently verify the frequencies identified (Fig. 4).

*Evolutionary spectral analysis.* Before analysis, the greyscale data were interpolated to 0.1 years (from an average spacing of 0.041 years) and the XRF data were interpolated to 0.4 years (from an average time step of 0.44 years), using piecewise-linear interpolation. EHA using the Thomson multi-taper method to determine power spectra was performed in the R package 'Astrochron'⁷⁰ using both the XRF and CT greyscale data. Outliers were removed from the series using the 'Trim' function in Astrochron, which uses a box plot algorithm with a coefficient of 1.5 to identify values greater than or less than 1.5 times the interquartile range from quartile 3 and quartile 1, respectively. For EHA on the CT greyscale data, a multi-taper method time-bandwidth product of 4, window width of 100 years and step size of 20 years were used. For the lower-resolution XRF data, a multi-taper method time-bandwidth product of 3, window width of 70 years and step size of 10 year were used. Resulting spectra were seen to be relatively insensitive to window width and step size, and time series analysis on other XRF productivity proxies (Ba/Ti, Si/Al) yielded similar results. In all datasets analysed, power was normalized so that the maximum power in each window is unity.

The manual laminae counts, binned into 10-year intervals, were then overlaid on the EHA results and show consistent centennial-scale shifts in the power of the 2–7 year frequency bands. This indicates that binned laminae frequencies are representative of the EHA results and are able to capture higher-frequency variations in bloom events.

XRF. XRF data were measured using an AVAATECH XRF core scanner at the Joint Oceanographic Institutions for Deep Earth Sampling (JOIDES) Resolution Science Operator XRF facility, located at the Gulf Coast Repository at Texas A&M

University Research Park. Measurements were undertaken at 0.5 cm resolution (where possible) with a 5 mm slit size using generator settings of 10 kV and current of 0.8 mA. The sampling time was set at 45 s, and scanning took place directly at the split core surface of the archive half. The split core surface was covered with a 4-micron thin SPEXCerti Prep Ultralene1 foil to avoid contamination of the XRF measurement unit and desiccation of the sediment.

Biogenic silica concentrations in sediments (%BSi) are commonly used as an indicator of past diatom and radiolarian productivity in high-latitude marine sediments (see, for example, refs. ²¹,²²). Silicon (Si) is the main component of biogenic opal, and Si-based ratios are commonly used as %BSi proxies³. Estimating %BSi from Si content or Si-based ratios obtained by XRF scanner requires site-specific calibration, but comparison with the Si/Ti ratio shows distribution almost parallel with %BSi records as a function of depth (see, for example, ref. ⁷³). Nevertheless, use of Si as a productivity proxy should be applied with caution, because Si can also be controlled by siliciclastic material during low-productivity periods, even in polar regions⁷⁴, and light elements, such as Si or Al, have low detectability by XRF scanner measurements when present in low concentrations⁷⁵.

Site U1357B is a laminated diatom ooze. Diatom content estimated from smear slides has a mean of 91% (ref. ¹⁷). %BSi content in this study ranges from 30% to 63% with an average of 48 %BSi for the late Holocene. Si detection by the XRF scanner is not an issue, as the average Si peak area is ~200,000 counts. In any case, the high opal content masks Si input related to siliciclastic material. To correct dilution effect and obtain a first-order discrimination between biogenic and detrital Si, we normalized Si to Ti. This normalization assumes that Ti is a conservative element associated only with the terrigenous fraction and Si/Ti ratio of the terrigenous matter remains almost constant over the period studied. We use the obtained Si/Ti ratio as a semi-quantitative record of the siliceous productivity, in agreement with previous studies that use Si/Ti or equivalent ratios as a productivity proxy in both marine⁷⁶,⁷⁷ and lacustrine records⁷⁸,⁷⁹.

*Correlation analysis between laminae counts, biogenic MAR and sand percent.* Laminae counts, biogenic MAR and sand percent were linearly interpolated to a common 100-year step. Regression statistics were calculated from 10,050 BP onwards, as the glaciated environment before this time is not representative of current relationships.

## Data availability

The raw greyscale data, light laminae depths, light laminae sand percent, XRF silicon, XRF titanium and HBI diene data for IODP Site U1357B are available at https://doi.org/10.1594/PANGAEA.933380.

## References

51. Yamane, M. et al. Compound-specific ¹⁴C dating of IODP expedition 318 core U1357A obtained off the Wilkes Land coast, Antarctica. *Radiocarbon* **56**, 1009–1017 (2014).

52. Blaauw, M. & Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. *Bayesian Anal.* **6**, 457–474 (2011).

53. Canuel, E. A. & Martens, C. S. Reactivity of recently deposited organic matter: degradation of lipid compounds near the sediment-water interface. *Geochim. Cosmochim. Acta* **60**, 1793–1806 (1996).

54. Ohkouchi, N., Kawamura, K. & Taira, A. Fluctuations of terrestrial and marine biomarkers in the western tropical Pacific during the last 23,300 years. *Paleoceanography* **12**, 623–630 (1997).

55. Yokoyama, Y. et al. Widespread collapse of the Ross Ice Shelf during the late Holocene. *Proc. Natl Acad. Sci. USA* **113**, 2354–2359 (2016).

56. Prothro, L. O. et al. Timing and pathways of East Antarctic Ice Sheet retreat. *Quat. Sci. Rev.* **230**, 106166 (2020).

57. Yamane, M. Small- to ultra-small-scale radiocarbon measurements using newly installed single-stage AMS at the University of Tokyo. *Nucl. Instrum. Methods Phys. Res. B* **455**, 238–243 (2019).

58. Berkman, P. A. & Forman, S. L. Pre-bomb radiocarbon and the reservoir correction for calcareous marine species in the Southern Ocean. *Geophys. Res. Lett.* **23**, 363–366 (1996).

59. Denis, D. et al. Holocene glacier and deep water dynamics, Adélie Land region, East Antarctica. *Quat. Sci. Rev.* **28**, 1291–1303 (2009).

60. Mortlock, R. A. & Froelich, P. N. A simple method for the rapid determination of biogenic opal in pelagic marine sediments. *Deep Sea Res. A.* **36**, 1415–1426 (1989).

61. St-Onge, G. & Long, B. F. CAT-scan analysis of sedimentary sequences: an ultrahigh-resolution paleoclimatic tool. *Eng. Geol.* **103**, 127–133 (2009).

62. Boespflug, X., Long, B. F. N. & Occhietti, S. CAT-scan in marine stratigraphy: a quantitative approach. *Mar. Geol.* **122**, 281–301 (1995).

63. *Horos Project* (Horos, 2017).

64. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. *Nat. Methods* **9**, 676–682 (2012).

65. Belt, S. T. et al. A novel chemical fossil of palaeo sea ice: IP25. *Org. Geochem.* **38**, 16–27 (2007).
66. Massé, G. et al. Highly branched isoprenoids as proxies for variable sea ice conditions in the Southern Ocean. *Antarct. Sci.* **23**, 487–498 (2011).

67. Jensen, S., Renberg, L. & Reutergårdh, L. Residue analysis of sediment and sewage sludge for organochlorines in the presence of elemental sulfur. *Anal. Chem.* **49**, 316–318 (1977).

68. Riis, V. & Babel, W. Removal of sulfur interfering in the analysis of organochlorines by GC-ECD. *Analyst* **124**, 1771–1773 (1999).

69. Johns, L. et al. Identification of a C₂₅ highly branched isoprenoid (HBI) diene in Antarctic sediments, Antarctic sea-ice diatoms and cultured diatoms. *Org. Geochem.* **30**, 1471–1475 (1999).

70. Meyers, S. *Astrochron: An R Package for Astrochronology* (CRAN, 2014).

71. Ragueneau, O. et al. A review of the Si cycle in the modern ocean: recent progress and missing gaps in the application of biogenic opal as a paleoproductivity proxy. *Glob. Planet. Change* **26**, 317–365 (2000).

72. Iwasaki, S., Takahashi, K., Ogawa, Y., Uehara, S. & Vogt, C. Alkaline leaching characteristics of biogenic opal in Eocene sediments from the central Arctic Ocean: a case study in the ACEX cores. *J. Oceanogr.* **70**, 241–249 (2014).

73. Brown, E. T. in *Micro-XRF Studies of Sediment Cores. Developments in Paleoenvironmental Research* (eds. Coudace, I. W. & Rothwell, R. G.) 267–277 (Springer, 2015).

74. Jimenez-Espejo, F. J. et al. Changes in detrital input, ventilation and productivity in the central Okhotsk Sea during the marine isotope stage 5e, penultimate interglacial period. *J. Asian Earth Sci.* **156**, 189–200 (2018).

75. Rothwell, R. G. & Coudace, I. W. in *Micro-XRF Studies of Sediment Cores. Developments in Paleoenvironmental Research* (eds. Coudace, I. W. & Rothwell, R. G.) 25–102 (Springer, 2015).

76. Agnihotri, R., Altabet, M. A., Herbert, T. D. & Tierney, J. E. Subdecadally resolved paleoceanography of the Peru margin during the last two millennia. *Geochem. Geophys. Geosyst.* **9**, Q05013 (2008).

77. Dickson, A. J., Leng, M. J., Maslin, M. A. & Röhl, U. Oceanic, atmospheric and ice-sheet forcing of South East Atlantic Ocean productivity and South African monsoon intensity during MIS-12 to 10. *Quat. Sci. Rev.* **29**, 3936–3947 (2010).

78. Martin-Puertas, C., Brauer, A., Dulski, P. & Brademann, B. Testing climate-proxy stationarity throughout the Holocene: an example from the varved sediments of Lake Meerfelder Maar (Germany). *Quat. Sci. Rev.* **58**, 56–65 (2012).

79. Melles, M. et al. 2.8 Million years of arctic climate change from Lake Elgygytgyn, NE Russia. *Science* **337**, 315–320 (2012).

80. Killick, R., Fearnhead, P. & Eckley, I. A. Optimal detection of changepoints with a linear computational cost. *J. Am. Stat. Assoc.* **107**, 1590–1598 (2012).

81. Lavielle, M. Using penalized contrasts for the change-point problem. *Signal Process.* **85**, 1501–1510 (2005).

## Acknowledgements

This research used samples and data provided by IODP expedition 318, sponsored by the US National Science Foundation (NSF) and participating countries under the management of the Consortium for Ocean Leadership, including the Australian and New

Zealand International Ocean Discovery Program Consortium. Funding was provided by Royal Society Te Apārangi Marsden Fund (18-VUW-089 to R.M.M. and 15-VUW-131 to N.A.N.B.) and the New Zealand Ministry of Business, Innovation and Employment through the Antarctic Science Platform (ANTA1801). Funding was also provided by the New Zealand Ministry of Business, Innovation and Employment Strategic Science Investment Fund (SSIF) through GNS Science (grant 540GCT32). We acknowledge funding from the Dumont d'Urville NZ-France Science and Technology Programme, MARICE project (Marine and Ice core reconstruction of East Antarctic sea ice variability over the past 2,000 years) (project nos. 45455NF and 19-VUW-047-DDU Catalyst Fund, RSNZ). J.E. and X.C. acknowledge funding by the ERC StG ICEPROXY (203441), the ANR CLIMICE and the FP7 Past4Future (243908) projects. F.J.J.-E. was funded by project 201830I092 (Spanish Research Council). C.E. and F.J.J.-E acknowledge funding by the Spanish Ministry of Science and Innovation (grant CTM2017-89711-C2-1-P), co-funded by the European Union through FEDER funds. C.R.R. was funded by a University of Otago research grant and a L'Oréal-UNESCO For Women in Science Australia and New Zealand Fellowship. The Natural Environment Research Council funded K.E.A. (CENTA PhD; NE/L002493/1) and J.B. (standard grant Ne/I00646X/1). Y.Y. was funded by the Japan Society for Promotion of Science (JSPS) grant no. JP20H00193. S.F.P. was supported by National Science Foundation grant OPP-0732796. R.B.D. was supported by National Science Foundation grants PLR-1644118 and OCE-1129101. The authors acknowledge the Norwegian Polar Institute's Quantarctica package, and the use of imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov/), part of the NASA Earth Observing System Data and Information System.

## Author contributions

K.M.J., R.M.M. and N.A.N.B designed the study and wrote the paper with input from all authors. K.M.J., R.M.M. and H.J.H. analysed the X-ray CT data. R.M.M. and A.A. conducted the grain size analyses. J.E. produced the HBI data. F.J.J.-E. produced the XRF geochemical data. C.R.R. conducted the opal (%BSi) measurements. M.Y. and Y.Y. analysed and provided the compound-specific ¹⁴C ages. R.B.D. and C.E. were lead proponents of the ancillary IODP expedition 318 proposal to core IODP Site 1357. All authors contributed to the interpretations of data and finalization of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41561-021-00816-y.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-021-00816-y.

**Correspondence and requests for materials** should be addressed to Katelyn M. Johnson.

**Peer review information** *Nature Geoscience* thanks Kaarina Weckström and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor(s): James Super.

**Reprints and permissions information** is available at www.nature.com/reprints.
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 1** | Corrected and calibrated radiocarbon age models for U1357B and MD03-2601. Age Model for U1357B (blue symbols show calibrated and corrected ages, and blue line of calculated Bayesian age-depth curve, using BACON R package; dotted lines are 2 sigma uncertainty) and MD03-2601 (red symbols and lines as for U1357B) show that sediment advection is a regional signal as sedimentation rates covary. Compound specific ages from U1357A are shown in green with 2 sigma uncertainties. U1357B is a longer core with 87 <sup>14</sup>C dates, leading to a much higher resolution age model. U1357B also includes the last 1,000 years, which is lacking in MD03-2601.
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 2** | Grain size analysis of light laminae. Light laminae in three intervals (Interval 1: ~4.5 to present; Interval 2: ~8-4.5 ka; Interval 3: ~11.4-8 ka;) have a dominant silt-fine sand mode and are less than 125 µm. Sand percent provides a measure of this coarse silt to fine sand mode and is used as a proxy to capture the upper values of relative current strength.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 3** | Correlations of Adélie Drift Proxies. (a,b) Overlay plot and correlation between laminae counts (frequency of biological blooms) and biogenic mass accumulation rates (advection of biological material by wind driven current strength). (c,d) Overlay plot and correlation between sand percent (grain size proxy for wind driven current strength) and biogenic mass accumulation rates (advection of biological material by wind driven current strength). (d,e) Overlay plot and correlation between sand percent (grain size proxy for wind driven current strength) and laminae counts (frequency of biological blooms).
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 4 | Image Analysis example.** From left to right: Line-scan core photo, CT image, raw greyscale curve with light laminae picked in orange, and XRF titanium data. The CT image visually enhances the laminations, and provides a sub-mm resolved greyscale curve that better captures the rapid shifts in sedimentation. Frequencies from the manually picked laminae (Fig. 3d) match those extracted from the greyscale curve (Fig. 4).
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 5** | Greyscale data in relation sedimentological changes at the site. Greyscale data compared to various sedimentological proxies of the core. (a) GRA Bulk Density (b) Natural Gamma Radiation (NGR) (c) XRF peak area of silicon (d) XRF peak area of titanium (e) XRF peak area of silicon/titanium ratio used to indicate long-term changes in biological productivity (f) Greyscale curve from the CT images. XRF, GRA Bulk Density, and Greyscale interpolated to 5 cm. Black curves indicate a robust LOESS smoothing using 2% of the data points for all data sets. GRA Bulk Density and NGR data from ref. ¹⁷.
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 6** | Holocene Mean Changepoints. Calculation of shifts in the mean of various Holocene records using statistical changepoint analysis. This method identifies the point at which the mean of a timeseries changes most significantly by finding the point at which the total residual error from the mean of each section is minimized (ref. ⁸⁰,⁸¹). Shifts in Adélie Land proxies (**a-d**) occur between ~4.3 and ~4.7 ka.
<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 7 | IBRD in U1357B.** Section 19H5 (169.98-170.67 m) from U1357B is the only portion of core with substantial IBRD. The top of this section is dated to -11.365ka.