Reconstructing Holocene temperature and salinity variations in the western Baltic Sea region: a multi-proxy comparison from the Little Belt (IODP Expedition 347, Site M0059)

Abstract

Sediment records recovered from the Baltic Sea during Integrated Ocean Drilling Program Expedition 347 provide a unique opportunity to study paleoenvironmental and climate change in central and northern Europe. Such studies contribute to a better understanding of how environmental parameters change in continental shelf seas and enclosed basins. Here we present a multi-proxy-based reconstruction of paleotemperature (both marine and terrestrial), paleosalinity, and paleoecosystem changes from the Little Belt (Site M0059) over the past $\sim$ 8000 years and evaluate the applicability of inorganic- and organic-based proxies in this particular setting.

All salinity proxies (diatoms, aquatic palynomorphs, ostracods, diol index) show that lacustrine conditions occurred in the Little Belt until$\sim$ 7400 cal yr BP. A connection to the Kattegat at this time can thus be excluded, but a direct connection to the Baltic Proper may have existed. The transition to the brackish–marine conditions of the Littorina Sea stage (more saline and warmer) occurred within $\sim$ 200 years when the connection to the Kattegat became established after $\sim$ 7400 cal yr BP. The different salinity proxies used here generally show similar trends in relative changes in salinity, but often do not allow quantitative estimates of salinity.

The reconstruction of water temperatures is associated with particularly large uncertainties and variations in absolute values by up to 8 ${}^{\circ }$C for bottom waters and up to 16 ${}^{\circ }$C for surface waters. Concerning the reconstruction of temperature using foraminiferal Mg $/$ Ca ratios, contamination by authigenic coatings in the deeper intervals may have led to an overestimation of temperatures. Differences in results based on the lipid paleothermometers (long chain diol index and TEX${}_{\mathrm{86}}^{\mathrm{L}}$) can partly be explained by the application of modern-day proxy calibrations to intervals that experienced significant changes in depositional settings: in the case of our study, the change from freshwater to marine conditions. Our study shows that particular caution has to be taken when applying and interpreting proxies in coastal environments and marginal seas, where water mass conditions can experience more rapid and larger changes than in open ocean settings. Approaches using a multitude of independent proxies may thus allow a more robust paleoenvironmental assessment.

Introduction

Recent climate change is showing a persistent trend to a more positive state of the North Atlantic Oscillation, resulting in wetter conditions, especially over northern Europe (Hurrell, 1995; Visbeck et al., 2001). Further increased westerly winds together with a changed barotropic pressure gradient may result in increased inflows of open marine waters and therewith more saline deep waters in the Baltic Sea in the future (see Meier et al., 2006; Meier, 2015). Concomitantly, the increase in continental run-off due to increased precipitation may result in a freshening of the surface waters in the Baltic Sea (Matthäus and Schinke, 1999; Gustafsson and Westman, 2002).

Both temperature and salinity changes will thus have important consequences for the Baltic Sea environment (e.g. Meier et al., 2012). To improve our understanding of the impact and magnitude of future environmental changes in the Baltic Sea region, it is essential to generate high-resolution paleoreconstructions in order to investigate how salinity and temperature varied in the past (e.g. Zillén et al., 2008; Andrén et al., 2015a; Ning et al., 2017). A multi-proxy approach, comprising proxies representative for bottom water, surface water, and air (terrestrial) conditions, is used here to reconstruct a wide array of environmental change (in particular temperature and salinity) from the same or closely neighbouring samples and thus how conditions simultaneously changed within the same age constraints.

The occurrence of specific species of ostracods and foraminifera is an important indicator of bottom water parameters in the case of the Baltic Sea primarily salinity (Lutze, 1965, 1974; Murray, 2006; Frenzel and Boomer, 2005; Frenzel et al., 2010; Viehberg et al., 2008). Additionally, stable oxygen and possibly carbon isotopes in foraminifera are commonly related to changes in salinity, though stable oxygen isotopes are also influenced by temperature and global ice volume, while productivity is another control factor for carbon isotopes (Kristensen and Knudsen, 2006; Filipsson et al., 2016). Bottom water temperatures can be reconstructed using Mg $/$ Ca in benthic foraminifera (Raitzsch et al., 2008), but large gradients in salinity (and low absolute salinity) as well as parameters related to the carbonate system (e.g. alkalinity and carbonate saturation state) in the Baltic Sea may significantly influence the application of Mg $/$ Ca ratios (Groeneveld and Filipsson, 2013).

Due to the low salinity and overall shallowness of the Baltic Sea, planktonic foraminifera are absent from the Baltic Sea region (Seidenkrantz, 1993; Andrén et al., 2015a). Diatoms and marine palynomorphs such as dinoflagellates, however, are common and often very good indicators of changes in salinity (Andrén et al., 2000; Snoeijs and Weckström, 2010; Ning et al., 2017).

A proxy, which appears to be only dependent on temperature and may therefore be ideal for the Baltic Sea, is clumped isotopes on molluscs (Henkes et al., 2013; Wacker et al., 2014). This proxy has not yet been applied to samples from the Baltic Sea. “Clumped isotopes” (or, more correctly, “clumped isotopologues”) are isotopic molecules that contain more than one rare isotope (e.g. the carbonate isotopologue Ca${}^{\mathrm{13}}$C${}^{\mathrm{18}}$O${}^{\mathrm{16}}$O${}_{\mathrm{2}}$). The abundance of these molecules in carbonate is determined by the formation temperature of the mineral.

The lipid paleothermometer TEX${}_{\mathrm{86}}$ (tetraether index of tetraethers consisting of 86 carbon atoms), based on the sedimentary distribution of isoprenoid glycerol dialkyl glycerol tetraethers (GDGTs) produced by pelagic Thaumarchaeota, has been successfully used in the reconstruction of seawater temperatures in numerous paleoceanographic settings (Schouten et al., 2002, 2013b). A derivative of this proxy is the TEX${}_{\mathrm{86}}^{\mathrm{L}}$, which in low temperature environments (generally below 15 ${}^{\circ }$C) shows a high correlation with sea surface temperature (SST; Kim et al., 2010). In the Baltic Sea, the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ has been shown to be best correlated with late summer SSTs and used to investigate climate-induced variations in water temperature over the past $\sim$ 1000 years (Kabel et al., 2012), while complete TEX${}_{\mathrm{86}}^{\mathrm{L}}$ records of Holocene climate change are still missing. A more recently introduced proxy that shows promise in paleoenvironmental studies is the long chain diol index (LDI; Rampen et al., 2012). This proxy makes use of the sedimentary distribution of long chain diols synthesized by eustigmatophytes and has been suggested to either reflect an annual (Rampen et al., 2012; de Bar et al., 2016) or a summer SST signal (Lopes dos Santos et al., 2013) in marine settings. In addition to SST, however, salinity has been observed to also affect the distribution of long chain diols in aquatic environments and consequently the diol index (DI) has been proposed as a tool to qualitatively assess relative changes in freshwater input to marine environments (Versteegh et al., 1997). An alternative means of investigating past climate dynamics is provided by terrestrial palynomorphs such as pollen grains. The modern analogues technique (MAT) and similar approaches further allow the calculation of quantitative climate data (e.g. temperature, precipitation, seasonality) from pollen assemblages (e.g. Guiot, 1990; Kotthoff et al., 2008, 2011).

Despite a large amount of proxy-based research performed in the Baltic Proper (e.g. Zillén et al., 2008, and references therein; Kotilainen et al., 2014), there are no continuous records from the Baltic Sea and its connections to the Kattegat (Little Belt, Great Belt, and Öresund) enabling multi-proxy studies in high temporal resolution for the entire Holocene (Andrén et al., 2015a). In particular, the Holocene history of the Little Belt has not been the focus of research (Bennike and Jensen, 2011). Continuous high-resolution records would not only allow correlation between the different basins in the Baltic Sea region but also provide a connection to terrestrial archives.

In the framework of International Ocean Drilling Program (IODP) Expedition 347, a series of long, continuous sediment records that allow the reconstruction of environmental change in central and northern Europe in an unprecedented manner were recovered from the Baltic Sea region (Andrén et al., 2015a, b). Here we present a centennial-resolution multi-proxy study of IODP Expedition 347 Site M0059 from the Little Belt, one of the connections between the Baltic Sea and the North Sea, using a variety of different proxies that are commonly used in paleoceanography to reconstruct past seawater temperature and salinity. We investigate how these proxies perform in a proximal setting like the Little Belt and how they link to seasonality. We selected the last $\sim$ 8000 years at Site M0059 (top $\sim$ 53 m of sediment), which includes the transition from freshwater to brackish–marine conditions. We show that a multi-proxy approach allows the unravelling of the different factors that have influenced past conditions in the Little Belt, suggesting that its history may have been partly different from that of the Baltic Proper.

Regional setting and methods

Regional setting

The Baltic Sea is connected to the North Sea by the Belt Seas, including the Öresund, the Great Belt, and the Little Belt, through the Kattegat and Skagerrak. It is one of the world's largest brackish water bodies, with a surface water salinity ranging from close to zero in the innermost parts of the Gulf of Finland, 7 in the Baltic Proper, and between 8 and 24 in the Belt Sea. A halocline, varying between 20 and 80 m water depth in the different subbasins, separates the upper water masses from more saline bottom waters. The mean water depth is 54 m, varying between 25 and 200 m with the deepest basin, the Landsort Deep (central Baltic Sea) having a maximum depth of 459 m (e.g. Leppäranta and Myrberg, 2009). The Baltic Sea basin is surrounded by two biomes, from temperate forest with mixed coniferous and broadleaved trees in the south to boreal forest with taiga-like conditions in the north.

Currents in the Little Belt consist of either relatively low-salinity (14–17) surface water flowing out of the Baltic Sea into the Kattegat or higher-salinity ($>$ 24) currents flowing into the Baltic Sea. Depending on the volume and mixing rate, the denser inflowing water sinks after passing the main sills and can therefore be encountered at different water depths in the Baltic Proper. Because the Baltic Sea has a surplus of freshwater input, i.e. the input of water by precipitation and rivers exceeds evaporation, a positive barotropic pressure gradient exists between the Baltic Proper and the Kattegat, so that under “normal” conditions surface waters flow from the Baltic Proper into the Kattegat (HELCOM, 1986; Jakobsen and Ottavi, 1997). The contribution of water mass transport through the Little Belt is minor ($<$ 10%) in comparison to the Great Belt and the Öresund (HELCOM, 1986). Occasionally the surface flow reverses and more saline water enters over the sills into the straits (Matthäus and Schinke, 1999; Mohrholz et al., 2015). The higher salinity causes the inflowing water to sink and continue as subsurface water into the Baltic Proper. Because the Little Belt is very narrow, with high current velocity and turbulence, a major marine inflow from the Kattegat can lead to mixing of the water column (Jakobsen and Ottavi, 1997).

Mean surface water salinity in the Little Belt was $\mathrm{15.9}\pm \mathrm{0.4}$ for a 10-year period (2004–2014), varying between $\mathrm{14.2}\pm \mathrm{0.9}$ in spring and $\mathrm{17.6}\pm \mathrm{1.5}$ in autumn, while mean bottom water salinity was $\mathrm{23.6}\pm \mathrm{1.0}$, varying between $\mathrm{21.5}\pm \mathrm{0.9}$ in winter and $\mathrm{25.2}\pm \mathrm{1.1}$ in summer (ICES, 2017). During the same time period, the mean annual SST was $\mathrm{10.6}\pm \mathrm{0.6}$ ${}^{\circ }$C (ICES, 2017). The mean winter SST was $\mathrm{2.5}\pm \mathrm{1.5}$ ${}^{\circ }$C, while the mean summer SST was $\mathrm{17.6}\pm \mathrm{0.7}$ ${}^{\circ }$C with maximum water temperatures of 21–22 ${}^{\circ }$C observed from late July to mid-August (ICES, 2017). Bottom water temperatures averaged $\mathrm{7.6}\pm \mathrm{0.6}$ ${}^{\circ }$C and varied between $\mathrm{4.3}\pm \mathrm{0.7}$ ${}^{\circ }$C in winter and $\mathrm{9.3}\pm \mathrm{1.3}$ ${}^{\circ }$C in summer (ICES, 2017).

Map of the southern Baltic Sea (redrawn from Groeneveld and Filipsson, 2013). The positions of IODP Site M0059 (red star; 55${}^{\circ }$0.29${}^{\prime }$ N, 10${}^{\circ }$6.49${}^{\prime }$ E; Andrén et al., 2015a), Lake Belau (Dörfler et al., 2012), and Lake Flarken (Seppä et al., 2005) are indicated. Continuous (surface) and dashed (bottom) lines indicate dominant current directions.

[Figure omitted. See PDF]

