Disentangling influences of climate variability and lake-system evolution on climate proxies derived from isoprenoid and branched glycerol dialkyl glycerol tetraethers (GDGTs): the 250 kyr Lake Chala record

High-resolution paleoclimate records from tropical continental settings are greatly needed to advance understanding of global climate dynamics. The International Continental Scientific Drilling Program (ICDP) project DeepCHALLA recovered a 214.8 m long sediment sequence from Lake Chala, a deep and permanently stratified (meromictic) crater lake in eastern equatorial Africa, covering the past ca. 250 000 years (250 kyr) of continuous lacustrine deposition since the earliest phase of lake-basin development. Lipid biomarker analyses on the sediments of Lake Chala can provide quantitative records of past variation in temperature and moisture balance from this poorly documented region. However, the degree to which climate proxies derived from aquatically produced biomarkers are affected by aspects of lake developmental history is rarely considered, even though it may critically influence their ability to consistently register a particular climate variable through time. Modern-system studies of Lake Chala revealed crucial information about the mechanisms underpinning relationships between proxies based on isoprenoid (iso-) and branched (br-) glycerol dialkyl glycerol tetraethers (GDGTs) and the targeted climate variables, but the persistence of these relationships in the past remains unclear. Here we assess the reliability of long-term climate signals registered in the sediments of Lake Chala by comparing downcore variations in GDGT distributions with major phases in lake-system evolution as reflected by independent proxies of lake depth, mixing regime and nutrient dynamics: seismic reflection data, lithology and fossil diatom assemblages. Together, these records suggest that during early lake history (before ca. 180–200 ka) the distinct mixing-related depth zones with which specific GDGT producers are associated in the modern-day lake were not yet formed, likely due to more open lake hydrology and absence of chemical water-column stratification. Consequently absolute GDGT concentrations dating to this period are relatively low, proxies sensitive to water-column stratification (e.g., branched versus isoprenoid tetraether (BIT) index) display highly irregular temporal variability, and correlations between proxies are dissimilar to expectations based on modern-system understanding. A sequence of lake-system changes between ca. 180–200 and ca. 80 ka first established and then strengthened the chemical density gradient, promoting meromictic conditions despite the overall decrease in lake depth due to the basin gradually being filled up with sediments. From ca. 180 ka onward some GDGTs and derived proxies (e.g., crenarchaeol concentration, BIT index and IR_6Me) display strong ∼ 23 kyr periodicity, likely reflecting the predominantly precession-driven insolation forcing of Quaternary climate variability in low-latitude regions. Our results suggest that GDGT-based temperature and moisture-balance proxies in Lake Chala sediments reflect the climate history of eastern equatorial Africa from at least ca. 160 ka onwards, i.e., covering the complete last glacial–interglacial cycle and the penultimate glacial maximum. This work confirms the potential of lacustrine GDGTs for elucidating the climate history of tropical regions at Quaternary timescales, provided they are applied to suitably high-quality sediment archives. Additionally, their interpretation should incorporate a broader understanding of the extent to which lake-system evolution limits the extrapolation back in time of proxy-climate relationships established in the modern system.

1 Introduction

Reliable methods to reconstruct past climate variability during both glacial and interglacial phases of the Quaternary are needed to more precisely model global climate dynamics and to correctly project future climate change due to anthropogenic global warming at the regional scale. However, whereas mid- to high-latitude continental regions are well represented by suitably long and high-quality climate reconstructions (e.g., Petit et al., 1999; S. Barker et al., 2011; Melles et al., 2012; Cheng et al., 2016; Wagner et al., 2019), long records of past climate variability from the tropics are still scarce. On the African continent, the largest landmass straddling the Equator, sediments accumulating without interruption in long-lived lakes are the principal natural archive of Quaternary climate history across its full range of variability (e.g., Lake Malawi: Scholz et al., 2007; Cohen et al., 2007; Stone et al., 2011; Johnson et al., 2016; Lake Tanganyika: Scholz et al., 2007; Tierney et al., 2008, 2010c; Lake Bosumtwi: Scholz et al., 2007; Miller et al., 2016) and provide a diverse array of climate information registered through climate-controlled biological and geological processes in the lake and surrounding terrestrial environments (Cohen, 2003).

An increasingly important source of information on past climate change extracted from lake sediments, both in Africa and elsewhere, is isoprenoid (iso-) and branched (br-) glycerol dialkyl glycerol tetraethers (GDGTs), membrane lipids produced by species of archaea and bacteria, respectively. These organic biomarkers are useful for paleoclimate reconstruction owing to their ubiquitous presence in natural settings, resilience to degradation, and strong response to environmental parameters such as temperature and pH (Schouten et al., 2013). IsoGDGTs consist of two ether-bound C₄₀ isoprenoid alkyl chains that can have varying numbers (0 to 8) of cyclopentyl moieties (i.e., isoGDGT-0 to isoGDGT-8; see GDGT molecular structures in Fig. S1 in the Supplement; De Rosa and Gambacorta, 1988). Crenarchaeol (as well as its isomer cren^′), an isoGDGT with 4 cyclopentyl moieties and 1 cyclohexyl moiety (Sinninghe Damsté et al., 2002; Holzheimer et al., 2021), is only known to be produced by chemolithotrophic, ammonia-oxidizing Thaumarchaeota (e.g., Sinninghe Damsté et al., 2002, 2018; Schouten et al., 2013; Elling et al., 2017; Bale et al., 2019). By contrast, isoGDGT-0 is synthesized by Thaumarchaeota (e.g., Sinninghe Damsté et al., 2012b; Schouten et al., 2013; Elling et al., 2017; Bale et al., 2019) as well as anaerobic methane-oxidizing archaea (e.g., Pancost et al., 2001; Schouten et al., 2001) and methanogenic Euryarchaeota (Schouten et al., 2013, and references therein), and isoGDGT-1 to isoGDGT-3 have been demonstrated to occur in Thaum-, Eury- and Crenarchaeota (Schouten et al., 2013, and references therein). Empirical observations from marine surface sediments suggesting that ring formation in isoGDGTs is controlled by (sub-)surface temperature led to the development of the paleothermometer of the TetraEther indeX of 86 carbon atoms (TEX₈₆; Table 1) to reconstruct past sea surface temperature (SST; Schouten et al., 2002; Kim et al., 2010). TEX₈₆ has also been used to reconstruct lake surface temperature (LST) from isoGDGTs in the sediments of mainly large lakes (e.g., Powers et al., 2005, 2011; Tierney et al., 2008, 2010a; Woltering et al., 2011; Blaga et al., 2013; Sun et al., 2020). However, the use of TEX₈₆ in lakes may be complicated by contributions of isoGDGTs from methanotrophs, methanogens and other archaea. Moreover, the position of the oxycline in the water column appears to strongly influence the niche available to Thaumarchaeota and hence the in situ TEX₈₆ signal (e.g., Zhang et al., 2016; Cao et al., 2020; Baxter et al., 2021; Sinninghe Damsté et al., 2022). The strong influence of lake size and depth on oxycline formation may also imply that small and shallow lakes are less suited for application of the TEX₈₆ proxy (Powers et al., 2010; Baxter et al., 2021; Sinninghe Damsté et al., 2022).

Schouten et al. (2002)Hopmans et al. (2004)Baxter et al. (2021)Sinninghe Damsté et al. (2012a)De Jonge et al. (2015)De Jonge et al. (2014)De Jonge et al. (2014)Sinninghe Damsté (2016)Baxter et al. (2019)Pearson et al. (2011)

Table 1 Formulas of GDGT-based environmental proxies employed in this study, with 6-Me brGDGTs indicated by a prime symbol. GDGT identifiers within square brackets refer to fractional abundances within the respective group (iso- or brGDGTs), and the others refer to absolute abundances (i.e., integrated peak area).

BrGDGTs contain two linear C₂₈ alkyl chains methylated at C-13 and C-16 that most likely formed from the tail-to-tail linkage of two C₁₅ iso fatty acids (Sinninghe Damsté et al., 2000; Fig. S1). The basic tetramethylated brGDGT is usually accompanied by penta- or hexamethylated forms, where additional methyl group(s) occur at the C-5 (Sinninghe Damsté et al., 2000; Weijers et al., 2006a) or C-6 (De Jonge et al., 2013, 2014) position. Cyclic brGDGTs contain 1–2 cyclopentane moieties, formed by cyclization involving the mid-chain methyl groups (Weijers et al., 2006a). Combined lipid–16S rRNA and culture studies suggest that brGDGTs are likely produced by Acidobacteria (Weijers et al., 2009; Sinninghe Damsté et al., 2011, 2014, 2018; Chen et al., 2022; Halamka et al., 2023), although other bacterial phyla are likely also capable of producing these lipids (e.g., Sinninghe Damsté et al., 2011, 2018; Weber et al., 2018; De Jonge et al., 2019; van Bree et al., 2020; Sahonero-Canavesi et al., 2022; Halamka et al., 2023).

As brGDGTs are particularly abundant in soils (Weijers et al., 2006b), their abundance relative to that of aquatically produced crenarchaeol, quantified in the branched versus isoprenoid tetraether (BIT) index (Hopmans et al., 2004; Table 1), was initially used to determine the input of soil material to coastal marine sediments as well as lakes (Sinninghe Damsté et al., 2009; Blaga et al., 2009). However, it is now established that also the production of brGDGTs within lakes is significant, sometimes even dominant (Tierney and Russell, 2009; Sinninghe Damsté et al., 2009; Woltering et al., 2012; Weber et al., 2015, 2018; van Bree et al., 2020). In Lake Chala and possibly also other deep stratifying lakes, the BIT index may instead reflect long-term changes in lake depth and then serve as a proxy for climatic moisture balance rather than rainfall amount per se (Baxter et al., 2021).

Besides the BIT index, several other climate proxies based on brGDGT distributions have been developed. In nearly all studied settings the distribution of these lipids displays a strong correlation to temperature, which is best reflected in the degree of methylation of 5-Me brGDGTs, as captured by the methylation of branched tetraether (${\text{MBT}}_{\text{5Me}}^{\prime }$) index (Table 1; De Jonge et al., 2014). However, despite strong correlation between ${\text{MBT}}_{\text{5Me}}^{\prime }$ in lacustrine surface sediments and temperature (Russell et al., 2018; Martínez-Sosa et al., 2021; Raberg et al., 2021), only a few downcore applications of lake-based temperature calibrations have proved successful (Feakins et al., 2019; Stockhecke et al., 2021; Zhao et al., 2021; Zhang et al., 2021; Garelick et al., 2021; Ramos-Roman et al., 2022; Parish et al., 2023), partly due to continued uncertainty about the exact source(s) of brGDGTs in lakes. In addition, several other environmental factors may influence brGDGT distributions in lakes, such as the lake's depth (Tierney et al., 2010b; Loomis et al., 2014a), trophic level (Loomis et al., 2014b; Martínez-Sosa and Tierney, 2019), conductivity (Shanahan et al., 2013; Raberg et al., 2021), dissolved oxygen content (Loomis et al., 2014a, b; Martínez-Sosa and Tierney, 2019; van Bree et al., 2020; Yao et al., 2020; Wu et al., 2021) and redox conditions. The latter, in particular, has been shown to substantially influence the concentration and distribution of brGDGTs in lacustrine sediments (Loomis et al., 2014a; Wu et al., 2021), as well as their distribution within the water column of stratifying lakes (Weber et al., 2018; van Bree et al., 2020; Yao et al., 2020). As a result, some studies use modified GDGT indices when these appear better suited to the particular study site or reconstruction (e.g., Bittner et al., 2022; Baxter et al., 2023).

Given that GDGT-based climate reconstructions from lake sediment records are based on space-for-time substitution of empirical proxy-climate relationships among a suitably large number of present-day lakes situated along regional- to global-scale climate gradients, the investigation of the specific lake system and depositional environment “hosting” the reconstruction is crucial to identify the influence of confounding factors on the exact relationship between specific GDGTs and temperature or moisture balance. However, despite the substantial effort involved in such modern-system studies, monitoring of proxy variation during multiple seasons or even multiple years does not necessarily suffice to explain proxy variation at the much longer timescale of climate reconstruction from sedimentary GDGT distributions. Long-lived lakes are by nature dynamic systems experiencing large-scale physical, chemical and biological changes related to the geological and tectonic evolution of the lake basin since its formation; its gradual infilling with sediments; and changes in the lake's hydrology or connection to a regional hydrographic network. The influence of these long-term lake-basin changes on local aquatic microbial communities may significantly impact the reliability of GDGT-based climate proxies, but regrettably climate reconstructions rarely take this important source of uncertainty into consideration.

