Holocene environmental change in Rotsee and its impact on sedimentary carbon storage

To assess the long-term impact of climate change and human influence on lakes and their sedimentary carbon storage, paleo-environmental approaches using well-dated lake sediment cores can be employed. Here, we reconstruct carbon mass accumulation rates for organic and inorganic carbon since 13 ka BP in Rotsee, a perialpine lake near the Swiss Alps, using a 12m sediment core. A multiproxy approach (XRF, carbon and nitrogen isotopes, organic macromolecule chemical compositions, aDNA) was used to explore changes in the lake system that affect sedimentary carbon storage. The Early Holocene (11.8 to 7 cal ka BP) was characterized by a mixed phytoplankton and watershed-derived provenance of organic matter, and the deposition of inorganic and organic sedimentary carbon. Warming during the Holocene Thermal Maximum (9.8 to 8.8 cal ka BP) increased sedimentary carbon storage. In the mid-to-late Holocene (7 to 1 cal ka BP), the sedimentary record indicates an increased influx of allochthonous, vascular plant-derived organic matter, and low production or conservation of phytoplankton-derived carbon. Organic carbon storage increased, while inorganic carbon became negligible. Larger deforestation events, potentially during Neolithic times (around 4 ka BP), but especially during Roman times (2 ka BP), coincided with further increased organic carbon MARs. Recent sediments, influenced by eutrophication in the last century, show higher carbon accumulation rates compared to earlier Holocene periods. Rotsee serves as a case study of how climate warming and human land use changes have influenced lake development and sedimentary carbon storage, with broader implications for understanding carbon dynamics in high-altitude lakes and their future carbon balance.

Methods
Core collection and on-site subsampling

Three short cores and two long cores were collected at a water depth of 5.5 m (47°04’27.81”N, 8°19’25.7”E WGS 84; 667230/214087 LV95) between 03/10/2021 and 05/10/2021. The short cores (40-60 cm long) were collected from a vessel using a gravity corer with clear plastic liners (UWITEC; inner liner diameter 90 mm). The two long cores (sections ROT21-1-1 to ROT21-1-5 and ROT21-1-6 to ROT21-1-9) were taken 4 m apart from each other, using a shoreline moored platform using a piston coring system with a manual hammer, without a re-entry cone (UWITEC; inner liner diameter: 59.5 mm). Long cores were taken in sections of 3 m (except for 1 section that was 2 m). The recovery of two parallel long cores, vertically offset by one meter, was necessary to obtain a high-quality, complete sedimentary sequence. All short and long cores were brought onshore for sampling immediately after core recovery. Short cores were maintained in vertical position and sampled by extrusion, whereas long cores were first accessed horizontally through ‘windows’ that were cut into the core liner. From each depth sampled, sediments for determination of porosity and bulk density, DNA analyses and macromolecular organic matter analyses were collected with sterile cut-off syringes. Samples for DNA extraction were immediately frozen in liquid N2, before storage at -80 °C, whereas samples for organic matter analysis were frozen and stored at -20 °C. Afterwards, the core sections were split lengthwise before subsampling 2 cm slots for bulk carbon and nitrogen analyses using solvent cleaned spatulas. One core half (the so-called archive half) remained intact for imaging and XRF scanning.

XRF scanning

Elemental compositions were measured at 5 mm resolution using a µXRF core scanner (Avaatech XRF) with an Oxford 100 W X-Ray tube and a rhodium anode equipped with a Canberra X - PIPS/DSA 1000 (MCA) detector. Prior to analysis, core surfaces were flattened to ensure uniform scanning and covered with 4 µm Ultralene foil. Elemental groups with lower energy levels were measured at 10 kV (1500 A, no filter, 15 s exposure), while mid-energy groups were measured at 30 kV (2000 A, Pd thin filter, 40 s exposure). Prior to determining the variability in XRF parameters (excluding Mo, Ar and coherent and incoherent scatter) using a principal component analysis (Supp. Fig. 1), the cps counts were transformed using a centered log-ratio transformation (Bertrand et al., 2024) and scaled. Based on untransformed cps counts (Supp. Fig. 2), selected XRF log-ratios were calculated.