The Baltic Sea region has experienced several climate-driven hydrological changes during the Holocene. Its history is highly dynamic and governed by the regional isostatic rebound and global sea level changes, resulting in alternating freshwater and brackish phases and complex shoreline development (e.g. Björck, 2008; Knudsen et al., 2011; Andrén et al., 2011). The time interval investigated in this study comprises the transition from freshwater to the brackish–marine conditions of the Littorina Sea stage to the present brackish Baltic Sea. In the Baltic Proper, the Baltic Ice Lake ($\sim$ 16 000 to $\sim$ 11 700 cal yr BP) was followed by the partly brackish Yoldia Sea stage ($\sim$ 11 700 to $\sim$ 10 700 cal yr BP), the lacustrine Ancylus Lake stage ($\sim$ 10 700 to 10 200 cal yr BP), a transitional low-salinity phase ($\sim$ 10 200 to $\sim$ 7400 cal yr BP), and the Littorina stage ($\sim$ 7400 cal yr BP; e.g. Sohlenius et al., 1996; Andrén et al., 2000, 2011; Björck, 2008).

Material and methods

During IODP Expedition 347, cores were recovered from five holes (A to E) at Site M0059 (Little Belt; Fig. 1), at $\sim$ 37 m water depth (55${}^{\circ }$0.29${}^{\prime }$ N, 10${}^{\circ }$6.49${}^{\prime }$ E). The results of our multi-proxy analyses are based on sediments from Holes A, C, D, and E. For most proxies, $\sim$ 35 samples have been analysed.

Accelerator mass spectrometry radiocarbon dates from Expedition 347 Site M0059; for more detailed information compare Table S1 and Van Helmond et al. (2017).

Laboratory ID	Mean composite	Material dated	${}^{\mathrm{14}}$C Age	Error	Calibrated age
(LuSNo)	depth (m)		(yr)	(yr)	(yr; median)
11289	2.84	Abra alba	700	35	405
11476	4.97	Abra alba	1265	35	899
11291	8.79	Abra alba	1585	35	1213
11292	14.24	Macoma balthica	2740	40	2556
11477	21.14	Macoma balthica	3780	35	3829
11295	23.49	Macoma balthica	4035	40	4187
11479	26.60	Macoma balthica	4395	40	4664
11297	27.03	Macoma balthica	4500	40	4785
11480	28.20	Macoma balthica	4540	40	4845
11300	29.51	Macoma balthica	4780	40	5147
11481	31.80	Macoma balthica	4890	40	5313
11482	34.05	Abra alba	5155	40	5603
11301	35.39	Abra alba	5415	40	5885
11483	38.04	Abra alba	5770	40	6289
11303	39.75	Macoma balthica	6130	40	6662
11304	48.64	Bivalve fragments	6845	45	7437

Sedimentology and age model

The sedimentology of Holes M0059A, D, and E was assessed in the framework of the expedition's onshore scientific sampling party (Andrén et al., 2015a, b). The sedimentary sequence of IODP Site M0059 has been divided into seven lithostratigraphic units based on visual core descriptions and smear slide analyses (Andrén et al., 2015a), from which subunits Ia and Ib have been investigated in this study. Subunit Ia (0–49.37 m composite depth – m c.d.) is composed of mostly homogeneous black to greenish black clay, with some millimetre-scale laminations. It was deposited under brackish–marine conditions, as indicated by diatom assemblages (Andrén et al., 2015a; Fig. 2). Subunit Ib (49.37–53.57 m c.d.) is a continuation of Subunit Ia. Its base consists of 10–15 cm of black, laminated clays. The remainder of Subunit Ib comprises greenish to gray silty clay, intercalated by centimetre-scale pale green laminae. The presence of freshwater diatoms indicates lacustrine conditions (Andrén et al., 2015a, Fig. 2). Magnetic susceptibility data (supported by natural gamma ray and gamma ray attenuation density data) collected from Holes M0059A–E were used to correlate between each hole and to construct a composite splice section for IODP Site M0059 including a metres composite depth (m c.d.) scale (Andrén et al., 2015a), which was used in this study.

The age model for IODP Site M0059 is based on 16 radiocarbon dates derived from bivalve fragments and intact bivalve specimens (Abra alba, Macoma balthica; Table 1, and Table S1 in the Supplement). Age–depth modelling was performed with clam (version 2.2; Blaauw, 2010) with 2000 iterations using the Marine13 calibration data set (Reimer et al., 2013) and with a deviation ($\mathrm{\Delta }R$) of $-$90 $\pm$ 53 from the Marine13 reservoir age. A detailed description of the age reconstruction is given by Van Helmond et al. (2017). Below 48.64 m c.d., the ages have been extrapolated linearly, assuming a constant sedimentation rate, because no material suited for dating purposes was recovered from the deeper sections of the sediment sequence. In order to confirm the precision of the ${}^{\mathrm{14}}$C-based age model, we have biostratigraphically compared the marine pollen record from IODP Site M0059 (Fig. 2) with a pollen record from varved sediments in Lake Belau in northern Germany (Fig. 1; Dörfler et al., 2012).

Comparison of the pollen record from IODP Site M0059 with the varved pollen record of Lake Belau (northern Germany; Dörfler et al., 2012). Sedimentological units, lithology (Andrén et al., 2015a), and position of ${}^{\mathrm{14}}$C ages (Van Helmond et al., 2017) plotted vs. depth and plotted vs. age (extrapolated below $\sim$ 48.64 m c.d.). (a) Environmental zones (EZs; IODP Site M0059) and Firbas pollen zones (Lake Belau). (b) Selected pollen taxa (with horizontal bars: IODP Site M0059; filled: Lake Belau). (c) Total organic carbon (TOC; white: shipboard data; triangles: Hole A; rhombi: Hole C; pentagons: Hole D; circles: Hole E; grey squares: onshore data) for IODP Site M0059. (d) Climate phases and chronozones.

[Figure omitted. See PDF]

Geochemistry

In order to measure the total organic carbon (TOC) content, sediment samples from Holes C and E were freeze-dried, powdered, and homogenized using an agate mortar and pestle. About 0.3 g of powdered sediment per sample was decalcified using 1 M HCl. Subsequently, the samples were washed repeatedly with Milli-Q water, after which they were dried for 72 h at 60 ${}^{\circ }$C. Finally, the samples were powdered and homogenized again, after which they were measured using a Fisons Instruments NA 1500 NCS analyser at Utrecht University, the Netherlands. Results were normalized to international standards and TOC was calculated upon correction for weight loss by decalcification. The average analytical uncertainty of 0.07 wt % was calculated based on duplicate analysis of sediment samples (Van Santvoort et al., 2002).

Palynology and pollen-based quantitative climate reconstructions

For the analysis of terrestrial and marine palynomorph assemblages, 36 samples of $\sim$ 1 cm thickness were taken from 0 to $\sim$ 53 m c.d. Per sample, 1 to 6 g of sediment was processed using standard palynological techniques (HCl and HF treatment and sieving with 7 $\mathrm{\mu }\mathrm{m}$ mesh) for marine sediments (e.g. Kotthoff et al., 2008). Lycopodium marker spores were added to the samples in order to calculate palynomorph concentrations (Stockmarr, 1971; Figs. 2, 3 in this paper). When possible, $>$ 200 terrestrial palynomorphs (excluding spores and bisaccate pollen grains) were determined and counted per sample under 400 to 1000 $\times$ magnification. In addition, organic-walled dinoflagellate cysts (dinocysts, as indicators for marine conditions; e.g. de Vernal and Marret, 2007) and freshwater algae (as indicators for freshwater influence; e.g. Mudie et al., 2002) were determined and counted when occurring. We did not aim at counting a certain sum of dinocysts per sample since their absolute amount varied greatly over the analysed interval. On average, 50 dinocysts were identified per sample. This relatively low value is explained by the general rarity of dinocysts in some samples and relatively low dinocyst diversity in the Baltic Sea (e.g. Ning et al., 2016).

Quantitative climate data were calculated from pollen data using the MAT (e.g. Guiot, 1990). For this technique numerical methods are used to identify recent analogues for fossil pollen assemblages (e.g. Guiot, 1990; Kotthoff et al., 2008, 2011). Climate conditions for the best modern analogues are considered to be most similar to the conditions during which the fossil pollen assemblage was deposited. The MAT reconstructions are based on a database with $>$ 3500 modern pollen spectra from Europe (Bordon et al., 2009) and the Mediterranean area (Dormoy et al., 2009). Here we used the 10 modern assemblages with the smallest chord distances for the reconstructions of the climate parameters. Bisaccate pollen was removed from the evaluated samples and the database due to its over-representation in marine pollen records (e.g. Kotthoff et al., 2008).

Sedimentological units, lithology (Andrén et al., 2015a), and position of ${}^{\mathrm{14}}$C ages (Van Helmond et al., 2017) and IODP Site M0059 precipitation and salinity proxies plotted vs. depth and age (extrapolated below $\sim$ 48.64 m c.d.). (a) Environmental zones (EZs) based on diatom data. (b) Pollen-based precipitation reconstructions based on the modern analogues technique. (c) Diatom salinity-based affinities and life forms (benthic vs. pelagic) following Snoeijs (1993), Snoeijs and Vilbaste (1994), Snoeijs and Potapova (1995), Snoeijs and Kasperoviciene (1996), and Snoeijs and Balashova (1998) with species assignments provided in Table S2. (d) Palynomorph $/$ pollen ratios based on total amount of terrestrial pollen, freshwater algae comprising Botryococcus and Pediastrum species, dinocysts predominately comprising cysts of Protoceratium reticulatum (Operculodinium cysts), Lingulodinium, and Spiniferites species (Gymnodinium excluded). (e) Ostracod shallow water and deeper-living marine taxa plotted from right to left. (f) Total abundance of benthic foraminifera per cubic centimetre of dried sediment with the relative abundance of the most common species; benthic foraminiferal taxa (Elphidium incertum, E. selseyense, E. magellanicum, and Ammonia beccarii) shown as a percentage of total benthic foraminiferal fauna; total abundance of benthic foraminifera shown as specimens per cubic centimetre of dried sediment; stable carbon and oxygen isotope data analysed on E. selseyense. (g) Diol index indicating relative changes in surface water salinity.

[Figure omitted. See PDF]

Diatoms

Sediment samples ($n=\mathrm{36}$; thickness: 1 cm) for diatom analysis were freeze-dried. Subsequently, a known mass of sediment, between 0.03 and 0.07 g dry weight, was subsampled from each sample and treated according to Warnock and Scherer (2015). Briefly, sediments were cleaned with dilute HCl to remove carbonates and 10 % H${}_{\mathrm{2}}$O${}_{\mathrm{2}}$ to remove organic carbon. The resultant sediment slurry was dispersed in a 2 L square cross section beaker containing distilled water and a small glass table with known area including a coverslip, allowing for ultimate abundance determination. Coverslips were dried and permanently mounted to glass slides with the mounting media Naphrax (refractive index $=$ 1.65).

Each slide was counted to a minimum of 300 vegetative diatom valves. Absolute diatom abundance (ADA) can be calculated by counting the number of fields of view needed to reach 300 diatoms and scaling the number of diatoms counted in that known area relative to the total amount of surface area the sample settled onto and the original dry weight of the sediment sample (Warnock and Scherer, 2015). Diatoms were identified to species level according to Snoeijs (1993), Snoeijs and Vilbaste (1994), Snoeijs and Potapova (1995), Snoeijs and Kasperoviciene (1996), Snoeijs and Balashova (1998), Krammer and Lange-Bertalot (1986, 1988, 1991a, b), and Patrick and Reimer (1966). The salinity preferences and life forms (benthic or pelagic) of diatoms were classified according to Snoeijs (1993), Snoeijs and Vilbaste (1994), Snoeijs and Potapova (1995), Snoeijs and Kasperoviciene (1996), and Snoeijs and Balashova (1998) (compare Table S2). All slides were analysed using differential interference contrast at 1000 $\times$ magnification and oil immersion. Chaetoceros resting spores (CRSs) were counted but not identified to the species level and were not included in the ADA calculations for this study. ANOVA was used to evaluate the differences between the means of these data, accompanied by a Tukey–Kramer test to evaluate pairwise relationships (for all statistical tests $\mathit{\alpha }=\mathrm{0.05}$). Because data derived from sediment cores are inherently chronological, the results of pairwise tests are only reported for adjacent environment zones (EZs). All statistical relationships were evaluated using PAST v. 3.10 (Hammer et al., 2001).

Benthic foraminifera