The present study aims to address this issue in the context of a GDGT-based climate reconstruction at Lake Chala in eastern equatorial Africa, a 90 m deep volcanic crater lake from where the International Continental Scientific Drilling Program (ICDP) project DeepCHALLA recovered a 214.8 m long sediment sequence covering ca. 250 kyr of continuous lacustrine deposition since shortly after lake-basin formation. This is done through detailed examination of the concentrations and distributions of isoGDGTs, brGDGTs and associated proxies in relation to major phases in the history of the lake's limnology, ecology and sedimentation dynamics as revealed by independent paleoenvironmental proxies derived from seismic reflection data, sediment lithology and fossil diatom assemblages. Lake Chala may provide a particularly valuable natural archive of regional climate history because its large relative depth, characteristic for crater lakes, has promoted the formation and persistence of a permanently stratified, oxygen-deprived lower water column allowing continuous undisturbed deposition of finely laminated sediments that are often rich in organic matter (Verschuren, 2003; Zolitschka, 2007). Situated in steep-sided basins, crater lakes also have a restricted catchment area lacking distinct stream inflows, so their hydrological setting is relatively simple, and past changes in lake water budget can be expected to be tied strongly to changes in the climate-controlled balance between precipitation and evaporation (e.g., Jones et al., 2001). Our interpretation of sedimentary GDGT data from Lake Chala also builds on a good understanding of GDGT proxy-climate relationships informed by diverse studies of the modern-day lake system and less-ancient sediment records from this location (Sinninghe Damsté et al., 2009; Buckles et al., 2013, 2014, 2016; van Bree et al., 2020; Baxter et al., 2021). This integrated analysis allows us to disentangle the influences of lake-basin development and climate variability on GDGT-based climate proxies extracted from the long and continuous sediment archive of Lake Chala, which may provide a unique view of climate and landscape history in eastern equatorial Africa spanning two complete glacial–interglacial cycles.

2 The modern system and history of Lake Chala

2.1 Setting of the study site

Lake Chala (3°19^′ S, 37°42^′ E) is a relatively large (4.2 km²), deep (ca. 90 m in 2016) volcanic crater lake situated at ca. 880 m above sea level in the southeastern foothills of Mt. Kilimanjaro. Lake-surface evaporation (1735 mm yr⁻¹; Payne, 1970) greatly exceeds mean annual rainfall (565 mm yr⁻¹; De Wispelaere et al., 2017; Griepentrog et al., 2019). Therefore, besides rainfall on the lake and on the steep inner slopes of the crater basin, substantial subsurface inflow is required to balance the lake's water budget (Payne, 1970). This subsurface inflow is derived from rainfall on the forested and subalpine zones of Mt. Kilimanjaro (Hemp, 2006; Bodé et al., 2020) that reaches the lake 3–4 months later (P. A. Barker et al., 2011). The modern-day Lake Chala is a fresh, slightly alkaline (surface-water pH 8.4–9.3) and unproductive tropical lake with high concentrations of silica but low concentrations of phosphorus and nitrogen in the mixed surface layer, although these nutrients accumulate in the hypolimnion (Wolff et al., 2014). The lake has a typical crater-lake morphology, with steep crater walls up to 170 m above the lake's surface and steep underwater slopes down to ∼ 60–70 m which level off to form a flat central lake bottom (Moernaut et al., 2010) and a total catchment area of 5.6 km², only 30 % larger than the surface area of the lake itself (Fig. 1). From about 10 m above the 2016 lake level, i.e., the upper limit of shallow caves formed by wave erosion during past high stands, more significant outflow is possible through the porous upper crater walls.

Figure 1 Bathymetry of Lake Chala (Moernaut et al., 2010; depth contours in m) surrounded by its steep-sided crater catchment, demarcated by the solid bold line. The 2005 CHALLACEA (25 ka to present; grey triangle) and 2016 DeepCHALLA (ca. 250 ka to present; yellow diamond) drilling sites are indicated, as well as the fixed sampling locations of suspended particulate matter (SPM; dark grey square) and settling particles (dark grey circle) in long-term monitoring studies of the modern system (Sinninghe Damsté et al., 2009; Buckles et al., 2014; van Bree et al., 2020; Baxter et al., 2021).

In the semi-arid tropical climate regime characterizing eastern equatorial Africa, mean monthly air temperatures are highest in February–March (night and daytime temperature of 21 and 33 °C) and lowest in July–August (18 and 28 °C; Buckles et al., 2014). Orographically isolated from moisture sourced from the Atlantic or Congo Basin (Sepulchre et al., 2006), Lake Chala is situated east of the Congo Air Boundary (CAB) year-round (Verschuren et al., 2009; Tierney et al., 2013) and thus part of the so-called greater Horn of Africa region which depends entirely on relatively modest rainfall from the Indian Ocean (Wainwright et al., 2019). The region's climate is characterized by a strongly bimodal pattern of seasonal rainfall associated with the shifting latitudinal position of the intertropical convergence zone (ITCZ). At the latitude of Lake Chala, rain seasons occur from late October to December (“short rains”) and from March to May (“long rains”), separated by the main dry season from June to September, i.e., during Southern Hemisphere (SH) winter, and a short dry season in January–February.

2.2 Water-column depth zones and mixing regime

Following Buckles et al. (2014), the water column of Lake Chala can be separated into six distinct zones (Fig. 2), differentiated by their frequency of mixing as reflected in physical and chemical properties. Zone 1 represents the daily mixed layer and is fully oxygenated with uniform temperature and pH. Zone 2 is characterized by oxic to sub-oxic conditions and positioned from immediately below the principal thermocline to the oxycline which demarcates its base. Zone 3 is thus anoxic with pH falling to ca. 7.2 and sharply increasing concentration of dissolved methane (Baxter et al., 2021). Zones 1–3 together constitute the mixolimnion, i.e., the portion of the water column that mixes at least once each year (Lewis, 1983; De Crop and Verschuren, 2021). Zones 4–6 together constitute the monimolimnion and are defined by a stable temperature of 22.3 °C and permanent anoxia, except that rare deep-mixing events reaching into Zone 4 may occasionally inject oxygen that is however quickly consumed by bacterial activity. Across Zone 4, pH decreases further to ca. 7.0, and the dissolved-ion concentration (measured as specific conductance) increases with depth from ca. 350 to ca. 450 µS cm⁻¹ (Barker et al., 2013; Wolff et al., 2014), creating a chemical density gradient across Zone 4 which largely prevents temperature-driven convective mixing (and oxygen injection) beyond the Zone 3–4 boundary (De Crop and Verschuren, 2021). Stable pH and dissolved-ion concentrations throughout Zone 5 indicate lack of mixing even on multi-annual timescales. This also applies to Zone 6, but being positioned directly above the profundal lake bottom the local water chemistry and redox conditions are affected by diffusion out of the uncompacted surficial sediments subject to diagenesis.

Figure 2 Schematic representation of the spatial distribution of GDGTs and their producers in the water column of modern-day Lake Chala during the two seasonal extremes of (a) highly stratified (shallow oxycline) conditions associated with the rainy season and (b) enhanced upper-water-column mixing associated with windy and dry conditions, as shown by analysis of SPM and settling particles (Buckles et al., 2013, 2014; van Bree et al., 2020; Baxter et al., 2021). The water column is divided in six distinct mixing zones with different mixing frequencies and chemical properties, following Buckles et al. (2014). Also shown is the location of DeepCHALLA drilling through the complete package of Quaternary sediments underlying Lake Chala, as documented by seismic profiling (Moernaut et al., 2010), and defining five major depositional stages (I–V) in its ca. 250 kyr history (Maitituerdi et al., 2022). See Sect. 2.3 for further description.

Lake Chala is characterized by a strong seasonal mixing pattern relating to the oscillation between windy dry and calm wet seasons, with substantial variability in the expression of these seasons between successive years. From September until May, i.e., the period encompassing the long and short rain seasons and the intervening warm dry season, high lake-surface temperatures and/or lower wind speeds result in reduced mixing of the upper water column, promoting stronger temperature and chemical stratification (Fig. 2a). Except for the daily mixed layer (Zone 1) oxygen renewal is diminished, and due to heterotrophic bacterial activity that is being promoted by high water temperatures, the depth range of sub-oxic transitional conditions (Zone 2) is greatly reduced or even eliminated, thereby shifting the oxycline (top of Zone 3) to a shallower position (∼ 10 to 15 m), with most intensely stratified conditions occurring during SH summer (Wolff et al., 2014; van Bree et al., 2018, 2020). This expansion of the anoxic (but seasonally mixing) Zone 3 increases the overall volume of anoxic water (Zones 3–6) relative to the oxygen-rich surface layers (Zones 1–2).

During the main dry season, lower air temperatures and higher wind speeds cause deep turbulent and convective mixing of the upper water column down to 42–46 m (Wolff et al., 2014; Buckles et al., 2014; van Bree et al., 2020; Fig. 2b). During this deep-mixing period, which normally starts at the end of May and finishes by mid-September (Wolff et al., 2014; van Bree et al., 2018), nutrient-rich deep water is brought up to the nutrient-starved epilimnion, promoting phytoplankton productivity (for example causing a pronounced diatom bloom; Wolff et al., 2014; van Bree et al., 2018). Simultaneous replenishment of the upper water column with oxygen causes a dramatic expansion of Zone 2 and shrinking of Zone 3, such that the oxygenated part of the water column (Zones 1–2) increases relative to the anoxic part (Zones 3–6). A second period of deep mixing interrupts the long period of upper-water-column stratification during the short dry season in January–February, but hampered by high surface-water temperature it is limited to the uppermost 20–25 m (Wolff et al., 2014; van Bree et al., 2018, 2020). To generalize, there are two seasonal extremes of mixing states in Lake Chala, experienced during the mostly calm–wet season and the windy–dry season and associated with shallow or deep oxycline conditions, respectively.

2.3 Depositional history inferred from seismic stratigraphy

Seismic profiling of the Chala crater basin revealed that the lake overlies at least ca. 210 m of near-continuous lacustrine sedimentation (Moernaut et al., 2010). Extrapolation of mean sedimentation rate in the previously cored and dated upper portion of the sequence (Blaauw et al., 2011) to the base of the seismic profiles indicated that the complete sequence may encompass 250 kyr of deposition, i.e., the last two glacial–interglacial cycles (Verschuren et al., 2013). The seismic-stratigraphic sequence is characterized by a sequence of either basin-wide (“draped”) sedimentation deposited under mostly high-lake-level conditions or basin-focused (“ponded”) sedimentation reflecting periods of low lake level (Moernaut et al., 2010), defining five successive phases in the basin's evolution (depositional stages V–I) characterized by pronounced changes in lake depth. Importantly, because (i) ponded sedimentation reflects the greater sediment focusing which occurs when turbulent mixing reaches the profundal lake bottom (Maitituerdi et al., 2022) and (ii) the steep-sided crater basin entails constancy through time in the minimum water depth allowing undisturbed accumulation of soft organic sediments (Håkanson and Jansson, 1983), seismic stratigraphy serves as a lake-level proxy tied to absolute water-column depth rather than to water-column structure or mixing regime. Also, over the complete 250 kyr history of Lake Chala, the surface level of the lake is a reflection of both climate-driven variability in lake water balance and changes in basin hydrology as the crater basin progressively filled with sediments. During Stage I (ca. 248–207 ka), lacustrine sedimentation was limited to a ring-shaped depositional area surrounding the at that time still exposed cones of volcanic tuff in the center of the basin and characterized by thick mass-wasting deposits in the steeply sloping basin periphery (Maitituerdi et al., 2022). Lake depth was probably low, but this could be due either to dry climate conditions or to the still leaky nature of the crater basin during the period shortly after lake formation. Stage II (ca. 207–113 ka) is defined by the complete burial of the central tuff cones and the start of basin-wide sedimentation. During the first half of this stage (ca. 207–147 ka) a gradual transition to more draped sedimentation takes place, indicating that the water column became progressively taller. During the second half (ca. 147–113 ka) Lake Chala may have attained the greatest depth of its entire 250 kyr history, albeit partly due to the still relatively thin underlying sediment package. Stage III (ca. 113–99 ka) represents a period of significantly reduced lake level implying severe climatic drought. During this stage Lake Chala developed the flat central lake floor which still exists today, and consequently from this stage onwards the total depositional area of the crater basin has remained fairly constant through time. Stage IV (ca. 99–19 ka) was a period of mostly high-lake-level conditions, except for a short-lived low stand at ca. 60 ka. Lastly, Stage V (19 ka to present) represents a period of fluctuating lake level, as inferred from a succession of ponded lenses reflecting intermittent sediment focusing sandwiched between units of draped sediments (Maitituerdi et al., 2022).

2.4 Ecological history inferred from fossil diatom assemblages

