Past climate variations recorded in needle-like aragonites correlate with organic carbon burial efficiency as revealed by lake sediments in Croatia

The drivers of organic carbon (OC) burial efficiency are still poorly understood despite their key role in reliable projections of future climate trends. Here, we provide insights on this issue by presenting a paleoclimate time series of sediments, including the OC contents, from Lake Veliko jezero, Croatia. The Sr/Ca ratios of the bulk sediment are mainly derived from the strontium (Sr) and calcium (Ca) concentrations of needle-like aragonite in Core M1-A and used as paleotemperature and paleohydrology indicators. Four major and six minor cold and dry events were detected in the interval from 8.3 to 2.6 calibrated kilo anno before present (cal ka BP). The combined assessment of Sr/Ca ratios, OC content, carbon/nitrogen (C/N) ratios, stable carbon isotope (δ13C) ratios, and modeled geochemical proxies for paleoredox conditions and aeolian input revealed that cold and dry climate states promoted anoxic conditions in the lake, thereby enhancing organic matter preservation and increasing the OC burial efficiency. Our study shows that the projected future increase in temperature might play an important role in the OC burial efficiency of meromictic lakes.

The drivers of organic carbon (OC) burial efficiency are still poorly understood despite their key role in reliable projections of future climate trends. Here, we provide insights on this issue by presenting a paleoclimate time series of sediments, including the OC contents, from Lake Veliko jezero, Croatia. The Sr/Ca ratios of the bulk sediment are mainly derived from the strontium (Sr) and calcium (Ca) concentrations of needle-like aragonite in Core M1-A and used as paleotemperature and paleohydrology indicators. Four major and six minor cold and dry events were detected in the interval from 8.3 to 2.6 calibrated kilo anno before present (cal ka BP). The combined assessment of Sr/Ca ratios, OC content, carbon/nitrogen (C/N) ratios, stable carbon isotope (δ 13 C) ratios, and modeled geochemical proxies for paleoredox conditions and aeolian input revealed that cold and dry climate states promoted anoxic conditions in the lake, thereby enhancing organic matter preservation and increasing the OC burial efficiency. Our study shows that the projected future increase in temperature might play an important role in the OC burial efficiency of meromictic lakes.
Lakes have disproportionally large annual amounts of buried organic carbon (OC) compared to oceans 1 and are of great importance for the global carbon budget and cycle; thus, they might also have a vast impact on climate change in the future. However, the influence of climate on OC burial efficiency is still not precisely understood. The influence of temperature on the OC burial efficiency has been studied extensively, mainly for lakes at higher northern latitudes [2][3][4][5][6] . Although a large number of studies have focused on this topic, a consensus has not been reached on the possible influence of temperature on the OC burial efficiency. Moreover, the main driving mechanisms that are directly responsible for the observed changes in OC content of lake sediments have not been identified. Studies have inferred that higher temperatures enhance OC burial mainly as a consequence of denser vegetation cover 3,7-10 , with an increase in temperature positively correlated with OC mineralization 2,9,11 . In more recent studies, anthropogenic influences have been reported, primarily through the role of reactive nitrogen and phosphorus on OC burial efficiency 4,6,12 , and the temperature effect has been negated. Furthermore, oxygen exposure time and redox conditions in the water may also play a prominent role in the OC burial efficiency 5,9,13,14 .
The majority of previous studies on OC burial efficiency have focused on recent and subrecent lake sediments, where climate effects can easily be concealed by anthropogenic influences. Additionally, in latitudinal studies, the temperature effect can be masked by other variables, such as changes in vegetation cover and/or precipitation rate. To infer the climate influence on OC burial efficiency, we studied sediments of Lake Veliko jezero, Croatia ( Fig. 1), for the period 8.3 to 2.6 cal ka BP. In detail, we interpreted the OC content in the context of high-resolution relative paleoclimate (log ratio of Sr/Ca) and paleoredox proxies (log ratio of Mo/detrital elements) and as an indicator for past aeolian activity (log ratio of Zr over Al). Our approach for correlating paleoclimate and paleoredox proxies was performed by considering the properties of compositional data 15 , which enabled more precise interpretations. The Sr/Ca ratios of the bulk sediment are mainly derived from Sr and Ca concentrations of needle-like aragonite in Core M1-A and used as paleoclimate, temperature proxies, with higher ratios www.nature.com/scientificreports/ white and dark laminas mainly composed of aragonite and organic matter. The top part of that unit consists of a few centimeters of oxidized lake sediment. The following unit from 241 to 343 cm is marsh to shallow lake sediment. Finally, the bottom unit (343-417 cm) is terra rossa-type soil interbedded with thick tephra. Core M2 resembles Core M1-A with the difference that unit boundaries are at different depths. Therefore, the studied interval spans from 127 to 266 cm in this core.