For foraminiferal analyses, 36 samples of 2 cm sediment thickness (20 cm${}^{\mathrm{3}}$ sediment) were prepared and wet-sieved through sieves with 63, 100, and 1000 $\mathrm{\mu }\mathrm{m}$ mesh sizes. Each fraction was subsequently dried at 40 ${}^{\circ }$C and weighed. The 100–1000 $\mathrm{\mu }\mathrm{m}$ fraction was used for the foraminiferal assemblage analysis. In samples with a high concentration of mineral grains, the 100–1000 $\mathrm{\mu }\mathrm{m}$ fraction of the samples was subjected to heavy liquid treatment using the heavy liquid tetrachloroethylene (C${}_{\mathrm{2}}$Cl${}_{\mathrm{4}}$) with a specific gravity of 1.6 g cm${}^{-\mathrm{3}}$ to separate foraminiferal tests from mineral grains prior to analysis. The residual fraction after heavy liquid separation was checked for foraminifera; only very few specimens were found and thus only the 100–1000 $\mathrm{\mu }\mathrm{m}$ light fraction was used for foraminiferal analyses. To ensure statistical validity of the assemblages, a minimum of 300 individual foraminifera were counted from each sample, when possible. Benthic foraminifera were counted and identified to species level to characterize the faunal assemblages. The benthic foraminiferal concentrations were calculated as number of specimens per cubic centimetre of sediment.

Ostracods

Ostracods were divided into predominant ecological groups and their ecological data were used to reconstruct paleosalinity. A total of 75 samples (each comprising 30 cm${}^{\mathrm{3}}$) were processed for ostracod analysis. Since ostracod material was not abundant enough to produce reliable percentage data, we used the total abundance of ostracods per 30 cm${}^{\mathrm{3}}$. Samples were freeze-dried and washed over a 63 $\mathrm{\mu }\mathrm{m}$ sieve and subsequently oven dried in paper filters at 40–50 ${}^{\circ }$C. Ostracods were picked from the entire sample residues and their valves identified and counted.

Two ecological classes of ostracods in relation to salinity were distinguished based on ecological data and distribution of modern taxa in the Baltic Sea and adjacent areas (Frenzel et al., 2010): freshwater and brackish marine. Among brackish–marine ostracods, very shallow and deeper water species were distinguished based on ecological data (Frenzel et al., 2010). Predominance of a very shallow water group implies salinity above 6–10 and an environment such as a lagoon, estuary, or very shallow water open sea. The deeper-living ecological group includes species that can be found in the open sea environments at both very shallow and deeper locations with salinity $>$ 7–14.

Stable isotopes in foraminifera

Stable oxygen and carbon isotope (${\mathit{\delta }}^{\mathrm{18}}$O, ${\mathit{\delta }}^{\mathrm{13}}$C) measurements for 36 samples were performed on 10–20 specimens (100–1000 $\mathrm{\mu }\mathrm{m}$ fraction) of Elphidium selseyense and Elphidium incertum. The analyses were performed on a Finnigan MAT 251 gas isotope ratio mass spectrometer equipped with a Kiel I automated carbonate preparation device at MARUM, University of Bremen, Germany. The stable isotopic data were calibrated relative to the Vienna Pee Dee Belemnite (VPDB) using the NBS19 standard. The standard deviation of the house standard (Solnhofen limestone) over the measurement period was 0.03 ‰ for ${\mathit{\delta }}^{\mathrm{13}}$C and 0.06 ‰ for ${\mathit{\delta }}^{\mathrm{18}}$O.

Trace metal $/$ calcium in foraminifera

Up to 40 specimens ($>$ 150 $\mathrm{\mu }\mathrm{m}$) of E. selseyense and 25 specimens of E. incertum ($>$ 150 $\mathrm{\mu }\mathrm{m}$) were selected for trace metal $/$ calcium analysis. If the number of specimens was not sufficient, both species were combined. Due to poor preservation, not enough specimens were present for the analysis of Mg $/$ Ca ratios in depths $<$ 5 and $>$ 30 m c.d. Trace metal $/$ Ca in foraminifera from the Nordic and Baltic seas has been shown to be occasionally affected by diagenetic coatings (Groeneveld and Filipsson, 2013; Ezat et al., 2016). Standard cleaning procedures for trace metal $/$ Ca in foraminifera after clay removal either involve only an oxidation step to remove organic matter (Barker et al., 2003) or also include a reduction step before the oxidation (Martin and Lea, 2002) to remove (oxy)hydroxide coatings. Both methods, however, triggered a reaction on samples of Site M0059, which turned the foraminiferal fragments (dark) brownish. This suggests that an additional contaminating phase, which was not removed by reduction but does respond to oxidation, was present on the foraminiferal tests. A likely source for this are (Fe) sulfides, which are easily formed in Baltic sediments including the Little Belt and adsorb cations (Raiswell and Canfield, 1998; Hardisty et al., 2016; Van Helmond et al., 2017). Oxidation would then turn the (Fe) sulfide into an (Fe) hydroxide, which has a brownish colour. To remove this phase from the foraminiferal tests, the standard combination of reduction and oxidation was reversed. First, oxidation was performed to form the hydroxide, followed by a reduction step to remove this again. As a reduction agent 0.1 M of hydroxylamine-hydrochloride buffered in a 1 M Na acetate solution was used (Shen et al., 2001; Steinke et al., 2010). An amount of 250 $\mathrm{\mu }\mathrm{L}$ reduction agent was added to the samples; the samples were heated for 30 min at 70${}^{\circ }$ C and rinsed three times with Seralpur water. After cleaning, the samples were dissolved and centrifuged for 10 min (6000 rpm) to exclude any remaining insoluble particles from the analyses. The samples were diluted with Seralpur water before analysis with inductively coupled plasma optical emission spectrometry (Agilent Technologies, 700 Series with autosampler ASX-520 CETAC and micro-nebulizer) at MARUM, University of Bremen, Germany. Instrumental precision was monitored after every five samples using analysis of an in-house standard solution with a Mg $/$ Ca of 2.93 mmol mol${}^{-\mathrm{1}}$ (standard deviation of 0.020 mmol mol${}^{-\mathrm{1}}$ or 0.67 %). A limestone standard (ECRM752-1, reported Mg $/$ Ca of 3.75 mmol mol${}^{-\mathrm{1}}$) was analysed to allow inter-laboratory comparison (Greaves et al., 2008; Groeneveld and Filipsson, 2013). The long-term average of the ECRM752-1 standard, which is routinely analysed twice after every 50 samples in every session, is $\mathrm{3.78}\pm \mathrm{0.073}$ mmol mol${}^{-\mathrm{1}}$. Analytical precision for Elphidium spp. was 0.08 % for Mg $/$ Ca. Reproducibility using replicates of the same samples but cleaned separately was $\pm$0.17 mmol mol${}^{-\mathrm{1}}$ ($n=\mathrm{18}$). Potential contamination was monitored using Al $/$ Ca (clays), Mn $/$ Ca, and Fe $/$ Ca (both for diagenetic coatings) (see Table S3 and Fig. S1 for details).

Mg $/$ Ca ratios reach numbers that are comparable to modern core top analyses of the same species (Groeneveld and Filipsson, unpublished data). Mg $/$ Ca ratios before 4000 cal yr BP, however, are still relatively high, suggesting that the cleaning did not fully remove the contamination. Although FeS and FeS${}_{\mathrm{2}}$ are present in the studied sediments of Site M0059 (Van Helmond et al., 2017), it remains unclear what specific kind of sulfide would have been present on the foraminiferal fragments as Fe $/$ Ca did not vary between the different methods (Fig. S1). Mean Mn $/$ Ca ratios were much lower for the samples in the oxidation and reduction series (10.09 vs. 14.81 mmol mol${}^{-\mathrm{1}}$; Table S3 and Fig. S1), but these are still values that are much higher than commonly accepted for uncontaminated foraminiferal calcite (Barker et al., 2003). It is possible that the high Mn $/$ Ca ratios are truly part of the primary foraminiferal calcite and relate as such to redox conditions in the sediment. Groeneveld and Filipsson (2013) showed Mn $/$ Ca ratios as high as 10.58 mmol mol${}^{-\mathrm{1}}$ in living specimens of Globobulimina turgida from the hypoxic zone in the Gullmarn Fjord. For future work it may be required to perform a second reduction step instead of just the one step performed here.

Clumped isotopes in molluscs

The analyte for clumped isotope thermometry in the case of carbonate is CO${}_{\mathrm{2}}$ liberated by acid digestion, where ${\mathrm{\Delta }}_{\mathrm{47}}$ is the measure of the number of mass-47 isotopologues of CO${}_{\mathrm{2}}$ (primarily ${}^{\mathrm{13}}$C${}^{\mathrm{18}}$O${}^{\mathrm{16}}$O) relative to that predicted by a random arrangement of atoms (Dennis et al., 2011; Eiler, 2011). A robust relationship between ${\mathrm{\Delta }}_{\mathrm{47}}$ and the formation temperature of marine molluscs has been determined (Henkes et al., 2013). We measured the clumped isotope composition of mollusc material from Holes M0059A and D and report ${\mathrm{\Delta }}_{\mathrm{47}}$ inferred temperatures for 11 samples. Mollusc material was separated from sediment samples ($\sim$ 3 cm thickness) by wet sieving. A cold, 10 min soak of dilute H${}_{\mathrm{2}}$O${}_{\mathrm{2}}$ was used for samples in which organic matter was particularly difficult to remove. All samples were triple-rinsed in distilled water and hand-cleansed with brushes under a microscope. After species identification, samples were ground in a mortar and pestle and analysed at Johns Hopkins University (USA) on a Thermo Scientific MAT-253 isotope ratio mass spectrometer coupled to an automated acid digestion and CO${}_{\mathrm{2}}$ purification system as described in Henkes et al. (2013). Samples were analysed (three to seven replicates, as material permitted) alongside carbonate standards of varying composition and CO${}_{\mathrm{2}}$ gases approaching equilibrium at 1000 and 30 ${}^{\circ }$C, from January to March 2016. Measurements of the temperature-equilibrated CO${}_{\mathrm{2}}$ gas samples were used to create an absolute reference frame for normalization of ${\mathrm{\Delta }}_{\mathrm{47}}$ values (Dennis et al., 2011), with the calculations performed using a MATLAB^® script that accounts for temporal drift in analytical conditions (Passey et al., 2010). Finally, we used the following empirically determined relationship between mollusc ${\mathrm{\Delta }}_{\mathrm{47}}$ and temperature (Henkes et al., 2013) to convert our measured ${\mathrm{\Delta }}_{\mathrm{47}}$ to paleotemperatures: ${\mathrm{\Delta }}_{\mathrm{47}}=\mathrm{0.0327}\times {\mathrm{10}}^{\mathrm{6}}/T^{\mathrm{2}}+\mathrm{0.3286}$.

Biomarkers

A total of 40 samples collected from Holes M0059A and D were analysed for their lipid biomarker content. For this, $\sim$ 1–5 g of sediment was freeze-dried and extracted using a modified Bligh–Dyer technique (Rütters et al., 2002). An aliquot of each Bligh–Dyer extract ($\sim$ 5 mg) was separated into an apolar and polar fraction using column chromatography with activated Al${}_{\mathrm{2}}$O${}_{\mathrm{3}}$ as stationary phase and $n$-hexane $:$ dichloromethane (DCM) ($\mathrm{9}:\mathrm{1}$, $v/v$) and DCM $:$ MeOH ($\mathrm{1}:\mathrm{1}$, $v/v$) as respective eluents. The polar fractions were subsequently dried under a gentle stream of N${}_{\mathrm{2}}$ and separated into two aliquots for analyses of GDGTs and long chain diols.