The diatom community of modern-day Lake Chala consists mainly of the open-water (pelagic) planktonic species Afrocymbella barkeri and Nitzschia fabiennejansseniana. Analysis of fossil diatom assemblages indicate that Lake Chala has been a true freshwater lake throughout its 250 kyr history (Tanttu, 2021). Diatom species associated with littoral (near-shore) and benthic (bottom) habitats are common only during early lake history (ca. 248–221 ka), indicating that Lake Chala was at that time a fairly shallow and well-mixed lake with adequate nutrient availability. After that time fossil diatom assemblages are entirely dominated by N. fabiennejansseniana and other needle-shaped Nitzschia species, indicating a purely pelagic environment with deep water column, weak upper-water-column turbulence and nutrient-starved conditions. A. barkeri makes its first appearance at ca. 144 ka, and from then on the diatom community had a composition similar to that of the modern-day lake, except for variation in the dominance of A. barkeri versus one or more needle-like Nitzschia spp. Dominance of the former is interpreted to reflect episodes with greater mixed-layer turbulence and more efficient nutrient cycling, whereas dominance of the latter reflects a shallower mixed upper layer and therefore more extreme nutrient depletion. To the extent that alternation between these conditions was climate-controlled, A. barkeri dominance during the periods of ca. 108–96 and 28–13 ka may therefore indicate reduced lake depth under a drier climate regime (Tanttu, 2021).

2.5 GDGT biogeochemistry

The distribution and sources of specific GDGTs in Lake Chala, in surrounding catchment soils and in recent sediments has been the focus of extensive modern-system studies, including monthly monitoring of physical limnology, suspended particulate matter (SPM) and settling particles over multiple years, and analysis of GDGT distributions in profundal surface sediments throughout the basin and in the mid-lake sediment record of the last 25 kyr drilled by the CHALLACEA project (Sinninghe Damsté et al., 2009; Buckles et al., 2013, 2014, 2016; van Bree et al., 2020; Baxter et al., 2021). The results of these studies provide important context for interpretation of GDGT distributions in the complete, ca. 250 kyr sediment record recovered by the DeepCHALLA project.

The BIT-index record of Lake Chala over the last 25 kyr showed good agreement with a first-order lake-level reconstruction based on seismic stratigraphy (Verschuren et al., 2009). As soils in the hills surrounding Mt. Kilimanjaro contain high amounts of brGDGTs (Sinninghe Damsté et al., 2008), temporal variation in the sedimentary BIT index was initially interpreted to reflect varying input of soil-derived brGDGTs to the lake due to varying precipitation and consequent soil erosion. Following the discovery that brGDGTs in Lake Chala are abundantly produced within the lake itself (Buckles et al., 2013, 2014), additional modern-system studies aimed to elucidate the exact nature of the relationship between the Chala BIT index and hydroclimate (van Bree et al., 2020; Baxter et al., 2021). These established that brGDGTs are primarily produced in the anoxic zone of the water column (Zones 4–6) and that their depth range follows the seasonal cycle of lake mixing and stratification. Specifically, during the deep-mixing season between June and September they are restricted deeper in the water column, while under conditions of strong upper-water-column stratification during the rest of the year, the oxycline moves upwards, thereby expanding the brGDGT production zone (Fig. 2; van Bree et al., 2020). On the other hand, Group I.1a Thaumarchaeota, which are the main producers of crenarchaeol in Lake Chala (with secondary contributions from Group I.1b; Buckles et al., 2013; Baxter et al., 2021), are primarily restricted to the (sub-)oxic zone between the principal thermocline and the oxycline (Zone 2), where the degree of sunlight is much less intense than in the uppermost layer and where deep-mixing events provide nitrogen in the form of ammonium (Buckles et al., 2013; Baxter et al., 2021; Fig. 2a). During periods of prolonged shallow oxycline conditions, the depth niche of Group I.1a Thaumarchaeota (Zone 2) is eliminated, and their annual “bloom” is suppressed (Buckles et al., 2013, 2014; Baxter et al., 2021; Fig. 2b). SPM sampled throughout the water column during a year when exceptionally shallow oxycline conditions prevailed yielded only gene copies of Group I.1b Thaumarchaeota, while crenarchaeol concentrations were several orders of magnitude lower than recorded previously (Buckles et al., 2016; Baxter et al., 2021). With substantially lower amounts of crenarchaeol settling on the lake bottom during such intervals, the accumulating sediments attain higher BIT-index values (Baxter et al., 2021). During periods of sustained deep mixing, the reverse situation of prolonged oxygenation of the upper water column (Zones 1–2) promotes development of Thaumarchaeota, thus increasing crenarchaeol production and lowering BIT-index values. Therefore, the BIT index tracks changes in the relative size of the anoxic and oxygenated zones in the water column. Within a single year, the oxycline position is controlled by the timing and duration of seasonal deep mixing. Heavy rainfall and low wind speeds when the ITCZ is overhead cause stratification of the upper water column and hence shallow oxycline conditions, while high wind speeds when the ITCZ is located north or south of the Lake Chala region promote deep mixing, which pushes the oxycline down (van Bree et al., 2020; Baxter et al., 2021). On the long timescales of paleoclimate reconstruction, the relative proportion of the anoxic and oxic zones is also influenced by changes in overall lake depth. High-stand episodes of greater lake depth will be associated with an overall taller anoxic zone, whereas during low stands the anoxic zone will shrink, increasing the relative volume of the upper mixed layer (Verschuren, 1999, 2001). Hence, on the long timescales registered in Lake Chala sediments, the BIT index reflects temporal variation in hydrological moisture balance (Baxter et al., 2021), which in the rather dry tropical region of eastern equatorial Africa is chiefly determined by changes in the strength of the Indian Monsoon and temperature effects on continental evaporation (Baxter et al., 2021, 2023).

The ephemeral nature of Zone 2 in the water column of Lake Chala, where Group I.1a Thaumarchaeota are most abundant, has a major influence on sedimentary proxies based on isoGDGTs (Baxter et al., 2021). During periods of exceptionally shallow oxycline and thus Thaumarchaeotal bloom suppression, greater contributions to the isoGDGT pool from methanotrophs, methanogens and other anaerobic archaea render TEX₈₆-based temperature estimates untrustworthy (Sinninghe Damsté et al., 2012a; Baxter et al., 2021). Hence, in line with the results of other modern-system studies (e.g., Zhang et al., 2016; Cao et al., 2020; Dang et al., 2016; Sinninghe Damsté et al., 2022), temporal variation in water-column stratification is a crucial factor controlling the sedimentary TEX₈₆ signal in Lake Chala, potentially equally important as temperature variation itself (Baxter et al., 2021). As methanogens produce relatively high amounts of isoGDGT-0, the ratio between isoGDGT-0 and crenarchaeol (isoGDGT-0 / cren) can be used to assess the contribution of methanogens to the sedimentary isoGDGT pool (e.g., Blaga et al., 2009; Bechtel et al., 2010). Like brGDGTs, isoGDGT-0 is produced most abundantly in the anoxic lower water column of Lake Chala (Buckles et al., 2013; Baxter et al., 2021). Therefore, the isoGDGT-0 / cren ratio reflects changes in the relative volume of the anoxic and oxic portions of the water column and is relatively higher during highly stratified lake conditions and lower during periods of deep mixing (Fig. 2). A greater relative abundance of the crenarchaeol isomer may reflect periods during which the contribution of Group I.1b Thaumarchaeota to the isoGDGT pool is increased, as these archaea produce a greater amount (typically 14 %–29 %) of the isomer than Group I.1a Thaumarchaeota (typically only 0 %–3 %; Pitcher et al., 2010, 2011; Kim et al., 2012; Sinninghe Damsté et al., 2012b; Elling et al., 2017; Bale et al., 2019). Notably, Group I.1b Thaumarchaeota do not produce isoGDGTs with the same temperature dependency of ring formation as Group I.1a (e.g., Elling et al., 2017). In Lake Chala, both crenarchaeol and its isomer most abundantly occur in Zones 1–2 (Fig. 2). In the 98-month data set of settling particles, higher than average f[CREN^′] values (a measure of the contribution of the crenarchaeol isomer; Table 1) were recorded during a period in 2013 when only Group I.1b Thaumarchaeota gene copies were detected (Baxter et al., 2021). Hence, it appears that the two groups of Thaumarchaeota are differently impacted by lake stratification, with Group I.1a being severely diminished when the (sub-)oxic Zone 2 is eliminated or reduced (Fig. 2b). The f[CREN^′] proxy can therefore be used as an indicator of prolonged shallow oxycline conditions in Lake Chala (Baxter et al., 2021).

The 5-Me and 6-Me isomers of brGDGTs, from their side, occupy spatially distinct zones within the anoxic lower water column of Lake Chala (van Bree et al., 2020): the 5-Me brGDGTs are produced mainly in the anoxic but seasonally variable Zone 3, whereas 6-Me brGDGTs are produced most abundantly in the equally anoxic but permanently stratified Zones 4–6 (Fig. 2). During periods of deep lake mixing Zone 3 is greatly reduced, hence limiting the growth of 5-Me brGDGT producers and increasing the relative contribution of 6-Me brGDGT producers. The isomer ratio (IR_6Me; Table 1) captures this relative contribution of 6-Me to 5-Me brGDGTs. Indeed, low IR_6Me values in settling particle data (van Bree et al., 2020) correspond to trends in other proxies (BIT index, isoGDGT-0 / cren) indicative of an unusually shallow oxycline (Baxter et al., 2021). Like in the case of the BIT index, on the long timescales registered by sedimentary records, changes in the IR_6Me ratio will be predominantly controlled by changes in lake depth because the increased inputs of fresh water which cause the lake level to rise lead to the expansion of Zone 3, where 5-Me brGDGT producers proliferate, hence causing lower IR_6Me values in the sediment (Baxter et al., 2023).

3 Materials and methods

3.1 Construction of the sediment sequence, lithofacies description and age model

The 214.8 m long sediment sequence from Lake Chala was recovered from five drill holes (A–E) at a single location in the eastern depocenter of Lake Chala (Fig. 1), with overlapping 3 m sections achieving complete (100 %) recovery in the upper 123 m (ca. 160 ka to present) of the sediment sequence and near-complete (∼ 85 %) recovery of the lower 92 m. Core splitting, imaging and preliminary lithological description were carried out at the U. S. National Lacustrine Core Facility (LacCore) in Minneapolis (USA), as described by Baxter et al. (2023) with respect to the upper part of the sediment sequence. The entire drilled sequence consists of fine-grained and diatom-rich organic muds with visually clear millimeter-scale lamination or centimeter-scale banding, allowing cross-correlation of core sections from different drill holes with millimeter-scale precision. Before extraction of sediment samples for organic biomarkers and other proxy analyses, all event deposits (turbidites) with thickness > 2 cm were excluded from the composite depth scale to obtain a provisional “event-free” depth scale and to ensure that samples reflect so-called “matrix” sediments of primary lacustrine deposition at the drill site. Samples were extracted from the work halves of core sections at predetermined constant depth intervals on the event-free depth scale such that proxy time series have a more or less uniform temporal resolution throughout the sediment record. Absolute dating efforts of the DeepCHALLA sequence are ongoing. Considering the focus of the present study on long-term lake-basin development rather than paleoclimate reconstruction, we use a preliminary sediment chronology based on the transfer of the high-resolution ¹⁴C-based age model for the last 25 kyr at the CHALLACEA site (Blaauw et al., 2011) to the DeepCHALLA site; links between the seismic stratigraphy of Chala basin at both sites and known near-global climate events back to 140 ka (Moernaut et al., 2010); and linear extrapolation of the average sedimentation rate over this 140 kyr interval to the base of the DeepCHALLA core sequence at 214.8 m below the lake floor (Martin-Jones et al., 2020).

During sampling of the DeepCHALLA sequence in June 2017 the general appearance of the sediment in each 2 cm thick sampled depth interval was noted, with reference to the preliminary lithological description done at LacCore. Matrix sediments were classified into one of two primary lithofacies, namely millimeter-scale (varve-like) laminations and centimeter-scale (banded) sediments. A third lithofacies type is used to describe core sections where these two facies alternate frequently (Baxter et al., 2023). Millimeter-scale laminated sediments are interpreted to have been deposited under stable stratification and a permanently anoxic lower water column as exists today, whereas centimeter-scale banding reflects post-depositional disturbance of the uppermost few centimeters of originally millimeter-scale laminated muds due to bottom currents associated with occasional complete water-column mixing. Although such events may have injected some oxygen to the near-bottom environment, almost certainly this must have been consumed rapidly (in days rather than weeks) by bacterial activity (Lewis, 1987; De Crop and Verschuren, 2019) such that for all practical purposes the lower water column would still have been permanently anoxic. Nevertheless, sections of millimeter-scale lamination can be considered to represent periods of stable meromixis, whereas centimeter-banded sections represent periods when complete water-column mixing occurred at least occasionally on a decadal timescale.

3.2 GDGT analysis and calculation of derived proxies