Dating and age model

The top 50 cm of a short core was sectioned at 1 cm resolution and used for 210Pb and 137Cs dating (Fig. 1A; Supp. Table 1A). 137Cs peaks were determined to date the sediment layers deposited in 1987 and 1963 Anno Domini (AD). In addition, radiocarbon dating on 19 macrofossils, including 12 seeds, leaf remains and twigs of terrestrial plants, and 7 macrodetrital remains of aquatic macrophytes was performed (Fig. 1B, Supp. Table 1B). After wet sieving, macrofossils were subjected to an acid-alkali-acid treatment at room temperature (Norris et al. 2020) to remove carbonates, acid soluble humic material, and humic acids. At two depths, bulk sediments were acidified using fumigation (described in Haas et al. 2019) and weighed in for 14C dating, with the aim of constraining the reservoir age during the Younger Dryas (Supp. Table 1C). The reservoir age was used to correct the uncalibrated 14C ages measured on the aquatic macrophytes. 14C measurements were carried out on an Accelerator Mass Spectrometer (AMS) with an Elemental Analyzer unit (EA) using the MIni CArbon DAting System (MICADAS) at the Laboratory for Ion Beam Physics of the Eidgenössische Technische Hochschule (ETH) in Zurich. An age-depth model was created using rplum, which allows unsupported 210Pb values, 137Cs ages and uncalibrated 14C ages to be combined based on Bayesian statistics (Blaauw and Christen 2011). Radon measurements are available as estimates of supported 210Pb, assumed constant by the model. A memory strength of 10 and memory mean of 0.5 is used. In this model, 14C ages are converted to calendar ages using the INtCal20 calibration curve (Reimer et al. 2020).

Bulk nitrogen and carbon content and isotopes

Sediments were freeze-dried and homogenized using an agate mortar and pestle. Total nitrogen (%; TN) and δ15N-TN values were determined on between 3 to 200 mg of unacidified sediments using an EA-IRMS system composed of a Vario Pyro Cube coupled to a Isoprime 100 IRMS (Elementar, Germany). Repeated measurement of reference materials Acetanilide #1 (Schimmelmann, USA, δ15N = +1.18 ± 0.02) and High Organic Sediment Standard (HEKAtech, Germany, δ15N = +4.32 ± 0.20) were used to determine the instrument precision, which was determined to be generally below 0.05 ‰ δ15N for the Acetanilide standard, and below 0.07 ‰ δ15N for the sediment standard. Offsets between the measured and known δ15N values of the Acetanilide standard (average offset = 0.19 ± 0.08) were used to correct the δ15N-TN values of the bulk sediments. The contents of total carbon (TC), total organic (TOC) and total inorganic (TIC) carbon were measured on 50 mg of sample weighed into a ceramic crucible, on the SoliTOC® Cube (Elementar Analysensysteme, Germany). The reported TOC is the summed amounts of TOC400 and refractory organic carbon (ROC), with TOC400 determined at 400 °C and ROC between 400 °C and 600 °C, and TIC between 600 °C and 900 °C. Low (B2152) and high organic carbon content standards (B2150) from Elemental Microanalysis (United Kingdom) were run with each batch to determine the instrument accuracy, which was determined to be 98.9 ± 0.6%. δ13C-TOC of bulk sediments was measured on an EA-IRMS system, EA Vario Pyro Cube (Elementar Analysensysteme, Germany) and Isoprime IRMS (GV Instruments, UK), after acidification. For acidification, samples with inorganic carbon were subjected to a 1N HCl treatment until no more visible reaction occurred. To calibrate the instrument an Acetanilide #1 (Schimmelmann, USA, δ13C = -29.52 ± 0.01) standard was used, as well as a High Organic Content Sediment (SA990894; δ13C = -28.85 ± 0.10) and Low Organic Soil (SA33802152; δ13C = -22.88 ± 0.10) standards from Hekatech Analytics. In general, instrument precision during the runs was below 0.06 ‰ δ13C for Acetanilide and below 0.16 ‰ δ13C for the sediment and soil standard, and an accuracy better than  0.02 for δ13C for Acetanilide and 0.1 ‰ δ13C for the sediment and soil standard. No corrections of the δ13C values were performed.