To analyse the GDGT content, aliquots of the polar fractions were dissolved in $n$-hexane:2-propanol ($\mathrm{99}:\mathrm{1}$; $v/v$) and passed through a 0.4 $\mathrm{\mu }\mathrm{m}$ polytetrafluoroethylene filter. Analysis of GDGTs was performed using an Alliance 2695 (Waters, UK) high-performance liquid chromatograph coupled to a ZQ (Micromass, UK) single quadrupole mass spectrometer (MS) as detailed in Zink et al. (2016) following the analytical protocols of Hopmans et al. (2000) and Liu et al. (2012). Isoprenoid and branched GDGTs were detected using single ion recording of their protonated molecules [M$+$H]${}^{+}$ as outlined in Schouten et al. (2013a). The TEX${}_{\mathrm{86}}$, TEX${}_{\mathrm{86}}^{\mathrm{L}}$, and TEX${}_{\mathrm{86}}^{\mathrm{H}}$ values were calculated using the equations given by Schouten et al. (2002) and Kim et al. (2010): $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\text{TEX}}_{\mathrm{86}}=\left(\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right],+,\left[{\text{Cren}}^{\prime }\right]\right)/\\ & {\displaystyle }& {\displaystyle }\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right],+,\left[{\text{Cren}}^{\prime }\right]\right){\displaystyle }& {\displaystyle }{\text{TEX}}_{\mathrm{86}}^{\mathrm{L}}=\log \left(\left[\text{GDGT-2}\right],/,\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right]\right.\right.\\ & {\displaystyle }& {\displaystyle }\left..\left..+,\left[\text{GDGT-3}\right]\right)\right){\displaystyle }& {\displaystyle }{\text{TEX}}_{\mathrm{86}}^{\mathrm{H}}=\log \left(\left(\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right],\mathit{\dots },+,\left[{\text{Cren}}^{\prime }\right]\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right],+,\left[{\text{Cren}}^{\prime }\right]\right)\right).\end{array}$$ TEX${}_{\mathrm{86}}^{\mathrm{L}}$ values were transferred to absolute temperatures in the brackish–marine sediments of EZ2–4 using the Baltic Sea surface sediment calibration (Eq. 4), which is best correlated with late summer SSTs (Kabel et al., 2012). In addition, SSTs (Eq. 5) and subsurface water temperatures (Eq. 6) were calculated using the calibrations of Kim et al. (2010, 2012), respectively. In the freshwater interval of EZ1, we used the summer lake surface temperature (SLST) calibration (Eq. 7) of Powers et al. (2010) to convert TEX${}_{\mathrm{86}}$ values into water temperatures. $$\begin{array}{}& {\displaystyle }& {\displaystyle }\text{SST}=\mathrm{34.0}\cdot {\text{TEX}}_{\mathrm{86}}^{\mathrm{L}}+\mathrm{36.7},& {\displaystyle }& {\displaystyle }\text{SST}=\mathrm{68.4}\cdot {\text{TEX}}_{\mathrm{86}}^{\mathrm{H}}+\mathrm{38.6},& {\displaystyle }& {\displaystyle }\text{SST}=\mathrm{54.7}\cdot {\text{TEX}}_{\mathrm{86}}^{\mathrm{H}}+\mathrm{30.7},& {\displaystyle }& {\displaystyle }\text{SLST}=\mathrm{46.6}\cdot {\text{TEX}}_{\mathrm{86}}-\mathrm{5.6}.\end{array}$$ In order to assess whether the reconstruction of water temperatures using TEX${}_{\mathrm{86}}^{\mathrm{L}}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$ is affected by a contribution of terrestrially derived GDGTs, we calculated the branched isoprenoid tetraether (BIT) index (Hopmans et al., 2004): $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\text{BIT}=\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right]\right)/\\ {\displaystyle }& {\displaystyle }\left(\left[\text{Crenarchaeol}\right]\right.+\left[\text{GDGT-1}\right]+\left[\text{GDGT-2}\right]\\ & {\displaystyle }& {\displaystyle }\left..+,\left[\text{GDGT-3}\right]\right).\end{array}$$ For the analysis of long chain diols, an aliquot of each polar fraction was dissolved in DCM at a concentration of 2 mg mL${}^{-\mathrm{1}}$ and the mixture was silylated by the addition of N,O-bis(trimethylsilyl)trifluoroacetamide and pyridine and heating at 80 ${}^{\circ }$C for 2 h. Gas chromatography coupled to mass spectrometry (GC–MS) of long chain diols was performed using an Agilent 7890A GC instrument coupled to an Agilent 5975B MS instrument following the method described in Rampen et al. (2012). Long chain diols were quantified using selected ion monitoring of the masses $m/z$ 313 (C${}_{\mathrm{28}}$ 1,13-diols and C${}_{\mathrm{30}}$ 1,15-diols) and $m/z$ 341 (C${}_{\mathrm{30}}$ 1,13-diols and C${}_{\mathrm{32}}$ 1,15-diols). The LDI was calculated and converted to SST for sediments deposited in EZ2–4 using the equations provided by Rampen et al. (2012): $$\begin{array}{cc}{\displaystyle }\text{LDI}& {\displaystyle }=\left[{\text{C}}_{\mathrm{30}},\mathrm{1},,,\mathrm{15},\text{-diol}\right]/\left(\left[{\text{C}}_{\mathrm{28}},\mathrm{1},,,\mathrm{13},\text{-},\text{diol}\right]\right.\\ & {\displaystyle }& {\displaystyle }+\left..\left[{\text{C}}_{\mathrm{30}},\mathrm{1},,,\mathrm{13},\text{-diol}\right],+,\left[{\text{C}}_{\mathrm{30}},\mathrm{1},,,\mathrm{15},\text{-diol}\right]\right),& {\displaystyle }\text{LDI}& {\displaystyle }=\mathrm{0.033}\cdot \text{SST}+\mathrm{0.095}.\end{array}$$ To transfer relative abundances (RAs) of long chain diols found in the freshwater interval of EZ1 into surface water temperatures (SWTs), the lake calibration of Rampen et al. (2014) was applied: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\text{SWT}=\mathrm{26.8}-\mathrm{25.9}\cdot {\text{RA C}}_{\mathrm{28}}\mathrm{1},\mathrm{13}\text{diol}\\ & {\displaystyle }& {\displaystyle }-\mathrm{54.3}\cdot {\text{RA C}}_{\mathrm{30}}\mathrm{1},\mathrm{13}\text{diol}+\mathrm{7.4}\cdot {\text{RA C}}_{\mathrm{30}}\mathrm{1},\mathrm{15}\text{diol}.\end{array}$$

The DI was calculated according to Versteegh et al. (1997): $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\text{DI}=\mathrm{100}\cdot \left[{\text{C}}_{\mathrm{30}},\mathrm{1},,,\mathrm{15},\text{-diol}\right]/\\ & {\displaystyle }& {\displaystyle }(\left[{\text{C}}_{\mathrm{30}},\mathrm{1},,,\mathrm{15},\text{-diol}\right]+\left[{\text{C}}_{\mathrm{32}},\mathrm{1},,,\mathrm{15},\text{-diol}\right]).\end{array}$$ The methane index was determined as described by Zhang et al. (2011): $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\text{MI}=\left(\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right]\right),/\right.\\ {\displaystyle }& {\displaystyle }\left(\left[\text{GDGT-1}\right],+,\left[\text{GDGT-2}\right],+,\left[\text{GDGT-3}\right]\right.\\ & {\displaystyle }& {\displaystyle }\left..\left..+,\left[\text{Crenarchaeol}\right],+,\left[{\text{Cren}}^{\prime }\right]\right)\right).\end{array}$$

Results

Comparison of ${}^{\mathrm{14}}$C age model with palynology-based chronology

The ${}^{\mathrm{14}}$C-based age model for Site M0059 is discussed in detail by Van Helmond et al. (2017). The mean sedimentation rate was $\sim$ 6.5 mm yr${}^{-\mathrm{1}}$ and there are no obvious signs of hiatuses in Unit I. A comparison of the pollen record from Site M0059 with that of the varved high-resolution pollen record from Lake Belau in northern Germany (Dörfler et al., 2012; Figs. 1, 2), situated $\sim$ 100 km south of our study site, reveals a close congruency between the pollen signals in the lacustrine and marine records. Particularly clear signals occurring in both Lake Belau and M0059 pollen records are the arrival of Fagus (beech) in the catchment area $\sim$ 6000 cal yr BP, a characteristic decline in Ulmus (elm) around $\sim$ 5600 cal yr BP after preceding values between $\sim$ 5 and $\sim$ 10 %, maximum percentages of Alnus (alder) at $\sim$ 5200 cal yr BP, and maximum percentages of Fagus between $\sim$ 1800 and $\sim$ 600 cal yr BP (Fig. 2 and Dörfler et al., 2012). Grains of cultivated Poaceae taxa, i.e. cereals (Fig. 2), particularly Triticum (wheat) and Secale (rye), are consistently present after $\sim$ 1100 cal yr BP and increase significantly after $\sim$ 800 cal yr BP. This signal is also congruent with the Lake Belau record, although the cereal pollen signal at Lake Belau is stronger (Dörfler et al., 2012), probably due to higher agricultural activity close to the lake. The general palynomorph preservation is excellent over the analysed interval, and pollen concentration varies around 500 000 grains per gram sediment for most of the analysed samples. This is in accordance with the high TOC values ($>$ 2 wt %) encountered for most samples (Fig. 2). Only in the uppermost part ($<$ 100 cal yr BP), does pollen concentration decrease to values around 200 000 grains per gram.

Pollen-based climate reconstruction

The pollen-based climate reconstructions reveal a general increasing trend from $\sim$ 7800 to $\sim$ 2000 cal yr BP in both annual precipitation and annual temperature (TANN; Figs. 3, 4). The increase in the latter from $\sim$ 5 to almost 9 ${}^{\circ }$C is tied to an increase in coldest month temperatures (MTCOs; $\sim$ $-$8 to almost 0 ${}^{\circ }$C), while warmest month temperatures (MTWAs) show a decreasing trend starting at $\sim$ 7000 cal yr BP ($\sim$ 19 ${}^{\circ }$C), which ends at ca. 2000 cal yr BP ($\sim$ 17 ${}^{\circ }$C). Between $\sim$ 2000 and $\sim$ 1000 cal yr BP, MTWAs increase to $>$ 19 ${}^{\circ }$C, and subsequently decrease again to $<$ 17 ${}^{\circ }$C. The MTWA maximum around $\sim$ 1000 cal yr BP is coeval with relatively low annual precipitation values ($>$ 640 mm compared to $>$ 800 mm earlier and subsequently).

Environmental zones of the Holocene record at Site M0059.

Zone	Depth (m c.d.)	Age (cal yr BP)	Characteristics
EZ4	$\sim$ 6 to top	$\sim$ 1000 to present	Increase in brackish diatoms, Gymnodinium cysts, decrease in diatom abundance
EZ3	$\sim$ 22 to $\sim$ 6	$\sim$ 4100 to $\sim$ 1000	Increase in diatom and Chaetoceros spore abundance,
			decreased indicators of marine influence
EZ2	$\sim$ 49 to $\sim$ 22	$\sim$ 7400 to $\sim$ 4100	Increased abundance of indicators of marine influence
EZ1	$\sim$ 53 to $\sim$ 49	$\sim$ 7800 to $\sim$ 7400	High absolute and relative values of freshwater proxies

Sedimentological units, lithology (Andrén et al., 2015a), and position of ${}^{\mathrm{14}}$C ages (Van Helmond et al., 2017) vs. depth and IODP Site M0059 temperature proxies plotted vs. age (extrapolated below $\sim$ 48.64 m c.d.). (a) Environmental zones (EZs) are based on diatom data. (b) Pollen-inferred climate reconstructions are based on the modern analogues technique, with the shaded area reflecting the range of the best 10 analogues. (c) SSTs are reconstructed with the long chain diol index (LDI) and TEX${}_{\mathrm{86}}^{\mathrm{L}}$, both reflecting a summer signal. (d) ${\mathrm{\Delta }}_{\mathrm{47}}$-inferred temperatures are based on the clumped isotope composition of mollusc shells. (e) Mg $/$ Ca-based bottom water temperature (BWT) is reconstructed using the benthic foraminifera E. selseyense and E. incertum; stable oxygen isotopes based on the benthic foraminifer E. selseyense. (f) Climate phases and chronozones.

[Figure omitted. See PDF]

Marine ecosystem changes

The Holocene record of Site M0059 was divided into four EZs (Table 2) based on statistically assessed changes in diatom assemblages (Sect. 3.3.1) and congruent signals of other marine proxies (benthic foraminifera – Sect. 3.3.2, ostracods – Sect. 3.3.4, and aquatic palynomorphs – Sect. 3.3.5). While EZ1 is characterized by freshwater conditions, EZ2–4 indicate a brackish–marine depositional environment that persisted in the Baltic Sea during the Littorina Sea stage. EZ1 corresponds with lithological Subunit Ib (49.37–53.57 m c.d.), while EZ2–4 correspond with Subunit Ia.

Diatoms

Relative abundances of salinity-based diatom assemblages, as well as ADA plotted as the number of valves per gram dry weight of sediment ($v/$ dgw), absolute abundance of CRSs ($v/$ gdw), and the ratio of benthic to planktonic diatoms are presented in Fig. 3. The division into four EZs as described above was based on salinity preference, ADA, CRS abundance, and the ratio of benthic to planktonic diatoms (B $:$ P). These are in general agreement with signals reflected by other proxies (see below).

EZ1 ($\sim$ 7800 to $\sim$ 7400 cal yr BP) is based on three diatom samples characterized by dominance of freshwater taxa, rare CRS abundance, (average $\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{5}}$), low ADA (average $\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{7}}$), and low B $:$ P (0.48). EZ2 ($\sim$ 7400 to $\sim$ 4100 cal yr BP) is defined by a distinct increase in brackish and marine species at the expense of freshwater diatoms. In addition, there is a clear increase in the B $:$ P ratio (average 4.17), ADA (average $\mathrm{4.65}\times {\mathrm{10}}^{\mathrm{7}}$), and CRS absolute abundance (average $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{7}}$). EZ2 contains both the highest B $:$ P and the greatest range in values of B $:$ P. EZ3 ($\sim$ 4100 to $\sim$ 1000 cal yr BP) shows an increase in marine, brackish–marine, and freshwater species at the expense of brackish species. The B $:$ P ratio (average 2.65) declines some in EZ3; however, ADA (average $\mathrm{5.1}\times {\mathrm{10}}^{\mathrm{7}}$) and CRS absolute abundance (average $\mathrm{2.2}\times {\mathrm{10}}^{\mathrm{7}}$) are highest in this interval compared to the rest of the record. Marine and brackish–marine species increase in abundance from the base of this zone to $\sim$ 1100 cal yr BP concurrent with a decrease in brackish species. Both diatom and CRS absolute abundance increase from the base of this zone to $\sim$ 1200 cal yr BP. B $:$ P increases from the base of the zone to $\sim$ 1000 cal yr BP.