In the present study a total of 949 sediment horizons from throughout the DeepCHALLA sequence are analyzed, each 2 cm thick and extracted at a regular interval of 16 cm in matrix sediments (i.e., skipping turbidites > 2 cm thick). Detailed inventory of all turbidites (Swai, 2018) revealed that 73 sediment horizons extracted and analyzed for GDGTs (7.5 % of the total) partly consist of thin turbidites (< 2 cm thick). Here only the 33 samples containing > 25 % of turbidite material (3.3 % of the total) were excluded from the final proxy time series, which hence consist of 916 sediment horizons spanning the past ca. 250 kyr. Methods of sample preparation and GDGT analysis have been described previously (Baxter et al., 2023). In short, freeze-dried and powdered sediments (0.3–1.2 g dry weight) were extracted with a Dionex accelerated solvent extraction (ASE) system using a 9:1 $v/v$ mixture of dichloromethane (DCM) and methanol, and 1 µg of internal standard (synthetic C46 glycerol trialkyl glycerol tetraether; GTGT) was added to the total lipid abstract (TLE; Huguet et al., 2006). TLEs were dissolved in DCM, passed through a Na₂SO₄ column and dried under N₂ gas before being separated into apolar, ketone and polar fractions using eluents of hexane$/$DCM (9:1, $v/v$), hexane$/$DCM (1:1, $v/v$) and DCM$/$methanol (1:1, $v/v$), respectively, and passing through an Al₂O₃ column. The fractions were dried under N₂ gas and the polar fractions, containing the GDGTs, were redissolved in hexane$/$isopropanol (99:1, $v/v$) prior to being filtered using a PTFE 0.45 µm filter. Measurement of GDGTs was carried out using an Agilent 1260 Infinity ultrahigh-performance liquid chromatography (UHPLC) system coupled to an Agilent 6130 single quadrupole mass detector following the method of Hopmans et al. (2016). GDGTs were identified by $[\mathrm{M}+\mathrm{H}{]}^{+}$ ion detection in selected ion monitoring (SIM) mode for $m/z$ 1018.0, 1020.0, 1022.0, 1032.0, 1034.0, 1036.0, 1046.0, 1048.0 and 1050.0 (brGDGTs); $m/z$ 1292.3, 1294.3, 1296.3, 1298.3, 1300.3 and 1302.3 (isoGDGTs); and $m/z$ 743.8 (internal standard) with a mass window of 1.0. Peak area integration of the peaks representing GDGTs in the $[\mathrm{M}+\mathrm{H}{]}^{+}$ mass chromatograms was done using Agilent MassHunter software. A peak area of 3×10³ units was used as the detection threshold, with peaks below this threshold being excluded for proxy calculation.

Absolute concentrations of isoGDGTs and brGDGTs were normalized to the organic carbon (C_org) content of the sampled intervals and hence expressed in µg g⁻¹ C_org. For C_org determination, ca. 1.0 g wet sediment from the same 2 cm core increments as the GDGT samples was weighted immediately after extraction and again after freeze-drying to measure the loss in mass as an estimate of water content (%H₂O). The freeze-dried samples were homogenized, split in two and analyzed using a Primacs carbon analyzer at the University of Haifa (Israel), which determines total carbon content by combusting the sample at 1050 °C and measuring the evolved carbon dioxide. One subsample was first treated with phosphoric acid to remove any inorganic carbon present and thus measure C_org only (Maitituerdi, 2023). The GDGT concentration time series comprises 909 sediment horizons due to missing dry sample weight or %C_org values for seven GDGT samples. The TEX₈₆ index was calculated according to Schouten et al. (2002) (Table 1). The BIT index was calculated according to Hopmans et al. (2004), modified to explicitly show the inclusion of both the 5- and 6-Me brGDGTs (De Jonge et al., 2014). IsoGDGT-0 / crenarchaeol, f[CREN]^′ and %isoGDGT-2 were calculated to investigate the producers contributing to the isoGDGT pool (Sinninghe Damsté et al., 2012a; Baxter et al., 2021). In addition, the relative abundance of 6-Me versus 5-Me isomers (IR_6Me) was calculated according to De Jonge et al. (2015), the methylation of 5-Me branched tetraether index (${\text{MBT}}_{\text{5Me}}^{\prime }$) and cyclization of branched tetraether index (CBT^′) according to De Jonge et al. (2014), and the degree of cyclization (DC) of brGDGTs according to Sinninghe Damsté (2016) and Baxter et al. (2021). For temperature reconstruction we applied the global lake calibration of Pearson et al. (2011), which represents mean summer temperature (MST) and was determined by Baxter et al. (2023) to be the brGDGT-based paleotemperature inference model best suited to the setting of Lake Chala at the intended timescale. In Table 1, the original calibration is rewritten to highlight the inclusion of both 5-Me and 6-Me isomers, with GDGTs in square brackets referring to the fractional abundances. The resulting 250 kyr MST record was then rescaled to the mean temperature range of an ensemble temperature reconstruction for the last 25 kyr based on seven independent GDGT-based temperature reconstructions from other East African lakes, according to Baxter et al. (2023).

3.3 Numerical, statistical and periodicity analysis

The relationships between temporal variation in individual GDGT compounds and selected proxies throughout the DeepCHALLA sediment sequence were explored using univariate and multivariate analyses, mostly in the R statistical package FactoMineR (Lê et al., 2008). Univariate analyses were performed on organic-matter normalized absolute concentrations, whereas multivariate principal component analyses (PCAs) were performed on the fractional abundances of individual GDGTs, relative to either the full suite of iso- and brGDGTs (22 compounds) or only the brGDGTs (15 compounds), isoGDGTs (7 compounds) or the sub-set of four isoGDGTs used in TEX₈₆ calculation (isoGDGT-1, isoGDGT-2, isoGDGT-3 and cren^′). To assess the influence of lake-basin development and water-column mixing regime on GDGT distributions, and by extension their sensitivity to climate variability, the 916 sediment horizons (or 909, in the case of GDGT concentrations) were grouped according to seismic lake-history stage or lithofacies to explore trends in GDGT concentrations, PCA scores, proxies and possible relationships to changes in lake properties through time. A Tukey multiple comparison of means with a 95 % family-wise confidence interval was used to test if the GDGT concentrations and proxy values of these groups are significantly different; means are considered significantly different at the 5 % level of significance.

The time series of GDGT concentrations and derived proxies were subjected to periodicity analysis using Acycle 2.3.1 software (Li et al., 2019), after standardization (mean = 0, standard deviation = 1) and third-order detrending. For wavelet (Morlet) analysis, the standardized and detrended GDGT time series were first resampled at the median temporal resolution (between 206 and 200 years depending on the selected time interval). Gaussian bandpass filtering was also performed on select GDGT proxies using bandwidths that targeted periodicities compatible with obliquity (41 kyr) and precession (23 kyr) orbital insolation forcing, as revealed by wavelet analysis.

4 Results

4.1 Trends in isoGDGT and brGDGT concentrations

The overall composition of GDGTs in the drilled Lake Chala sediments is dominated by crenarchaeol and isoGDGT-0, with mean fractional abundances of 0.30 and 0.28, respectively, followed by tetramethylated (0.14), 6-Me (0.10) and 5-Me brGDGTs (0.07). GDGT concentrations display a high degree of high-frequency variability with major swings between adjacent samples (Fig. 3). Smoothing of the time series (black curves indicate a 10-point rolling mean) reveals that substantial changes also occurred with regard to longer-term trends. The concentrations of crenarchaeol (Fig. 3c; range = 0–2740 µg g⁻¹ C_org, average = 201 µg g⁻¹ C_org) and the less abundant isoGDGTs used to calculate TEX₈₆ (“minor” isoGDGTs; Fig. 3d; range = 0.2–900 µg g⁻¹ C_org, average = 62 µg g⁻¹ C_org) are highly correlated (R=0.99, p<0.001; Fig. S3), as they both fluctuate in the same semi-regular pattern throughout the record with the exception of two longer periods within Stage IV (ca. 80–70 and ca. 50–30 ka) during which these compounds are often nearly absent. Concentrations of isoGDGT-0 (Fig. 3e; range = 11–740 µg g⁻¹ C_org, average = 142 g⁻¹ C_org) and the summed brGDGTs (Fig. 3f; range = 6–810 µg g⁻¹ C_org, average = 164 µg g⁻¹ C_org) are also highly correlated with one another (R=0.83, p<0.001; Fig. S3). The concentrations of all GDGTs are notably lower in the lowermost, oldest portion of the record compared to later stages (Fig. 3). For example, average summed concentrations of iso- and brGDGTs during Stage I are, respectively, 149 and 60 µg g⁻¹ C_org, only 27 % of their average concentrations in the full record (405 and 164 µg g⁻¹ C_org). All GDGT concentrations generally increase from the oldest horizon until a pronounced maximum halfway through Stage II (dated to ca. 153 ka) and display notably similar trends until near the end of Stage II when all GDGT concentrations experience a pronounced minimum dated to ca. 125 ka. Using the diatom-based lake phases as reference, correlations between the concentrations of crenarchaeol, minor isoGDGTs (isoGDGT-1, isoGDGT-2, isoGDGT-3 and cren^′), isoGDGT-0 and brGDGTs are universally high in the period before the establishment of the modern diatom community at ca. 144 ka (R=0.69–0.99, p<0.001; Fig. S5), while since then (Fig. S7) the concentrations of crenarchaeol and minor isoGDGTs are strongly correlated (R=0.99, p<0.001) and likewise those of isoGDGT-0 and the brGDGTs (R=0.84, p<0.001) but correlations across these two groups are severely reduced (R=0.25–0.52) and no longer significant (p>0.05), as is also the case in the present-day lake (see Sect. 2.5). The broad interval of ca. 125–75 ka is characterized by lower-than-average isoGDGT-0 and brGDGT concentrations, after which concentrations increase again and generally remain high throughout the upper part of the record.

Figure 3 Downcore profiles of the organic-specific concentration of selected (groups of) GDGTs in the DeepCHALLA sediment sequence in relation to Lake Chala system evolution over the past 250 kyr. Indicated on top are the timing of three major phases in lake ecology as registered in fossil diatom assemblages (Tanttu, 2021), the lithofacies category of each sediment horizon (colored bar) and the depositional stages (V–I) based on seismic stratigraphy (Maitituerdi et al., 2022). Subsequent panels show (a) lake level reconstruction based on seismic stratigraphy (Maitituerdi et al., 2022), (b) sedimentary organic carbon content (%TOC) and TOC-normalized concentrations of (c) crenarchaeol, (d) minor isoGDGTs (isoGDGT-1, isoGDGT-2, isoGDGT-3 and the crenarchaeol isomer), (e) isoGDGT-0, and (f) brGDGTs. Black curves represent the 10-point rolling averages of the full data series, shown in grey. The box plots compare GDGT concentration values in samples grouped by lithofacies category or depositional stage (n = number of samples), with median (black line) and average (solid or open symbol) values, the first and third quantiles (lower and upper hinges), whiskers (thin lines) extending to 1.5 × the interquartile range from the hinges, and data points beyond this range treated as outliers (open black circles). Letters a–d near the box plots indicate statistically significant groups of data determined using Tukey multiple comparison of means with a 95 % family-wise confidence interval. Means that do not share letters are significantly different at the 5 % level of significance. In panels (b–f) the y axis has been stretched to improve readability of the 10-point rolling average trends, but such that many individual data points with very high values, and corresponding box plot outliers, are off-scale. At the bottom the timing of the marine isotope stages (MISs; Lisiecki and Raymo, 2005) is shown for reference.

Samples extracted from millimeter-scale laminated and centimeter-scale banded sediments contain significantly different GDGT concentrations (boxplots in Fig. 3). Millimeter-scale laminated sediments contain higher amounts of isoGDGT-0 (on average 146 µg g⁻¹ C_org) and brGDGTs (206 µg g⁻¹ C_org) and lower amounts of crenarchaeol (120 µg g⁻¹ C_org) and minor isoGDGTs (36 µg g⁻¹ C_org) compared to centimeter-scale banded sediments (averages of, respectively, 112, 103, 319 and 98 µg g⁻¹ C_org, p<0.05 on all differences). Concentrations of these four classes of GDGTs in samples containing a mixture of the two lithofacies are intermediate, with differences also being statistically significant.

4.2 Trends in the distribution of individual GDGTs

Principal component analysis (PCA) was performed on the fractional abundances of all 22 measured GDGTs together and separately on the brGDGTs, isoGDGTs and minor isoGDGTs (Figs. 4 and S8). First considering the full suite of GDGTs, the first principal component (PC1) accounts for 91.3 % of the variance and is chiefly related to the strongly opposing behavior of crenarchaeol versus isoGDGT-0 and the brGDGTs, in particular brGDGTs Ia, IIa^′ and IIa (Fig. 4a). PC2 accounts for only 6.9 % of the variance and mainly separates the isoGDGTs from the brGDGTs. In the PCA of only the isoGDGTs (Fig. 4b), PC1 accounts for a remarkable 99.7 % of the variance, and as seen in the PCA on the full suite of GDGTs, it relates overwhelmingly to the opposing behavior of crenarchaeol and isoGDGT-0. PC2 accounts for only 0.2 % of the variance and mainly separates isoGDGT-1 from isoGDGT-2 and isoGDGT-3. In the PCA of only the minor isoGDGTs contributing to the TEX₈₆ proxy (i.e., isoGDGT-1, isoGDGT-2, isoGDGT-3 and cren^′; Fig. 4d), PC1 accounts for 74.3 % of the variance and is mainly controlled by the different loadings of isoGDGT-1 and isoGDGT-2. PC2 accounts for 24.3 % of the variance, with isoGDGT-3 plotting strongly on the negative side along with the weak negative position of cren^′ versus the positive position of isoGDGT-1 and isoGDGT-3. Finally, in the PCA highlighting the differing distributions of individual brGDGTs (Fig. 4c), PC1 accounts for 65.1 % of the variance, mainly separating the 6-Me brGDGTs (in order of contribution: IIa^′, IIIa^′ and IIb^′), which plot on the positive side, from the acyclic tetramethylated brGDGT Ia and the pentamethylated 5-Me brGDGT IIa, which plot on the negative side. PC2 accounts for 15.4 % of the variance, with as the largest contributors the opposing groups of IIIa^′ and IIa versus Ia and IIb^′.

