Marine sediment record
The Calypso Square gravity core AMD15-CASQ1 (77°15.035′ N, 74°25.500′ W, 692 m water depth) and accompanying box core (BC; same location) were retrieved aboard the CCGS Amundsen during the ArcticNet 2015 Leg 4a expedition in 2015, in accordance with relevant permits and local laws. The CASQ corer recovered a sequence 543 cm long, while the box core was 40 cm long. Sediment material from these cores is stored at the Geological Survey of Denmark and Greenland and available upon reasonable request to the first and corresponding author (SRI).
Computed Tomography (CT) scanning of the core was performed using a Siemens SOMATOM Definition AS + 128 at the Institut National de la Recherche Scientifique (INRS), Quebec, Canada. The tomograms were converted into digital DICOM format using a standard Hounsfield scale (HU scale) from −1024 to 3071, where −1024 corresponds to the density of air, 0 to the density of water and 2500 to the density of calcite.
The age control on the marine sediment record was provided by 11 accelerator mass spectrometry (AMS) radiocarbon dates on mollusc shells (Supplementary. Table 1) at the Keck Carbon Cycle AMS Facility, University of California, Irvine, US, and 210Pb/137Cs measurements conducted on 20 samples at the Gamma Dating Center, Copenhagen University, Denmark. In the box core, the content of unsupported 210Pb showed a clear exponential decline with depth (Supplementary Fig. 1). A clear 137Cs peak was not detected, but the 210Pb-based chronology dates the earliest sample with 137Cs to 1969 ± 2 years, which is close to the expected date, 1963, for the global 137Cs peak induced by nuclear weapons testing in the atmosphere. This, and the very uniform exponential decline in unsupported 210Pb with depth, gives confidence in the calculated chronology. A mixed age-depth model, using both 210Pb and 14C dates, was constructed using BACON, an open-source package of ‘R’54. This Bayesian accumulation model code allows for greater flexibility in sedimentation rates between dated intervals than traditional linear age-depth models54. The AMS radiocarbon dates were calibrated with the Marine13 IntCal1355, and the regional marine reservoir offset was estimated based on existing 14C data from marine specimens collected before the mid-1950s. Distinct regional offset values have been proposed for Arctic Canada, but do not include the Smith Sound region56. Existing data from NW Greenland show local reservoir correction (ΔR) values ranging from -40 years in the Inglefield Fjord to +320 years in Ellesmere Island (the latter consistent with the proposed 335 ± 85 years for the Canadian Arctic Archipelago56). However, these samples have been retrieved from shallow sites (<40 m water depth), which are unlikely to reflect influence from the West Greenland Current. Data from the only deeper site in the NOW region are based on measurements of the mollusc Astarte montagui and indicate an offset of 140 ± 60 years57. This value is consistent with measurements from the Disko Bay region, also under the influence of the West Greenland current58. We have therefore chosen to use an offset value of 140 ± 60 years. Ages are reported in calibrated years before 2000 (cal yrs b2k) and year CE (for the box-core record presented in Fig. 7).
TOC measurements were carried out at the Geological Survey of Denmark and Greenland, at 5 cm intervals in the CASQ core and at 2 cm intervals in the box-core. Dried sediment samples (~0.5 g) were powdered (<250 micron) and subjected to Rock-Eval type bulk flow pyrolysis using a HAWK instrument (Wildcat Technologies, Texas). Sets of one control-standard (in-house standard) and one blank were run every 10 samples to ensure instrument stability.
Sea ice biomarkers
A small number of common pan-Arctic diatoms belonging to the Haslea and Pleurosigma genera are known to produce the sea ice biomarker IP25—a mono-unsaturated highly branched isoprenoid (HBI) alkene, biosynthesized in the sea-ice matrix21,22 and deposited in marine sediments following ice melt. Due to its source-specificity and good preservation potential in marine sediments, IP25 constitutes direct evidence for seasonal sea ice. In marine settings, high sedimentary IP25 content generally reflects increasing seasonal sea ice concentrations, whereas the absence of IP25 can either indicate perennial sea-ice cover or open-water23. As such, downcore changes in sedimentary IP25 fluxes from a given location can be interpreted to reflect temporal fluctuations in sea ice conditions. A related lipid biomarker HBI III is produced by diatoms blooming in the often ice-loaded and relatively fresher and cooler surface waters typical of the marginal ice zone (e.g. ref. 59). It has recently been shown that the relative abundances of IP25 and other HBIs remain essentially unaltered in trophic food webs and faecal pellets60. These findings suggest that source HBI distributions remain unaltered following grazing, which implies that changes in grazing efficiency do not have a significant impact on the sedimentary signature of HBIs.
Prior to analytical treatment, an internal standard (7-hexylnonadecane) was added to 0.5 g of freeze-dried and homogenised sediment. Total lipids were ultrasonically extracted (×3) using a mixture of dichloromethane (DCM: CH2Cl2) and methanol (MeOH) (2:1, v/v). Extracts were pooled together, and the solvent was removed by evaporation under a slow stream of nitrogen. The total extract was subsequently resuspended in hexane and purified through open column chromatography (SiO2). Hydrocarbons (including IP25 and HBI III) were eluted using hexane (8 mL). Procedural blanks and standard sediments were analysed every 15 samples. Hydrocarbon fractions were analysed using an Agilent 7890 gas chromatograph (GC) fitted with 30 m fused silica Agilent J&C GC columns (0.25 mm i.d. and 0.25 µm phase thickness) and coupled to an Agilent 5975 C Series Mass Selective Detector (MSD). The following oven temperature programme was used: 40–300 °C at 10 °C min−1, followed by an isothermal interval at 300 °C for 10 min. The data were collected using Chemstation and analysed using the MassHunter quantification software. IP25 and HBI III were identified on the basis of retention time and comparison of mass spectra with authenticated standards. Abundances were obtained by comparison between individual GC–MS responses relative to those of the internal standard. Biomarker data presented here are reported as fluxes to account for changes in sedimentation rate.
Changes in marine primary production
Diatom fluxes were used to infer changes in marine primary production. For diatom quantification, sediment samples were treated with hydrogen peroxide (H2O2, 30%) and hydrochloric acid (HCl, 10%) to remove the organic material and carbonate, respectively. Residues were then rinsed several times with distilled water. A known volume of the final residue, homogenised in suspension, was added to a coverslip. Once the samples were completely dried, microscopy slides were mounted in Naphrax® for observation. Quantification of diatom valves was done using an optical microscope (Olympus BX43) with phase contrast optics at a magnification of 1000x. Concentrations were calculated based on the surface area of the slide that was analysed. Fluxes were calculated by combining diatom concentrations (ind. g−1) with mass accumulation rates (g cm−2 yr−1).
Lake sediment record
Sediment core NOW25c was collected from a lake at Annikitsoq on the Cape York Peninsula (76°2.100′ N, 67°36.540 W, 8.1 m a.s.l.) on July 30th 2015. The sediment core was recovered from 34 m water depth using a highly portable piston corer specially adapted for remote location use. The 177cm-long core was kept upright and drained of water, the sediment surface was secured by packing the core top with a rigid foam block (Oasis), and the core was kept cool and dark before and during transport from Greenland to Denmark. Fieldwork was conducted in accordance with local laws and permits. Sediment material is deposited at Aarhus University and available upon reasonable request from the last author (TAD).
The lake was not stratified at the time of sampling (July 30–August 2 2015) and was partially ice-covered on the day of arrival. The surface waters were 4 °C and oxygen saturation was over 100% all the way to the lakebed.
The age control of NOW25c was attained by 10 accelerator mass spectrometry (AMS) radiocarbon dates at Aarhus AMS Centre (AARAMS), Aarhus University, 9 on terrestrial moss remains and 1 on humic extraction of a bulk sample (Supplementary Table 1). The radiocarbon ages of the samples were converted into calendar years using the IntCal13 calibration curve55. The age model was calculated using the R routine BACON54. Ages reported and used in the figures are median modelled ages converted to calibrated years before 2000 (cal yrs b2k).
The Loss on Ignition (LOI) technique61 was used to determine the organic matter content at a 1 cm resolution. The sediment was dried to calculate water content and then heated to 550 °C for two hours and reweighed to calculate the percentage organic matter.
The lake sediment core was split along its length then placed in an ITRAX core scanner to obtain high-resolution pictures and measure micro-XRF. The XRF scans were made at the Aarhus University core scanning facility with a molybdenum tube set at 30 kV and 30 mA with a dwell time of 4 s. Prior to analysis, the sediment surface was flattened and covered with a 4 µm ultralene film. A step size of 0.1 mm was selected to capture possible elemental variations even in small laminations. Count readings for less abundant elements, such as Ti, maybe too low with a 4 s dwell time, so counts were summed to at least 1 mm for analysis and presentation.
Tracking the presence, absence, and relative abundance of little auks
We used a combination of δ15N (stable isotope of nitrogen), the ratio of Cadmium (Cd) and Titanium (Ti) concentrations in the lake sediments, concentrations and fluxes of cholesterol and β-sitosterol and diatom assemblage composition changes to assess the presence/absence and relative size of the adjacent little auk colony through time. δ15N differs markedly between marine and freshwater systems and has been shown to provide an unequivocal signal of marine-derived nutrients (MDN)13. Cd is also more concentrated in the marine system and thus can be used to trace seabird influence62. Cd is abundant in the seabird excrement, but it may also be present in the lake’s catchment and so the ratio Cd:Ti is used, as Ti is a proxy of catchment input to the lake. The development of soil, peat and permafrost in the catchment stores C, N and perhaps to a lesser extent Cd, reducing the quantities delivered to the lake and affecting isotope fractionation. The extremely high values of δ15N and cholesterol at the time of bird arrival likely reflect the absence of soil in the catchment and the unimpeded flow of bird-derived compounds into the lake. In the event that the inputs of C, N and Cd are reduced due to reduced supply as bird numbers fall, the catchment has the potential to act as a source. Thus, whilst the isotopes of N give a clear signal of bird arrival, they have the potential to be less reliable in tracking absence as a result of their storage and subsequent release. This may depend on the residence time of the lake and the degree of flushing by non-bird impacted waters. The combination of total cholesterol flux and the ratio of cholesterol to cholesterol plus β-sitosterol has been used to identify marine bird influence in freshwater systems in the Arctic (e.g. ref. 63).
Marine zooplankton, upon which the little auks feed, are rich in cholesterol but contain virtually no β-sitosterol. In contrast, terrestrial and freshwater primary producers contain a high proportion of β-sitosterol64. We found very high cholesterol concentrations (1497 µg g−1) and low β-sitosterol concentrations (15 µg g−1) in little auk excrement, contrasting with goose excrement collected in the same region (16 µg g−1 of cholesterol and 159 µg g−1 of β-sitosterol), a species which largely grazes on terrestrial vegetation. Thus, when seabirds are abundant, cholesterol values in the lake sediments can be expected to be high, and the ratio of cholesterol to cholesterol plus β-sitosterol should also be high. In the NOW region, the ratio is 10-fold higher in little auk excrement (0.99) compared with goose excrement (0.09). This ratio of cholesterol to cholesterol plus β-sitosterol has been previously used as an index of seabird influence, with a value of around 0.6 indicating that the majority of the cholesterol originates from seabirds64. This index was derived from studies of fulmar colonies in relatively oligotrophic systems, which is not the case at Annikitsoq, where the lake is situated in extremely eutrophic systems (due to the little auk colony) with lush vegetation in the catchment and freshwater algae thriving in the water column, both of which are sources of β-sitosterol. Here, the arrival of little auks, and with them large quantities of MDN, transformed nutrient availability, prompting a period of exceptionally high terrestrial and aquatic biological productivity. This is demonstrated by the fact that 2 m of peat accumulated in the catchment adjacent to the lake in just 1000 years, following little auk colonisation65, and the extremely high LOI values recorded in the lake (75% instead of <10% as expected for a High Arctic lake). Therefore, in addition to the large input of marine-derived cholesterol, there was also the input of cholesterol of terrestrial and freshwater origin, which has a higher proportion of β-sitosterol. Thus, a slight drop in bird input combined with increased terrestrial production, as nutrient levels remain sufficient, would result in a lower ratio, even when birds may still be present. Thus, whilst the ratio is still a useful indicator (especially when <0.4), we have not used a cut-off value signalling the dominance of seabird influence but instead simply present the index value over time. Bird absence is most likely when both cholesterol and the ratio of cholesterol to cholesterol plus β-sitosterol are low.
Sterol analysis
Sterol analysis followed standard protocols66. Specifically, Androstanol (0.1 mg mL−1) was added as an internal standard to each sample of approximately 0.5 g of dried, homogenised sediment. Lipid compounds were extracted with solvents (DCM:MeOH, 3:1) using Microwave-Assisted Extraction, saponified and separated into neutral and acid fractions using aminopropyl SPE columns. The neutral fraction of each sample was then separated using silica gel column chromatography. The sterol fraction was trimethylsilylated using N,O-bis(trimethylsilyl)trifluoroacetamide (BSTFA)/trimethylchlorosilane (TMCS) (99:1 v/v) and heated at 70 °C overnight. Excess BSTFA-TMCS was removed by drying gently under nitrogen. Samples were dissolved in 50–100 µl of ethyl acetate prior to gas chromatography-flame ionisation detection (GC-FID) and gas chromatography–mass spectrometry (GC–MS) analysis. GC–MS analyses were performed on an Agilent 7890B GC injector (280 °C) linked to an Agilent 5977B MSD in full scan mode (50–600 amu s−1). Separation was performed on an Agilent fused silica capillary column (HP-5, 60 m, 0.25 mm ID, 0.25 um df) with Helium as a carrier gas. Sterol derivatives were analysed using the following temperature programme: 50 °C (held for 2 min) to 200 °C at 10 °C min−1, then to 300 °C at 4 °C min−1 and held for 20 min. GC–MS peaks were identified through comparisons with known mass spectra (NIST08) and standards where possible. Analytes were quantified based on internal standards. For 13 of the samples, lipids were extracted from the dry sediment using a chloroform:methanol 2:1 mixture and then sonicated for 10 min, after which 0.75 mL of distilled water was added. Lipids were fractionated into neutral lipids (NLs; including sterols), glycolipids, and phospholipids (PLs) using a Bond Elut (0.5 mg) silica cartridge. First, the resin of the cartridges was conditioned using 5 mL of chloroform. Subsequently, the total lipids (1 mL) were applied to the resin, rinsed using chloroform, and then the NLs (including sterols) were collected under vacuum using 10 mL of chloroform. Sterols from the NL fraction were silylated with BSTFA, TMCS, and pyridine at 70 °C for 1 h. Trimethylsilyl (TMS) derivatives of sterols were analysed with GC–MS (Shimadzu) and GC-FID (Shimadzu) equipped with a Phenomenex (USA) ZB-5 Guardian column (30 m × 0.25 mm × 0.25 μm). Cholesterol and β-sitosterol were identified using characteristic ions of GC–MS runs64 and quantified with GC-FID using authentic standard solutions of plant sterol mixture from Larodan (including 53% β-sitosterol, 7% stigmasterol, 26% campesterol, 13% brassicasterol), and cholesterol from Sigma-Aldrich. The recovery percentage of the sterol samples was calculated using 5-α-cholestane (Sigma-Aldrich) as an internal standard.
Stable isotope analyses
For stable isotope analysis, samples were taken at 2 cm intervals, freeze-dried for 48 h and ground into fine powder. The samples were packed into tin cups and analysed at UC Davis Stable Isotope Facilities, California, USA. Here, carbon (13C) and nitrogen (15N) isotope analyses were conducted using an elemental analyser and a continuous flow isotope ratio mass spectrometer (IRMS). Specifically, an Elementar Vario EL Cube (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced with an Isoprime VisION IRMS (Elementar UK Ltd, Cheadle, UK). Samples were combusted at 1080 °C in a reactor packed with chromium oxide and silvered copper oxide. After combustion, a reduction reactor trap removed oxides and a helium carrier then flowed through a water trap (magnesium perchlorate and phosphorous pentoxide). CO2 was held in an adsorption trap until the N2 peak was analysed; after which the CO2 was released by heating to the IRMS. Reference materials included: IAEA-600, USGS-40, USGS-41, USGS-42, USGS-43, USGS-61, USGS-64, and USGS-65. A sample’s isotope ratio is expressed relative to a reference gas peak analysed with each sample. These provisional values are finalised by correcting the values for the entire batch based on the known values of the included laboratory reference materials. The long-term standard deviation was 0.2 per mil for 13C and 0.3 per mil for 15N. The delta values are expressed relative to international standards VPDB (Vienna Pee Dee Belemnite) and Air for carbon and nitrogen, respectively.
Diatom analyses of the lake record
Diatom analyses were carried out on 77 samples at a resolution of 1–3 cm, covering the entire lake record. Sediment samples were cleaned using H2O2 and HCl, and permanent slides were prepared using Naphrax®. A minimum number of 400 valves were counted per sample, and relative species abundances were calculated as percentages of the total counts in each sample.
Modern sea ice concentration in the NOW region
We assessed mean sea-ice concentration during the months May–September (MJJAS) in two regions of interest—NOW and Ice Arch (Fig. 7). We used the gridded sea ice concentration and extent product based on satellite observations during 1979–2015, available from the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration67.
Numerical methods
Monte Carlo simulations and principal components analysis
To derive principal components of the marine and lake records, we first performed Monte Carlo simulations of the constituent time series from each record (Supplementary Figs. 4–6). We performed 10,000 simulations of both the marine and lake records that acknowledged both measurement and dating uncertainties of the constituent time series. The 10,000 Monte Carlo simulations yield 2D frequency histograms depicting the probability density of a given measurement at a given time for each time series. We took the median (50th percentile) as the ‘most likely’ time series for each variable and used the 5th and 95th percentiles to define a 90% confidence envelope around this ‘most likely’ time series. For the marine record, each of the three constituent time series analysed (Diatom flux, Chaetoceros spores flux, and HBI III) were assumed to have a measurement uncertainty of ±7.5% (5–10%) and the depth-dependent dating uncertainty shown in Fig. 2, assumed to represent 2-sigma uncertainty. For the lake record, the four constituent time series were assumed to have measurement uncertainties of 2% for Cd:Ti, 0.5 per mil for δ15N, and 0.001 mg cm−2 yr for fractional sterol, with depth-dependent dating uncertainty (as shown in Fig. 3) assumed to represent 3-sigma uncertainty. To ensure that sterol was not over-represented as two independent variables in the subsequent principal component analysis, we employed fractional sterol, which is sterol flux times sterol ratio. Similarly, for the fourth constituent series, lake diatoms, we used a weighed sum to combine the two principle components (resulting from analyses carried out on centred and standardised, square root transformed relative abundance data) into a single index. This weighed sum is proportional to the per cent variance explained by each principal component (34% for PC1 and 27% for PC2) and assumed an uncertainty equivalent to 10%.
We performed separate principal component analyses on the median time series of the three constituent marine time series and the four constituent lake time series and only considered the resulting first principal component. We propagated the Monte Carlo uncertainty envelopes as fractional mean uncertainties bounding the first principal components of both the marine and lake records (Supplementary Table 2). In the marine record, Diatoms and Chaetoceros spores both load positively on the first principal component (PC1 loadings of 0.63 and 0.62, respectively), while HBI III loads negatively (PC1 loading of −0.47). We interpret this as suggesting that the first marine principal component is positively correlated with primary production. In the lake record, all-time series (Cd:Ti, δ15N, sterols and diatoms) load positively on the first principal component (Supplementary Table 2). We interpret this as suggesting that the first lake principal component is positively correlated with little auk abundance. We do not consider principal components beyond the first but acknowledge that they may also contain climate and/or population signals.
When plotting both the first marine and lake principal components (Fig. 4), we included the respective Monte Carlo uncertainty envelopes, which have been propagated as the fractional mean uncertainties of the constituent time series underlying each principal component. The Y-envelope corresponds to uncertainty in principal component magnitude, while the X-envelope corresponds to uncertainty in time. The temporal uncertainties we discuss are associated with these latter horizontal uncertainty envelopes.
General additive models for detection of significant changes
To identify points of significant change in our proxy data over time, we used general additive models (GAM), and took a flexible approach to the degree of smoothing employed in modelling the response of both the marine and lake record proxies. Where possible, the degree of smoothing was estimated using restricted maximum likelihood (REML), and the model contained a continuous-time first-order autoregressive process (CAR(1)) to account for temporal autocorrelation. This solution was optimal for the marine record proxies Diatoms, Chaetoceros spores, IP25 and HBI III, and for δ15N and the diatom PC1 and PC2 axes scores from the lake record. However, for Cd:Ti and cholesterol from the lake record, this approach produced an over-smoothed model, and so in these cases, we used a generalised cross validation (GCV) approach, which optimises predictive accuracy and allows for heteroscedasticity in the data. After the model was fitted (Supplementary Figs. 7–15), posterior simulation was run with 20 random draws from the posterior distributions of the fitted GAMs to reflect the degree of certainty of the model at a given time (Supplementary Figs. 7–15), and to allow for the calculation of confidence intervals. Finally, first derivatives (black lines) and associated 95% confidence intervals were estimated for the GAM trend of each proxy. Using this approach, a significant change is demonstrated when the confidence interval of the first derivative does not include 0 (Supplementary Figs. 7–15)68. GAM analyses were carried out in R version 4.0.069, using the package mgcv70,71 and additional code outlined in ref. 68.
Source: Ecology - nature.com