EZ4 ($\sim$ 1000 cal yr BP to present) begins with an increase in brackish species. The relative abundance of marine and brackish–marine taxa, as well as B $:$ P ratio (average 2.07), declines rather gradually throughout this interval, whereas ADA (average $\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{7}}$) and CRS abundance (average $\mathrm{5.3}\times {\mathrm{10}}^{\mathrm{6}}$) declines rapidly at $\sim$ 1000 cal yr BP and remains at low levels until the core top.

The significance of the differences between the four identified EZs was assessed in regards to the percent abundance of all five salinity affinities, as well as ADA, CRS abundance, and B $:$ P. ANOVA reveals a significant difference between the means of all of these measures (Table S4). The Tukey–Kramer pairwise test was then used to evaluate the differences between individual EZs. Statistical differences between adjacent zones, representing change through time, are reported in Table S5. Pairwise statistical analysis reveals significant differences between EZ1 and EZ2 in terms of all ecological metrics, with the exception of ADA. EZ2 and EZ3 can be distinguished statistically on the basis of the percent of diatoms with a brackish affinity, which decreases from EZ2 to EZ3, as well as CRS, which increases from EZ2 to EZ3. EZ3 can be significantly distinguished from EZ4 in terms of ADA and CRS abundance, both of which decrease. Finally, EZ3 does not have a significant change in salinity with respect to EZ4; however, there is a shift in the species present between these zones.

Benthic foraminiferal assemblages

Benthic foraminifera are found in sediments covering the interval since $\sim$ 7400 cal yr BP until today. The foraminiferal assemblages comprise very few species, with Elphidium selseyense (formerly named E. excavatum forma selseyensis; see Darling et al., 2016) being dominant. Elphidium incertum is also found throughout the entire core, while all other species occur only in specific intervals. In the lower part of EZ2, from $\sim$ 7400 to $\sim$ 6600 cal yr BP, the foraminiferal assemblages are characterized by very high frequencies of E. incertum and, around $\sim$ 6900 cal yr BP, relatively higher abundance of Ammonia beccarii ($\sim$ 15% compared to $<$ 1 % in samples above and below; Fig. 3). After $\sim$ 6600 cal yr BP the relative abundance of A. beccarii is strongly reduced, while E. incertum continues to dominate and increases in abundance until $\sim$ 5700 cal yr BP. During most of the interval since $\sim$ 5700 cal yr BP (top zone EZ3 to zone EZ1), only E. selseyense and E. incertum are found consistently, with only minor, short-term occurrences of Elphidium albiumbilicatum and Elphidium magellanicum at $\sim$ 5900 cal yr BP and $\sim$ 4850 cal yr BP, as well as a short-term peak of A. beccarii at $\sim$ 3800 cal yr BP. The samples representing the past $\sim$ 460 cal yr BP contain no or only very few foraminifera (E. selseyense and a few specimens of E. incertum and E. albiumbilicatum; Fig. 3).

Benthic foraminifer isotopes and trace metal $/$ calcium

The difference in average Mg $/$ Ca between E. selseyense and E. incertum was not significant (1.52 vs. 1.53 mmol mol${}^{-\mathrm{1}}$; Fig. 1S); thus, for the remainder only results of E. selseyense are included, as they are of higher resolution. Mg $/$ Ca values decrease from 2.96 mmol mol${}^{-\mathrm{1}}$ at $\sim$ 4000 cal yr BP to values between 0.62 and 1.10 mmol mol${}^{-\mathrm{1}}$ between $\sim$ 2200 and 800 cal yr BP (Fig. 4). As a species-specific calibration is absent for Elphidium spp., we employed the Mg $/$ Ca temperature calibration for Melonis barleaanum, which has been shown to give reasonable temperatures for the Baltic Sea area (Kristjánsdóttir et al., 2007; Anjar et al., 2012). Calculated bottom water temperatures decrease from 10.9 ${}^{\circ }$C at $\sim$ 4000 cal yr BP to 0–4 ${}^{\circ }$C between $\sim$ 2200 and $\sim$ 800 cal yr BP.

${\mathit{\delta }}^{\mathrm{18}}$O values show an initial increase from $-$1.2 to 0.4 ‰ from $\sim$ 7400 to $\sim$ 6700 cal yr BP. ${\mathit{\delta }}^{\mathrm{18}}$O values average $-$0.2 ‰ for the interval $\sim$ 7400 to $\sim$ 3700 cal yr BP and increase to an average of 0.4 ‰ afterwards, suggesting either a more gradual decrease in bottom water temperature, an increase in bottom water salinity, or both (Figs. 4, S1). ${\mathit{\delta }}^{\mathrm{13}}$C values increase from $-$3.2 ‰ at $\sim$ 7000 cal yr BP to $-$2 ‰ at $\sim$ 5000 to $\sim$ 4500 cal yr BP (Fig. 3). Afterwards, values oscillate around an average of $-$1.3 ‰, suggesting improved ventilation of the water column.

Ostracods

Ostracod abundance was relatively low ($\sim$ 5 valves per 30 cm${}^{\mathrm{3}}$, 23 samples barren), but preservation was good in most samples. EZ1 is entirely comprised of freshwater species. Abundant freshwater taxa Candona spp. and Cytherissa lacustris were identified in the interval between $\sim$ 7800 and $\sim$ 7400 cal yr BP. The rest of the record contains marine and brackish–water ostracod taxa. The ostracod assemblages imply a subdivision of EZ2 (Fig. 3). EZ2a ($\sim$ 7400 to $\sim$ 6800 cal yr BP) is mainly represented with two shallow water taxa, Palmoconcha spp. and Hirschmannia viridis. The abundance peak of shallow water taxa of 6 valves per 30 cm${}^{\mathrm{3}}$ is observed at $\sim$ 7000 cal yr BP. EZ2b ($\sim$ 6800 to $\sim$ 3600 cal yr BP) contains both shallow water and deeper-living taxa, with deeper-living taxa being slightly more abundant and showing two abundance peaks in this interval. In the lower part ($\sim$ 6,800 to $\sim$ 5600 cal yr BP), a peak of 7 valves per 30 cm${}^{\mathrm{3}}$ is observed for deeper-living taxa at $\sim$ 6500 cal yr BP. A short interval ($\sim$ 5600 to $\sim$ 5300 cal yr BP) is barren of ostracods. The upper part of the interval, $\sim$ 5300 to $\sim$ 3600 cal yr BP, includes a peak in deeper-living taxa of 5 valves per 30 cm${}^{\mathrm{3}}$ at $\sim$ 5000 cal yr BP. Sarsicytheridea bradii, Cytheropteron latissimum, Elofsonella concinna, Robertsonites tuberculatus, and Paraciprideis sp. are the dominant taxa in this interval. These are typical representatives of Arctic and/or Boreal shallow water shelf faunas (Stepanova et al., 2007). An open sea environment and a salinity of 14–16 or higher is expected for this group of species (Frenzel et al., 2010). Shallower water species like Palmoconcha spp. and Leptocythere spp. are the most abundant among shallow water taxa.

EZ3 ($\sim$ 3600 to $\sim$ 1000 cal yr BP) is characterized by an overall increase in abundance and predominance of shallow water taxa (reaching a peak of 15 valves per 30 cm${}^{\mathrm{3}}$ at $\sim$ 3500 cal yr BP). The most abundant taxa here are the shallow water genera Leptocythere spp. and Palmoconcha spp. These taxa are present in the underlying EZ2, but become dominant here. These genera inhabit lagoons and open sea environments with a salinity above 7–12. A short interval between $\sim$ 2600 and $\sim$ 2200 cal yr BP is barren of ostracods.

EZ4 (covering the past $\sim$ 1000 yr) is marked by a decrease in the taxonomic diversity; it is almost solely represented by Sarsicytheridea bradii, a taxon that can be found in open sea environments and at a salinity as low as 7 (Frenzel et al., 2010).

Aquatic palynomorphs

Percentages of dinocysts and freshwater algae (Fig. 3) are based on the total pollen sum (thus representing a dinocyst $/$ pollen or a freshwater algae $/$ pollen ratio, respectively). The freshwater algae Pediastrum and Botryococcus are present in the entire record, but values $>$ 3 % are only found $>$ 7700 cal yr BP (Fig. 3). We have combined the dinocysts of phototrophic taxa belonging to gonyaulacoids into one group (“G cysts”; predominately comprising Protoceratium reticulatum $/$ Operculodinium centrocarpum, Spiniferites spp., and Lingulodinium machaerophorum), and microreticulate, gymnodinioid cysts of Gymnodinium in another. Cysts of heterotrophic taxa are rare and thus not depicted in Fig. 3. While not all Gymnodinium cysts were identified to species level, most of the encountered Gymnodinium cysts probably belong to G. nolleri (referred to as Gymnodinium cf. nolleri in the following). The G cyst percentages are $>$ 3 % between $\sim$ 7800 and $\sim$ 7400 cal yr BP (EZ1) but increase to $>$ 5 % around $\sim$ 7400 cal yr BP. Thereafter, values of this proxy stay rather constant at $\sim$ 10 % until $\sim$ 4000 cal yr BP (EZ2). For the latest $\sim$ 4000 cal yr BP, G cysts percentages vary around 5 %, with particularly low values around 500 cal yr BP (EZ4). Gymnodinium cf. nolleri cysts are almost completely absent until $\sim$ 7000 cal yr BP, increase rapidly at $\sim$ 6800 cal yr BP, but in the following years decline and are almost absent at $\sim$ 4000 cal yr BP. A second interval of high Gymnodinium cf. nolleri abundances occurs between $\sim$ 3000 and $\sim$ 300 cal yr BP, with maximum values at $\sim$ 1000 and $\sim$ 200 cal yr BP.

Mollusc clumped isotopes

The analysed mollusc shell fragments were assigned to the species Abra alba, Arctica islandica, Corbula gibba, and Macoma balthica (see Table S6). The oldest analysed mollusc is from EZ2 ($\sim$ 6800 cal yr BP) and the youngest analysed shell fragment is from EZ4 ($\sim$ 140 cal yr BP). The temperatures inferred from clumped isotopes range from $\sim$ 2 to $\sim$ 12 ${}^{\circ }$C (Fig. 4). In EZ2, temperatures vary between 5 and 11.8 ${}^{\circ }$C, with a majority of temperatures around $\sim$ 8 ${}^{\circ }$C. The two temperatures obtained in EZ3 are 2.1 and 4.8 ${}^{\circ }$C, and the single temperature obtained in EZ4 is 7.7 ${}^{\circ }$C (Fig. 4).

Biomarkers

The DI shows consistently low values of $\mathrm{55.2}\pm \mathrm{2.7}$ in EZ1 but rapidly increases to a value of 90.5 with the establishment of brackish–marine conditions at the start of the Littorina Sea phase (Fig. 3). Thereafter, DI values gradually decline and show a minimum of 54.3 around $\sim$ 4000 cal yr BP, after which values gradually increase again until $\sim$ 3000 cal yr BP, when they peak at 80.1. In the following, the DI remains rather constant, averaging $\mathrm{77.0}\pm \mathrm{3.0}$ for the remainder of the record with an exception found around $\sim$ 700 cal yr BP, at which the DI decreases to values as low as 61.7.

TEX${}_{\mathrm{86}}$-based SLSTs generally evidence a warming trend throughout EZ1 with temperatures increasing from 12.9 to 14.8 ${}^{\circ }$C, except in the uppermost sample of this interval in which a drop in temperature to 12.5 ${}^{\circ }$C is observed (Fig. 4). At the transition to the marine-brackish phase of EZ2, TEX${}_{\mathrm{86}}^{\mathrm{L}}$-inferred temperatures rapidly increase to 19.5 ${}^{\circ }$C and thereafter show a gradual decline to yield minimum SSTs of $\sim$ 14 ${}^{\circ }$C in the topmost part of EZ2 ($\sim$ 5000 cal yr BP). In the following, TEX${}_{\mathrm{86}}^{\mathrm{L}}$ SSTs gradually increase and cumulate in a peak SST of 26 ${}^{\circ }$C around $\sim$ 1000 cal yr BP, after which SSTs decline again and are at a minimum of 14 ${}^{\circ }$C around $\sim$ 300 cal yr BP. The uppermost part of the record ($<$ 300 cal yr BP) is characterized by increasing SSTs that maximize at 21.1 ${}^{\circ }$C in the sample closest to the core top.

