A high-resolution perspective on climate drivers of lake stratification and phototrophic community dynamics in Late Glacial Central Europe

Predicting the trajectory of aquatic deoxygenation under global warming requires a mechanistic understanding of lacustrine responses to rapid climate shifts. We investigated how climate-driven changes in catchment vegetation and local iron-rich lithology regulated lake stratification and ecosystem resilience in the maar lake Holzmaar (Central Europe). We focused on the Late Glacial, specifically on transitions during Dansgaard-Oeschger Event 1 (DOE-1; ca. 14 690–11 700 cal yr BP), a period of rapid natural warming and cooling that serves as an analogue for future high amplitude climate variation and for modern Arctic lakes undergoing rapid climate-driven transitions. Combining non-destructive hyperspectral imaging (HSI) of sedimentary pigments with high-resolution XRF geochemistry, we resolved parts of the ecosystem trajectory during DOE-1.

The primary producer community shifted from an oligotrophic cyanobacterial and low-light Pleniglacial assemblage to a stable, stratified Allerød assemblage, characterized by the planktonic diatom Stephanodiscus minutulus and anoxygenic purple sulphur bacteria (PSB) in the photic zone. While regional warming (mean summer temperature increased ∼ 2.8 °C) provided the physical potential for lake stratification, our data suggest that intense anoxia was primarily triggered by the expansion of Betula in the watershed. This afforestation stabilized the water column through wind shielding. The termination of the anoxic phase coincided with the onset of the Younger Dryas cooling and increased aridity, which effectively destabilized the existing stratification. While the shift from Betula to Pinus forest may have caused a change in the terrestrial-aquatic linkage, the primary driver of the transition was the physical forcing (lake mixing) of the climatic shift (cooling).

Geochemically, the lake exhibited remarkable resilience. Unlike carbonate-dominated systems prone to internal phosphorus loading, Holzmaar efficiently sequesters nutrients via a dual mechanism of reactive iron binding (authigenic vivianite) and stable mineral burial. The phosphorus trap prevents nutrient release by permanently sequestering P in the sediment, allowing rapid ecosystem recovery without delay once the specific climate and vegetation drivers shift. Our findings demonstrate that in volcanic maar lakes, catchment vegetation characteristics and local lithology can modulate, and even override, the direct effects of climate warming on aquatic anoxia.

1 Introduction

Hypoxia and anoxia in freshwater ecosystems represent a critical global environmental challenge. The depletion of dissolved oxygen in bottom waters threatens aquatic biodiversity, alters nutrient cycling, and compromises water quality (Jenny et al., 2016; Woolway and Merchant, 2019). While often driven by anthropogenic nutrient loading in the Anthropocene, anoxia is also a natural feature of many deep lake systems, fundamentally linked to physical mixing regimes and climate forcing. However, predicting the trajectory of these systems under future warming remains difficult because separating the direct effects of temperature from the indirect effects of catchment change (land cover, nutrients) is challenging in modern, human-impacted landscapes.

The development of hypolimnetic anoxia typically follows a predictable cascade of physical and biogeochemical processes. Rising atmospheric temperatures increase the thermal stability of the water column, extending the duration of seasonal stratification and isolating deep waters from atmospheric oxygen exchange (Kienel et al., 2017). Simultaneously, warmer surface waters and enhanced nutrient fluxes stimulate primary production, increasing the export of organic matter (OM) to the hypolimnion. The microbial decomposition of this OM consumes the limited oxygen pool, eventually driving the system toward anoxia once physical re-aeration is inhibited. This coupling has been well-documented in high-resolution studies of varved sediments across Central Europe. For instance, in Lake Żabińskie and Lake Jaczno (NE Poland), seasonal anoxia is tightly coupled to productivity cycles and physical mixing thresholds (Zander et al., 2021a; Żarczyński et al., 2019; Butz et al., 2017). Similarly, studies from Swiss lakes have demonstrated how changes in mixing regimes directly influence sediment phosphorus retention (Tu et al., 2020). In many modern temperate lakes, this cycle is seasonal: anoxia builds up during the stratified summer and is disrupted when surface cooling and wind stress induce deep convective mixing in late autumn/winter or spring (Schwefel et al., 2016).

The Late Glacial period (Dansgaard-Oeschger event 1: Greenland Interstadial 1, GI-1, ca. 14 690–12 90 cal yr BP; Greenland Stadial 1, GS-1, 12 900–11 700 cal yr BP) offers a unique natural laboratory to test these mechanisms over centennial timescales. We hypothesize that Late Glacial lakes in Central Europe serve as valuable analogues for modern high-latitude lakes, oscillating between distinct thermal states as climate shifts alter aquatic dynamics. During the rapid Bølling-Allerød warming, lakes entered a long-term phase of stable stratification and anoxia, driven not just by rising temperatures, but also by the expansion of catchment forests that shielded the water column from wind mixing (Makri et al., 2020; Tu et al., 2021). Conversely, the Younger Dryas cooling mimics a transition back to a polar mixing regime. Despite a shift in catchment vegetation from Betula to Pinus, decreased mean annual temperatures and humidity, along with a shorter growing season, likely reduced thermal gradients and promoted more frequent lake mixing to reoxygenate the hypolimnion. However, the resilience of these lake ecosystems to such shifts, specifically, whether they recover immediately or exhibit a delayed response, depends heavily on internal biogeochemical feedbacks (Schmidt et al., 2002; Lone and Balakrishna, 2023; Schouten et al., 2026).

A key uncertainty lies in how catchment characteristics (lithology, topography, land cover and soils) modulates this response. While physical thresholds for anoxia have been established, such as the ca 80 % arboreal pollen threshold for wind shielding at Soppensee, Moossee and Żabińskie (Makri et al., 2020; Zander et al., 2021b; Schouten et al., 2026), a critical gap remains regarding the biological and geochemical response. Previous studies at Holzmaar have shown that during Dansgaard-Oeschger Event 1 (ca. 14 690–11 700 cal yr BP), a rapid climate-driven vegetation succession occurred from steppe tundra to shrub tundra, and finally to Betula and Pinus forests during the Allerød (Litt and Stebich, 1999). Coincident with this catchment stabilization, the lake transitioned from ultra-oligotrophic to mesotrophic conditions, marked by the replacement of benthic Staurosira construens with planktonic Stephanodiscus minutulus and an increased accumulation of organic carbon, primary production (based on δ¹³C) and biogenic silica (Zolitschka et al., 2000; García et al., 2022; Lücke et al., 2003). Although the subsequent Younger Dryas cooling (ca. 12 750–11 500 cal yr BP) temporarily opened the forest canopy and reduced aquatic production (Brauer et al., 2001), the high-resolution dynamics of these primary producer communities, specific algal turnovers, and discrete elemental nutrient cycles remain completely unresolved.

Here, we investigate Holzmaar (Westeifel Volcanic Field, Germany) to resolve the coupling between climate forcing, catchment vegetation, and the internal drivers of eutrophication, stratification and aquatic anoxia. The Eifel region is uniquely suited for this, as recent quantitative temperature reconstructions using branched GDGTs provide a specific local thermal history (Zander et al., 2024), allowing us to discriminate local climate forcing from hemispheric trends (NGRIP). We combine high-resolution hyperspectral imaging (HSI) and HPLC-derived pigment concentrations with XRF geochemistry and sequential extractions of P, Mn, and Fe. These biogeochemical data are further supported by pollen, diatom, and temperature records to address the following key questions:

1. How do phototrophic communities change in response to periods of rapid climate shifts (both warming and cooling) during Dansgaard-Oeschger Event 1 (14 690–11 700 yr cal BP)?

2. How do internal biogeochemical mechanisms (such as iron-mediated nutrient retention) and catchment vegetation modulate the resilience of the lake and ecosystem reversibility when these climatic drivers shift?

By comparing the trajectory of Holzmaar to regional records (Tu et al., 2021; Makri et al., 2020), we hypothesize the specific roles of catchment vegetation and iron availability in regulating lake resilience.

2 Material and methods

2.1 Study site

Holzmaar (50°7^′ N, 6°53^′ E; 425 m a.s.l.) is a maar lake located in the Westeifel Volcanic Field, Germany (Fig. 1A, Zolitschka et al., 2000). The catchment is characterized by iron-rich Devonian bedrock and known for its iron deposits associated with volcanism (Schmincke, 2007; Lottermoser et al., 1997; Schmitt et al., 2023). The lake basin exhibits a circular morphology with a small surface area of 58 000 m² (5.8 ha) relative to its maximum depth of 20 m (Negendank and Zolitschka, 1993). Under modern conditions, due to this great depth and steep basin slopes (Fig. 1B), deep-water mixing is typically inhibited below 10–15 m (Brauer et al., 2001). While this specific crater morphometry naturally predisposes the lake to meromixis and a high sensitivity to climatic forcing (cf. Lami et al., 1994), Holzmaar operates primarily as a monomictic meso- to eutrophic lake, or dimictic in years with winter ice cover (Lücke, 1998; Moschen et al., 2009). It experiences seasonal thermal stratification from spring (April/May) to late autumn or winter (October/December) (Lücke, 1998; Messyasz et al., 2012), with a steep thermocline reaching a maximum depth of approximately 6 m (Raubitschek et al., 1999). The resulting physical isolation of the hypolimnion drives oxygen depletion and solute accumulation, establishing a distinct metalimnetic chemocline roughly between 6 and 8 m depth (Messyasz et al., 2012). This seasonal hypolimnetic anoxia ultimately favors the preservation of the annually laminated (varved) sediments utilized in this study (Zolitschka, 1998; García et al., 2022).