Figure 4 Principal component analysis (PCA) biplots of the fractional abundances of (a) all GDGTs, (b) isoGDGTs, (c) brGDGTs and (d) minor isoGDGTs (i.e., isoGDGT-1, isoGDGT-2, isoGDGT-3 and the crenarchaeol isomer) in the 250 kyr DeepCHALLA sediment sequence, with the loadings of individual GDGTs contributing at least 1 % to variability among samples on either PC1 or PC2 (or both) shown as arrows. Sediment horizons from depositional stages I–V are presented by different symbols and colored according to lithofacies. Red numerals I–V represent the average PC scores of samples from those respective depositional stages. Note the different scales of PC axes for individual samples (sediment horizons) and variables (GDGTs).

In all PCAs described above the PC scores of samples originating from each of the three lithofacies as defined above are noticeable different (color codes in Fig. 4). Generally, samples from centimeter-scale banded sediments show greater contributions of crenarchaeol, the crenarchaeol isomer, 6-Me brGDGTs and isoGDGT-3, whereas samples from millimeter-scale laminated sediments have a greater contribution of isoGDGT-0 and the 5-Me brGDGTs. Further, in the PCA of the brGDGTs, ∼ 30 samples from centimeter-scale banded and centimeter–millimeter alternating sediments deposited during Stage V are separated from all other samples (Fig. 4c; upper right-hand corner) by unusually high fractional abundances of IIIa^′.

4.3 Trends in GDGT-based proxies

4.3.1 The BIT index

As presented above, PCA on the full suite of GDGTs strongly separates the brGDGTs from crenarchaeol along PC1. Moreover, variation in the BIT index (Fig. 5b) is highly correlated (R=0.99, p<0.001) to the PC1 scores of the respective sediment horizons (Figs. 6a and S3), demonstrating that this proxy captures the most significant variability in GDGT distributions throughout the 250 kyr sediment sequence. Lake Chala sediments cover nearly the full range of possible BIT-index values (range = 0.06–1; Fig. S2) with an overall average value of 0.5. From the oldest sediments up to ca. 160 ka, variation in the BIT index is erratic, switching rapidly between high (> 0.8) and low (< 0.3) values without discernible long-term trends. From ca. 160 ka onwards changes in the BIT index become less erratic, and from ca. 138 to ca. 84 ka sustained periods of very low values (< 0.2) are interrupted by periods with predominantly high BIT values. BIT-index values are consistently low (< 0.2) during ca. 114–96 ka, largely overlapping with Stage III, which is identified as a pronounced lake low-stand episode (Maitituerdi et al., 2022; Fig. 5a and b). The lowest BIT-index values of the entire record (0.06) occur at the very start of this interval. From ca. 80 ka onward, corresponding almost exactly with a sustained lithofacies shift from centimeter-scale banded to millimeter-scale laminated sediments, a major change is observed in the nature of BIT-index variability, with values approaching unity that are sustained for longer periods. BIT-index values are continuously high until ca. 24 ka (roughly at the Stage IV–V transition), except for a ∼ 10 kyr long episode from ca. 62–52 ka and a brief interruption dated to ca. 37 ka. Sustained low BIT-index values between 20 and 14 ka again correspond to a sustained period of centimeter-scale banded sedimentation interpreted as a pronounced lake low stand (Maitituerdi et al., 2022). Thereafter, BIT-index values rise steadily before experiencing a brief reversal to low BIT values at 13–11 ka, followed by a period with values close to unity during 11–9 ka. Over the last 6 kyr, the Chala BIT index has fluctuated around 0.6.

Figure 5 Downcore profiles of selected GDGT-derived proxies in the DeepCHALLA sediment sequence in relation to Lake Chala system evolution over the past 250 kyr. Indicated on top are the timing of three major phases in lake ecology as registered in fossil diatom assemblages (Tanttu, 2021), the lithofacies category of each sediment horizon (colored bar) and the depositional stages (V–I) based on seismic stratigraphy (Maitituerdi et al., 2022). Subsequent panels show (a) lake level reconstruction based on seismic stratigraphy (Maitituerdi et al., 2022), (b) the BIT index, (c) isoGDGT-0 / crenarchaeol ratio, (d) f[CREN^′], (e) %isoGDGT-2, (f) degree of cyclization (DC) and (g) IR_6Me. Data points with a red highlight indicate values of BIT index>0.8, isoGDGT-0 / cren>0.7, $\mathrm{f}\left[{\text{CREN}}^{\prime }\right]>\mathrm{0.04}$ and %isoGDGT-2>45 % (see “Material and methods”). The associated box plots are as described in Fig. 3. Also shown is the timing of the marine isotope stages (MISs; Lisiecki and Raymo, 2005), for reference.

Figure 6 Scatter plots comparing GDGT-based proxies with the PC scores of individual sediment horizons in the 250 kyr DeepCHALLA sequence. (a) BIT index versus PC1 from the PCA on all GDGTs. (b) IR_6Me versus PC1 from the PCA on brGDGTs. (c) ${\text{MBT}}_{\text{5Me}}^{\prime }$ versus PC2 from the PCA on brGDGTs. (d) IsoGDGT-0 / cren versus PC1 from the PCA on isoGDGTs. (e) TEX₈₆ versus PC1 scores from the PCA on minor isoGDGTs (isoGDGT-1 to isoGDGT-3 and cren^′). Sediment horizons from depositional stages I–V are presented by different symbols and colored according to lithofacies. Red numerals I–V represent the average PC scores of samples from those respective depositional stages.

A comparison of the BIT-index time series with those of the absolute concentrations of crenarchaeol and summed brGDGTs reveals that BIT-index variation is predominantly influenced by crenarchaeol abundance and considerably less by the contribution of brGDGTs (Fig. 7a and b), as is also the case in settling particle data (Baxter et al., 2021). However, sediment horizons with BIT-index values < 0.2 generally have brGDGT concentrations < 300 µg g⁻¹ C_org, suggesting that brGDGT variability does affect BIT-index variation to a small extent. Also the average BIT index of the three lithofacies differs significantly (boxplots in Fig. 5b), with centimeter-scale banded sediments most often showing low values (average = 0.2, including a handful of distinct outliers), whereas millimeter-scale laminated sediments have much higher values (average = 0.7), and samples from the millimeter–centimeter-scale alternating lithofacies have intermediate values (average = 0.4).

Figure 7 Scatter plots comparing variation in selected GDGTs and GDGT-based proxies. (a) BIT index versus crenarchaeol concentration, (b) BIT index versus summed brGDGT concentration, (c) IR_6Me versus BIT index, (d) IR_6Me versus summed 5-Me brGDGT concentration, (e) IR_6Me versus summed 6-Me brGDGT concentration, (f) isoGDGT-0 / cren versus BIT index and (g) isoGDGT-0 / cren versus f[CREN]. Sediment horizons from depositional stages I–V are indicated by different symbols and colored according to lithofacies. Red numerals I–V represent the average PC scores of samples from those respective depositional stages.

4.3.2 IsoGDGT-derived proxies

The isoGDGT-derived proxies isoGDGT-0 / cren, f[CREN^′] and %isoGDGT-2 (Table 1) allow investigation of changes through time in Thaumarchaeotal production (Sinninghe Damsté et al., 2012a; Baxter et al., 2021). In the 250 kyr record of Lake Chala (Fig. 5c–e), the isoGDGT-0 / cren ratio is strongly correlated to PC1 of the PCA on all GDGTs (R=0.59, p<0.001; ρ=0.99, p<0.001; Fig. S3) and PC1 of the isoGDGTs (ρ=1.0, p<0.001; Figs. 6d and S3). Hence isoGDGT-0 / cren is also strongly correlated with the BIT index (R=0.61, p<0.001; ρ=0.95, p<0.001; Fig. 7f), as expected because both of them are controlled by crenarchaeol abundance, and isoGDGT-0 and brGDGTs have overlapping ecological niches in the anoxic lower water column of Lake Chala (Fig. 2). IsoGDGT-0 / cren variability is highly erratic in the lower part of the DeepCHALLA sequence (ca. 250–160 ka) and becomes more structured thereafter, mirroring the trends displayed by the BIT index. f[CREN^′] is moderately correlated with PC1 scores of all GDGTs (R=0.51, p<0.001) and of the isoGDGTs separately (R=0.58, p<0.001) and more strongly correlated to isoGDGT-0 / cren (R=0.70, p<0.001; Figs. 7g and S3). Through most of the DeepCHALLA sequence values of f[CREN^′] are low (< 0.025; Fig. 5d), presumably reflecting the often limited presence of Group I.1b Thaumarchaeota (Baxter et al., 2021). The very few instances of higher f[CREN^′] values before 60 ka mostly consist of single sediment horizons. Sustained periods of high f[CREN^′] values (> 0.04) occurred from ca. 60 to 25 ka in Stage IV and from 11 to 9 ka in Stage V, mostly in millimeter-scale laminated sediments. The %isoGDGT-2 proxy (Fig. 5e), used by Sinninghe Damsté et al. (2012a) to determine whether the isoGDGTs in Lake Chala are predominantly produced by Thaumarchaeota, is strongly inversely correlated to PC1 of the minor isoGDGTs ($R=-\mathrm{0.93}$, p<0.001). Before ca. 80 ka its value hovers stably around 30 %, with only two instances of more elevated values. After this time its baseline shifts to higher percentages. From ca. 55 to 24 ka in Stage IV, %isoGDGT-2 is highly variable, frequently exceeding 50 %. Two other periods of high %isoGDGT-2 values occurred between 20–11 and 9–0 ka.

TEX₈₆ is strongly correlated to PC1 of the minor isoGDGTs ($R=-\mathrm{0.95}$, p<0.001; Figs. 6e and S3). From the oldest sediment horizon until ca. 180 ka, the index is relatively stable around 0.45, after which it increases to ∼ 0.55 by ca. 165 ka, remaining at that level until ca. 80 ka (Fig. 8b). The period ca. 80–20 ka is characterized by highly variable TEX₈₆ values ranging between 0.43 and 0.77. At 20 ka the index rises dramatically, reaching peak values of 0.85 at ∼ 5 ka and remaining high until the top of the sequence.

Figure 8 Downcore profiles of selected GDGT-based paleotemperature proxies in the DeepCHALLA sediment sequence in relation to Lake Chala system evolution over the past 250 kyr. Indicated on top are the timing of three major phases in lake ecology as registered in fossil diatom assemblages (Tanttu, 2021), the lithofacies category of each sediment horizon (colored bar) and the depositional stages (V–I) based on seismic stratigraphy (Maitituerdi et al., 2022). Subsequent panels show (a) lake level reconstruction based on seismic stratigraphy (Maitituerdi et al., 2022), (b) TEX₈₆ values with red circles indicating samples excluded from the time series using the filtering method of Baxter et al. (2021), (c) ${\text{MBT}}_{\text{5Me}}^{\prime }$ index, (d) mean daily insolation at the Equator in September (blue line) and March (grey line), and (e) mean summer temperature (MST), calculated according to Pearson et al. (2011) and rescaled using an ensemble reconstruction of post-glacial (20–0 ka) temperature history from seven other East African lakes following Baxter et al. (2023). The associated box plots are as described in Fig. 3. Also shown is the timing of the marine isotope stages (MISs; Lisiecki and Raymo, 2005), for reference.

Millimeter-scale laminated and centimeter-scale banded sediments again differ significantly in average isoGDGT-0 / cren ratio, f[CREN^′] and %isoGDGT-2 values (boxplots in Fig. 5c and d), with millimeter-scale laminated sediments having higher average values (39.6 %, 0.04 % and 37.1 %) than centimeter-scale banded sediments (1.1 %, 0.01 % and 30.2 %), and millimeter–centimeter-scale alternating sediments having values intermediate (isoGDGT-0 / cren ratio, %isoGDGT-2) or near-identical to centimeter-scale banded sediments (f[CREN^′]). In contrast, average TEX₈₆ values of the three lithofacies are not significantly different from one another (boxplots in Fig. 8b).

4.3.3 BrGDGT-derived proxies