M1-
High-resolution XRF scanning and age-depth model. The Sr/Ca records were obtained by highresolution XRF core scanning at split core surfaces at 1 cm and 2 mm resolutions. We compared the Sr/Ca data and ratios between the two investigated cores (M1-A and M2) to confirm the robustness of our record and exclude any potential analytical artifacts. M1-A core chronology is based on four tephra layers and four C-14 dates 25 . Briefly, the studied interval chronology is very well constrained since in ̴ 1.25 m, there are three tephra layers, which serve as depth-age model anchor points: two charcoal C-14 dates free of reservoir effects and one marine shell C-14 date with known reservoir effects (Fig. 2). The geochronological age of Core M2 is based on one plant remaining sample, which was dated by radiocarbon methodology. Three additional data points were included in the age-depth model of this core. Two of these ages are based on visible tephra layers correlated to the known volcanic eruptions previously recognized in Core M1-A 25 . The third date corresponds to a time of marine intrusion into the lake, which is marked by a sharp boundary between the top and underlying lithological units ( Fig. 2) in both cores. It was radiometrically dated in Core M1-A 25 . Individual radiocarbon dates were calibrated using the R package rbacon 34 (the age-depth model for core M2 is presented in the Supplementary file, Figure  S4). Despite the lower resolution of the Core M2 record, it is evident that major peaks and troughs in Sr/Ca are well presented in both datasets (Fig. 3).
Sr/Ca ratio as paleoclimate proxy. Our working hypothesis was that the Sr/Ca ratio in bulk sediment almost exclusively reflects the Sr/Ca ratio of needle-like aragonite in the sediment (Fig. 2). The temperature and Mg/Ca ratio are the two most important variables for the precipitation of inorganic aragonite 35 . The main source of Mg for the aragonites from the Veliko jezero lake could be seawater, and the second source could be from the surrounding carbonate rocks, which are mainly dolomites. Although the Mg/Ca ratio does not affect the crystallization of aragonites directly, it prevents the formation of calcite and the transformation of aragonite into calcite [36][37][38] . Additionally, the Sr/Ca ratio in aragonite is not affected by variations in salinity, sulfate or carbon dioxide content 37,39,40 . Needle-like aragonite currently precipitates during late spring and summer in the adjacent Lake Malo jezero 41 . Since the Sr/Ca ratio of needle-like aragonite is largely temperature dependent 17,18 , we assume that the Sr/Ca ratio of the studied sediment could be used as a relative paleotemperature proxy. Mineralogical analysis performed with X-ray diffraction (XRD) and scanning electron microscopy coupled with www.nature.com/scientificreports/ energy dispersive spectroscopy (SEM-EDS) proved that inorganic needle-like aragonite is the main mineral and nearly the only carbonate phase in our investigated samples (Fig. 2). The only exception is an interval from 201 to 214 cm, where inorganic rhombohedral Mg-calcite is the main carbonate phase, thus confirming previously published results 33 . According to the age model of Core M1-A, the occurrence of Mg-calcite coincides with known time intervals of wet climate during which lake levels were high and generally correspond to pluvial periods observed in the wider Mediterranean region (Figs. 3 and 4, from ca. 7.6 to 7 ka BP) [28][29][30][42][43][44][45][46][47] . We propose that during this period, increased freshwater input lowered the Mg/Ca ratio in the lake water, leading to Mgcalcite precipitation and hindering the precipitation of aragonite. This consequently led to a Sr/Ca decrease in our records because Mg-calcite incorporates Sr in the crystal lattice less effectively than aragonite 48 ; thus, the Sr/ Ca ratio of bulk sediment cannot be interpreted as a relative paleotemperature proxy in this interval. However, the predominance of Mg-calcite over aragonite points to wetter climate conditions, which are also observed in the wider region during this time period [28][29][30][42][43][44][45][46][47] .
To utilize the Sr/Ca ratio as a relative paleotemperature proxy, we have proven that detrital Sr and Ca inputs are negligible. XRD analyses revealed aragonite as the main mineral, with minor quartz content limited only to the oldest portion of the studied period. Furthermore, large variances of the centered log ratio transformed variables of Sr and Ca relative to detrital elements, such as titanium (Ti) and aluminum (Al) (Supplementary file, Figure S2) combined with small logratio varaince of Sr and Ca relative to the inorganic carbon (INC) (Supplementary file, Figure S2) additionally confirmed that Sr and Ca represent carbonate components and are not of detrital origin.
Finally, the Sr budget of the Lake Veliko jezero water and consequently the Sr/Ca ratio of the bulk sediment can also be influenced by hydrological variability. Changes in hydrological regime would theoretically affect the relative marine influence at this location because of a limited connection of the lake to the Adriatic Sea through permeable karst: ocean water is characterized by higher Sr concentrations than freshwater 49 . Two possible scenarios emerge if hydrologically induced Sr availability was the limiting factor for the Sr/Ca ratio of the bulk sediment. MD 90-917 28  www.nature.com/scientificreports/ First, during cold periods, Sr/Ca in the lake water would be lower because of the reduced evaporation rate, i.e., decreased marine influence, which would finally result in a relative decrease in the Sr concentration of the lake water and consequently the Sr/Ca of the lake sediments. The opposite would be the case for warmer climate conditions. This scenario, however, can be discarded based upon the good correlation of the maxima in our Sr/ Ca record with cold events that were previously observed in multiple paleoclimate archives in the Mediterranean region (Fig. 3). In Fig. 3, the locally estimated scatterplot smoothed (LOESS) Sr/Ca curve, with a smoothing factor of 0.09, displays four distinct peaks centered at 8.3, 6.0, 4.25 and 2.9 cal ka BP. The first Sr/Ca maxima in our record, which are centered at 8.3 ka, are coeval with cold events described in the pollen record from the Adriatic Sea 28 and a decrease in the sea surface temperature 29 recognized in the same core. In the Alboran Sea, a drop in the temperature during that time period was recorded as well 26,50 , while a minor drop in the temperature was also observed in Gemini Lake 31 . These events correlate with the North Atlantic cold event (NAC5) 51 . Following the pluvial period (7.6 to 7 cal ka BP) characterized by Mg-calcite deposition instead of aragonite, another maximum in the Sr/Ca record from the M1-A core, which is centered at approximately 6 cal ka BP, can be correlated with cold spells recognized in the Adriatic 27,28 and Tyrrhenian and less clearly in the Alboran Sea 26 as well as the NAC4 event in the North Atlantic 51 . Temperature reconstructions based on chironomid communities from Gemini and Verdarolo lakes also indicate cold conditions at approximately 6 cal ka BP 31 . A Sr/Ca maximum at approximately 4.25 cal ka BP is correlated with a decrease in the temperature in the Adriatic Sea and Italian lakes 28,31 and in the north 51 . Finally, an increase in the Sr/Ca centered at approximately 2.9 cal ka BP correlates well with the cold spells already recognized in the Adriatic, Alboran and Tyrrhenian seas [26][27][28]50 as well as in the North Atlantic 51 . A less pronounced temperature decrease during this time interval was observed in Lake Gemini, Italy 31 .
The second scenario implies that the environment during cold events was also dry, while during warmer events, it was wetter, which would lead to a decreased relative marine influence and the dilution of the Sr budget in the lake water. Such changes would finally result in lower Sr/Ca during warm periods and higher values during cold periods. To further investigate the role of wet versus dry conditions as potential drivers of Sr/Ca variability, we studied additional chemical elements acquired by XRF scanning, such as Br, Rb, Si, K, Na, Al, Fe, Mn, Zr, Ti and Mo, which were primarily validated by wet-chemical analysis on discrete samples, e.g., through inductively coupled plasma mass spectrometry (ICP-MS). We used all of these elements to model proxies called balances; for paleotemperature, aeolian input and paleoredox conditions via orthonormal log ratio transformation (OLR) 52 , which enabled firm statistical parameters to be established for reliable interpretations of the involved processes 53 . Sr/Ca (detrended) 12 14 16 July air temperature(°C) 10 15 20 C/N  www.nature.com/scientificreports/ The methodology is fully explained in a previous study 53 and described briefly in the Methods section. Ultimately, five balances were modeled, each of which is a proxy for a certain process: b3 and b4 are proxies for paleoredox conditions, b5 is a proxy for aeolian activity, and b2 is a proxy for paleoclimate Sr/Ca (rationale for the balance construction is shown in the Supplementary file, text S1, and sign matrix table used for OLR transformation is  shown in the Supplementary file, table S1). The data analyses show that an increase in Sr/Ca correlates with increases in anoxic conditions (r (b2-b3) = 0.55, p (0.05) = 0.00001) and aeolian activity (r (b2-b5) = 0.56, p (0.05) = 0.00001) (Supplementary file, table S2) (Fig. 4). A comparison between b2 and Sr/Ca is presented in the Supplementary file, figure S3). Furthermore, a speleothembased paleohydrological reconstruction from Corchia cave in Italy 30 implies that during the time corresponding to the Sr/Ca maxima in our record, climate conditions were generally drier (Fig. 4). However, it is also evident that the Sr/Ca ratios in our record do not exhibit exceptionally high values during the most widespread dry event in the Mediterranean and wider region at 4.2 ka 27,45,[54][55][56][57] . This finding implies that hydrological conditions probably played less of a role in driving Sr/Ca ratios than the temperature effect. Following all lines of evidence, we propose that the Sr/Ca ratio of Lake Veliko jezero bulk sediment represents the Sr/Ca ratio of inorganic needle-like aragonite, which mainly reflects relative paleotemperature changes, while hydrological variability likely plays a secondary role.
One of the main advantages of XRF core scanning in addition to being nondestructive is the high resolution. If Sr/Ca ratios are a temperature indicator, then they may also be applicable for short events; for example, we observe several relatively brief cold events at 8.0, 6.6, 5.4, 4.8, 4.0 and 3.3 cal ka BP. From these, the events at 8.0, 6.6, 5.4 and 4.0 cal ka BP were also recorded in Lake Gemini while the event at 3.3 ka BP was recognized in Lake Verdarolo 31 . Most of these events are characterized by an increase in anoxic conditions and aeolian input, indicating not only cooler but also drier climate conditions at our sites (Fig. 4). These events can be detected due to the combined effects of the limited lake size and detrital influence on the studied lake(s) and the high resolution of the data. The small size of this lake (surface area of 1.44 km 2 ) has very limited detrital influence, i.e., small effects of internal and landscape filters 58 result in increased sensitivity in recording smaller-scale climate events.