Branched isoprenoid tetraether index values are generally high in the freshwater interval of EZ1 with an average of $\mathrm{0.93}\pm \mathrm{0.04}$ (Fig. 4). With the establishment of brackish–marine conditions (base of EZ2), the branched isoprenoid tetraether index rapidly declines to values of $\mathrm{0.29}\pm \mathrm{0.04}$ and shows only little variation in the organic-rich deposits of the Littorina Sea stage. Likewise, values of the methane index are, with an average of $\mathrm{0.42}\pm \mathrm{0.02}$, comparatively high in EZ1 but at the transition to EZ2 decrease substantially to an average of $\mathrm{0.12}\pm \mathrm{0.05}$ and stay low throughout the remainder of the record.

Long chain diol distributions in EZ1 yield SWTs varying from 18.5 to 21.3 ${}^{\circ }$C (average $\mathrm{19.7}\pm \mathrm{1.3}$ ${}^{\circ }$C). At the transition to EZ2, SSTs based on the LDI rapidly increase to a maximum of 24.3 ${}^{\circ }$C. Subsequently, reconstructed SSTs gradually decline and show a minimum of 15.9 ${}^{\circ }$C at the transition from EZ3 to EZ4 ($\sim$ 1000 cal yr BP) that is followed by a short-lived warming event with peak temperatures of 20.8 ${}^{\circ }$C around $\sim$ 900 cal yr BP. A second minimum in LDI-based SSTs is observed between $\sim$ 800 and $\sim$ 400 cal yr BP, after which temperatures again increase to 18.3 ${}^{\circ }$C close to the core top.

Discussion

The complex climatic development during the Holocene in the Baltic Sea region (Björck, 2008; Zillén et al., 2008; Andrén et al., 2011) enables us to compare and validate several temperature reconstruction methods, with pollen reflecting terrestrial conditions, the LDI and the TEX${}_{\mathrm{86}}$ with its different derivatives indicating changes in the surface water, and mollusc clumped isotopes and foraminiferal Mg $/$ Ca and ${\mathit{\delta }}^{\mathrm{18}}$O allowing the reconstruction of bottom water conditions. Of the various salinity proxies used in this study, diatoms, cysts of phototrophic dinoflagellates, and the DI mainly reflect surface water conditions, while foraminifera and ostracods indicate changes in the bottom water.

Salinity and productivity changes

All biogenic proxies indicate a phase of low salinity during EZ1 (before $\sim$ 7400 cal yr BP); diatom and ostracod data suggest freshwater conditions, and diatom data additionally indicate low productivity in the lowermost part of EZ1 (Fig. 3), though this is not fully consistent with relatively high TOC values in the lowermost part of EZ1 (Fig. 2) and results of Van Helmond et al. (2017). Among the aquatic palynomorphs, remains of the green algae Pediastrum sp. and Botryococcus sp., together with low occurrences of marine dinoflagellate cysts compared to pollen grains, also indicate a significant freshwater influence. Likewise, the DI expresses the lowest values in this interval that are on the same order of magnitude reported from other lacustrine environments (Versteegh et al., 1997). Comparison with the varved pollen record from Lake Belau (Fig. 2; Dörfler et al., 2012) suggests that EZ1 is probably not older than $\sim$ 9000 cal yr BP, as indicated by the presence of Alnus and Ulmus ($\sim$ 5 % each), for example, in the lowermost samples. Our findings are thus congruent with those of Bennike and Jensen (2011) and Van Helmond et al. (2017), who suggested that during the early Holocene a large lake developed in the southern part of the Little Belt region.

Fully marine conditions in the southern Kattegat were already established by $\sim$ 9300 cal yr BP (Bendixen et al., 2016), but not sooner than $\sim$ 8500 to $\sim$ 7500 cal yr BP in the central Baltic Proper (Björck, 2008). This implies that during EZ1 the Little Belt may have been connected to the Baltic Proper but not to the Kattegat. Bennike and Jensen (2011) suggested that the transition to brackish conditions may have started at $\sim$ 8500 cal yr BP in the Little Belt, although the oldest dated marine shell from the Little Belt is dated to $\sim$ 7700 cal yr BP (Bennike and Jensen, 2011), consistent with our findings.

The change to more marine conditions in EZ2 ($\sim$ 7400 to $\sim$ 4100 cal yr BP) occurred quickly, as already implied by Andrén et al. (2015b), since most proxies show a rapid decrease or disappearance of freshwater indicators (e.g. diatoms, green algae, ostracods). At the same time, the applied salinity proxies indicate a shift to more saline conditions, e.g. increasing CRS abundance, as most Chaetoceros require brackish to marine conditions (Snoeijs, 1993; Snoeijs and Vilbaste, 1994; Snoeijs and Potapova, 1995; Snoeijs and Kasperoviciene, 1996; Snoeijs and Balashova, 1998), and a significant rise in the DI values. Maximum abundances of dinocysts of phototrophic taxa occurred in EZ2 and are coeval with abundance maxima at the southern and southeastern coast of Sweden found by Ning et al. (2017) or Yu and Berglund (2007) and in the Gotland and the Fårö basins (Brenner, 2005; Willumsen et al., 2013).

We suggest that maxima in Gymnodinium cf. nolleri cysts between $\sim$ 6800 and $\sim$ 5200 cal yr BP may indicate more saline conditions, though they may also be related to increased temperature (e.g. Thorsen et al., 1995; Thorsen and Dale, 1997; see below), evidenced here by high SSTs from TEX${}_{\mathrm{86}}^{\mathrm{L}}$ and the LDI in this interval. The Gymnodinium peak at $\sim$ 6800 cal yr BP has, to our knowledge, not yet been described from the Baltic Sea and does not seem to have occurred in the Baltic Proper (Brenner, 2005; Yu and Berglund, 2007; Ning et al., 2017). Thorsen et al. (1995) and Harland and Nordberg (2011) describe later Gymnodinium mass occurrences from the Skagerrak and the Kattegat, but sediments of equivalent age containing our Gymnodinium peak at $\sim$ 6800 cal yr BP were not recovered in the framework of those studies.

Increasing content of foraminifera per gram of sediment within a short interval also points to a salinity increase, as well as the near disappearance of the foraminiferal species A. beccarii ($<$ 1 %) at $\sim$ 6600 cal yr BP and the relatively high frequencies of E. incertum until $\sim$ 5700 cal yr BP. Ammonia beccarii is a euryhaline species (Murray, 2006), but it increases in frequency in lower-salinity environments in the Kattegat–Baltic Sea region (e.g. Lutze, 1974; Seidenkrantz, 1993). Elphidium incertum is a shallow infaunal species that favours sandy substrate in brackish, inner shelf areas (salinity $>$ 25) where it is particularly frequent just below the halocline in stratified waters (Lutze, 1974; Darling et al., 2016). Foraminiferal ${\mathit{\delta }}^{\mathrm{18}}$O values are also in line with increasing salinity, increasing from $\sim$ $-$1.2 to $\sim$ $+$0.4 ‰ (Fig. 3). The following drop in E. incertum abundance and consequent increase in E. selseyense in the upper part of EZ2 may potentially be ascribed to somewhat reduced bottom water salinities, as E. selseyense is an opportunistic species that is widespread in tidal to shelf areas with relatively large variations in temperature and salinity (Murray, 2006; Darling et al., 2016). As stable ${\mathit{\delta }}^{\mathrm{18}}$O values do not suggest a strong change in salinity between $\sim$ 6600 and $\sim$ 3700 cal yr BP, the shift in benthic foraminiferal fauna may also be linked to a reduction in bottom water oxygenation, as E. incertum is also reported to require relatively high oxygen concentrations (Lutze, 1974; Murray, 2006). This, however, is not corroborated by the results of Van Helmond et al. (2017), who observed a significant reduction in bottom water oxygenation only at the onset of this interval.

In terms of ostracod fauna, EZ2 can be subdivided into two subzones, with the lower interval $\sim$ 7400 to $\sim$ 6800 cal yr BP representing a transitional early marine stage environment with a salinity above 10–14 and the second interval $\sim$ 6800 to $\sim$ 3900 cal yr BP with a salinity above 14–16. Both ostracod intervals within EZ2 reflect the onset of the Littorina Sea and the salinity increase pattern as evidenced by the highest concentration of euryhaline taxa. At the same time, the transition from lower to higher salinities must have been gradual, with the lower part of EZ2 containing higher numbers of shallower water meso- and polyhaline taxa decreasing upwards. This implies that bottom water salinity change was not abrupt but occurred more gradually. There was only a little change in diatom-based primary productivity compared to EZ1 and likely a deepening of the photic zone or shallowing of the water column as indicated by the sharp increase in B $:$ P. Decreasing DI values and Gymnodinium percentages during the second half of EZ2 imply slightly decreasing surface water salinity between $\sim$ 5500 and $\sim$ 4100 cal yr BP. This is also corroborated by a decline in the marine and brackish–marine diatom species. The interval represented by EZ2 ($\sim$ 7400 to $\sim$ 4100 cal yr BP) is thus characterized by brackish to marine conditions with higher salinity at the onset of the Littorina Sea stage and slowly decreasing surface water salinity afterwards, while salinity slightly increased in the bottom water. Gustafsson and Westman (2002) regard the interval between $\sim$ 6500 and $\sim$ 5000 cal yr BP as the most marine period in the Baltic Sea in general and suggest that low precipitation could be the main cause for the particularly high salinity during this interval compared to the rest of the record. While our data point to an earlier maximum in salinity (Fig. 3), a connection between low precipitation and high salinity is in accordance with the generally low precipitation values reconstructed via the MAT ($<$ 600 mm yr${}^{-\mathrm{1}}$ excluding one sample until 4000 cal yr BP) for EZ2 compared to values from the late 20th century, which are generally $>$ 600 mm yr${}^{-\mathrm{1}}$ and often reach $\sim$ 800 mm yr${}^{-\mathrm{1}}$ (Omstedt et al., 1997). In EZ3 ($\sim$ 4100 to $\sim$ 1000 cal yr BP), however, high salinity seems not to be coupled with low precipitation (Fig. 3, see below).

At the transition between EZ2 and EZ3 ($\sim$ 4100 to $\sim$ 1000 cal yr BP), salinity was probably lower for a short interval, as indicated by the low values of the DI and a short increase in the abundance of A. beccarii, which appears to be a relevant signal even though it is only reflected in one sample. In EZ3, the DI, diatom and benthic foraminiferal assemblages, benthic ${\mathit{\delta }}^{\mathrm{18}}$O, and to some degree ostracod occurrences (with particularly high values at $\sim$ 3500 cal yr BP and the complete disappearance of the very shallow water ostracod taxa) indicate another increase in both surface and bottom water salinities between $\sim$ 4100 and $\sim$ 3000 cal yr BP. The composition of the diatom assemblages implies that EZ3 includes another particularly marine phase of the Little Belt region. Peak marine conditions, as indicated by diatom species' salinity affinities, occur from $\sim$ 1200 to $\sim$ 1000 cal yr BP, at the transition to EZ4. This is evident in a statistically significant decrease in brackish diatom species, while the proportions of marine and brackish–marine diatom species both increase throughout this interval. The DI also shows high values at the transition between EZ3 and EZ4 compared to most of the record, though relatively high values had already been established since $\sim$ 3000 cal yr BP. In this case, there is no clear coupling between increased salinity and decreased precipitation, but several temperature proxies indicate high temperatures around $\sim$ 1000 cal yr BP, pointing to higher evaporation.

EZ3 shows diatom B $:$ P ratios that are statistically similar to EZ2, implying a deep photic zone or shallow water column (supported by a peak in shallow water ostracod taxa in EZ3). In addition, this zone is characterized by particularly high primary productivity, as indicated by ADA, CRS, and somewhat increased benthic foraminifera and ostracod abundances. Precipitation was also increased, likely causing increased freshwater run-off indicated by the pollen assemblages (Fig. 3).