The IR_6Me of penta- and hexamethylated brGDGTs varies substantially throughout the DeepCHALLA sequence (0.24–0.89, average = 0.61; Fig. 5g, note the reversed y axis). There is strong correlation between IR_6Me and PC1 of the PCA on brGDGTs (R=0.95, p<0.001; Fig. 6b), meaning that this ratio reflects an important aspect of variability in the distribution of brGDGTs in Lake Chala sediments. IR_6Me values are relatively low and stable in the early part of the time series, becoming more variable from ca. 170 ka onwards. Temporal variability in IR_6Me shows moderately negative correlation with the BIT index ($R=-\mathrm{0.49}$, p<0.001; Fig. 7c) and isoGDGT-0 / cren ($R=-\mathrm{0.33}$, p<0.001). At ca. 142 ka a sharp transition occurs from low to high IR_6Me values, almost coeval with the sustained transition to centimeter-scale banded sedimentation and establishment of the modern-day diatom community. From this time onwards similarity between trends in IR_6Me and BIT index (and also isoGDGT-0 / cren) is enhanced, resulting in stronger negative correlation between IR_6Me and the BIT index after 144 ka ($R=-\mathrm{0.64}$, p<0.001; Fig. S7). Comparing variation in IR_6Me to the concentrations of 5-Me and 6-Me brGDGTs indicates that both groups of brGDGTs have a comparable influence on the IR_6Me signal (Fig. 7d and e). For example, sediments with the highest IR_6Me values (> 0.7) also have the highest concentrations of 6-Me brGDGTs, while those with IR_6Me below 0.7 generally have greater amounts of 5-Me brGDGTs.

The CBT^′ index is partly affected by the degree of cyclization of the brGDGTs, but as shown by its strong positive correlation with IR_6Me throughout the DeepCHALLA sequence (R=0.95, p<0.001; Fig. S3), in Lake Chala this index is mainly controlled by variation in the relative abundance of 5-Me and 6-Me brGDGTs. Accordingly, the CBT^′ index also shows strong correlation with PC1 of the PCA on the brGDGTs, which separates the 5-Me and 6-Me brGDGTs (R=0.99, p<0.001; Fig. S3). On the other hand, the degree of cyclization of the brGDGTs (expressed with the DC ratio) shows only weak positive correlation with IR_6Me (R=0.35, p<0.001; Fig. S3) and is also only weakly correlated with either PC1 and PC2 of brGDGT distribution (R=0.31, p<0.001; $R=-\mathrm{0.30}$, p<0.001). This suggests that DC is not a strong measure of the variance in the distribution of brGDGTs in Lake Chala sediments. Accordingly, its time series does not show large changes, remaining mostly around 0.2, except during the period between ca. 140 and ca. 100 ka and in the last 10 kyr when larger temporal variability is observed (Fig. 5f).

The 250 kyr time series of ${\text{MBT}}_{\text{5Me}}^{\prime }$ is strongly correlated to both PC1 (R=0.66, p<0.001; Figs. 6c and S3) and PC2 ($R=-\mathrm{0.68}$, p<0.001) of the brGDGTs and, hence, naturally also to IR_6Me (R=0.84, p<0.001) and CBT^′ (R=0.68, p<0.001). This index ranges from 0.48–0.85 and besides a sharp rise at the base of the record is relatively stable in sediments deposited before ca. 175 ka (Fig. 8c). Thereafter, ${\text{MBT}}_{\text{5Me}}^{\prime }$ first increases to near-maximum values (∼ 0.8) at ca. 160 ka and then drops abruptly to a sustained minimum (∼ 0.6) dated to ca. 150–140 ka which terminates as sharply as it started. The period from ca. 140 until ca. 80 ka is characterized by variable but generally high ${\text{MBT}}_{\text{5Me}}^{\prime }$ values (frequently above 0.8), again ending in an abrupt drop to values below 0.7. ${\text{MBT}}_{\text{5Me}}^{\prime }$ reaches its lowest values just after ca. 60 ka, then gradually increases until 14 ka when a short-lived spike occurs, lasting until 11 ka and followed by a gently rising then falling trend over the last 10 kyr.

The time series of MST spanning the complete DeepCHALLA sequence is strongly correlated to PC2 of the brGDGTs ($R=-\mathrm{0.65}$, p<0.001, Fig. S3) and to the DC (R=0.76, p<0.001). Rescaled MST values (Baxter et al., 2023) range between 20.3 and 25.6 °C. As reconstructed, rescaled MST is rather stable (generally 22–23 °C) from ca. 250 until 180 ka, at which time it gently increases to ∼ 24 °C then again decreases to a pronounced minimum (∼ 21.5 °C) dated to ca. 145 ka (Fig. 8e). After this, MST attains peak values at ca. 134 and 121–108 ka, separated by a minimum bottoming out at ca. 126 ka. From ca. 108 to 24 ka, rescaled MST generally hovers between 22 and 23.5 °C, with the exception of an inferred cooler episode (∼ 21.5 °C) centered at ca. 60 ka. Sustained low MST values also occur from 21 to 14 ka. Thereafter, MST rises sharply to a mid-Holocene maximum (8.5–4.5 ka) followed by slightly lower values in sediments deposited during the last 4 kyr.

As is the case with most GDGTs, also average values of the proxies IR_6Me, DC, ${\text{MBT}}_{\text{5Me}}^{\prime }$ and rescaled MST (°C) are significantly different between samples extracted from millimeter-scale laminated and centimeter-scale banded sediments, now with centimeter-scale banded sediments producing higher values on average (respectively, 0.71, 0.30, 0.72 and 23.0; see boxplots in Figs. 5 and 7) than millimeter-scale laminated sediments (0.54, 0.32, 0.66 and 22.7).

4.4 Periodicities in GDGT concentrations and proxies

Periodicity analysis on downcore profiles of selected GDGT concentrations and proxies (Figs. 9 and S9–S11) show strong cyclicities likely related to orbital precession (23 kyr) and perhaps obliquity (41 kyr). However, periodicities of 22.2 and 44.4 kyr in the concentration of crenarchaeol become apparent only from ca. 180 and ca. 110 ka onwards (see wavelet analysis; Fig. 9c). The lack (or weakness) of these cycles in the older portion of the sediment record is also evident in the varying amplitude of bandpass filters with those periodicities (Fig. 9b), which is very modest in Stage I sediments and increases in the course of Stage II. Moreover, in the global wavelet spectrum restricted to the period 180–0 ka, the signatures of precession- and obliquity-scale cycles are enhanced compared to that covering the complete record, and the cycle lengths shift closer to the astronomical solutions (Fig. 9d and e). Also the IR_6Me and BIT-index time series display periodicities close to 23 and 41 kyr (Figs. S9 and S10). Although less clear as in the case of crenarchaeol, also the wavelet spectra and bandpass-filtered time series of these two proxies indicate weak expression of the astronomical cycles during Stage I and the first half of Stage II, becoming more pronounced when analysis is restricted to the last 180 kyr.

Figure 9 Periodicity analysis of crenarchaeol concentration in the 250 kyr DeepCHALLA sediment sequence. Indicated on top are the timing of three major phases in lake ecology as registered in fossil diatom assemblages (Tanttu, 2021), the lithofacies category of each sediment horizon (colored bar) and the depositional stages (V–I) based on seismic stratigraphy (Maitituerdi et al., 2022). Subsequent panels show (a) lake level reconstruction based on seismic stratigraphy (Maitituerdi et al., 2022), (b) variation in crenarchaeol concentration in light grey with green and pink curves representing the time series after bandpass filtering with periods similar to those of orbital precession and obliquity, as revealed by wavelet analysis, and (c) visual representation of wavelet analysis using a Morlet function, with warm and cold colors reflecting high and low values of the power spectrum. Also shown is the timing of the marine isotope stages (MISs; Lisiecki and Raymo, 2005), for reference. Panels (d) and (e) summarize the wavelet spectra of crenarchaeol concentration in the complete DeepCHALLA sequence (250–0 ka) and in the section 180–0 ka only to highlight more pronounced precession and obliquity cycles in the latter. The red stippled line in (d) and (e) represents the 95 % confidence interval; dominant frequencies are labeled at the top of the maxima.

5 Discussion

5.1 Influence of lake-basin evolution on niches of GDGT producers

The low concentrations of all GDGTs during the first ca. 50–70 kyr of lacustrine deposition in the Lake Chala basin (Fig. 3) might be interpreted to indicate poorer preservation of the aquatically produced GDGTs due to longer exposure to oxygen while settling through the water column and/or in unconsolidated surficial sediments after settling but prior to permanent burial. However, GDGTs are found abundantly in well-oxygenated environments such as aerated soils and are generally considered to be resilient to degradation at least on the timescale of this study (Schouten et al., 2013). Alternatively, the GDGT concentrations may only appear low because of normalization against higher C_org values. However, the C_org content of Lake Chala sediments deposited during this early period is not systematically higher than in later stages (Fig. 3b). We therefore interpret the low concentrations of all GDGTs in this section of the record (Fig. 3) to reflect lower GDGT production because a shallow and/or well-mixed and oxygenated water column was less favorable to either brGDGT producers or Thaumarchaeota (Fig. 10a). Seismic stratigraphy does indicate that lake depth during Stage I (ca. 248–207 ka) was markedly lower than today (Maitituerdi et al., 2022), and also the common occurrence of benthic diatoms in the first ca. 30 kyr of lake history (Tanttu, 2021) confirms the proximity of a shallow-water habitat to the drill site. The first meaningful rise in GDGT concentrations at ca. 215 ka (Fig. 3) is broadly coeval with the disappearance of benthic diatoms, indicative of a substantial increase in water depth; and with the first consistent deposition of millimeter-scale laminated sediments (Fig. 3), suggesting that this deepening promoted greater persistence of bottom-water anoxia and thus of the sub-oxic and anoxic niches required by aquatic GDGT producers.

Figure 10 Schematic representation of the ecological niches and distribution of GDGTs and their producers in the water column of Lake Chala during major successive stages of lake-basin history over the past 250 kyr, based on GDGT time series in the DeepCHALLA sediment sequence (red line; this study) and inferences from SPM and settling particle data in the modern system (Sinninghe Damsté et al., 2009; Buckles et al., 2014; van Bree et al., 2020; Baxter et al., 2021). The gradual infill of the crater basin with sediments is based on seismic stratigraphy (Maitituerdi et al., 2022), with darker grey shading reflecting low-stand conditions.

Besides lower-than-average GDGT concentrations, several other aspects of GDGT variability were different before ca. 200–180 ka compared to later on. First, concentrations of all GDGTs (crenarchaeol, isoGDGT-0, minor isoGDGTs and summed brGDGTs) are positively correlated with each other (see Fig. S4 in comparison to correlations during 180–0 ka in Fig. S6), in conflict with their niche partitioning in the modern system which predicts distinct temporal trends in the two groups of GDGTs associated with either the upper mixed layer (crenarchaeol and minor isoGDGTs) or the anoxic zones (isoGDGT-0 and brGDGTs). Only when the GDGT time series are limited to the last ca. 144 kyr do correlations between these GDGTs reflect their different association with distinct water-column zones in the modern lake. Notably, the last ca. 144 kyr is the period for which fossil diatom assemblages attest that nutrient budgets and cycling in Lake Chala, and thus presumably also its water-column structure, have resembled the present-day situation. Similar temporal shifts occur in the correlations between GDGT-based proxies. For example, inverse correlation between the BIT index and IR_6Me improves from $R=-\mathrm{0.31}$ (p<0.001) in the section before ca. 144 ka to $R=-\mathrm{0.64}$ (p<0.001) after ca. 144 ka, reflecting the strong connection between them based on the niches of GDGT producers identified in the modern system. Also, temporal variation in some GDGT proxies influenced by water-column mixing and stratification (BIT index, isoGDGT-0 / cren) is highly erratic before ca. 170 ka, and others (IR_6Me, f[CREN^′]) seem largely unresponsive before ca. 200 ka. In summary, these proxies do not appear to provide a cohesive account of changes in lake depth and mixing regime during the early part of lake history. During Stage II (ca. 207–113 ka), seismic stratigraphy indicates that total water depth increased steadily, i.e., the surface level of Lake Chala rose faster than the rate of sediment accumulation filling the basin. However, this lake deepening was accompanied by only a very modest increase in lake surface area due to the near-vertically sloping crater walls (Maitituerdi et al., 2022). Comparative analysis of 60 crater lakes in western Uganda showed that the depth of water-column mixing in these systems is most strongly related to lake surface area, not the height of the crater rim that might be thought to provide wind shelter (De Crop and Verschuren, 2021). Therefore in Lake Chala, where lake surface area is nearly constant across a large range of lake depths (Maitituerdi et al., 2022), changes in lake depth alone cannot be responsible for changes in the depth of the mixed layer. In the absence of > 100 m of sediments deposited later on, the crater basin during Stage II was substantially deeper than today, and the water column of Lake Chala may have attained its overall greatest height ever (> 200 m). Under these conditions, mixing of the entire water column might be expected to have been much less likely than today. However, given the thinner sediment infill at that time, the crater bottom and side walls probably were still fairly porous, allowing greater subsurface outflow and thus removal of the dissolved solids that might otherwise accumulate in the lower water column by decomposition of organic matter or dissolution of photosynthetically precipitated calcite. This more “leaky” nature of the crater basin prevented development of the chemical density gradient which promotes permanent stratification (biogenic meromixis, sensu Hutchinson, 1937) and thus allowed at least occasional mixing of the deeper water column at multi-annual timescales: under the right weather conditions (for example, strong lake-surface cooling during an exceptionally windy dry season), this deep but chemically dilute water column could experience complete turnover, allowing oxygen to penetrate to the profundal lake bottom. Therefore, not only Stage I but also the deep early phase of Stage II were likely characterized by sporadic events of complete mixing, contrasting with the present-day situation of a permanently stratified, and permanently anoxic, lower water column. We suggest that the erratic variation in GDGT-based proxies such as the BIT index, isoGDGT-0 / cren and f[CREN^′] prior to ca. 170 ka may be explained by the occasional occurrence of such deep-mixing events, with intermittent peak values indicating that multi-annual episodes of stratification and deep-water anoxia did occur but inconsistently (Fig. 5). In support of this inference, frequent interruptions of the millimeter-scale laminated sediments by alternating millimeter–centimeter-scale lamination imply that events of complete mixing did indeed occur. Most importantly, the lack of chemical stratification during Stage II meant that meromictic conditions were maintained (at least most of the time) by the temperature gradient alone (thermogenic meromixis, sensu Katsev et al., 2010) and thus that the water column of Lake Chala did not comprise today's six well-defined depth zones (Fig. 10a). Zones 4–5, which at present constitute the permanent anoxic zone and have higher dissolved-ion content and lower pH than mixing Zones 1–3, probably did not yet exist (Fig. 10b). Consequently the depth of oxygen penetration during seasonal deep-mixing events must have been highly variable between years because it was not restricted by a static chemocline (De Crop and Verschuren, 2021). Moreover, due to the more dilute water mass in addition to low biological productivity typical of deep lakes, perhaps also the average depth of visible light and UV penetration was greater than today (Secchi disk depth varies seasonally between 1 and 8 m; van Bree et al., 2018; Fig. 10a and b). As Thaumarchaeota are photosensitive (Merbt et al., 2012; Horak et al., 2018), this may have influenced their ability to grow in the upper water column and, hence, their presence and/or production. Given these marked differences in aquatic environment, it seems unsurprising that distributional relationships among different GDGTs and GDGT-based proxies from before ca. 170 ka do not match those observed in the modern system and younger sediments. As the association of the different GDGTs with particular niches in the modern-day water column of Lake Chala underpins understanding of proxy-climate relationships, simple extrapolation of this understanding to the earliest phases of lake history is unwarranted.