Disentangling the drivers of OC burial efficiency.
To decipher the potential drivers of OC burial efficiency, we analyzed the OC, inorganic carbon (INC) and organic nitrogen (N) and explored their relationship with enhanced anoxic episodes, cold spells and aeolian input. Additionally, to better characterize the provenance of organic matter (OM), we analyzed δ 13 C throughout the same interval.
Higher C/N ratios (> 10) indicate potential mixing of land-derived and autochthonous organic matter 59,60 , while more negative δ 13 C values could also be related to the increase in the land-derived component of organic matter 61 . Our results demonstrate a slight long-term decreasing trend in the C/N record that coincides with a much stronger increasing trend of δ 13 C values. This anticorrelated pattern of the two proxies indicates that land-derived organic matter is partly decreasing, which is consistent with the overall detrital influence during our studied interval (Fig. 4). Additional evidence comes from the occurrence of quartz, which is only found at the base of the studied interval. A combination of two factors is most likely responsible for the described trends. First, lake deepening caused by the Holocene sea level rise moved the shoreline away from the core site. More specifically, the distance between Core M1-A and the Pomena doline, which is a small terra rossa soil patch adjacent to Lake Veliko jezero (Pomena field, Fig. 1), increased. Second, due to the Holocene sea level rise, a gradual increase in marine influence through permeable karst occurred. Sea level rise changed the lake biota 33 and consequently the organic carbon content and composition.
During the studied interval, only two different pollen zones occurred, one with Juniperus and Phillyrea ca. 8 to 6.5 ka BP and another with Quercus ilex from 6.5 ka BP to the present 62 . This finding implies minor changes in terrestrial vegetation with negligible impacts on carbon content, composition and variability.
The correlation of the OC content with the Sr/Ca record reflecting temperature variations on the millennial time scale suggests that the OC content increased during cold exposure (Fig. 4). The observed OC increase during cold events is in line with a previous study 2 , where a temperature decrease leads to low mineralization of OC. However, the temperature effect through OC mineralization was not a substantial factor for OC preservation in the Lake Veliko jezero sediments. If temperature is a driver of OC preservation on millennial time scales, then cold climate conditions would decrease the degradation rate of algal (autochthonous) organic matter 63 , which would result in lower C/N values, higher OC and more positive δ 13 C values, which we did not observe in our data. Indeed, we observe slightly higher C/N and more negative δ 13 C values during cold spells, which would imply that land-sourced organic matter increased during cold spells as a result of enhanced aeolian input as confirmed by correlation analyses r(b2-b5) (Fig. 4). However, we argue that an increase in land-sourced organic matter is not the main mechanism underlying the overall OC increase because a more than double increase in OC during cold spells would have caused a substantial increase in C/N, which is not observed in our record. An exception is the 4.2 event, when maximum C/N ratios occurred.
The OC amount and variability also correlate with the paleoredox proxy (Fig. 4), i.e., an increase in anoxic conditions corresponds to an increase in OC content. Based on this finding, we propose that an increase in anoxic conditions is the main factor that led to OC preservation at our site. A decrease in temperature and possibly drier conditions during cold events would shift the redox zone boundary and thermocline closer to the surface 64 , which would prevent mixing of the water throughout the water body, thereby causing the anoxic boundary to move upward and leaving the majority of the water column under anoxic conditions. This supposition is also confirmed by the high correlation between paleotemperature (b2) and paleoredox (b3) proxies. Although the depth of a thermocline depends on a number of factors, such as the lake size, dissolved organic content, temperature, wind activity, etc. 65 . We believe that the temperature was the key factor controlling the thermocline and redox zone boundary depth in Lake Veliko jezero. These findings are underpinned by a study of modern processes in Lake www.nature.com/scientificreports/ Veliko jezero lake, where the thermocline occurs only during summer months 66 . Finally, a shallow thermocline/ redox zone boundary would cause most of the organic matter produced in and transported into the lake to be prevented from decomposing, resulting in higher OC values during cold spells.