Bulk compound classes and hydrocarbons

To determine the OM macro-molecular composition, pyrolysis gas chromatography mass spectrometry was used, following the set-up as described in Gajendra et al. 2023. Between 100 - 500 mg of freeze-dried sediments were weighed into autosampler containers (Eco-cup SF, Frontier Laboratories, Japan) and pyrolyzed at 450 °C and 650 °C, according to Tolu et al. (2015). The pyrolizer was connected to an autosampler (PY-2020iD and AS-1020E, FrontierLabs, Japan) and to a GC/MS system (Trace 1310, Thermo Scientific and ISQ 7000, Thermo Scientific) equipped with a DB-5MS capillary column (30 m x 0.25 mm i.d., 0.25 mm film thickness; J&W, Agilent Technologies AB, Sweden). Chromatograms were analyzed in R (version 2.15.2, 64 bits) based on a pipeline that identifies common mass fragments as pyrolysis products (Tolu et al. 2015). Pyrolysis products were then classified and annotated according to Tolu et al. (2015) and Ninnes et al. (2017). On average 27% of the total peak area occurred in peaks that didn’t provide conclusive structural information. Areas of individual compounds within each compound class were summed up (Supp. Table 2 for identity of individual compounds), and compound classes expressed as relative abundances (peak area of each compound class relative to total characterizable Py-GC/MS peak area). To summarize the main trends in compositional variability, a PCA was performed based on the standardized fractional abundance of the compound classes (Supp. Fig. 3).

Mass accumulation rates

Dry bulk density values, the mass (weight) of the dry solids divided by the total volume of the wet sample, were measured on 7 mL of sediments sampled with a cut-off syringe, based on weights before and after drying (n = 68). Using the linterp command from the astrochron package (Meyers 2014), the bulk dry density values were afterwards interpolated at a 1 cm resolution. Mass accumulation rates (MAR) were then calculated by multiplying the interpolated dry bulk density with measured weight percentages of TOC, TIC, and normalized per year, using a smoothed sedimentation rate (autoplot, smoothing with a smoothing width of 800, using the astrochron package; Meyers, 2014). Supp. Fig. 4 shows the variability of measured and interpolated parameters that are used to calculate the MAR values through time.

aDNA analysis

Sedimentary DNA was extracted according to Lysis Protocol I of Lever et al. 2015. Briefly, 0.2 g of sediment was added to 2-mL screw-cap tubes filled to 15 % with 0.1 mm Zr beads. For the vast majority of samples, 100 µl of 10 mM adenosine triphosphate (ATP; prepared by dissolving adenosine 5’-triphosphate disodium trihydrate in molecular grade water) solution was added to reduce DNA sorption. The only exceptions were deep glacial clay layers, in which recovery was significantly enhanced by increasing the ATP concentration to 300 mM. 0.5 ml of lysis solution I was added to all samples (Lever et al. 2015). Samples were then shaken for one hour at 50 °C (600 rpm; ThermoShaker, Eppendorf), and subsequently washed twice with cold 24:1 chloroform-isoamyl alcohol and precipitated with linear polyacrylamide (20 µg ml-1), 5 M sodium chloride and ethanol. The pellets were dried using a vacuum centrifuge (50 °C; Thermo Fisher Scientific, USA), resuspended in molecular grade water and purified with the CleanAll DNA/RNA Clean-up and Concentration Micro Kit (Norgen Biotek Corp., Canada; Protocol A). All extracts of samples and extraction negative controls (extraction blanks) were stored at -80 °C. Eukaryotic 18S rRNA and rbcL genes were quantified by real-time PCR (Lightcycler® 480; Roche) with SYBR®Green as dye. Eukaryotic 18S rRNA genes were amplified using the All18S primer pair (Deng et al. 2020), whereas chloroplast genes encoding the large subunit of ribulose-1,5bisphosphate carboxylase (rbcL) were quantified using published assays for vascular plants (Willerslev et al. 2003), Chlorophyta and Ochrophyta (both Han et al. 2022).