Climate influence on organic carbon burial efficiency as revealed from a multi-proxy lake sediment record in Croatia

The drivers of the efficiency in organic carbon (OC) burial are still poorly understood despite their key role in reliable projections of future climate trends. Here we shed new light on this question by presenting paleoclimate time series including OC content in sediments from Lake Veliko jezero, Croatia. We identified the Sr/Ca ratio in bulk sediments as (comparative) palaeotemperature and palaeohydrology indicator. Four major and six minor cold and dry events were detected in the 8.3 to 2.6 cal ka BP interval. The combined assessment of Sr/Ca ratios, OC content, C/N ratios, δ13C data, and modelled proxies for palaeoredox conditions and aeolian input reveals that cold and dry climate state promoted anoxic conditions in the lake enhancing preservation of organic matter and leading to increased OC burial efficiency. Our study contributes to that projected future increase of temperature might play an important role in OC burial efficiency of meromictic lakes.


Introduction
Lakes, with disproportionally large amounts of buried organic carbon (OC) per year if compared to the oceans 1 are of great importance for the global carbon budget and cycle and thus might have a vast impact on climate changes also in the future. However, the climate influence on OC burial efficiency is still not precisely understood. Temperature influence on the OC burial efficiency has been studied extensively, mainly for the lakes of higher northern latitudes [2][3][4][5][6] . Despite extensive efforts as represented by the large number of studies on this topic, there is still no overall supposition on a possible temperature influence on the OC burial efficiency. Nor are the main driving mechanisms recognized, which are likely and directly responsible for the observed changes in OC content of the lake sediments. Some studies infer that higher temperatures are enhancing OC burial mainly as a consequence of a denser vegetation cover 3,[7][8][9][10] , while at the same time an increase in temperature is positively correlated to OC mineralization 2,9,11 . In more recent studies the anthropogenic influence is preferred, primarily through the role of the reactive nitrogen and phosphorus on OC burial efficiency 4,6,12 , and a temperature effect 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 of OC burial efficiency were focussed on recent and sub-recent lake sediments, where climate effects can easily be concealed by anthropogenic influence. 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 have studied sediments of Lake Veliko jezero, Croatia ( Figure 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 palaeoclimate and palaeoredox proxies as well as an indicator for past aeolian activity. During the studied interval Lake Veliko jezero was a brackish, meromictic lake 15 , which is now submerged because of the Holocene sea level rise.
Our approach for correlating palaeoclimate and palaeoredox proxies was made taking into account compositional nature of the geochemical data 16 , that enabled more precise interpretations. The Sr/Ca ratios of the bulk sediment mainly derives from Sr and Ca concentrations of needle like aragonite in Core M1-A and was used as a paleoclimate, temperature proxy with higher ratio indicating cooler conditions. The Sr/Ca ratio of aragonite is typically used as a proxy for sea surface temperature (SST) variability 17 because incorporation of Sr into aragonite corral skeletons as well as in inorganic aragonite is temperature dependent, i.e. an increase in temperature lowers the Sr/Ca ratio in inorganic aragonite 18,19 and in coral skeletons 20 . However, this simplistic interpretation has been challenged because it has been proven that incorporation of Sr into aragonite corral skeletons is also affected by vital effects 21,22 , variability of the Sr/Ca ratio in the oceans 23,24 and unwanted effects caused by algal symbionts 25 . In contrast to numerous palaeoclimate studies of biogenic aragonite, palaeoclimate potential of inorganic aragonite in marine/lake sediments yet needs to be explored. This is mainly because sedimentary environments where inorganic aragonite precipitates are rare and studies show that Sr/Ca ratio of inorganic aragonite is much less sensitive to temperature changes when compared Sr/Ca ratio of coral aragonite skeletons 22 . Despite this obstacle, the main advantage of inorganic aragonite is lack of the vital effects problem. We validated the reliability of this novel proxy by also acquiring the same kind of data for Core M-2 that was taken from a different basin of the same lake (Age-Depth model for core M2 is presented in Supplementary file, Figure S4). An age-depth model of Core M1-A ( Figure 2) is previously described 26 . The significance of Sr/Ca ratio as a palaeotemperature proxy is also confirmed by correlation with other previously published studies including the main Holocene climatic events from the wider Mediterranean region 27-32 .

Results and discussion
Core M1-A recovered from Lake Veliko jezero is 417 cm long, four lithological units were distinguished (Figure 2). The first unit (417 to 343 cm) is a terra rossa type soil, the second unit (343 to 241 cm) is marsh to shallow lake sediment, the third unit (241 to 121 cm) consists of deep lake sediment characterized by alternations of white and dark laminas mainly composed of aragonite and organic matter. The last unit (121 to 0 cm) is made up of gray homogenous marine sediment with a few centimeters of oxidized lake sediment in the lowermost part. The third unit and the oxidized part of the unit four lake sediment were in the focus of this study. Core M2 resembles Core M1-A with the difference that unit boundaries are in different depths. Therefore, the studied interval spans from 127 to 266 cm in this core.

2.1.High-resolution XRF scanning and age-depth model
The Sr/Ca records were obtained by high-resolution XRF core scanning at split core surfaces in 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 robustness of our record and to exclude any potential analytical artefacts. The geochronological age of Core M2 is based on one plant remain sample, which was dated by radiocarbon methodology. Additional three datums 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 26 . The third date corresponds to a time of marine intrusion into the lake, which is marked by a sharp boundary between lithological units three and four ( Figure 2) in both cores. It was radiometrically dated in Core M1-A 26 . Individual radiocarbon dates were calibrated using the R package rbacon 33 . Despite lower resolution of the Core M2 record it is evident that major peaks and throughs in Sr/Ca are well presented in both datasets ( Figure 3).

2.2.Sr/Ca ratio as palaeoclimate 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 ( Figure 2). Needle-like aragonite precipitates nowadays during late spring and summer in adjacent lake Malo jezero 34 . Since Sr/Ca ratio of needle-like aragonite is largely temperature dependant 18,19 we assume that the Sr/Ca ratio of the studied sediment could be used as a relative palaeotemperature proxy.
Mineralogical analysis performed with a X-ray diffraction (XRD) and scanning electron microscope coupled with 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 ( Figure 2). The only exception is an interval from 201 to 214 cm where inorganic rhombohedral Mg-calcite is the main carbonate phase confirming previously published results 15 . 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 corresponding to pluvial periods observed in the wider Mediterranean region ( Figure   3 and 4, from ca. 7.6 to 7 ka BP) [29][30][31][35][36][37][38][39][40] . We propose that during this period increased freshwater input lowered the Mg/Ca ratio in the lake water leading to the Mg-calcite precipitation and hindering the precipitation of aragonite. This consequently led to Sr/Ca decrease in our records because Mg-Calcite incorporates Sr in the crystal lattice less effectively compared to the aragonite 41 thus the Sr/Ca ratio of bulk sediment cannot be interpreted as a relative palaeotemperature 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 29-31,35-40 . In order to be able to utilize the Sr/Ca ratio as a relative paleotemperature proxy we have proven that detrital Sr and Ca input is negligible. XRD analyses revealed aragonite as the main mineral, with minor quartz content limited only to the oldest portion of the studied period. 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 relative marine influence at this location because of a limited connection of the lake to the Adriatic Sea through permeable karst: the ocean water is characterized by higher Sr concentration compared to the freshwater 42 . Two possible scenarios emerge if hydrologically induced Sr availability was the limiting factor for Sr/Ca ratio of the bulk sediment.
Firstly, during cold periods, Sr/Ca in the lake water would be lower because of reduced evaporation rate i.e. decreased marine influence. This would finally result in a relative decrease in Sr concentration of the lake water and consequently Sr/Ca of the lake sediments. The opposite would be the case for warmer climate conditions. This scenario however can be discarded based upon good correlation of maxima in our Sr/Ca record with cold events which were previously observed in multiple paleoclimate archives in the Mediterranean region ( Figure 3). In the figure 3 locally estimated scatterplot smoothed (LOESS) Sr/Ca curve, with smoothing factor of 0.09, displays four distinct peaks centred at 8.3, 6.0, 4.25 and 2.9 cal ka BP. First Sr/Ca maxima in our record, centred at 8.3 ka is coeval to cold event described in pollen record from the Adriatic Sea by 29 and drop in the sea surface temperature 30 recognized in the same core. In Alboran Sea drop in the temperature at that time period was recorded as well 27,43 while minor drop in the temperature was also observed in the Gemini lake 32 . These events correlate with the North Atlantic cold event (NAC5) 44 . Following the pluvial period Second scenario implies that the environment during cold events was also dry, while during warmer events it was more wet. This would lead to a decreased relative marine influence and the dilution of the Sr budget in the lake water, which would finally result in lower Sr/Ca during warm periods and higher values during the cold periods. To further investigate the role of wet versus dry conditions potentially driving the Sr/Ca variability we studied additional chemical elements acquired by XRF scanning like Br, Rb, Si, K, Na, Al, Fe, Mn, Zr, Ti and Mo, most of them validated by wet-chemical analysis on discrete samples, e.g. through ICP-MS. We used all of these elements to model proxies called balances for palaeotemperature, aeolian input and palaeoredox conditions. This was done through orthonormal log ratio transformation 45 which enabled firm statistical parameter to be established for reliable interpretation of involved processes 46 . The methodology is fully explained in 46 and is given briefly in method section. In the end five balances were modelled, each is a proxy for certain process; b3 and b4 are proxies for palaeoredox conditions, b5 is a proxy for aeolian activity while b2 is a palaeoclimate proxy equivalent to Sr/Ca (rationale for balance construction is shown in Supplementary file, text S1, sign matrix table used for OLR transformation is shown in Supplementary file, table S1).
The data analyses show that an increase in Sr/Ca correlates to 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) (Figure 4). Comparison between b2 and Sr/Ca is presented in Supplementary file, figure S3). Furthermore, a speleothem based palaeohydrological reconstruction from Corchia cave in Italy 31 implies that during the time of Sr/Ca maxima in our record climate conditions were generally drier ( Figure 4). However, it is also evident that 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 28,38,[47][48][49][50] . This implies that hydrological conditions probably played, compared to the temperature effect, less of a role in driving Sr/Ca ratios. Following all lines of evidences, we propose that the Sr/Ca ratio of Lake Veliko jezero bulk sediment is representing Sr/Ca ratio of inorganic needle-like aragonite that mainly reflects relative palaeotemperature changes while hydrological variability likely playing a secondary role.
One of the main advantages of XRF core scanning despite being non-destructive is the high resolution. If Sr/Ca ratios is a temperature indicator, then it is maybe also 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 3.3 ka BP event was recognized in Lake Verdarolo 32 . 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 ( Figure 4). These events can be detected due to the combined effects of the limited size and detrital influence on the studied lake(s) and high-resolution data. The small size of this lake (surface area of 1.44 km 2 ) with very limited detrital influence i.e. small effects of internal and landscape filters 51 results into increased sensitivity in recording smaller-scale climate events.