Figure 1 (A) Regional topographic map (derived from WorldClim v2.1, Fick and Hijmans, 2017) showing the study area in western Germany, highlighting the location of Holzmaar (red triangle). The inset map displays the broader location within Germany, with administrative boundaries sourced from Natural Earth (https://www.naturalearthdata.com/, last access: June 2026). (B) Bathymetric map of Holzmaar (modified after Zolitschka, 1998). The red markers denote the locations of the four adjacent sediment cores (HZM19-7, HZM19-8, HZM19-10, and HZM19-11) retrieved from the deep central basin in 2019.

2.2 Coring, physical and chemical properties

Four adjacent sediment cores (Fig. 1B, HZM19-07, HZM19-08, HZM19-10, HZM19-11) were recovered in August 2019 from the centre of Holzmaar at 18–19 m water depth using a UWITEC piston corer (García et al., 2022) and compiled to a composite sediment record (Birlo et al., 2023). The X-ray fluorescence (XRF) core scanning for elemental composition with a Cr X-ray tube, high-resolution logging of magnetic susceptibility, and determination of total organic carbon (TOC), total nitrogen (TN), and biogenic silica (BSi) data at a relatively low resolution of 16 cm (García et al., 2022) were complemented by hyperspectral scanning data acquired at the University of Bern and higher resolution diatom stratigraphy (every 4 cm) established at the University of Bremen.

Sedimentary pigments were quantified non-destructively using hyperspectral imaging with a Specim PFD-CL-65-V10E line scan camera operating in the visible and near-infrared range (VNIR, 400–1000 nm). The camera set-up (frame rate: 8 Hz, exposure: 120 ms) delivered data with a spectral resolution of 1.56 nm and a spatial resolution of 80 µm (pixel size). The hyperspectral data were normalized using a white reference (BaSO₄) and dark reference (closed camera aperture).

Hyperspectral data post processing was performed using napari-sediment (Witz and michaelh00, 2025), an open-source Python plugin for the napari image viewer, specifically designed for interactive analysis of hyperspectral sediment core imagery. napari-sediment provides an integrated workflow that combines visualization, preprocessing, dimensionality reduction, and spectral analysis in a single graphical user interface following the data processing proposed in Butz et al. (2015) and further described in the Appendix A. Relative Absorption Band Depth (RABD) indices were calculated to quantify absorption troughs associated with specific photopigments. Following established protocols by Butz et al. (2015), Schneider et al. (2018) and Zander et al. (2022), absorption trough minima were identified at the following wavelengths:

• 619 nm – Phycocyanin from cyanobacteria (Wienhues et al., 2025)

• 670 nm – Total chlorophyll a (Tchl) from oxygenic phototrophs (Rein and Sirocko, 2002)

• 715 nm – bacteriopheophytin e (bphe e) from green sulphur bacteria (GSB) (Zander et al., 2023)

• 845 nm – bacteriopheophytin a (bphe a) from purple sulphur bacteria (PSB) (Butz et al., 2015)

RABD index values were calculated using the formula of Butz et al. (2015):

$$\begin{array}{}\text{(1)}& {\text{RABD}}_{\mathit{\lambda }}={\displaystyle \frac{\frac{X_{\text{right}}\cdot R_{\text{left}}+X_{\text{left}}\cdot R_{\text{right}}}{X_{\text{right}}+X_{\text{left}}}}{R_{\mathit{\lambda }}}}\end{array}$$

where R_λ is the reflectance at the absorption minimum (λ), X stands for the number of bands from the left side of the trough to λ (X_left), or from the right side to the λ (X_right), and R_left and R_right are the reflectance values at start and end of the reflectance trough.

RABD index values from hyperspectral imaging were calibrated to absolute pigment concentrations measured by spectrophotometry using linear regression models. Calibration samples were selected to represent the full range of pigment concentrations present in the core following a normal distribution.

2.3 Chronology

The chronology for the Holzmaar composite profile is based on the Bayesian age-depth model established by Birlo et al. (2023). This model integrates the original high-resolution varve chronology of Zolitschka et al. (2000) with radiometric dating to optimize the age transfer to the new HZM19 sediment cores. The age model is stratigraphically anchored by the Laacher See Tephra (LST) (Reinig et al., 2021), a distinct isochron visible in the lithology (Cluster IV, HClust, dark red in Fig. 2), which serves as a critical time marker for the Allerød. This varve-based framework provides annual to sub-decadal temporal resolution, allowing for precise correlation of biogeochemical changes during the Late Glacial period (ca. 14 690–11 700 cal yr BP).

The biostratigraphy for the Holzmaar record follows the regional pollen stratigraphy (Lang et al., 2023; Litt and Stebich, 1999). In this study, we adopt the term Pleniglacial (Birlo et al., 2023; García et al., 2022) to describe the period corresponding to the Oldest Dryas in Swiss and South German biostratigraphies, e.g., Lotter and Kienast (1992). Similarly, the period referred elsewhere as the Meiendorf is here termed the Bølling. The subsequent sequence follows the regional standard: the Older Dryas stadial, the Allerød interstadial with the Laacher See Tephra, the Younger Dryas stadial. The exact dates of their stratigraphic boundaries are from Birlo et al. (2023).

2.4 Sedimentary pigment quantification

To combine the advantages of rapid, non-destructive, and ultra-high-resolution hyperspectral imaging (HSI) with the precise, compound-specific quantification of traditional wet-chemistry methods, we calibrated the continuous HSI indices against discrete pigment concentrations. These independent pigment measurements were performed using UV-visible absorption spectrophotometry (Shimadzu UV-1800) and high-performance liquid chromatography (HPLC). Pigment extractions were performed on 0.5 to 1.5 g wet and homogenized sediment subsamples using 100 % HPLC-grade acetone following the protocol of Lami et al. (1994); Pniewski (2020); Strickland and Parsons (1972). In short, we extracted pigments by adding 5 to 13 mL acetone in several steps (2–3 mL per step), keeping the samples in dark and at 4 °C conditions over night, or 30 min in −4 °C between the extraction steps, until the supernatant was visually colourless. The samples were always centrifuged (3500 rpm, 10 min), the supernatant was decanted and filtered through hydrophobic PTFE 0.22 µm filters (13 mm syringe filters) and stored at −20 °C in darkness until analysis.

For spectrophotometric analysis, pigment extracts were measured using a UV-visible spectrophotometer (wavelength range 350–900 nm, spectral resolution 0.1 nm) with 1 cm path length PP cuvettes. Total chlorophyll a and derivatives (pheophytins, pheophorbides) were quantified spectrophotometrically using absorption at 663 and 665 nm and translated to concentration using the Lambert-Beer law and the extinction coefficient for bulk green pigments (chlorophylls, pheophytins and pheophorbides) in 90:10 Acetone:Water solution ($\mathit{ϵ}=\mathrm{80.8}\cdot {\mathrm{10}}^{-\mathrm{3}}$ L mg⁻¹ cm⁻¹, adjusted after Jeffrey and Humphrey, 1975). Total bacteriopheophytin a concentration was quantified using the absorption peak at 745 nm in 90 % acetone (±10 % of water originating from the wet sediment), and applying the specific absorption coefficient ($\mathit{ϵ}=\mathrm{52.9}\cdot {\mathrm{10}}^{-\mathrm{3}}$ L mg⁻¹ cm⁻¹,  Fiedor et al., 2002).

Pigment extracts (100 µL aliquots) were analysed using a Dionex Ultimate 3000 series reverse-phase HPLC system equipped with a C18-ODS column (Agilent Omnisphere-5) and a Diode Array Detector (DAD-3000RS), following the method of Lami et al. (1994) described in detail in Schouten et al. (2025). Chromatogram peaks were integrated using Chromeleon 7.2 (ThermoFischer Scientific^®) with a Gaussian fit approach. Individual pigments were identified based on retention times and absorption spectra, and quantified using pre-determined linear regression coefficients for authentic standards (Schouten et al., 2025).

To calibrate hyperspectral imaging (HSI) data with pigment concentrations measured by spectrophotometry, we employed a proxy-proxy linear regression model. After testing for normality using Shapiro-Wilk tests and screening for outliers using Rosner's test, we established linear relationships between RABD indices and pigment concentrations. For total chlorophyll (TChl), we fitted the model:

$$\begin{array}{}\text{(2)}& \mathrm{TChl}\left(\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\right)=\mathrm{467.6}\times {\mathrm{RABD}}_{\mathrm{670}}-\mathrm{467.1}\end{array}$$

with R²=0.92, $p=\mathrm{4.6}\cdot {\mathrm{10}}^{-\mathrm{19}}$, RMSEP = 3.25 µg g⁻¹, (Fig. A1).

For bacteriopheophytin a (Bphe a), we fitted:

$$\begin{array}{}\text{(3)}& \mathrm{Bphe}a\left(\mathrm{\mu }\mathrm{g}{\mathrm{g}}^{-\mathrm{1}}\right)=\mathrm{131.3}\times {\mathrm{RABD}}_{\mathrm{845}}-\mathrm{129.2}\end{array}$$

with R²=0.82 $p=\mathrm{4.6}\cdot {\mathrm{10}}^{-\mathrm{13}}$, RMSEP = 0.97 µg g⁻¹ (Fig. A2). These calibration equations were then applied to the entire HSI dataset to generate spatially-resolved pigment concentration maps.

2.5 Diatoms

Diatom samples were prepared from 50 freeze-dried subsamples collected at discrete 4 cm intervals along the composite sediment profile. An aliquot of each sample was oxidized with 30 % H₂O₂ and heated in a water bath at approximately 80 °C for 2–5 min to eliminate organic matter, following standard methods (Battarbee and Berglund, 1986). Permanent slides were mounted using Naphrax mounting medium. A minimum of 400 diatom valves per slide were counted using light microscopy (Olympus CX40 with Plan Ach 100X/1.25 immersion objective and Zeiss Axioplan with Plan Neofluar 100X/1.30 immersion objective equipped with differential interference contrast) to calculate relative abundances and absolute concentrations (valves per gram dry sediment) using a modified (square cover slides) evaporation tray method (Battarbee and Berglund, 1986). Diatom identification followed standard European floras, and diatom accumulation rates were calculated by multiplying total diatom abundance by sediment accumulation rates derived from the varve chronology.

2.6 Sequential extraction of P, Mn, and Fe

Sequential extraction was performed using a modified protocol from Lukkari et al. (2007) and Scholtysik et al. (2022) to differentiate sedimentary P, Fe, and Mn into chemically distinct pools. The extraction scheme defined four operational fractions: (1) F1 (labile/porewater), extracted with 0.46 M NH₄Cl; (2) F2 (redox-sensitive), extracted with 0.11 M Bicarbonate-Dithionite (BD); (3) F3 (Al-bound/organic/hydroxides), extracted with 1 M NaOH; and (4) F4 (calcium-bound/residual), extracted with 0.5 M HCl. The extracts were analysed using ICP-MS (Agilent Series 7600) in a 1 % HNO₃ matrix-matched solution at the University of Bern. The concentrations of Fe, Mn and P were calculated based on external calibration curves from IPC-MS multi-element standard.

We acknowledge that strict anaerobic integrity could not be maintained during core transport and subsampling. Consequently, labile reduced species (e.g., porewater Fe²⁺, authigenic vivianite) susceptible to rapid oxidation may have shifted from their target fractions (F1 or F3) into the reducible fraction (F2). Therefore, we adopt a conservative operational definition, where we interpret Fraction 2 not as specific oxides, but as the total reactive pool, encompassing both original oxides and secondary oxidation products of authigenic minerals. To infer the original mineralogical hosts of phosphorus, we evaluate the molar Fe:P stoichiometry of this reactive pool, rather than relying on the operational separation of F1 and F2 alone.

2.7 Pollen and external data

To find the link of our in-lake biogeochemical data with the external drivers, we compiled existing pollen data of Holzmaar from the literature and databases (Litt and Stebich, 1999; Litt, 2001; Litt et al., 2009; Litt, 2010, 2019). To synchronize these records with our 2019 chronology, individual core sections from the 1992 campaign (HZM92 1a–6u, 2b–5o, 2b–5u) were correlated to the 2019 composite depth using visual lithostratigraphic tie-points and then stacked with the 1996 composite record (see Appendix B for full methodological details). The Holzmaar high-resolution pollen record covers only the period from 15 000 to 13 250 cal yr BP, and thus the trends for Betula, Pinus, and non-arboreal pollen (NAP) from 13 250 to 11 750 cal yr BP were modelled using the high-resolution record from the nearby Meerfelder Maar (Litt and Stebich, 1999, 2001). Missing values were reconstructed by applying scaling factors derived from the linear relationship between the two lakes during their overlapping stratigraphy (code is available as a part of the workflow at GitHub/Renku (Zahajská, 2026).

Moreover, we generated the annual, July and January insolation curve for the location of Holzmaar, with the resolution of 50 years from −50 to 15 950 cal yr BP using the model of Laskar et al. (2004). To further constrain the temperature proxy, we used the chironomid-based July temperature reconstruction for Gerzensee (Lotter et al., 2012) as well as for Egelsee (Larocque-Tobler et al., 2009; Kaufman et al., 2020). We also utilized the regional GDGT-based temperature reconstruction from the Eifel, including Holzmaar (Zander et al., 2024). Lastly, we included the North GRIP δ¹⁸O (NGRIP Members, 2004) and ice-core dust records (Ruth et al., 2003; Gkinis et al., 2014) into our set of external variables. These datasets demonstrate a clear seasonal pattern in dust concentrations that correlate with climatic variations.

2.8 Predictor variable selection and quality control

The first step was to evaluate relationships among external environmental variables and reduce dimensionality to minimize collinearity and identify the most informative drivers. Initial pre-selection of predictor variables was based on Pearson correlation matrices and Principal Component Analysis (PCA) of z-scored data (Fig. B2). Prior to redundancy analysis (RDA), we performed Variance Inflation Factor (VIF) assessment to quantify multicollinearity among the remaining candidate variables using the threshold VIF > 10 (indicating critical multicollinearity, Blanchet et al., 2008). This VIF-based filtering removed seven variables with high collinearity (annual insolation, January insolation, July insolation, Holzmaar Betula, Pinus, non-arboreal pollen, and arboreal pollen), retaining nine variables with VIF ≤ 10 for the final RDA model. We added passively on top of the RDA three additional variables (Holzmaar Betula, Pinus and January insolation), which did not pass the VIF just for the purpose of easier interpretation.

As a quality control measure, three GDGT-derived temperature samples with known analytical issues (ages 13 128, 13 340, 13 875 yr b2k) were excluded from all analyses. The final subset of retained environmental variables comprised: GDGT-derived temperature, chironomid-inferred July temperature from Gerzensee, chironomid-inferred July temperature from Egelsee, δ¹⁸O_NGRIP, NGRIP dust data, XRF-derived Ca and Ti variation, and regional pollen abundances (Juniperus, Salix from Holzmaar). This refined set of non-redundant variables was used for RDA analysis with pigment data.

2.9 Data analysis

High-resolution XRF core scanning data were normalized using centred log-ratio (CLR) transformation to address the closed-sum constraint inherent to compositional data and minimize matrix effects (Weltje and Tjallingii, 2008; Bertrand et al., 2024; Aitchison, 1982). Hyperspectral imaging indices were calibrated against discrete spectrophotometer bulk pigment measurements using linear regression, as described above, enabling the quantification of sedimentary pigment concentrations at high spatial resolution throughout the core. The XRF and hyperspectral imaging (HSI) datasets, acquired at native sub-millimetre resolution (HSI: 80 µm pixel size, see Sect. 2.2), were co-registered and resampled to a common 1 mm depth resolution prior to further analysis, whereas pollen, diatom, and external climate proxy records were resolved at coarser, decadal-to-centennial resolution. To facilitate multi-proxy comparison, including the external datasets, all datasets were harmonized and interpolated to a common (100 years) temporal resolution using Gaussian kernel smoothing.

Hierarchical cluster analysis (Ward's D2 method with Euclidean distance) and principal component analysis (PCA) were applied to the CLR-transformed XRF data to identify geochemical facies, with optimal cluster numbers determined through silhouette analysis. HPLC-derived pigments were grouped based on Pearson correlation coefficients of their downcore concentration profiles, and taxonomic assignments to algal and bacterial producer groups (e.g., cyanobacteria, green algae, diatoms, cryptophytes, purple sulphur bacteria) were made following established pigment biomarker literature (see Table A1).

Redundancy analysis (RDA) was performed to quantify the relationship between pigment community composition and environmental drivers. The RDA model was constructed using scaled CLR-transformed pigment data as response variables and the VIF-selected environmental variables as predictors. Statistical significance of the overall RDA model was assessed using permutation tests (999 unrestricted permutations) at a significance level of α=0.050. Permutation testing was chosen as it is more robust for multivariate data and does not assume normality. An ecological variable grouping approach was used to ensure that major environmental drivers (vegetation proxies: Betula and Pinus; winter insolation) were represented in the analysis even if they were excluded by strict VIF filtering (VIF > 10). The full workflow is available on Renku/GitHub (Zahajská, 2026) and fully linked to the dataset at Zenodo (https://doi.org/10.5281/zenodo.18429717, Zahajská et al., 2026).

3 Results

3.1 Lithology

The unconstrained hierarchical clustering of XRF data resulted in five distinctive clusters (Figs. 2 and B3) separating more allochthonous sediment (Ti, K, Zr, Al) from autochthonous sediment (Ca, S, $\mathrm{inc}/\mathrm{coh}$). In combination with the hyperspectral data, we interpret the clusters as follows (Fig. 2): Cluster I (lime-green) corresponds to clayey silts with low pigment abundance at the bottom of the core, followed by cluster II inter-layered with cluster III (orange, dark blue) representing the varved section of diatoms ($\mathrm{Si}/\mathrm{Ti}$) and organic matter ($\mathrm{inc}/\mathrm{coh}$, RABD620, the total chlorophylls, bacteriopheophytin a). Cluster IV (dark red) is specific for the LST followed by organic rich section of cluster III (dark blue). The core top is classified as cluster V (light blue), which carries similarities with cluster I and is mainly minerogenic.

Figure 2 Selection of data from non-destructive high-resolution methods presenting the Holzmaar paleoenvironmental record plotted against depth (left axis) and time (non-linear right axis). The panels display geochemical proxies for detrital input (Ti clr) and carbonate (Ca clr), relative diatom abundance ($\mathrm{Si}/\mathrm{Ti}$), total diatom concentration in green bars, and organic matter ($\mathrm{inc}/\mathrm{coh}$). Statistical analyses include hierarchical clustering (HClust and CONISS) defining stratigraphic zones, and Principal Component Analysis. PC1 (45.3 %) indicating autochthonous vs. allochthonous input; PC2 (11.5 %) is reflecting organic (redox) conditions vs. minerogenic input. An intensity colour map of hyperspectral indices is shown alongside pigment proxies for cyanobacteria (RABD620), total chlorophylls, and purple sulphur bacteria (bacteriopheophytin a). Vertical lines are representing the limit of quantification. Major biozones following Birlo et al. (2023): Pleniglacial, Bølling, Allerød, and Younger Dryas are indicated on the right and marked by red and blue horizontal lines. The gray horizontal band highlights the position of the Laacher See Tephra (LST).

Similar separation is observed in the PCA (Fig. B3), where PC1 (45.3 %, Figs. 2 and B3) separates the XRF data between allochthonous (positive) and autochthonous (negative) sediment sources. PC2 (11.5 %) represents the division between minerogenic (positive) vs organic matter and redox sensitive (negative) elements in the sediment record.

Lastly, CONISS clustering was used to explore the clusters constrained with depth. Similarly to the unconstrained hierarchical clustering, five clusters were defined, which coincide well with the previously published lithology (Zolitschka, 1998).

3.2 Primary producer communities

The sedimentary pigment data were examined using a heatmap, where individual pigments were clustered based on their time-series similarities using correlation (Fig. 3). This clustering resulted in three distinct primary producer assemblages that closely track the stratigraphic zones: (1) the Pleniglacial assemblage, (2) the transitional Bølling assemblage, and (3) the stratified Allerød assemblage (Fig. 4). The statistical separation of these communities and their distinct stratigraphic evolution, including the return to mixed conditions during the Younger Dryas, is further visualized in the PCA trajectory (Fig. B5).

The Pleniglacial assemblage (pigment group 1, Figs. 4, 3 and B4) is characterized by high abundances of pigments associated with cyanobacteria (echinenone, canthaxanthin) and cryptophytes (alloxanthin). This early assemblage also includes green algae (α-carotene, lutein) alongside diatoms and chrysophytes. Notably, this assemblage features strong signals of physiological stress markers, including diatoxanthin (a light-stress indicator in diatoms) and astaxanthin. The astaxanthin indicates contributions from freshwater chlorophytes such as Haematococcus pluvialis or dinoflagellates (e.g., the Peridinium or Glenodinium). The presence of benthic species of diatoms, such as Staurosira construens, supports low lake-level, well-mixed (holomictic) conditions. Additionally, the early presence of the planktonic diatom Pantocsekiella ocellata aligns with the detection of specific pigment for green sulphur bacteria (isorenieratene, Fig. 3), suggesting low-light conditions and anoxygenic niches. Further, we do not observe any significant formation of the redox-sensitive fractions of Mn (Fig. 4), P, or Fe (Fig. B6), indicating generally oxygenated bottom waters. This very diverse aquatic primary producer community occupied the lake during the relatively colder Pleniglacial period and was gradually replaced in the early Bølling by deeper water column community.

Figure 3 Stratigraphic heatmap and hierarchical clustering of sedimentary pigment assemblages. The heatmap displays the relative abundance (z-score normalized) of identified pigments across the sediment profile. Rows represent individual sediment samples (with ages in cal yr BP on the right), colour-coded by stratigraphic zone (left bar). Columns represent individual pigments, categorized based on their producer in the colour-coded top bar. The correlation-based clustering (top dendrogram) reveals three distinct ecological clusters further used as pigment groups. The left cluster contains pigments associated with the pigment group 1 (e.g., echinenone, canthaxanthin, alloxanthin and isorenieratene), which show high relative abundances in the Pleniglacial (bottom rows). The middle cluster contains pigments defining the pigment group 3 (okenone and fucoxanthin), which exhibit peak abundances during the Allerød (middle-upper rows). The right cluster represents pigment group 2 including green algae, cyanobacteria and green sulphur and purple non-sulphur bacteria, peaking in the early Bølling.

The transitional Bølling assemblage developed in the lake at the onset of the Bølling warming (14 443 cal yr BP, Birlo et al., 2023), indicated by a sharp peak in total oxygenic biomass (chlorophyll a and its derivatives, pigment group 2, Fig. 4) alongside cyanobacteria (zeaxanthin). This expansion was accompanied by the appearance of purple non-sulphur bacteria (OH-spheroidine, Fig. 3) and green sulphur bacteria (bacteriochlorophyll a). While the main diatom and chrysophyte community is represented by the continuous presence of diatoxanthin and diadinoxanthin, we also observe a singular, isolated peak of fucoxanthin. The phase is further characterized by a distinct turnover in the diatom community. Initially, the Pleniglacial benthic dominant, Staurosira construens, completely disappears. The benthic and periphytic niches abandoned by this species are rapidly occupied by a new shallow-water assemblage including Amphora indistincta, Pseudostaurosira brevistriata, Pseudostaurosira robusta, and Staurosirella lapponica (Fig. B7). As lake levels rose (Negendank and Zolitschka, 1993), A. indistincta gradually declined while the planktonic diatom Stephanodiscus minutulus rose in high abundances. By the mid-Bølling, the open-water assemblage further diversified with the appearance of Pantocsekiella ocellata. The specific composition of this transitional community suggests that the lake was partly nutrient limited, likely in nitrogen. However, the gradually increasing reactive manganese pool (F2 + F3) indicates that despite rapidly increasing primary production, the lake bottom waters remained at least seasonally oxygenated (Fig. B6).

The stratified Allerød assemblage is dominated by specific marker pigments for purple sulphur bacteria (okenone, bacteriopheophytin a) co-occurring with the continued high abundance of planktonic diatoms (S. minutulus and P. ocellata). The high concentrations of okenone are particularly significant, as they strictly indicate the development of stable, meromictic conditions where the chemical density gradient (chemocline) reached into the photic zone. This permanent stratification separated the actively circulating upper water column (mixolimnion) from the stagnant, euxinic deep water (monimolimnion). Crucially, this community structure remained remarkably stable throughout the Allerød, showing no substantial shifts in composition even following the deposition of the Laacher See Tephra (LST). Although the LST event introduced a significant amount of iron and silica, as well as abrupt changes in light penetration, the established meromictic stability buffered the ecosystem against this disturbance. The persistent euxinia below the chemocline was sufficient to maintain the robust bacterial plate of okenone-producing Chromatiaceae without interruption.

Figure 4 Stratigraphic variation of hyperspectral imaging data, selected sedimentary pigments, pigment communities derived from correlation (see Fig. 3), selected redox fractions (Mn), selected diatom species and pollen in Holzmaar (and Meerfelder Maar) (Litt and Stebich, 1999; Litt et al., 2009; Litt and Stebich, 2001; Litt, 2019, 2010). Downcore profiles of relative abundance of phycocyanin as RABD620, with solid line representing the mean value and dashed lines ±1σ, calibrated hyperspectral imaging bulk pigments data (expressed as [µg g⁻¹ of wet sediment]) with the limit of quantification as a solid vertical line, and HPLC-analysed pigment concentrations (expressed as [µg g⁻¹ of OM]) are plotted against age (cal yr BP). The pigment communities are plotted as z-scores. Benthic diatoms are plotted in dark red and planktonic diatoms are displayed in blue. Arboreal pollen is show in light green and non-arboreal pollen in orange (Litt and Stebich, 1999; Litt et al., 2009; Litt and Stebich, 2001; Litt, 2019, 2001). Finally, Holzmaar pollen is in solid colour, whereas Meerfelder Maar pollen is underlaying the Holzmaar pollen and is transparent. The horizontal gray bar marks the stratigraphic position of the Laacher See Tephra (LST) and the stratigraphy is following Lang et al. (2023); Litt and Stebich (1999) with the exact age boundaries following Birlo et al. (2023).

The Younger Dryas assemblage is defined by a sudden and profound restructuring of the phototrophic community, initiated by the sharp decline and eventual complete disappearance of the purple sulphur bacteria marker pigment okenone (Fig. 4). This total collapse of the anoxygenic bacterial plate is accompanied by a severe reduction in bacteriopheophytin a and total chlorophyll a concentrations. Diatom data for this interval are sparse, likely reflecting poor diatom preservation and a near-absence of siliceous primary production during this predominantly minerogenic, low-productivity interval, combined with the 4 cm sampling resolution of the diatom record; together these preclude a detailed reconstruction of the community turnover at this transition. For the subsequent re-establishment of the diatom community in the following Holocene, we refer readers to the higher-resolution diatom stratigraphy of Holzmaar reported by García et al. (2022). Despite the sparse diatom data during this phase, the high-resolution pigment stratigraphy (Fig. B4) demonstrates a clear multivariate departure away from the stable Allerød cluster, reflecting a return toward a lower-biomass, mixed-community space.

3.3 Redox sensitive fractions of Mn, P and Fe

The sequential extraction data differentiates the available/reactive pools of Mn, P and Fe under changing redox conditions and anoxia. All extracted fractions are shown in the supplementary data (Fig. B6). Based on these profiles, we focus on Fraction 2 (F2) to represent the total reactive pool of Mn, Fe and P available for remobilization from the sediment during lake stratification.

Overall, the most abundant fraction for all three elements was the residual/carbonate-bound fraction (F4), which is dissolved only at low pH. The second most abundant fraction was the reactive F2, which represents phases readily dissolvable under reducing conditions.

Total manganese concentrations were highest in the Pleniglacial sediments. With the onset of the Bølling and the increasing concentrations of okenone during the Allerød, the reactive fraction of Mn (F2) shows a decreasing trend (Fig. 4), consistent with release from the sediment to the anoxic water column (reductive dissolution). Only after the LST (LST) we observe a distinct increase in the Mn reactive pool sequestered in sediments (Fig. B6).

Reactive iron (F2) and phosphorus (F2) display strong stratigraphic covariance throughout the record. Fe is closely coupled with P within this reactive fraction, with the molar Fe:P ratio of F2 stabilizing at ∼ 1.5 throughout the Bølling and Allerød (Fig. B6). This value corresponds to the 3:2 stoichiometry of authigenic vivianite (Fe₃(PO₄)₂, Nriagu, 1972). In classic stratified systems, the onset of meromixis and deep-water anoxia typically drives internal nutrient loading: Fe(III)-(oxy)hydroxides undergo reductive dissolution, releasing their adsorbed phosphorus back into the water column. However, the stoichiometric coupling in our F2 pool indicates that in Holzmaar, the initial physical sequestration of P by meromixis was immediately followed by a post-depositional chemical reaction. The abundant availability of dissolved ferrous iron (Fe²⁺) under strict anoxia allowed this newly liberated phosphorus to be rapidly re-precipitated as authigenic vivianite. Therefore, it is this post-depositional mineral formation that permanently sequestered the phosphorus in the sediment, rather than mere adsorption onto iron oxides (which typically yields ratios >2, Slomp et al., 1996; Jensen et al., 1992). While vivianite is frequently below bulk XRD detection limits, this consistent stoichiometric coupling in the reactive pool is a widely accepted diagnostic tool for identifying its presence and role in phosphorus sequestration in lacustrine settings (O’Connell et al., 2015; Rothe et al., 2016).

Following the LST, the lake was geochemically enriched by the influx of volcanic material. Fe and Mn concentrations in the reactive fraction (F2) increase significantly while P remains stable, resulting in fluctuating Fe:P ratios and an excess of reactive iron relative to phosphorus. This mineral influx effectively transformed the lake into an even more efficient phosphorus trap. By locking phosphate into stable mineral phases such as vivianite, the sediment acts as an efficient trap that prevents nutrient release back into the water column. This iron-driven mechanism eliminates the internal loading effect common in carbonate-dominated or iron-poor lakes, ultimately preventing the ecosystem from being held in a prolonged eutrophic or anoxic state fueled by legacy phosphorus.

4 Discussion

4.1 Development of vegetation

The pollen records from Meerfelder Maar (MFM) and Holzmaar (HZM) reveal a transformation of the catchment during the transition from the Pleniglacial to the Late Glacial (Fig. 4). The onset of the Bølling (∼14 450 cal yr BP) is marked by an increase in Juniperus concurrent with a relative decline in Pinus. Only in the late Bølling, the gradual establishment of Betula is observed.

At the onset of the Allerød, Betula became the dominant tree taxon in the catchment, although total Arboreal Pollen (AP) abundances remained below the closed forest threshold (<80 %, Lang, 1994; Sugita et al., 1999). This expansion coincides stratigraphically with the onset of preservation of biochemical varves and the appearance of indicators of anoxia.

Progressing through the Allerød, the vegetation shifted from Betula to Pinus. The percentages of Betula pollen gradually declined, while Pinus pollen increased, eventually replacing Betula as the dominant tree species by the late Allerød. The termination of the anoxic phase (disappearance of okenone) at the onset of the Younger Dryas coincided with the decline of Betula. In contrast, regional AP values remained high (∼70 %, Litt and Stebich, 1999) due to the persistence of Pinus, indicating that the recovery of the lake's oxygenation state occurred despite sustained forest cover.

4.2 Shift of the climate-driven primary producer community

The high-resolution pigment and diatom stratigraphy of Holzmaar reveals distinct ecological shifts that mirror the physical state of the lake. Rather than representing an internal autogenic succession, these communities demonstrate highly sensitive, allogenic responses to climate-driven changes in temperature, wind shielding, and mixing regimes (Fig. 6).

• The Pleniglacial phase corresponds to a holomictic water column heavily influenced by cold temperatures and high detrital input (Ti). This signal reflects local loess deflation from the Rhine Valley, driven by high winds and limited vegetation cover, and cold mean annual temperatures (insolation, Figs. 5 and B2). The primary producer community during this phase is characterized by abundances of various benthic diatoms, such as Staurosira construens (García et al., 2022; Lotter and Bigler, 2000), alongside cryptophytes (alloxanthin) and cyanobacteria (echinenone and canthaxanthine, Fig. 4). Cryptophytes traced by the pigment alloxanthin are characteristic of cold, well-mixed, minerogenic glacial lakes adapted to low-light benthic conditions (Züllig, 1986). Further, the presence of these small fragilarioid diatoms infers a low-light environment and a highly unstable landscape. These small benthic taxa are well-documented low-light specialists possessing a distinct ecological advantage that allows them to thrive in highly turbid, minerogenic conditions or under prolonged seasonal ice cover (Kingsbury et al., 2012; Gushulak and Cumming, 2020). Consistent with these shallow-water conditions, the ratio of chlorophyll a to pheophytin a, a proxy for water-column photodegradation, shows a significant positive correlation with the abundance of S. construens ($r=+\mathrm{0.488}$, p=0.0014), supporting the hypothesized link between shallower water-column residence time and reduced pigment degradation (Leavitt and Hodgson, 2002).

  While overall primary production was low during the Pleniglacial, the community structure indicates that nitrogen was a key limiting nutrient. This inference of nitrogen limitation is supported by two concurrent signals: the relative abundance of cyanobacteria, which are highly efficient at operating in nitrogen-poor environments or utilizing nitrogen fixation (Guildford and Hecky, 2000; Engstrom et al., 2000), and the expansion of Juniperus in the catchment, a terrestrial pioneer taxon known to colonize young, nitrogen-poor soils (Iversen, 1954; Matthews, 1992).

  In the open water, the pelagic pigment profile further reveals a community operating under severe environmental and nutritional strain. The prominent signal of astaxanthin points to cold-tolerant microalgae like Haematococcus pluvialis or dinoflagellates. In H. pluvialis, astaxanthin accumulation is a well-documented physiological response to extreme environmental stress, including cold temperatures, high irradiance, and severe nutrient or nitrogen starvation, triggering the formation of protective resting cysts (Lemoine and Schoefs, 2010; Shah et al., 2016; Boussiba, 2000). This is further corroborated by the presence of diatoxanthin, a specific xanthophyll produced by diatoms undergoing light stress (Kuczynska et al., 2015). Together, these stress-induced pigments depict a harsh, ultra-oligotrophic water column.

  Despite these extremes, the presence of green sulphur bacteria (isorenieratene) highlights survival within highly specific low-light niches. Because green sulphur bacteria are obligate anaerobes and extreme low-light specialists, their co-occurrence with low-light tolerant diatoms suggests they occupied transient, shaded anoxic micro-niches. These habitats likely developed under seasonal ice cover or within the turbid cryolittoral zone, analogous to modern assemblages found in polar lakes (Antoniades et al., 2009) or the early developmental stages of Lake Jaczno (Makri et al., 2021).

• The Bølling phase reflects the immediate biological response to rapid climatic warming (δ¹⁸O_NGRIP, Figs. 5 and 7) and rising lake levels, resulting in a dynamic internal sequence of community shifts. The almost complete disappearance of S. construens marks the abrupt crossing of an ecological threshold, signaling the end of the harsh, shallow Pleniglacial environment. As lake levels initially rose, newly flooded, warmer littoral zones were rapidly colonized, driving the early peak in Amphora indistincta thriving on submerged macrophytes or mineral substrates in the littoral zone (Levkov, 2009). Its expansion, simultaneous with the sharp rise in cyanobacterial pigments (Fig. 4), likely reflects a rising lake level and the creation of new littoral habitats (Negendank and Zolitschka, 1993). However, as the basin continued to fill, these early shallow habitats were progressively submerged, leading to a gradual decline in the A. indistincta signal within the central basin. Simultaneously, the focus of primary production shifted to the newly expanding pelagic zone, marked by sharp rises in total chlorophyll a, cyanobacteria, and the early proliferation of the heavy planktonic diatom Stephanodiscus minutulus.

  Interestingly, the mid- to late-Bølling record reveals two distinct spikes of okenone, indicating brief episodes of photic zone euxinia and an early trend toward chemical stratification. This early trend toward stratification and deep-water oxygen depletion is independently supported by regional molecular lipid biomarker records (Zander et al., 2024). Specifically, the continuous presence of the branched GDGT isomer IIIa”, a compound whose fractional abundance strongly correlates with production under low-oxygen conditions, confirms the establishment of hypolimnetic anoxia in Holzmaar as early as 14.2 ka BP (Zander et al., 2024). Furthermore, high values of the %GDGT-0 index throughout the Late Glacial indicate substantial activity of methanogenic archaea within the lake (Zander et al., 2024). Together, these lipid biomarkers strongly corroborate our pigment-based evidence, confirming that a stable, highly reducing environment had already developed in the deep water well before the permanent photic-zone euxinia of the Allerød phase.

  However, the continued proliferation of S. minutulus requiring active turbulent mixing to overcome its natural sinking velocity (Huisman et al., 2002; Bradbury and Dieterich-Rurup, 1993) suggests that despite this incipient stratification, Holzmaar maintained a robust seasonal mixing regime (likely monomixis). This mixing was crucial for allowing nutrient regeneration and keeping the heavy diatoms in suspension. Alongside S. minutulus, the delayed appearance of P. ocellata by the mid-Bølling signifies the full establishment of a mature, deep, open-water habitat, seamlessly setting the physical and biological stage for the transition to stable meromixis in the subsequent Allerød.

  During the Bølling phase, the watershed was colonized by Juniperus (Fig. 4). As a light-demanding pioneer, its presence implies an open canopy structure before the establishment of forests. We propose that the interplay of rising lake levels (dilution) and developing littoral zones created dynamic niches that were rapidly filled by this secondary community. Geochemically, the gradual increase in bioavailable manganese (Mn F2 + F3) confirms that these changes were accompanied by efficient Mn sequestration during oxic conditions linked to regular lake mixing.

• The Allerød phase marks the transition to stable meromixis, driven by summer warming and the wind-shielding effect of Betula establishment in the watershed (Figs. 5 and 7). The dominance of okenone (Fig. 4) signals that the chemocline reached into the photic zone. This indicates the formation of a dense, self-shading bacterial plate comprised of purple (Chromatiaceae) and green (Chlorobiaceae) sulphur bacteria, analogous to the deep chlorophyll maxima observed in modern Lake Cadagno (Tonolla et al., 2005; Zander et al., 2023). The simultaneous deposition of these anaerobic bacterial pigments alongside planktonic diatoms (S. minutulus and P. ocellata) reflects a highly structured, compartmentalized water column. Our high-resolution HSI mapping (Fig. 6) captures the micro-stratigraphy, which reveals two distinct depositional scenarios: layers containing both diatom-associated chlorophylls and bacterial okenone (indicating simultaneous summer production), interbedded with layers where these pigments are strictly separated (indicating strict seasonal separation).

  This micro-stratigraphy is driven by a combination of physical mixing and biological adaptability, operating strictly above the permanent, euxinic deep water (the monimolimnion). Based on the HSI data, we propose two non-exclusive mechanisms for this community structure:

  • 1.

    Temporal separation (Seasonal dynamics): Because heavy diatoms like S. minutulus and P. ocellata rely on turbulent mixing to remain in suspension (Vossel et al., 2015; Kienel et al., 2017), their growth was actively supported during periods of robust spring and autumn circulation within the upper mixolimnion. This strict seasonal separation results in the interbedded layers of exclusively bacterial or algal pigments and siliceous layer.

  • 2.

    Spatial co-existence (Deep Chlorophyll Maximum): During the summer, the mixolimnion is itself thermally stratified. While mixing ceased in the upper epilimnion, the strong biological adaptability of S. minutulus (Fietz and Nicklisch, 2002; Gladyshev et al., 2010; Celewicz and Gołdyn, 2021; Estepp and Reavie, 2015) allowed it to maintain a competitive strategy in a light-limiting environment. These diatoms likely formed a Deep Chlorophyll Maximum (DCM) within the metalimnion (Spanbauer et al., 2014). Positioned deep in the water column, light availability was severely attenuated by the overlying epilimnetic production. Here, the diatoms successfully occupied a shaded niche vertically distinct from, and immediately above, the permanent self-shading bacterial plate at the chemocline (Guasch and Subater, 1995; Köhler et al., 2017; Pestryakova et al., 2018). The natural settling of these diatoms through the actively producing bacterial plate yields the sediment layers containing both pigments simultaneously.

  Ultimately, this meromictic structure proved to be highly resilient. The persistence of this community throughout the Allerød confirms that the lake entered a stable state of meromixis with persistent monimolimnetic anoxia, resilient to minor climate fluctuations, such as centennial-scale oscillations in warm-season temperatures, periodic shifts in humidity, and the short-lived cooling of the Intra-Allerød Cold Period (IACP) documented in the Eifel region (Zander et al., 2024; Litt and Stebich, 1999).

• The Younger Dryas marks the final major community shift in the Late Glacial sequence, demonstrating the clear and rapid reversibility of the Holzmaar ecosystem under intense physical forcing. Driven by sharp regional cooling and increased wind-driven wave action, the chemical and thermal density gradients that sustained the Allerød meromixis completely collapsed. The total disappearance of okenone signals the rapid breakdown of the deep bacterial plate at the chemocline. As the thermal gradient weakened, the lake transitioned from its highly stable meromictic state back into a holomictic, deeply circulating mixing regime.

  This complete seasonal water-column turnover re-oxygenated the previously stagnant deep waters, abruptly terminating the long-term photic zone euxinia of the monimolimnion. The trajectory shift in the continuous pigment record back toward a lower-biomass, mixed-community space (Fig. B5) highlights that the physical forcing of the climatic cooling rapidly overrode the internal biological and chemical buffers of the lake. Ultimately, while the meromictic structure proved resilient enough to withstand centennial-scale oscillations (such as the IACP) and acute volcanic impacts (such as the LST), it remained fundamentally vulnerable to the large-scale, persistent reorganization of the regional climate and mixing regime.

Figure 5 Redundancy Analysis (RDA) of sedimentary pigment assemblages and environmental drivers. The biplot displays the relationship between the primary producer community (pigments, response variables) and independent environmental proxies (explanatory variables). Individual samples are displayed as solid points, colour-coded according to their pollen zones. The identified pigments are colour-coded according to the three pigment groups assigned in Fig. 3. Black arrows represent the significant explanatory variables selected by VIF filtering (VIF ≤ 10, p<0.05), including catchment vegetation (Salix, Juniperus pollen), dust input (Dust_NGRIP), GDGT-derived lake temperature in months above freezing (T_GDGT, Zander et al. (2024)), chironomid-based July temperatures at Gerzensee, Switzerland (T_Ger July, Lotter et al. (2012)) and Egelsee (T_Egel July, Larocque-Tobler et al., 2009), allochthonous minerogenic input (Ti_XRF and Ca_XRF), and stable isotopes (δ¹⁸O_NGRIP). Green arrows indicate passively added ecological representatives added to ensure better data interpretability: Betula and Pinus vegetation proxies, and winter insolation. The first two axes explain 28.3 % and 19.8 % of the total variance, respectively.

Figure 6 Selected section of core HZM19_11_7 from 65.2–70.5 cm (∼ 13 015–13 102 cal yr BP) demonstrating hyperspectral imaging data (HSI) of varves and pigment distributions represented by RABDs immediately below the LST. (A) Green colour and lines (RABD671) correspond to chlorophyll a derivatives (green algae/diatoms), while purple/pink colour and lines (RABD845) indicate bacteriopheophytin a (purple sulphur bacteria). (B) The same section as in panel (A), where purple/pink colour and lines (RABD845) indicate bacteriopheophytin a (purple sulphur bacteria) and orange colour and lines (RABD715) show distribution of bacteriochlorophyll e (green sulphur bacteria). Both, the separation of layers as well as overlaps of the pigments are observed. Summer layers are highlighted in the plot area by red overlay and spring layers are highlighted in green. Dashed lines represent mean RABD values of the section.

4.3 Anoxia drivers

Disentangling the drivers of anoxia requires separating the regional temperature trends from local catchment forcing. While local warm-season temperatures increased moderately by $\sim +\mathrm{2.8}$ °C in the period from 14 642 to 12 846 cal yr BP (Zander et al., 2024), the primary driver of the presence of anoxia in the photic zone was likely the extension of the growth season (Bekaert et al., 2023). This shift prevented deep freezing and extended the thermal stratification period into warm-seasons, providing the physical boundary conditions for anoxia. Given the latitudinal position of Holmzaar (∼ 50° N), a permanent, perennial ice cover during the Pleniglacial is highly unlikely, and is further argued against by the continuous, undisturbed deposition of sediments preserved throughout both the Pleniglacial and the Younger Dryas. However, an extended seasonal or unusually severe winter ice cover cannot be excluded and would have reinforced deep-water isolation during these colder intervals. The specific timing of the euxinic phase, which lags the initial warming, suggests an influence of catchment-mediated processes, such as vegetation and soil development, to cross the anoxia threshold.

In addition to regional climate, the onset of anoxia was likely pre-conditioned by the specific basin morphometry of Holzmaar. With a relative depth (Z_r, a ratio of maximum depth to lake mean diameter) of 7.4 %, the volcanic maar basin is morphologically predisposed to stable thermal stratification by decoupling the hypolimnion from wind-driven mixing (Lehner, 2024; Magee and Wu, 2017). Global predictive models suggest that a high Z_r, combined with rising temperatures, moved the onset of spring stratification earlier in the season (Demers and Kalff, 1993).

While direct evidence for lake-level fluctuations at Holzmaar is limited, regional records indicate higher lake levels during the Bølling-Allerød across Central Europe (Magny, 2001). Such a deepening would have increased Z_r and physically primed the system for anoxia. Within this physically sensitive framework, the expansion of Betula in the watershed acted as a reinforcing factor. However, given the high crater rim ($<+\mathrm{35}$ m) and significant relative depth (−20 m), the threshold of arboreal pollen (AP) required to shield the lake from wind was likely much lower than in shallower basins. This implies that the lake was physically predisposed to stability regardless of subtle changes in forest density.

The onset of intense anoxia coincides precisely with the expansion of Betula in the watershed and a reduction in allochthonous sediment yield (from 16 to <2 t km⁻² yr⁻¹), as quantified by Zolitschka (1998). This confirms that vegetation stabilized the crater slopes and its catchment. We propose that Betula afforestation strongly reinforced lake stratification primarily through wind shielding. The forestation of the crater rim acted as a roughness filter, dampening the wind shear required for deep mixing (Dietrich and Seelos, 2010; Makri et al., 2020; Zander et al., 2021b; Wienhues et al., 2026).

At Holzmaar, anoxia disappeared precisely at the onset of the Younger Dryas (YD), although Pinus coverage increased (Fig. 7). However, because Pinus pollen is highly susceptible to long-distance aeolian transport, these regional AP values must be interpreted conservatively. This recovery highlights that climatic cooling and a regional shift towards less humid conditions were the decisive drivers for breaking stable chemical stratification. The temperature decline and increased aridity overrode the physical shielding of the crater rim by increasing water density and enhancing convective mixing. Furthermore, despite the highly arid regional climate, the development of deep seasonal soil frost and frozen ground conditions strongly restricted groundwater infiltration (Isarin, 1998; Renssen et al., 2001), which substantially amplified seasonal surface run-off (elevated Ti). This mechanism explains the sustained delivery of allochthonous minerogenic sediment into the basin despite reduced overall precipitation (Brauer et al., 1999). Consequently, the shift in forest composition (from Betula to Pinus) is interpreted as a climatically-driven consequence rather than the main cause of re-oxygenation. This abrupt physical and ecological transition is clearly reflected in the stratigraphic trajectory of pigment assemblages, which exhibits a rapid return toward the oligotrophic space (Fig. B5).

Figure 7 Multi-proxy summary of climatic, vegetational, and geochemical transitions at Holzmaar. The top panels display (from top to bottom): NGRIP δ¹⁸O and regional temperature reconstructions, including Gerzensee July T, the Alpine July T stack, and local GDGT-derived Temperatures of Months Above Freezing (GDGT TMAF). Mid-panels show total chlorophyll concentrations (a proxy for primary production) and Betula and Pinus pollen percentages (with modelled sections indicated). Lower panels illustrate the concentrations of the anoxia indicator bacteriopheophytin a (unsmoothed) alongside the relative abundance of phototrophic pigment communities (Pleniglacial, Early Bølling, and Allerød). The bottom panel presents reactive manganese (Mn) redox fractions F2 and F3. Vertical dashed lines delineate the major biozones. Data sources: NGRIP (NGRIP Members, 2004), Gerzensee July T (Lotter et al., 2012), July T stack (Heiri et al., 2014), T_GDGT (Zander et al., 2024), HZM Betula pollen (Litt, 2010, 2001).

4.4 Geochemical buffering and ecosystem resilience

The geochemical response of Holzmaar reveals a highly efficient nutrient retention mechanism that distinguishes it from purely iron-controlled (Søndergaard et al., 2003; Gächter and Müller, 2003) or carbonate-controlled lakes (Tu et al., 2020). Within the reactive pool (Fraction 2), the molar Fe:P ratio stabilizes at ∼1.5 throughout the anoxic phase (Fig. B6). This stoichiometry is the diagnostic signature of formation of authigenic ferrous phosphates such as vivianite (Fe₃(PO₄)₂⋅8H₂O) in iron-rich sediments (Nriagu, 1972; Rothe et al., 2016; Kleeberg et al., 2013). It indicates that the volcanic iron flux (Förster et al., 2019) provided a sufficient reactive buffer to sequester phosphorus in the sediment throughout the DOE-1, even under reducing conditions (O’Connell et al., 2015).

The sequential extraction data further identify the acid-soluble Fraction 4 (HCl, Fig. B6) as the dominant long-term sink for manganese, iron, and phosphorus. We interpret this stable pool as the result of a dual-retention mechanism: initially, the reactive iron pool (F2) traps porewater phosphorus as vivianite ((Fe₃(PO₄)₂⋅8H₂O), and subsequently, this phosphorus is buried within or converted into stable mineral phases in Fraction 4 that do not redissolve during redox shifts. By permanently locking P, Mn and Fe into this stable sink, the lake avoided the development of a chemical legacy or an internal phosphorus loading (Søndergaard et al., 2003). The biological consequence of this geochemical resilience is a high degree of community reversibility, as the pigment assemblages effectively return to the oligotrophic community once the climatic drivers of the Allerød (warming) were removed (Fig. B5).

This geochemical framework provides critical insight into lake ecosystem resilience and rapid high reversibility. We attribute this resilience to a specific decoupling of nutrient cycles under an iron-unlimited regime. While manganese followed the expected reductive dissolution-precipitation pattern, being recycled back into the water column during anoxia (depleted F2) and precipitated during the seasonal mixing (increase F2 + F3), phosphorus behaved differently. In many eutrophic or carbonate-dominated systems, such as Lake Burgäschi, Lake Baldeggersee or Lake Soppensee (Swiss Plateau, Tu et al., 2020; Lotter, 1998; Schouten et al., 2026), the release of legacy phosphorus creates a delayed recovery, sustaining anoxia long after external climate forcings are removed (Carpenter, 2005; Søndergaard et al., 2003). In those systems, sulfate reduction often strips iron from the sediments to form minerals like pyrite, releasing phosphorus back into the water column (Tu et al., 2020; Lotter, 1998; Schouten et al., 2026).

Consequently, regional cooling at the onset of the Younger Dryas allowed the system to re-oxygenate, evidenced by the decline of anoxia followed by the vegetation shift (Betula to Pinus) with a minimal time lag (Fig. 7). Ultimately, while climate and vegetation set the boundary conditions, the specific coupling of catchment characteristics and an iron-rich lithogenic buffer determined the resilience of the aquatic ecosystem.

5 Conclusions

This study reconstructs the high-resolution biological and geochemical responses of a volcanic maar lake to rapid natural warming (DOE-1), providing a mechanistic template for understanding climate-driven anoxia without anthropogenic interference.

We demonstrate that hyperspectral imaging can resolve micro-stratigraphic shifts in phototrophic communities within varved sediments. Holzmaar progressed from oligotrophic conditions to an eutrophic stage where the chemocline reached into the photic zone. This allowed purple sulphur bacteria to thrive alongside oxic low-light adapted planktonic diatom (Stephanodiscus minutulus), creating a complex biological structure analogous to modern meromictic systems where primary production is maintained across steep light and redox gradients.

Further, our multi-proxy approach reveals that atmospheric warming was a necessary driver for sustained anoxia, though its impact was modulated by local factors. The onset of stratification was likely pre-conditioned by the lake’s high relative depth (Z_r=7.4 %), a physical state that may have been amplified by climate-driven increases in lake level. Within this physically sensitive framework, the expansion of Betula in the watershed acted as a reinforcing factor for lake stratification through wind shielding, dampening the wind shear required for deep mixing. The rapid recovery of the lake was triggered primarily by the onset of the Younger Dryas cooling and increased aridity, which overrode this physical shielding and broke the stable stratification. The concurrent shift in catchment vegetation from Betula to Pinus is interpreted as a climatically-driven consequence of this cooling and aridification, rather than an active driver of the re-oxygenation.

Finally, we identify a distinct lithological control on aquatic resilience and ecosystem reversibility. In many iron-poor or carbonate-dominated lakes, recovery from anoxia is significantly delayed by the continuous feedback loop of internal phosphorus loading, the biogeochemical process by which sediment-bound nutrients are recycled back into the photic zone. Where sustained, this feedback loop is the proximate mechanism behind the broader phenomenon of ecosystem hysteresis: a delayed, non-linear recovery of the whole-lake ecosystem state that persists after the original climatic driver has been removed. Holzmaar, however, functioned as an exceptionally efficient nutrient trap. The volcanic catchment provided a continuous flux of reactive iron that permanently sequestered phosphorus into stable mineral phases like authigenic vivianite, preventing its release via reductive dissolution during anoxia. Due to the complete absence of this internal phosphorus loading feedback, the lake lacked the biogeochemical mechanism that produces hysteresis in other systems. Consequently, the biological community exhibited rapid ecosystem reversibility rather than a hysteretic response, breaking its meromictic state and returning to an oxygenated, mixed-community space immediately upon the onset of Younger Dryas climatic cooling and drying.

These findings suggest that while climate warming sets the physical potential for anoxia, the vulnerability and recovery of lake ecosystems are driven by the specific coupling between temperature, regional moisture balance, catchment vegetation, and their lithogenic redox buffer.

Appendix A: The napari-sediment workflow

napari-sediment workflow consists of three main components which follow the original workflow by Butz et al. (2015):

1. Image Loading and Preprocessing: Raw hyperspectral data stored in ENVI HDR format were imported and normalized using white and dark reference images. An interactive masking tool allowed manual annotation of corrupted regions (e.g., core breaks, edge effects, or moisture artifacts) that were excluded from subsequent analysis. Region-of-interest (ROI) selection tools enabled focusing analysis on specific core sections or depths.

2. Dimensionality Reduction: To reduce data volume and computational complexity while retaining diagnostic spectral features, two complementary dimensionality reduction approaches were applied. Spectral dimensionality reduction was performed using Minimum Noise Fraction (MNF) transformation, which maximizes the signal-to-noise ratio by separating signal from noise through two successive principal component analyses. The first transformation decorrelates and rescales the noise in the data, and the second decorrelates the noise-whitened data. MNF components are ordered by decreasing signal-to-noise ratio, with the first components containing coherent spectral information and later components dominated by noise. Spatial dimensionality reduction employed pixel purity index (PPI) algorithms to identify spectrally pure pixels representing distinct sediment facies. These pure pixels were subsequently clustered to identify representative spectral endmembers.

3. Spectral Index Calculation and Mapping: After endmember identification, spectral indices (RABD, RABA) were calculated based on characteristic absorption features. The software enabled interactive selection of wavelength ranges for baseline definition and absorption minimum identification. Index maps were generated showing the spatial distribution of each index across the core, and 1D projection profiles were extracted by averaging index values across user-defined ROI widths. These profiles were optionally smoothed using Savitzky-Golay filtering with adjustable window sizes. All processed data, including index maps, projection profiles, and endmember spectra, were exported for statistical analysis and calibration.

The interactive nature of napari-sediment provided several advantages: (i) real-time visualization of preprocessing effects enabled quality control at each processing step, (ii) manual masking allowed expert knowledge to guide data cleaning, (iii) interactive ROI selection facilitated targeted analysis of specific sedimentary features, and (iv) the ability to iteratively adjust spectral index parameters optimized signal detection for specific pigments. The software architecture, built on the napari viewer framework, enabled GPU-accelerated rendering of large hyperspectral datasets and seamless integration with Python scientific libraries (NumPy, SciPy, scikit-learn) for advanced data analysis.

Figure A1 Calibration and validation of the hyperspectral chlorophyll proxy. Linear regression model establishing the quantitative relationship between the spectral index (RABD₆₇₁; Relative Absorption Band Depth at 671 nm measured on wet sediment) and spectrophotometrically determined Total Chlorophyll concentrations (µg g${}_{\text{wet sed.}}^{-\mathrm{1}}$). The model demonstrates a strong linear correlation (R²=0.92, p<0.001, n=33), validating RABD₆₇₁ as a robust high-resolution proxy for primary production. The red line represents the linear fit, while the black dashed and green dotted lines indicate the 95 % confidence and prediction intervals, respectively. Cross-validation (RMSEP) indicates a prediction error of approximately 3.2 µg g⁻¹ (∼10 %).

Figure A2 Calibration and validation of the hyperspectral bacteriopheophytin a proxy. Linear regression model establishing the quantitative relationship between the spectral index (RABD₈₄₅; Relative Absorption Band Depth at 845 nm measured on wet sediment) and spectrophotometrically determined bacteriopheophytin a concentrations (µg g${}_{\text{wet sed.}}^{-\mathrm{1}}$). The model demonstrates a strong linear correlation (R²=0.82, p<0.001, n=33), validating RABD₈₄₅ as a robust high-resolution proxy for purple sulphur bacteria (PSB) abundance. The red line represents the linear fit, while the black dashed and green dotted lines indicate the 95 % confidence and prediction intervals, respectively. Cross-validation (RMSEP) indicates a prediction error of approximately 0.94 µg g⁻¹ (∼14.6 %), confirming the proxy's sensitivity for detecting anoxia-driven pigment layers.

Roy et al. (2011)Deshpande et al. (2014)Leavitt and Hodgson (2002)Leavitt and Hodgson (2002)Guilizzoni and Lami (2003)Lami et al. (1994)Bianchi and Canuel (2011)Roy et al. (2011)Puddick et al. (2023)Senger et al. (1993)Orosa et al. (2000)Takaichi and Mochimaru (2007)Deshpande et al. (2014)Roy et al. (2011)Puddick et al. (2023)Senger et al. (1993)Roy et al. (2011)Deshpande et al. (2014)Puddick et al. (2023)Roy et al. (2011)Bertrand (2010)Roy et al. (2011)Bertrand (2010)Roy et al. (2011)Lami et al. (2000)Senger et al. (1993)Orosa et al. (2000)Roy et al. (2011)Pennington et al. (1985)Lotter (2001)Züllig (1982)Roy et al. (2011)Senger et al. (1993)Orosa et al. (2000)Roy et al. (2011)Makri et al. (2018)Lemoine and Schoefs (2010)Shah et al. (2016)Oelze (1985)Smith et al. (2013)Coolen and Overmann (1998)Züllig (1986)Coolen and Overmann (1998)Züllig (1989)Vila et al. (1998)Biebl and Pfennig (1978)Simpson and Smith (1988)Chen et al. (2001)Glaeser et al. (2002)Zander et al. (2023)Albrecht et al. (1997)Davies et al. (1969)

Table A1 Sedimentary pigments identified in Holzmaar, their taxonomic affinities, and key references.

Appendix B: Construction of the Holzmaar Late Glacial stacked pollen record

To provide a continuous palynological context for the 2019 Holzmaar (HZM) sediment sequence, we integrated pollen datasets from Litt and Stebich (1999) and Litt et al. (2009). These records were derived from cores collected in 1992 (Litt, 2019) and 1996 (Litt, 2010), respectively. The integration into the 2019 master chronology followed a two-step alignment process:

1. Section-to-Core Correlation (1992 Core to 2019 Master)

The pollen data from the 1992 campaign Litt and Stebich (1999) were originally available only on a depth scale within individual core sections (specifically sections HZM92 1a–6u, 2b–5o, and 2b–5u). To be able to use these data, we needed to first transfer these data on our 2019 individual cores, and then we could assign our age-depth model (VT22, Birlo et al., 2023) to these data. As expected, the 1992 and 2019 coring locations differed, localized stratigraphic stretching and compression were observed.

• Correlation method: We performed a visual lithostratigraphic correlation using high-resolution core photographs. Distinctive marker layers and varve patterns served as tie-points to align these individual 1992 sections with the 2019 individual cores, and fuhrer to the HMZ19 composite depth.

• Transformation: This alignment allowed us to transfer the 2019 age-depth model to the 1992 sampling depths, synchronizing the pollen counts with our new high-resolution geochemical and pigment datasets as well as with the complementary, but lower resolution pollen data from Litt et al. (2009).

2. Multi-core stacking (1996 core)

The pollen data from the 1996 core (Litt et al., 2009) were available only on a calibrated age scale.

• Integration: Once the 1992 data were successfully migrated to the 2019 age scale (VT22), we stacked the two datasets (cores from 1992 and 1996) based on their age scale.

• Result: This produced a comprehensive stacked pollen record for Holzmaar, anchored to the 2019 chronology and depth scale. This approach minimizes stratigraphic gaps and provides a framework for interpreting vegetation-driven changes in lake stratification.

Figure B1 Correlation matrix of environmental variables used in RDA analysis. Pearson correlation coefficients between the 16 environmental predictor variables (T_GDGT_MAF = GDGT temperature of months above freezing, $\mathrm{JulyT}\mathit{\_}\mathrm{Ger}/\mathrm{Egel}$ = July temperatures from Gerzensee and Egelsee, ${\mathrm{insol}}_{\mathrm{A}}/\mathrm{July}/\mathrm{Jan}$ = annual/seasonal insolation, $\mathrm{18}\mathrm{O}\mathrm{\_}\mathrm{NGRIP}/\mathrm{Dust}\mathrm{\_}\mathrm{NGRIP}$ = stable isotope and dust proxy from NGRIP ice core, $\mathrm{Ca}/\mathrm{Ti}$ = geochemical elements from XRF analysis, pollen percentages = BetulaHZM, JuniperusHZM, PinusHZM, SalixHZM, NAP_HZM, AP_HZM). Circle size and color intensity indicate correlation strength, with red indicating positive correlations and blue indicating negative correlations. Absolute correlation values are displayed for all variable pairs. This matrix demonstrates substantial inter-relationships among climatic predictors, particularly between temperature metrics and orbital insolation variables, which justified the application of a VIF < 10 threshold to balance multicollinearity control while retaining ecologically meaningful temperature signals for the constrained ordination analysis.

Figure B2 The biplot displays the correlation structure among the independent proxies considered for the Redundancy Analysis (RDA), revealing the primary environmental gradients driving the system. The first axis (Dim1, 46.4 %) represents the dominant climatic and vegetational gradient. It strongly opposes indicators of warmth and forest development (positive loadings: Betula, summer insolation, regional summer temperatures) against markers of cold, open steppe environments (negative loadings: Non-Arboreal Pollen (NAP), Juniperus, winter insolation). The second axis (Dim2, 25.1 %) captures the sedimentological signal, distinguishing allochtonous detrital inputs (positive loadings: Ti, NGRIP Dust) from autochthonous carbonate precipitation (negative loading: Ca) and climate signals (δ¹⁸O_NGRIP). This analysis demonstrates the high collinearity between local vegetation and regional temperature records, justifying the VIF selection of a reduced subset of non-redundant variables for the final RDA model.

Figure B3 Principal Component Analysis (PCA) variable correlation plot of high-resolution XRF geochemical data. The biplot displays the loading vectors of elemental proxies to the first two principal components. The first principal component (PC1) explains 45.3 % of the total variance and represents the gradient between allochtonous catchment inputs and autochthonous lake production. It separates detrital/siliciclastic elements (strong positive loadings: Ti, K, Zr, Al, Si, Rb) from elements associated with endogenic carbonate precipitation (strong negative loadings: Ca, Sr). The second principal component (PC2) explains 11.5 % of the variance, primarily distinguishing minerogenic elements (Ti, K, Al, Ca, positive loading) from elements associated with organic matter binding and redox cycling (negative loadings: Mn, Fe, Cu, Zn, and coherent scattering ${\mathrm{Cr}}_{\mathrm{inc}/\mathrm{coh}}$). S shows a minor positive loading, likely being associated with diagenetic processes of organic matter under anoxic conditions,, forming minerals such as pyrite or vivianite. This analysis validates the use of Ti as a robust proxy for detrital/minerogenic material supply and Ca and Mn as indicators of in-lake geochemical processes. The ellipses and their colours correspond to the unconstrained cluster displayed in Fig. 2.

Figure B4 Stratigraphic profiles of absolute sedimentary pigment concentrations and degradation ratios at Holzmaar spanning Dansgaard-Oeschger Event 1 (ca. 14 690–11 700 cal yr BP). To effectively account for variations in sediment composition and matrix effects, concentrations for each identified HPLC pigment are presented using a dual-axis format: values normalized to total organic matter (µg g⁻¹ of OM) are shown on the top axes (darker bars), while absolute concentrations per gram of bulk wet sediment (µg g⁻¹ of bulk sediment) are shown on the bottom axes (lighter bars). The pigment profiles are color-coded to align with the primary producer groups discussed in the main text. The far-right panels display the Chlorophyll $a/$Total Pheophytins ratio ($\mathrm{Chl}/\mathrm{Pheo}$) and the Pheophorbides/Pheophytins ratio. While the Chl $a/$Pheophytin ratio is provided here to explore potential in-lake photodegradation dynamics related to changing water depths, these specific pigments degrade rapidly upon exposure to oxygen, light, and temperature changes during sediment coring and extraction. Consequently, this ratio may be susceptible to methodological artifacts and must be interpreted with caution.

Figure B5 Stratigraphic trajectory of sedimentary pigment assemblages in PCA space (PC1 vs. PC3). The scatterplot displays the temporal evolution of the phototrophic community, with sample scores connected in chronological order. Points are coloured by stratigraphic zone: Pleniglacial (green circle), Bølling (orange triangle), Allerød (purple square), and Younger Dryas (pink cross). The large red and green dots mark specific time markers (ages in cal yr BP) to indicate the direction of change. The trajectory reveals the reversibility of the primary producer community. The transition from the Bølling to the Allerød involves a shift toward negative PC1 and negative PC3 values, a region defined by the dominance of okenone (purple sulphur bacteria; see Fig. 3). The subsequent transition to the Younger Dryas marks a return toward the positive PC1 space and positive PC3 (cyanobacteria and low-light pigment producers), indicating the breakdown of stable euxinia and the recovery of the oxic/mixed community structure.

Figure B6 Fe, Mn and P fractions with the anoxic indicators: okenone and calibrated RABD845. Fe:P ratio of the fractions 2 are used to interpret what form of Fe and P we have present. Fe:P < 1, iron limitation, P in F2 can be present as organic bound P; if S present (okenone), Fe is bound as pyrite (FeS₂) and not available to bound with P. Fe:P between 1.5–2.5 is vivianite, Fe:P > 2.5, iron excess, likely to be bound as siderite or iron oxides.

Figure B7 Clustered heatmap of diatom taxa abundance across Late Glacial. Rows represent individual samples with their age (cal yr BP, right side), and columns represent individual diatom species. Taxa are categorized by ecological niche (planktonic vs. benthic + epiphytic) and clustered by similarity in occurrence patterns (correlation/euclidean distance). Left-side annotation (PolZ) indicates the pollen zones: Pleniglacial, Bølling, and Allerød. Colour intensity represents relative abundance (scaled).

Appendix C: VIF Threshold Justification

We applied a Variance Inflation Factor (VIF) threshold of 10 to select environmental variables for constrained ordination analysis, following standard statistical guidance (Kutner et al. 2005). This threshold balances the need to control multicollinearity while preserving ecologically meaningful predictors in paleoclimate data. Temperature-related variables (GDGT T months above freezing, July temperatures) show elevated VIF values (7.1–9.9) because they reflect inherently correlated seasonal warmth signals and co-vary with orbital insolation cycles; these relationships represent real ecological variability rather than statistical artifacts. A more conservative threshold (VIF < 5) would exclude critical temperature signals, sacrificing ecological interpretability for statistical stringency. Empirical analysis across thresholds (VIF < 3 to 10) confirmed that VIF < 10 optimally balances collinearity control (mean VIF = 6.2) while retaining 9 of 16 environmental variables, enabling robust multivariate inference. The correlation matrix (Fig. B1) demonstrates substantial but expected inter-relationships among climatic predictors, particularly between temperature metrics and insolation variables, validating the appropriateness of the VIF < 10 criterion for this paleoclimate dataset.