Following the peak marine interval, EZ4 (since $\sim$ 1000 cal yr BP) begins with increases in brackish diatom taxa and decreases in marine and brackish–marine diatom species to the core top. These changes between EZ3 and EZ4 imply that species shifts are more likely related to variations in primary productivity and nutrient conditions than salinity. However, the DI and Gymnodinium percentages also decrease during this interval. In the case of Gymnodinium, the decrease could be tied to temperature changes (according to Thorsen and Dale, 1997; see below). This might also be the case for the DI as a slight correlation between this proxy and water temperature has been reported previously (Rampen et al., 2012). Productivity decline from EZ3 to EZ4 is demonstrated by statistically significant decreases in ADA and CRS abundance. The uppermost diatom sample in this interval has the lowest recorded ADA as well as the lowest B $:$ P of any brackish water interval in this study. This sample ($\sim$ 5 cal yr BP) likely represents modern human influence to the Baltic Sea, i.e. eutrophication. Eutrophication would cause a decline in overall diatom abundance, as diatoms are often replaced by cyanobacteria (O'Neil et al., 2012; Michalak et al., 2013) and dinoflagellates (Wasmund and Uhlig, 2003) in cases of high levels of nutrient loading. In turn, this would lead to a shallowing of the photic zone resulting from sunlight being rapidly absorbed by cyanobacteria and dinoflagellates in the upper water column, as seen in the low B $:$ P ratio. The relative and absolute number of dinocysts decreases though in the uppermost palynology sample; thus, denser sampling will be needed to determine the anthropogenic impact on the Baltic Sea ecosystems in future studies. Human influence is, however, indicated by the significant increase in pollen of cultivated Poaceae taxa (Fig. 2), starting at $\sim$ 1000 cal yr BP, with persistent occurrences since $\sim$ 700 cal yr BP. The coeval decrease in pollen of tree taxa like Alnus and Fagus may be an anthropogenic deforestation signal due to increased agricultural activity (Fig. 2).

The Gymnodinium maxima encountered in EZ3 and EZ4 are probably temporally related to two mass occurrences of Gymnodinium in the eastern Kattegat as described by Harland and Nordberg (2011) for $\sim$ 2000 and $\sim$ 300 cal yr BP. The occurrence of Gymnodinium cysts in the sediments from the Skagerrak–Kattegat is tied to warm water conditions according to Thorsen and Dale (1997). This connection would be congruent with our findings that the maximum at the transition from EZ3 to EZ4 is coeval with the Medieval Warm Period (MWP, also called the Medieval Climate Anomaly) and with high temperatures. However, while Thorsen and Dale (1997) assign the cysts to G. catenatum, Harland and Nordberg (2011) claim that the Gymnodinium cysts found in Holocene sediments from the Skagerrak and Kattegat rather belong to G. nolleri (in accordance with our results) and that the known mass occurrences may not be tied to high temperatures. Considering our findings, Gymnodinium mass occurrences may as well be tied to particularly high salinity as implied by the good correlation with other salinity proxies, e.g. the DI (which varies around values of $\sim$ 80, Fig. 3), but higher-resolution studies beyond the scope of this study would be needed to be able to prove such a connection statistically. EZ4 is characterized by a taxonomic diversity decrease in ostracods; this assemblage is almost solely represented by one species, S. bradii. This change may imply an unfavourable habitat for ostracods probably related to shallowing of the site and low oxygen conditions.

Taken together, the different salinity proxies provide a consistent picture of how salinity in the water column changed in the Little Belt during the last $\sim$ 8000 cal yr BP, with freshwater conditions preceding $\sim$ 7400 cal yr BP and more saline conditions thereafter. Precipitation (reconstructed via pollen assemblages) may have been one of the factors influencing salinity, particularly between $\sim$ 8000 and $\sim$ 4000 cal yr BP. Slight discrepancies between the salinity proxies are likely explained by (i) differences between surface and bottom waters and (ii) other factors influencing microfossil assemblages such as productivity and temperature. While qualitative estimates are congruent for most proxies, quantitative salinity estimates could rarely be made in the framework of this study.

Temperature reconstructions

Inorganic- (e.g. Mg $/$ Ca and ${\mathit{\delta }}^{\mathrm{18}}$O of benthic foraminifera, clumped isotopes) and organic-based temperature proxies (e.g. TEX${}_{\mathrm{86}}$, LDI) have become indispensable tools in paleoenvironmental research as they provide quantitative information of past climate variations. However, proxy-based temperature estimates may vary significantly and/or show different trends even within the same sample set. This apparent mismatch between proxies severely complicates their application in paleoenvironmental and paleoclimate research. Ideally, multiple proxies should therefore be used simultaneously to provide independent paleotemperature records (e.g. Butruille et al., 2017; Jonas et al., 2017).

Pollen-based reconstruction of seasonal and mean annual air temperatures as well as estimates of surface water temperature variations using lipid paleothermometers were obtained continuously over the investigated sequence of Site M0059, but clumped isotope and Mg $/$ Ca data are missing, in particular from the freshwater interval reflected in EZ1 and in the lowermost section of EZ2 due to the lack of sufficient carbonaceous micro- and macrofossils. Despite the low resolution of some of the investigated parameters, there is a certain degree of similarity in temperature trends observed across the different proxies. The pollen-based temperatures reconstructed for the two uppermost samples from the interval between 100 cal yr BP and the present show an average value of $\sim$ 0.1 ${}^{\circ }$C for MTCO and of $\sim$ 16.4 ${}^{\circ }$C for MTWA, which is close to mean temperatures observed in Denmark during January ($\sim$ 1.5 ${}^{\circ }$C) and August ($\sim$ 17.2 ${}^{\circ }$C; Danish Meteorological Institute, 2017). The most complete temperature records are provided by the sedimentary distribution of isoprenoid GDGTs (TEX${}_{\mathrm{86}}$ and its derivatives) and long chain diols (in the form of the LDI).

Although absolute temperature estimates differ between the two proxies, trends are overall very similar in freshwater interval EZ1, showing a rather constant SWT of $\mathrm{13.5}\pm \mathrm{1.0}$ ${}^{\circ }$C (TEX${}_{\mathrm{86}}$) and $\mathrm{19.7}\pm \mathrm{1.3}$ ${}^{\circ }$C (long chain diols), followed by a rapid temperature increase in both lipid paleothermometers to maximum values of 19.5 ${}^{\circ }$C (TEX${}_{\mathrm{86}}^{\mathrm{L}}$) and 24.5 ${}^{\circ }$C (LDI) at the transition from EZ1 to EZ2 at $\sim$ 7400 cal yr BP (Fig. 4). An interval of generally higher mean annual air and water temperatures has previously been described from other records of the Baltic Sea region such as Lake Flarken, Sweden (Seppä et al., 2005) or the Skagerrak (Butruille et al., 2017; Krossa et al., 2017) and is consistent with the Holocene Climate Optimum (HCO, also Holocene Thermal Maximum; e.g. Seppä et al., 2005). The pollen-based temperature reconstructions imply that in the terrestrial realm, high summer temperatures were already established before $\sim$ 7400 cal yr BP, but that temperatures were very low during winter. In EZ2, the MTWA and SWTs based on the long chain diol index indicate a gradual and slow cooling trend by $\sim$ 1 to 4 ${}^{\circ }$C towards the top of EZ2. A similar cooling trend is also observed in the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ record, which despite generally high surface water temperatures, shows a decline in temperature by $\sim$ 6 ${}^{\circ }$C during EZ2. The general order of magnitude of the pollen-based temperatures for EZ2 ($\sim$ $-$8 to $\sim$ 1.5 ${}^{\circ }$C MTCO and $\sim$ 19 to $\sim$ 17.5 ${}^{\circ }$C MTWA) is in accordance with pollen-based climate reconstructions from central Germany, for example (Kühl and Moschen, 2012). Furthermore, the overall trends in all temperature reconstructions are congruent with trends reported by Davis et al. (2003) and Mauri et al. (2015) for western central and northern Europe based on pollen records. However, the presence of pollen of thermophilous species in Danish records from the Boreal and the Atlantic period ($\sim$ 9000 to 5000 cal yr BP) implies that winters have been mild enough for these taxa (e.g. Iversen, 1944).

In EZ3, the SWTs based on the LDI, pollen-inferred MTWA, and Mg $/$ Ca ratios suggest a continuation of the cooling that agrees with trends observed in other records (Seppä et al., 2005), while MTCO, TANN, and the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ indicate a gradual warming that cumulates in maximum air and surface water temperatures at the transition to EZ4 at $\sim$ 1000 cal yr BP and therewith in a time interval equivalent to the MWP. A congruent short-lived warming that interrupts the general cooling trend is also observed in the temperature record inferred from the LDI.

In EZ4, the SWTs based on the LDI and TEX${}_{\mathrm{86}}^{\mathrm{L}}$ show a rapid cooling, with minimum temperatures of 14 ${}^{\circ }$C (TEX${}_{\mathrm{86}}^{\mathrm{L}}$) and 17 ${}^{\circ }$C (LDI) observed at $\sim$ 400 cal yr BP. Such a cooling trend is also evident in all pollen-based temperature reconstructions and likely corresponds to the Little Ice Age. In contrast to the pollen record, however, the long chain diol index and TEX${}_{\mathrm{86}}^{\mathrm{L}}$ both show a warming trend in the uppermost part of the studied sequence, with SWTs close to 20 ${}^{\circ }$C in the sample closest to the core top. Surface water temperatures of similar magnitude are presently observed during late July to mid-August in the Little Belt (ICES, 2017) and imply that both lipid paleothermometers reflect a late summer signal in our setting. For the TEX${}_{\mathrm{86}}^{\mathrm{L}}$, this finding is in agreement with previous observations made by Kabel et al. (2012) in the Baltic Proper. Likewise, the LDI has been reported to reflect summer surface water temperatures in a variety of oceanographic regions (Lopes dos Santos et al., 2013; Jonas et al., 2017), although it has recently been advocated to reflect an annual water temperature signal in the Ångermanälven River estuary (Warnock et al., 2017).

Reconstructed absolute temperatures vary significantly between the different inorganic- and organic-based proxies (Fig. 4). This is – to some extent – explained by the different realms the proxies reflect, but even between the LDI and TEX${}_{\mathrm{86}}^{\mathrm{L}}$, which are both considered to reflect summer SWTs in the Baltic Sea, significant deviations are observed. Such discrepancies have previously been reported for proxies from both subpolar (Rodrigo-Gámiz et al., 2015) and tropical regions (Jonas et al., 2017) and were explained by differences in either the habitat depth and/or seasonality of the biological sources of both proxies. Long chain diols used to calculate the LDI are synthesized by eustigmatophyte algae (Volkman et al., 1992), which are oxygenic photoautotrophs, and therefore this organic temperature proxy should exclusively record SWTs. Indeed, SWT based on the LDI in the sample closest to the core top is 18.3 ${}^{\circ }$C and therewith matches the summer SWT of $\mathrm{17.6}\pm \mathrm{0.7}$ ${}^{\circ }$C observed in the Little Belt region (ICES, 2017). This finding is in accordance with the observation that the LDI reflects the late summer–autumn water temperature based on a comparison of long chain diol distribution patterns in globally distributed marine surface sediments and monthly satellite SST data (Rampen et al., 2012). In addition to the marine realm, long chain diols are also ubiquitously found in freshwater environments. Rampen et al. (2014) noted a linear correlation between the distribution of long chain diols in a set of globally distributed lake sediments, annual mean, and summer surface water temperatures, possibly allowing the reconstruction of past continental temperature changes from lacustrine systems. However, the riverine transport of such freshwater-derived long chain diols to marine environments may confound the robust application of the LDI in coastal systems as long chain diol distribution patterns in lacustrine and marine settings show different responses to temperature. Indeed, substantial offsets between SWTs based on the LDI and satellite-derived SWTs have been observed in front of the Tagus and Sado river mouths (Iberian Atlantic margin) with the LDI yielding temperatures that are up to 8 ${}^{\circ }$C lower than measured average or summer SWTs (de Bar et al., 2016). In our setting, the branched isoprenoid tetraether index, a relative measure for the loading of terrestrial organic matter to aquatic environments (Hopmans et al., 2004), is characterized by consistently low values in brackish–marine EZ2–4 that contradict an increased freshwater discharge that may have confounded the use of the LDI (see below).

Sedimentological units, lithology (Andrén et al., 2015a), and position of ${}^{\mathrm{14}}$C ages (Van Helmond et al., 2017) vs. depth and IODP Site M0059 temperature proxies and related datasets plotted vs. age (extrapolated below $\sim$ 48.64 m c.d.). (a) Environmental zones (EZs) based on diatom data. (b) Pollen-inferred climate reconstructions based on the modern analogues technique, with the shaded area reflecting the range of the best eight analogues. Distmin1 and Distmin2 are the squared chord distance between fossil sample and the best to 8th-best modern analogues, indicating similarity between fossil and modern samples. (c) Methane and branched isoprenoid tetraether (BIT) indices, TEX${}_{\mathrm{86}}^{\mathrm{H}}$-based subsurface water temperature, TEX${}_{\mathrm{86}}^{\mathrm{H}}$-based surface water temperature, TEX${}_{\mathrm{86}}^{\mathrm{L}}$-based surface water temperature from the Gotland Deep (Kabel et al., 2012) for comparison, TEX${}_{\mathrm{86}}^{\mathrm{L}}$-based surface water temperatures shown with lake calibrations for EZ1 (ALST: annual, SLST: summer, WLST: winter lake surface temperatures), and LDI-based sea surface temperatures with lake calibrations for EZ1. (d) Climate phases and chronozones.

[Figure omitted. See PDF]

Temperatures reconstructed using the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ vary from $\sim$ 12 to $\sim$ 26 ${}^{\circ }$C (average of $\mathrm{17.0}\pm \mathrm{3.5}$ ${}^{\circ }$C) with three distinct maxima observed at the base of EZ2 ($\sim$ 7500 to $\sim$ 7300 cal yr BP; $\mathrm{19.6}\pm \mathrm{0.1}$ ${}^{\circ }$C) as well as in the lowermost part of EZ4 ($\sim$ 900 cal yr BP; 25.9 ${}^{\circ }$C) and close to the core top (16 cal yr BP; 21.1 ${}^{\circ }$C). The latter TEX${}_{\mathrm{86}}^{\mathrm{L}}$-based water temperature is close to SWTs observed in the Little Belt region during late July to mid-August (ICES, 2017). This suggests that the TEX${}_{\mathrm{86}}^{\mathrm{L}}$, similar to the LDI, may reflect summer surface water temperatures in the Baltic Sea as previously also argued by Kabel et al. (2012), who observed the best correlation for isoprenoid GDGT distribution patterns in surface sediments of the central Baltic with average water temperatures from July to October. The TEX${}_{\mathrm{86}}^{\mathrm{H}}$ has also been shown to correlate best with summer surface water temperatures in the Baltic Sea but correlations between this proxy and water temperatures are less strong as compared to those derived from the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ (Kabel et al., 2012). However, the TEX${}_{\mathrm{86}}^{\mathrm{H}}$-reconstructed temperature for the sediment closest to the core top is 18.9 ${}^{\circ }$C and therewith it is closer to the average summer SWT of $\mathrm{17.6}\pm \mathrm{0.7}$ ${}^{\circ }$C measured in the modern Little Belt. In addition, absolute temperatures inferred from the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ are generally lower, ranging from 12.6 to 20.6 ${}^{\circ }$C, and thereby show a significantly smaller amplitude compared to those derived from the TEX${}_{\mathrm{86}}^{\mathrm{L}}$. The calculation of this proxy, however, requires the robust identification of the crenarchaeol regioisomer that in marine environments is ubiquitously distributed but in brackish and freshwater systems is generally less abundant (e.g. Schouten et al., 2013b). In the Little Belt, it could be detected with only little confidence in most samples, resulting in comparatively high scatter with several abrupt temperature shifts by 5 ${}^{\circ }$C or higher in TEX${}_{\mathrm{86}}^{\mathrm{H}}$-reconstructed SWTs. Therefore, the reconstruction of water temperatures based on the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ seems to be less robust compared to the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ in the Baltic Sea. Molecular analysis evidences that thaumarchaeotal abundance in the central Baltic Sea is highest at greater water depth associated with increased rates of nitrification in oxygen-deficient bottom waters (Berg et al., 2015). As isoprenoid GDGTs may thus also derive from deeper parts of the water column, we also calculated subsurface water temperatures using the calibration of Kim et al. (2012). Reconstructed temperatures varied from 9.0 to 16.3 ${}^{\circ }$C (average $\mathrm{12.7}\pm \mathrm{1.7}$ ${}^{\circ }$C) and therewith significantly exceeded average bottom water temperatures ($\mathrm{7.6}\pm \mathrm{0.6}$ ${}^{\circ }$C) as well as summer bottom water temperatures ($\mathrm{9.3}\pm \mathrm{1.3}$ ${}^{\circ }$C) in the Little Belt region (Fig. 5). It thus seems unlikely that deep-water-dwelling Thaumarchaeota constitute a major source of isoprenoid GDGTs at our site. Hence, in contrast to some recent studies that have argued for the TEX${}_{\mathrm{86}}$ and its derivatives to specifically record water temperatures at the depth of the thermocline (Lopes dos Santos et al., 2013; Jonas et al., 2017), the TEX${}_{\mathrm{86}}^{\mathrm{L}}$, at least in the western Baltic Sea, may reflect a surface water temperature signal. Other factors potentially complicating the reconstruction of TEX${}_{\mathrm{86}}^{\mathrm{L}}$-based SWTs are the input of large quantities of soil-derived isoprenoid GDGTs (Weijers et al., 2006; Schouten et al., 2013b) or admixtures of high abundances of GDGT derived from sedimentary Euryarchaeota involved in the anaerobic oxidation of methane (AOM; Weijers et al., 2011). A larger contribution of isoprenoid GDGTs synthesized by methanotrophic archaea can be excluded for our site as the methane index, a quantitative measure for the relative contribution of GDGTs derived from AOM-related archaea to the total GDGT pool (Zhang et al., 2011), generally does not exceed a value of 0.2 and thus argues against a strong impact of methanotrophic archaea in the brackish–marine sediments of EZ2–4 (Fig. 5). In order to determine the input of terrestrial organic matter to the Little Belt region, we employed the branched isoprenoid tetraether index. Low values of this proxy ($<$ 0.1) indicate open marine conditions, while high values are more representative for terrestrial settings ($<$ 0.9; Hopmans et al., 2004). Empirical studies suggest that the TEX${}_{\mathrm{86}}$ and its derivatives do not allow a reliable reconstruction of SWTs if branched isoprenoid tetraether index values exceed a threshold of $\sim$ 0.3 (Weijers et al., 2006; Zhu et al., 2011). Branched isoprenoid tetraether index values in the Little Belt have an average of $\mathrm{0.29}\pm \mathrm{0.04}$ in EZ2–4, close to but still below this threshold. The TEX${}_{\mathrm{86}}^{\mathrm{L}}$ does, therefore, not seem to be biased by large allochthonous contributions of isoprenoid GDGTs. This is also supported by an absence of correlation between the branched isoprenoid tetraether index and TEX${}_{\mathrm{86}}^{\mathrm{L}}$ values ($r^{\mathrm{2}}=\mathrm{0.04}$; $p>$ 0.001), which is only to be expected if the terrestrial influence on the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ is negligible (Schouten et al., 2013a). In fact, a TEX${}_{\mathrm{86}}^{\mathrm{L}}$ temperature record with similar trends, covering the last 1000 cal yr BP with temperature maxima during the MWP and the modern hypoxic period, has been reported from the Gotland Basin in the central Baltic Proper (Kabel et al., 2012), suggesting that the TEX${}_{\mathrm{86}}^{\mathrm{L}}$ indeed traces climate-driven variations in SWT in the brackish–marine Littorina Sea stage.

Temperatures inferred from carbonate-clumped isotopes and Mg $/$ Ca ratios of benthic foraminifera are largely similar to those observed during late summer in the Little Belt. Indeed, in temperate climate settings molluscs and benthic foraminifera precipitate carbonate only if environmental conditions permit; thus, shell geochemical proxies (such as carbonate clumped isotopes and Mg $/$ Ca) tend to be skewed towards the seasons when growth is favourable (Filipsson et al., 2004; Austin et al., 2006; Schöne, 2008; Skirbekk et al., 2016). Despite temperatures, which are within range of present-day variability, the cleaning experiments on the foraminiferal calcite show that the 8 ${}^{\circ }$C decrease in bottom water temperatures may not be realistic after all. It cannot be excluded that the foraminiferal tests were still contaminated by the presence of authigenic carbonates like rhodochrosite, which can also contain Mg, and which cannot be removed by the current cleaning methodologies. The occurrence of such authigenic carbonate deposits in the Baltic Sea is common (Huckriede and Meischner, 1996; Andrén et al., 2015b; Hardisty et al., 2016; Van Helmond et al., 2017) and likely also affects foraminiferal tests (Groeneveld and Filipsson, 2013). Remaining contamination is also suggested by the significant correlation coefficient between Mg $/$ Ca and Mn $/$ Ca ($r^{\mathrm{2}}=\mathrm{0.72}$) after the oxidation-reduction cleaning. This suggests that especially the older part of the Mg $/$ Ca record, in which both Mg $/$ Ca and Mn $/$ Ca display the highest ratios, is contaminated.

The presence of benthic foraminifera in EZ2–4 and their absence in EZ1 is related to salinity, i.e. foraminifera did not occur in the freshwater setting of EZ1. The absence of foraminifera in the upper metres of the record ($<$ 5 m c.d.), however, is more likely due to poor preservation of the foraminiferal calcite. Similarly, the foraminiferal abundance decreased in EZ2, such that enough specimens were present to perform stable isotope analyses, but there were too few to perform reliable Mg $/$ Ca analyses. Poor preservation may have been caused by the degradation of organic matter in the sediment, which led to lower pH values of the pore water such that dissolution of the foraminiferal tests took place.

Conclusions

In this study, we have investigated and compared how inorganic- and organic-based temperature as well as salinity proxies perform in the coastal setting of the Little Belt with highly variable environmental conditions. Sediments deposited since $\sim$ 7800 cal yr BP at IODP Site M0059 were analysed to capture the transition from freshwater to marine conditions and thus to study the effect of facies change on proxy applications.

Our study demonstrates that a multi-proxy approach allows the decipherment of the various factors that have influenced past oceanographic conditions in the Little Belt. Between $\sim$ 7800 and 7400 cal yr BP, a direct connection of the Little Belt to the Baltic Proper may have existed as salinity-specific proxies indicate lacustrine conditions, while a connection to the Kattegat can be excluded. Salinities and probably also water temperatures increased within 200 years after the onset of the Littorina Sea stage at $\sim$ 7400 cal yr BP. An interval of particularly saline conditions and high summer water and air temperatures followed, although both salinity and water temperatures declined with increasing precipitation until $\sim$ 4000 cal yr BP. After $\sim$ 4000 cal yr BP, the Little Belt witnessed decreasing temperatures during the warm season both in the marine and terrestrial realm, while salinity increased.

In addition, our study highlights the importance and value of a multi-proxy approach to reconstruct past oceanographic conditions. The different salinity proxies used here show generally similar trends in relative changes in salinity, but do not allow quantitative estimates of salinity (except for marine ostracods). In contrast, the reconstruction of temperatures is associated with particularly large uncertainties and variations in absolute values by up to 8 ${}^{\circ }$C for bottom waters and even up to 16 ${}^{\circ }$C for summer surface waters. For example, different cleaning techniques for Mg $/$ Ca in the foraminifera show different results, which partly correlate with indicators for contamination, especially in the deeper intervals studied. This suggests that the Mg $/$ Ca ratio of those samples is likely overestimated, so that the decreasing trend in Mg $/$ Ca ratios and thus calculated bottom water temperatures may not be as large. The differences in results based on the lipid proxies (long chain diol index and TEX${}_{\mathrm{86}}^{\mathrm{L}}$) can partly be explained with the application of modern-day proxy calibrations in areas that experienced significant changes in depositional settings (e.g. change from freshwater to marine conditions). Our study shows that particular caution has to be taken when applying and interpreting proxies in coastal environments, where water mass conditions can experience more rapid and larger changes than in open ocean settings. Approaches using a multitude of independent proxies may thus allow for more robust paleoenvironmental assessments.

Data not enclosed in the supplement will be uploaded to the PANGEA database (LINK). All ostracod data are accessible via PANGEA (10.1594/PANGAEA.873270, Stepanova, 2017).

The Supplement related to this article is available online at https://doi.org/10.5194/bg-14-5607-2017-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

This investigation used samples from IODP Expedition 347, and we are very grateful to IODP and ESSAC for funding and organizing the expedition. Stephen Obrochta, Sofia Ribeiro, Ingo Feeser, and Walter Doerfler are gratefully acknowledged for insightful discussions and comments on the age model, dinoflagellate cyst taxonomy and ecology, and the Lake Belau pollen record. We thank the crew and scientific team of IODP Expedition 347 as well as the technical staff at MARUM, Bremen, for their support during the expedition's off- and onshore phases. We are grateful for technical support from the different laboratories involved in the study. We thank Audrey Limoges, Filipa Naughton, and the anonymous referee for their detailed input to the discussion paper, and we thank Markus Kienast for editing. This research was supported by the German Research Foundation (grants KO3944/6-1 to Ulrich Kotthoff and BA3841/5-1 to Thorsten Bauersachs), the Carlsberg Foundation and Geocenter Denmark (grants to MSS and ORC), and US National Science Foundation (NSF grant no. T347A13). Thomas Andrén, Jeroen Groeneveld, Nadine Quintana Krupinski, Helena L. Filipsson thank the Swedish research councils VR and FORMAS, the LAM Foundation, the Crafoord Foundation and the Center for Environmental and Climate Research at Lund University for funding. Thomas Andrén and Elinor Andrén acknowledge funding from The Foundation for Baltic and East European Studies (grants 1562/3.1.1/2013 and 2207/3.1.1/2014). Jeanine L. Ash was supported by the National Science Foundation (NSF-DGE-1144087 and OCE-0652315). Caroline P. Slomp and Niels A. G. M. Van Helmond acknowledge funding from the European Research Council (ERC Starting Grant 278364) and the Netherlands Organisation for Scientific Research (NWO-Vici 865.13.005). Edited by: Markus Kienast Reviewed by: Audrey Limoges, Filipa Naughton, and an anonymous reviewer