Further support for a diminished influence of climate on the GDGT proxies measured in the older sediments can be gleaned from the periodicity analysis. Namely, the apparent influence of orbital insolation forcing on certain GDGT concentrations (e.g., crenarchaeol) and proxies (e.g., BIT index and IR_6Me) is strong in the younger portion of the time series but weak before ca. 180 ka in the case of precession and before ca. 90 ka in the case of obliquity (Figs. 9, S9 and S10). Precession is known to strongly influence monsoon dynamics and hydroclimate in tropical Africa, as it controls the amount and distribution of solar radiation received there (Singarayer and Burrough, 2015). Due to variation in orbital eccentricity modulating the amplitude of precessional insolation forcing (Scholz et al., 2007; Blome et al., 2012), the prominence of a precession signature in African climate history can be expected to vary through time. However the amplitude of precession was larger, on average, during ca. 250–180 ka (late MIS8 to early MIS6) than during the last glacial period (MIS4–MIS2; Fig. 8d), arguing against this being a feasible explanation for our data. Also, delayed expression of the obliquity cycle relative to precession could partly relate to the polar ice sheets reaching their greatest extent during the last glacial period (MIS4–MIS2), as this enhanced the influence of high-latitude climate dynamics on the tropics (Tjallingii et al., 2008), but this fails to explain the lack of an obliquity signature during the penultimate glacial period (MIS6). Moreover the timing of when an apparent orbital period appears in the GDGT time series is not consistent among proxies and most often began before MIS4 (Figs. 9, S9 and S10). In conclusion, the lack of orbital insolation signatures in GDGT time series from the older part of the sequence confirms that climate variability was probably not the dominant mechanism driving variability in sedimentary GDGT distributions during this period, unlike during more recent lake history.

Progressive accumulation of fine-grained sediments in the Lake Chala basin (in total ∼ 70 m by ca. 180 ka; Maitituerdi et al., 2022; Fig. 10b) gradually diminished the leaky nature of the basin floor such that more dissolved solids were retained, in turn promoting development of the chemical density stratification that could maintain biogenic meromixis. Thus mixing Zones 4 and 5 probably developed during the later phase of Stage II (Fig. 10b). Judging from our GDGT data, the water column of Lake Chala attained its modern-day structure of six mixing zones from ca. 180–170 ka onward. Nevertheless, the establishment of the modern-day diatom community at ca. 144 ka (Tanttu, 2021) implies a dramatic change in open-water nutrient budgets or cycling ∼ 30 kyr later. As the heavily silicified A. barkeri presumably requires more nutrients for population persistence than Nitzschia species (Tanttu, 2021), the establishment of the former taxon in Lake Chala at ca. 144 ka testifies to improved upcycling of nutrients from the hypolimnion. A coincident transition to mainly centimeter-scale banded sedimentation (Fig. 3a) indicates that deep mixing must indeed have occasionally reached the lake bottom, causing modest sediment disturbance. The occurrence of such events would have been promoted if an episode of drier climatic conditions reduced lake depth by lowering its surface level. However, as evidence of lake-level lowering at this time is lacking in seismic stratigraphy (Fig. 3a), we surmise that both the facies transition and the improved nutrient budget may have resulted from shallowing of the water column due to the progressive accumulation of sediments (∼ 100 m by that time; Maitituerdi et al., 2022). The combined evidence suggests that despite development of the six-zone water-column structure, the stability of chemical stratification was still relatively low, so it could be overcome by the deep convective mixing occurring during an exceptionally cool or windy dry season. Another ∼ 30 kyr later, Stage III started with a major drop in lake surface level (Moernaut et al., 2010; Maitituerdi et al., 2022), indicating that water loss by lake-surface evaporation markedly exceeded water inputs from catchment precipitation and subsurface inflow. The consequent reduction in lake volume must have led to an overall increase in dissolved-ion concentration such that by the end of Stage III the resistance to deep mixing approached conditions similar to today (Fig. 10c). Presumably, the transition to a more positive water balance which ended Stage III also strengthened chemical stratification across the Zone 3–4 boundary by adding dilute water to the annually mixed surface layer (Zones 1–3) such that by ca. 80 ka Lake Chala enjoyed stable meromixis, as evidenced by the almost uninterrupted deposition of millimeter-scale laminated sediments throughout the remainder of Stage IV. In summary, the sequence of events which first created the six depth zones in the water column of Lake Chala and then enhanced their distinctness in terms of ambient physical and chemical conditions allowed the niches of GDGT producers to become increasingly relatable to those observed in the modern-day system. Consequently between ca. 180 and ca. 80 ka there is increasing consistency between lake conditions as inferred from the distributions of GDGTs based on modern-system understanding, and the evidence from lithostratigraphy, fossil diatom assemblages and seismic reflection data which reflect the long-term history of changes in, respectively, water-column mixing (lithostratigraphy and diatoms) and lake depth (seismics).

Between ca. 140 and ca. 80 ka, a period when predominantly centimeter-scale banded sediments indicate that permanent stratification (meromixis) in Lake Chala was less stable than before and after, several long sections of low BIT-index and isoGDGT-0 / cren values occur in combination with high IR_6Me values, suggesting a relatively thicker oxygenated layer (Zones 1–2) and greatly reduced or absent anoxic lower mixed layer (Zone 3). This is also indicated by lower-than-average concentrations of the anoxically produced brGDGTs and isoGDGT-0. These data indicate that deep mixing was frequent but most often halted by the chemocline at the top of Zone 4, which in turn suggests that also total lake depth was lower during this broad period than before and after. In particular, the longest period of consistently low BIT-index and isoGDGT-0 / cren values from ca. 113–99 ka matches near-perfectly with the timing of the major Stage III low stand (Figs. 5 and 10d). Concentrations of all GDGTs and derived proxies point to an unprecedented change in the structure of the water column of Lake Chala at ca. 80 ka, when inferred loss of Thaumarchaeota and increasing absolute and relative proportions of all GDGTs produced in anoxic waters suggest a rather abrupt establishment of strongly stratified conditions, as also indicated by the equally abrupt transition from centimeter-scale banded to millimeter-scale laminated sediments. The seismic record suggests that lake level was consistently high throughout Stage IV (Fig. 10e), except for a brief interlude of ponded sedimentation (i.e., an inferred lake low stand) at ca. 62 ka, also registered as a brief interruption of millimeter-scale laminated sedimentation (Fig. 3a). Between ca. 80 and 25 ka, GDGT concentrations and proxies mostly suggest an expanded anoxic zone under tall-water-column conditions and more limited upper-water-column mixing: reduced amounts of crenarchaeol and the minor isoGDGTs relative to isoGDGT-0 and the brGDGTs translate to mostly high BIT-index and isoGDGT-0 / cren values and generally low IR_6Me. These patterns in stratification-associated GDGTs and proxies reversed during a ∼ 10 kyr period dated to between ca. 69 and ca. 60 ka, in agreement with the seismic and lithofacies evidence for a low stand at ca. 62 ka. Also at the end of Stage IV and after transition to Stage V, patterns in these GDGT proxies remain consistent with the seismic and lithofacies data, which indicate another period of reduced lake level and deep mixing between 24 and 14 ka (Fig. 10f), apparently second in amplitude only to the Stage III low stand (Fig. 5). However by the onset of Stage V the additional accumulation of ∼ 55 m of sediments had markedly reduced the magnitude of lake-level drop required to generate bottom disturbance at a comparable position of the lake's surface level. After high-stand conditions between 12 and 9 ka (Fig. 10g), in the last 9000 years Lake Chala has been characterized by a moderately deep water column, with certain GDGT proxies (e.g., the BIT index) suggesting a fluctuating lake level but high stability of water-column stratification inferred from GDGT proxies matching the predominantly millimeter-scale laminated sedimentation.

5.2 Environmental controls on the IR_6Me proxy in the Lake Chala record

In soils and peats (Weijers et al., 2007; Peterse et al., 2012; De Jonge et al., 2014; Xiao et al., 2015; B. Naafs et al., 2017; B. D. Naafs et al., 2017) strong correlations occur between pH and the degree of cyclization of brGDGTs or the relative contribution of 6-Me brGDGTs. Although similar relationships have been reported in surface sediment datasets from lakes (e.g., Raberg et al., 2021; Dang et al., 2016), they are not unambiguous (e.g., Russell et al., 2018; Martínez-Sosa et al., 2021), possibly due to unaccounted variability in pH conditions with depth and season in lakes. In the DeepCHALLA sequence, it appears that the relative contribution of 5-Me and 6-Me isomers has a strong influence on the CBT^′ index and that the DC is not an obviously important measure for understanding variability in brGDGT distributions. Notably, whereas SPM data (van Bree et al., 2020) do not clearly reveal the environmental significance of the DC and CBT^′ indices in modern-day Lake Chala, they revealed that IR_6Me reflects variation in the spatially distinct niches of 5-Me and 6-Me brGDGT producers in the water column (see Sect. 2.5; van Bree et al., 2020). Notably, the pH range in their respective niches (7.2–8 in Zone 3 and 6.9–7.1 in Zones 4–5; Buckles et al., 2014; Baxter et al., 2021) opposes the trend found in soils, where 6-Me isomers dominate in high-pH conditions (De Jonge et al., 2014). This ambiguity suggests that IR_6Me variation in the Lake Chala sediment sequence is likely related to another factor than pH.