2.3.Disentangling the drivers of OC burial efficiency
To decipher potential drivers of OC burial efficiency we analyzed OC, inorganic carbon (INC) and organic nitrogene (N) and explored its 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 the land-derived and autochthonous organic matter 52,53 , while more negative δ 13 C values could also be related to the increase of the land-derived component of the organic matter 54 . Our results demonstrate a slight long-term decreasing trend in C/N record coinciding 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 in line with the overall detrital influence during our studied interval ( Figure 4). Additional evidence comes from the occurrence of quartz only found at the base of the studied interval. A combination of two factors is most probably causing the described trends. First, lake deepening caused by the Holocene sea level rise, moved 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, Figure 1) was increased. Second, due to the Holocene sea level rise a gradual increase in marine influence through the permeable karst did occur. The sea level rise changed the lake biota 15 and consequently organic carbon content and composition.
During the studied interval just two different pollen zones occurred, one with Juniperus and Phillyrea at ca. 8 to 6.5 ka BP and another of Quercus ilex from 6.5 ka BP to present 55 .
This finding implies minor changes in the terrestrial vegetation with negligible impact on carbon content, composition and variability.
The correlation of the OC content with Sr/Ca record reflecting temperature variations on millennial time scale suggests that the OC content increased during cold interludes ( Figure   4). The observed OC increase during cold events is in line with previous study 2 where a temperature decrease leads to low mineralization of OC. However, the temperature effect through OC mineralization was not the substantial factor for OC preservation in the Lake Veliko jezero sediments. If temperature was driving the OC preservation on millennial time scales, cold climate conditions would decrease the degradation rate of algal (autochthonous) organic matter 56 . This would result in all lower C/N values, higher OC and more positive δ 13 C values, which we did not see in our data. Indeed, we observe slightly higher C/N and more negative δ 13 C values during cold spells. This would imply that land-sourced organic matter increased during cold spells as a result of enhanced aeolian input, confirmed by correlation analyses r(b2-b5) (Figure 4). Yet, we argue that an increase in land-sourced organic matter is not the main mechanism in overall OC increase. If this would have been the case, then one would expect that more than double increase in OC during the cold spells would cause 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 correlates with paleoredox proxy (Figure 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. Decrease in temperature and possibly drier conditions during cold events would shift the redox zone boundary and thermocline closer to the surface 57 . This would prevent mixing of the water throughout the water body, thus causing the anoxic boundary to move upward, leaving the majority of the water column under anoxic conditions. This is also confirmed by high correlation between palaeotemperature (b2) and palaeoredox (b3) proxies.
Although the depth of a thermocline depends on a number of factors such as lake size, dissolved organic content, temperature, wind activity etc. 58 we believe that the temperature was the key factor controlling thermocline and redox zone boundary depth in Lake Veliko jezero. These findings are underpinned by a study of modern processes in Lake Veliko jezero lake, where the thermocline occurs only during summer months 59 . Finally, a shallow thermocline/redox zone boundary would cause that most of the organic matter produced in and transported into the lake is prevented from decomposing resulting in higher OC values during cold spells.

Conclusions
Data presented demonstrate that temperature changes may have a significant impact on OC burial efficiency. Temperature decrease and likely drier climate conditions caused shifting of anoxic boundary towards the surface of the lake and thus prevented OC mineralization in an oxic environment. This study is, unlike many previous studies, unbiased in respect to the anthropogenic influence, changes in latitude 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 resulting into oscillations of 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 known in the wider Mediterranean area, but the unique high resolution of our data enabled to also identify a number of short-term cold and dry events throughout the 8.3 to 2.6 cal ka BP period that mainly have not been found before. Further high-resolution studies on additional archives would be beneficial for investigating their wider regional character.

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 from 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 colour using both Munsell colour 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. This was also described in 26,46 and is reported here as well. For this study we used intervals from 114 to 240 cm in the core M1-A and from 127 to 266 cm in the core M2.

High-resolution XRF scanning
Both cores M1-A and M2 were analyzed at XRF Core Scanner III (AVAATECH Serial  Figure S1). This methodology is also described in 26 where only trend of the K and Zr were shown while no data was published. For each core, RGB and LAB colour parameters were obtained in 68 μm resolution. The analysis was done within 24 hours after the core is opened.

Semi-quantitative X-ray diffraction and Scanning electron microscope
To explore if non-aragonite 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, a step size of 0.02° 2θ, and a time per step of 2 s. After 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 a Philips X'Pert High Score search and match function to identify peaks and determine qualitative mineral compositions 46 . Additionally, black and white laminas were examined throughout the cores with scanning electron microscopy (SEM) to study the morphology of aragonite minerals in order to test its inorganic origin. In the interval from 201 to 214 cm SEM information was used to infer about morphology and potential non-biogenic origin of Mg-calcite, which was proven by both XRD and SEM -EDS (Energy-dispersive X-ray spectroscopy).   60 . It is very hard to measure all elements; therefore, in reality, we are analysing a subcomposition, i.e., only some of all possible parts. Prior to statistical analyses one should represent data that are originally given as elements of a simplex space in log ratio coordinates 16 . Orthonormal log ratio coordinates (olr) were used 45,61 to construct geochemical proxies 46 . The orthonormal basis for the olr compositional biplots was constructed using CoDa pack 62 . The construction of the sample basis was conducted by performing a Sequential Binary Partition (SBP) of a compositional vector 63 . Prior to proxies construction dimensionality of the problem was reduced with an aid of a compositional biplot 64 which helped to discard redundant elements i.e. those which carry geochemically similar information (have small variation between clr transformed variables). Proxies constructed in that way are free of compositional data restrictions regarding the multivariate statistics and correlation analyses. It is important to stress out that high Mg-calcite interval (group High D) because of inherently different geochemical affiliations was removed from correlation analysis of constructed proxies reported in the main paper. Rationale for balance construction and compositional biplot interpretation is presented in Supplementary information.

Data availability statement
The datasets generated during and/or analysed during the current study are available in the PANGAEA repository, https://doi.org/10.1594/PANGAEA.924331