Conclusions
The data presented here demonstrate that temperature changes may have a significant impact on OC burial efficiency. The temperature decrease and likely drier climate conditions caused shifts in the anoxic boundary towards the surface of the lake and thus prevented OC mineralization in an oxic environment. Compared with many previous studies, this study is unbiased with respect to anthropogenic influences, latitude changes or significant vegetation changes, which might have an effect on OC burial efficiency. Our results demonstrate that climate variability was able to trigger mechanisms inherent to the lake, thereby resulting in oscillations in OC burial efficiency. The Sr/Ca ratio of bulk sediment reflects the formation of aragonite needles in this special lake setting and is a novel approach that can be utilized for paleoclimate reconstructions. We were able to identify several cooling events observed in the wider Mediterranean area, although the unique high resolution of our data enabled us to also identify a number of short-term cold and dry events throughout the 8.3 to 2.6 cal ka BP period that have mainly not been previously observed. Further high-resolution studies on additional archives would be beneficial for investigating their wider regional character.

Methods
Collection and extended description of sediment cores. Sediment cores were recovered from the deepest part of the Veliko jezero (M1-A at -45 m) and the deepest part of the second largest basin in Veliko jezero (M2 at -40 m) using a 3 m long piston corer (60 mm diameter UWITEC piston corer) deployed from a floating platform. The M1-A core had a total length of 459 cm, while the M2 core had a total length of 406 cm. Before being split lengthwise and photographed, the entire cores were subjected to magnetic susceptibility loop sensor measurements at 2 cm intervals with a Bartington MS2 magnetic susceptibility system. The split cores were logged for their lithology (smear slides), grain size, textures, structures, clast size and color using both the Munsell color chart and diffuse reflectance measurements (CIELAB-International Commission on Illumination L*a*b*) at 1 cm intervals using an X-rite DTP22 digital swatchbook spectrophotometer. Additionally, the magnetic susceptibility values of the split cores were determined using a Bartington MS2E system at 1 cm intervals. Working halves were subsampled at 1 cm intervals and stored at + 4 °C in plastic bags until further analysis, while the archived halves were stored in D-tubes in a cold chamber at the same temperature. For this study, we used intervals from 114 to 240 cm in core M1-A and from 127 to 266 cm in core M2.
High-resolution XRF scanning. Both cores M1-A and M2 were analyzed at XRF Core Scanner III (AVAATECH Serial No. 12) at the MARUM Center for Marine Environmental Sciences at the University of Bremen, Germany. Core M1-A was scanned at resolutions of 1 and 0.2 cm, while for Core M2, data were collected every 1 cm. In both cases, generator settings of 10 and 30 kV were applied. A current of 1 mA was used at 30 kV, and 0.2 mA and 0.5 mA at 10 kV were used at 1 cm and 0.2 cm intervals, respectively. The sampling times were 20 s and 30 s for the 1 cm and 0.2 cm intervals, respectively. Before scanning, the split core surface was covered with a 4 µm thin SPEXCerti Prep Ultralene foil to avoid contamination of the XRF measurement unit as well as desiccation of the sediment. The data reported here were acquired by a Canberra X-PIPS Silicon Drift Detector (SDD; Model SXD 15C-150-500) with a 150 eV X-ray resolution, a Canberra Digital Spectrum Analyzer DAS 1000 and an Oxford Instruments 100 W Neptune X-ray tube with rhodium (Rh) target material. Raw data spectra were processed by the analyzing the X-ray spectra by the Iterative Least square software (WIN AXIL) package from Canberra Eurisys. Reliable data (expressed in counts per second (cps)) were collected for the following elements: Br, Rb, Sr, Zr, Al, Si, S, K, Ca, Ti, Mn, and Fe; however, Mo values were considered unreliable because the concentration was close to or below the instruments' detection limit, and they were validated with an additional ICP-MS analysis on selected samples (Supplementary file, Figure S1). For each core, the RGB and LAB color parameters were obtained at 68 μm resolution. The analysis was performed within 24 h after the core was opened. Semiquantitative X-ray diffraction and scanning electron microscopy. To explore whether nonaragonite carbonate phases are present in the studied sediment samples, we performed XRD analyses on samples from Core M1-A. To cover the complete record from this core, we analyzed the mineral composition of 42 composite samples that were defined by homogenization of 3 cm long intervals. The analyses were performed on a PANalytical X'Pert Powder X-ray diffractometer equipped with Ni-filter CuKα radiation, a vertical goniometer with θ/θ geometry and a PIXcel detector. The scan conditions were as follows: 45 kV and 40 mA, ¼ divergence slit and anti-scatter slits, 0.02° 2θ step size, and 2 s time per step. After the samples were ground with mortar and pestle, they were sieved through a 0.4 mm sieve. To reduce the grain size of the material to < 5 μm, powdered samples were ground in McCrone micronizing mills. XRD digital scans were analyzed using the Philips X'Pert High Score search and match function to identify peaks and determine qualitative mineral compositions. Additionally, black and white laminas were examined throughout the cores with scanning electron microscopy (SEM) to study the morphology of aragonite minerals to test their inorganic origin. In the interval from 201 to 214 cm, SEM information was used to infer the morphology and potential nonbiogenic origin of Mg-calcite, which was proven by both XRD and SEM-EDS (energy-dispersive X-ray spectroscopy). and total nitrogen (TN) of the sediment samples were analyzed using a Thermo Scientific FLASH 2000 Series Nitrogen and Carbon analyzer of the Croatian Geological Survey at a 1 cm resolution in the studied interval. This method allows for direct measurements of the total carbon (TC) and TN. The addition of HCl removes the carbonate component and allows for the determination of TOC. Carbonate was removed by treating 1 g of sediment sample with 12 ml 4 M HCl and heating in a centrifuge tube sitting in a hot block for 2 h. The insoluble residue was washed with Milli-Q water and centrifuged (2x), freeze-dried and weighed. The carbon content of the insoluble residue after HCl treatment is the TOC. The difference between TC and TOC was used for the calculation of TIC, whereas the calcium carbonate (CaCO 3 ) content was calculated from the obtained TIC values. The C/N mass ratio was calculated by dividing the TOC and TN. XRD analyses of insoluble residue were used to confirm that carbonates were not present in the samples after HCl treatment. A split of the acid-washed sample was weighed into tin capsules optimized for stable carbon and nitrogen isotope ratio measurements in the Sample Weight Calculator of the Stable Isotope Facility (SIF) at the University of California, USA. The samples were shipped to SIF in sealed evaporation plates with 96 wells. The sediments at SIF were analyzed for δ 13 C and δ 15 N using an Elementar Vario EL Cube (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced with an Isoprime VisION IRMS (Elementar UK Ltd, Cheadle, UK). The isotope data are expressed relative to international standards VPDB (Vienna Pee Dee Belemnite). The long-term standard deviation reported by the SIF is 0.2‰ for δ 13 C and δ 0.3‰ for δ 15 N. Stable carbon isotopes were measured at a 10 cm resolution.
Statistical methods. Geochemical data obtained with XRF are compositional data, i.e., all components are parts of the same whole 67 . It is very hard to measure all elements; therefore, in reality, we analyzed a subcomposition, i.e., only some of all possible parts. Prior to the statistical analyses, one should represent data that are originally given as elements of a simplex space in log ratio coordinates 15 . Orthonormal log ratio coordinates (olr) were used 52,68 to construct geochemical proxies 53 . The orthonormal basis for the olr compositional biplots was constructed using the CoDa pack 69 . The construction of the sample basis was conducted by performing sequential binary partitioning (SBP) of a compositional vector 70 . Prior to the proxy construction, the dimensionality of the problem was reduced with the aid of a compositional biplot 71 , which helped in discarding redundant elements, i.e., those that carry geochemically similar information (have small variance between clr transformed variables). Proxies constructed via this method are free of compositional data restrictions regarding the multivariate statistics and correlation analyses. It is important to stress that the high Mg-calcite interval (group High D) was removed from the correlation analysis of constructed proxies reported in the main paper because of inherently different geochemical affiliations. The rationale for balance construction and compositional biplot interpretation is presented in the Supplementary information.

Data availability
The datasets generated during and/or analyzed during the current study are available in the PANGAEA repository, https:// doi. org/ 10. 1594/ PANGA EA. 924331