Mechanistic understanding of the modern-day Lake Chala system strongly couples variation in IR_6Me on seasonal and longer timescales with BIT-index variation. An increase in IR_6Me can be linked to either enhanced mixing or reduced lake depth shrinking Zone 3, compensated by the relative expansion of Zone 2 where Group I.1a Thaumarchaeota are produced, thus lowering the BIT index. The opposing trends in these two proxies is clearly apparent in the last ca. 160 kyr ($R=-\mathrm{0.61}$, p<0.001) and imply that both the BIT index and IR_6Me register changes in the lake's water balance through its influence on the relative size of the different mixing zones. However, as hydrological proxies they also display important differences, such as the sharpness of recorded transitions in water-column structure. For example, during the Stage III low-stand period from ca. 113–99 ka, the BIT index displays sustained low values (most often < 0.2), suggesting that the strong reduction in lake depth and enhanced upper-water-column mixing started and ended abruptly (Fig. 5b). The IR_6Me during this period begins at low values but steadily increases, inferring a similarly dramatic but more gradual decrease in lake depth (Fig. 5g). BIT-index shifts through time may be rather abrupt because it depends on the proliferation of Thaumarchaeota, which appear highly sensitive to changes in the upper mixed layer (i.e., they either bloom or fail to bloom depending on ambient conditions; Baxter et al., 2021). By contrast, the brGDGTs controlling variation in IR_6Me are produced in the anoxic lower water column, which makes their abundance less sensitive to environmental changes impacting the upper water column, and they thus display a more paced response to gradual expansion/shrinking of the anoxic zones during phases of lake level rise or decline. Consequently the BIT index may be more sensitive and/or respond quicker to climate-driven changes in lake water balance (as attested by its truthful registration of relatively modest drought events in the last 200 years; Buckles et al., 2016), even though overall lake depth is likely the most important factor controlling the BIT index on long timescales. Episodes when both the BIT index and IR_6Me infer abrupt changes in lake depth, for example, at ca. 140 and 80 ka, thus likely do represent drastic changes in climatic moisture balance within a relatively brief period of only a few hundred years. Additionally, under extreme shallow oxycline conditions (i.e., strong upper-water-column stratification) the BIT index has a sensitivity threshold defined by the near absence of Thaumarchaeota (BIT-index values approaching 1), beyond which variations in inferred wet climate conditions are no longer registered (Baxter et al., 2021). In the DeepCHALLA record the BIT index approaches 1 during several sustained periods between ca. 80 and 24 ka, largely overlapping with seismic Stage IV and suggesting persistent high lake level and extremely stratified water-column conditions (Fig. 5b). The IR_6Me proxy continues to vary in this section (Fig. 5g) and therefore provides valuable additional information on changes in lake depth or water-column structure during such intervals. Certainly, the exact relationship between 5-Me and 6-Me brGDGTs and the particular water-column zones which allows their relative proportion to be used as a qualitative proxy for lake depth in Lake Chala is not likely universal across different lake systems. Even in the similarly meromictic Lake Lugano (Switzerland) markedly different brGDGT distributions with depth occur (Weber et al., 2018), implying that comprehensive local water-column profiling of GDGT distributions is necessary prior to interpretation of downcore BIT-index or IR_6Me records.

5.3 Reliability of GDGT-based temperature proxies in Lake Chala

Multiple lines of evidence presented above (see Sect. 5.1) indicate that GDGT variability in the deepest, oldest portion of the DeepCHALLA sediment sequence is not predominantly controlled by past climate variation, prompting restriction of GDGT-based temperature inference (using TEX₈₆, ${\text{MBT}}_{\text{5Me}}^{\prime }$ or MST) to the last 180 kyr of the record. Previous studies (Sinninghe Damsté et al., 2012a; Buckles et al., 2013; Baxter et al., 2021) already showed that the TEX₈₆ index of Lake Chala sediments is highly sensitive to contributions from non-Thaumarchaeotal sources to the sedimentary isoGDGTs, leading to incorrect temperature reconstruction. To distinguish between trustworthy and problematic temperature estimates, the isoGDGT-0 / cren, f[CREN^′] and %isoGDGT-2 proxies have been used to detect TEX₈₆ values which are likely influenced by mixed isoGDGT sources and therefore do not accurately reflect past temperature (Sinninghe Damsté et al., 2012a). In Lake Chala the non-Thaumarchaeotal contribution can be assessed using threshold values informed by GDGT distributions in multi-year SPM and settling particle data (Baxter et al., 2021). Accordingly, sediment horizons with either BIT index>0.8, isoGDGT-0 / cren>0.7, $\mathrm{f}\left[{\text{CREN}}^{\prime }\right]>\mathrm{0.04}$ or %isoGDGT-2>45 % (407 out of 798 in the 180–0 ka section or 51 % of the total) indicate that the sedimentary isoGDGTs are not largely derived from (Group I.1a) Thaumarchaeota and, hence, that TEX₈₆ temperature estimates from these horizons should be rejected. The consequent rejection of the overwhelming majority of measurements on Stage IV sediments prevents TEX₈₆-based temperature reconstruction over most of the glacial period (MIS4 to MIS2) between ca. 85 and 20 ka (Fig. 8b). Moreover, with the exception of a TEX₈₆ increase from ∼ 0.45 to 0.6 after 20 ka (translating to a warming of > 10 °C using the calibration of Tierney et al., 2010a), presumed trustworthy TEX₈₆ values remain relatively stable at ∼ 0.45–0.55 throughout the last ca. 180 kyr, which can be considered unrealistic as available climate records indicate that also low-latitude regions including tropical Africa experienced significant temperature variation at glacial–interglacial timescales (Blome et al., 2012; Johnson et al., 2016). Thus even though the named thresholds provide some guidance for identifying trustworthy values, the exclusion of more than half of the measured sediment horizons and apparent lack of response of the remaining time series to global temperature trends compromise the reliability of the TEX₈₆ index as a paleothermometer in Lake Chala.

Variation in the ${\text{MBT}}_{\text{5Me}}^{\prime }$ index throughout the 180–0 ka section of the DeepCHALLA sequence implies a temporal pattern of temperature change which is likewise incompatible with global temperature variation at glacial–interglacial timescales. For example, while ${\text{MBT}}_{\text{5Me}}^{\prime }$ correctly estimates regional temperature to have been above-average during the last interglacial period (MIS5), it shows no clear pattern of alternating warmer and cooler conditions during interstadials MIS5e–MIS5a. Also the inferred coldest interval of the last glacial cycle is reconstructed to have occurred around 50 ka within MIS3 rather than during the Last Glacial Maximum in MIS2. Notably the abrupt shifts increase and decrease in ${\text{MBT}}_{\text{5Me}}^{\prime }$ at, respectively, ca. 140 and ca. 80 ka and occur simultaneously with major hydrological changes in Lake Chala as recorded by the BIT index and IR_6Me (Figs. 5 and 8). This suggests a strong imprint of climate-driven hydrological change on ${\text{MBT}}_{\text{5Me}}^{\prime }$ and a reason to disqualify its use for temperature reconstruction at this site. Periodicity analysis of ${\text{MBT}}_{\text{5Me}}^{\prime }$ over the last 180 kyr does reveal cycles of 24.8 and 37.6 kyr potentially reflecting a response to insolation variation due to precession and obliquity, but only the obliquity-band signal is significant using the 95 % confidence interval. However, as similar cyclicities occur in the hydrological proxies crenarchaeol, BIT index and IR_6Me (Figs. 9, S9 and S10), and as especially precession strongly influences tropical monsoon dynamics, this does not suffice to use ${\text{MBT}}_{\text{5Me}}^{\prime }$ as a paleotemperature proxy at Lake Chala.

Recent GDGT-based climate studies on African lake records (e.g., Bittner et al., 2022) indicate that in order to achieve more consistent temperature reconstructions 6-Me brGDGTs should be included in the calibration of brGDGT-based temperature proxies (Baxter et al., 2023). Application of the global lake calibration by Pearson et al. (2011), which estimates mean summer temperature (MST) based on the combined abundance of 5-Me and 6-Me brGDGTs, to the DeepCHALLA sequence produces a temperature reconstruction displaying peak temperatures of 25 °C being reached during the current interglacial period (the Holocene) and between ca. 140 and ca. 130 ka, which likely (i.e., accounting for chronological uncertainty in the present age model) represents MIS5e, the warmest known episode of the last interglacial period (Fig. 8e). Conversely, inferred cool episodes centered at ca. 150, ca. 60 and 15 ka correspond to known periods of extensive high-latitude glaciation during the penultimate ice age (MIS6) and the last glacial period (MIS4 and MIS2). Importantly, the MST time series does not show large shifts at ca. 140 and ca. 80 ka, together with those in the BIT index and IR_6Me, indicating that this temperature proxy is not affected by the reorganization of the water-column structure related to changes in lake depth and is thus likely a more trustworthy tracer of changes in local air temperature than ${\text{MBT}}_{\text{5Me}}^{\prime }$. The MST time series also does not display clear changes associated with low BIT-index and isoGDGT-0 / cren values during the Stage III low stand, supporting the notion that the GDGT drivers of the MST proxy vary independently from changes in lake mixing. We therefore suggest that the MST time series covering the last 180 kyr of the DeepCHALLA sequence may constitute a reliable record of past temperature variation in the Lake Chala region covering almost two glacial–interglacial cycles. This suggestion appears to be supported by clear periodicities of 23.0 and 40.1 kyr in these data over the last 160 kyr, presumably reflecting the influences of orbital precession and obliquity on tropical African paleoclimate (Fig. S11e).

6 Conclusions

The analysis of iso- and brGDGTs in the 250 kyr DeepCHALLA sediment sequence from Lake Chala shows that the first ca. 70 kyr of sedimentation is characterized by relatively low GDGT concentrations and erratic variation in the BIT index and isoGDGT-0 / cren ratio, suggesting a poorly stratified water column with highly unstable oxycline position. A comparison with independent measures of lake-basin evolution indicates that water-column structure during this early period was dissimilar to the present-day situation because the lake had a leaky hydrology which prevented the accumulation of dissolved solids, thereby hampering chemical stratification and the formation of distinct mixing zones. Hence the differentiated niches of various GDGT producers as occurring in the water column today were not yet established, resulting in GDGT proxy values that were not predominantly controlled by climate variability. In line with this, time series of crenarchaeol concentration, the BIT index and IR_6Me only start to display periodicities reflecting orbital insolation forcing of the region's climate from ca. 180 ka onwards, suggesting that from around this time climate rather than lake-basin evolution exerted the primary control on niches of GDGT producers and hence GDGT-derived proxies in the lake. The connection between GDGT proxies and climate as understood on the basis of modern-system studies solidified between ca. 180 and ca. 80 ka as the lake gradually developed the strong chemical gradient in the lower water column characterizing modern-day conditions, permitting increasingly more trustworthy quantitative inferences of past climate variability during those more recent periods.

The IR_6Me, which captures the relative proportion of 6-Me brGDGTs and 5-Me brGDGTs, in Lake Chala is related to past changes in lake depth altering the relative size of the distinct niches where these lipids are most abundantly produced. Hence IR_6Me can be an important proxy for investigating past changes in climatic moisture-balance changes in this system alongside the BIT index. Detailed consideration of alternative GDGT-based paleothermometers resulted in the rejection of the TEX₈₆ temperature proxy, as previously set filtering criteria (Baxter et al., 2021) indicate that half of the sediment horizons younger than ca. 180 ka contain a significant fraction of non-Thaumarchaeotal isoGDGTs. In particular, the strong influence of upper-water-column mixing on Thaumarchaeota niche space casts doubt on the application of this temperature proxy in Lake Chala and likely also other (tropical) lakes experiencing shallow oxycline conditions. The ${\text{MBT}}_{\text{5Me}}^{\prime }$ index, commonly assumed to best capture the temperature dependence of brGDGTs (Russell et al., 2018; Martínez-Sosa et al., 2021), results in a reconstruction that lacks a clear glacial–interglacial pattern and shows evidence for an overprint of lake mixing influences. Consistent with indications of the importance of 6-Me brGDGTs in temperature proxies applied to Lake Chala (Baxter et al., 2023), we find that MST reconstruction using the global lake calibration of Pearson et al. (2011) displays a strong and temporally feasible alternation between inferred warm and cool episodes over the last glacial–interglacial cycle and contains clear periodicities related to the long-term variation in solar insolation due to orbital precession and obliquity. We stress that the niche partitioning of GDGT producers within the distinct water-column mixing zones of Lake Chala underpinning the temperature proxies as presented here is potentially unique (as discussed previously in Baxter et al., 2021). Dimictic and monomictic lakes, as are common in cold- and warm-temperate regions, may be equivalent to Zones 1–2 or Zones 1–3 of Lake Chala depending on their trophic status (i.e., level of primary productivity), as this would control the occurrence and seasonal persistence of hypolimnetic anoxia. Unstratifying (truly shallow) lakes enjoy year-round oxygen injection and nutrient upcycling and thus present GDGT producers with entirely different niches. At the other extreme of the depth gradient, permanent stratification of the lower water column (meromixis) can be created and maintained by different processes and with different prevalence in temperate and tropical lakes (Gulati et al., 2017); consequently the occurrence of GDGT niches as in Zones 4–6 of Lake Chala is more likely site-specific (and as documented here, may change over long timescales). Thorough modern-system studies should ideally be conducted before applying such proxies to the sedimentary record of any type of lake system.

Importantly, the types of chemical and physical changes that characterize lake-system evolution at Lake Chala are not altogether unique, and similar processes are certainly involved in the history of most lake basins. To date these potential confounding factors are generally not considered when interpreting biomarker-based climate records from lakes. This work shows the necessity of applying a comprehensive approach which incorporates lake-basin information when interpreting downcore trends in sedimentary GDGT proxies to reconstruct past climate history, in particular when involving GDGTs that are produced in situ in the water column or sediments. Based on our findings, particular caution is recommended when interpreting proxy records that extend to the initial filling stage of lakes or include episodes when lake-system functioning and sedimentation were clearly different from today, regardless of the apparent continuity of lacustrine deposition.