Late Holocene climate anomaly concurrent with ﬁ re activity and ecosystem shifts in the eastern Australian Highlands

The alpine area of the Australian mainland is highly sensitive to climate and environmental change, and poten- tially vulnerabletoecosystemtippingpoints.OverthenexttwodecadestheAustralianalpineregionispredicted toexperiencetemperatureincreasesofatleast1°C,coupledwithasubstantialdecreaseinsnowcover.Extending the short instrumental record in these regions is imperative to put future change into context, and potentially provideanaloguesofwarming.Wereconstructedpasttemperatures,usingalipidbiomarkerpalaeothermometer technique and mercury ﬂ ux changes for the past 3500 years from the sediments of Club Lake, a high-altitude alpine tarn in the Snowy Mountains, southeasternAustralia. Using a multi-proxy framework, includingpollen and charcoal analyses, high-resolution geochemistry, and ancient microbial community composition, supported by high-resolution 210 Pb and AMS 14 C dating, we investigated local and regional ecological and environmental changes occurring in response to changes in temperature. We ﬁ nd the region experienced a general warming trendoverthelast3500years,withapronouncedclimateanomalyoccurringbetween1000and1600calyrs.BP. Shifts in vegetation took place during this warm period, characterised by a decline in alpine species and an increase in open woodland taxa which co-occurred with an increase in regional ﬁ re activity. Given the narrow al- titudinal band of Australian alpine vegetation, any future warminghas the potential to result inthe extinction of alpine species, including several endemic to the area, as treelines are driven to higher elevations. These ﬁ ndings suggest ongoing conservation efforts will be needed to protect the vulnerable alpine environments from the combined threats of climate changes, ﬁ re and invasive species.


H I G H L I G H T S
• Historical climate and environmental change not well constrained in Australian Alps • To investigate these changes we analysed proxy indicators including lipid biomarkers. • We find increased temperatures in the Australian Alps between 1000 and 1600 years ago. • This change was concurrent with increase in regional fire activity and vegetation shifts. • Future warming may result in the extinction of (endemic) alpine species.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Alpine areas in Australia occupy only a small fraction of the largely flat continent, effectively making them 'habitat islands'. The ecosystems of these areas are characterised by high biodiversity and endemism, with 212 species of vascular plants, 21 of which are endemic to the Snowy Mountains alpine area (Costin et al., 2000;Venn et al., 2017). The alpine area is only 135 km 2 (Green and Stein, 2015). Despite this importance, these ecosystems are considered vulnerable to widespread phase changes that could fundamentally alter their dynamics (Elmendorf et al., 2012;Laurance et al., 2011). In addition, the marginal alpine snowpack of the area, a critical water resource for south-eastern Australia, is increasingly under threat (Pepler et al., 2015). Alpine areas face uncertainty for their future longevity, with high vulnerability to anthropogenic and natural disturbances, and low rates of ecosystem resilience and recovery following disturbance (Hoffmann et al., 2019).
Arguably, the most serious threat to these endemic ecosystems is climate change: by 2050 the Australian Alps are predicted to warm by between 0.6 and 2.6°C, with a predicted decrease in snow cover of 22-85% (Hennessy et al., 2008). Many of the alpine species of the Snowy Mountains exist near their climatic limits, constrained by altitude, and are at risk of becoming regionally extinct if these climatic thresholds are exceeded (Pickering et al., 2008). Indeed, a study measuring the changes in phenology and morphology of Australia's native alpine plants found that they do not appear to be adapting to changed conditions (Sritharan et al., 2021). In addition, fire activity in alpine landscapes in response to warming climate is expected to increase, with higher landscape flammability related to an increase in woody vegetation, resulting in a positive feedback (Camac et al., 2017). The competitive advantage of a range of invasive species may also be enhanced as a result of increasing temperatures and decreasing snow cover (Pickering et al., 2004) which may be further facilitated by the predicted increase in the frequency of disturbances, such as bushfires and drought (Pickering et al., 2008). Other perils, such as soil erosion and habitat fragmentation from overgrazing by feral species, as well as invasive fauna, pests, and pathogens, are likely to further reduce the resilience of these endemic ecosystems (Green, 2016). These changes are not just of future concern; over the previous few decades, there has been a significant reduction in the amount of snow cover and an earlier reduction in the annual snowpack, linked to increased temperatures (Green, 2010;Green and Pickering, 2009).
To better understand the potential vulnerability and resilience of these alpine ecosystems, a longer temporal perspective than is possible with instrumental records, which span only a few decades, is required (Hennessy et al., 2008). Quantitative palaeoclimatic data from Australia (and the wider Southern Hemisphere) are also rare, and this is especially true for the reconstruction of temperature. The exceptions are generally from further afield, such as the EPICA δD record from Antarctica (Röthlisberger et al., 2005), or limited in length, such as the temperature-calibrated carbon isotope speleothem reconstruction from Yarrangobilly Caves, situated within the Australian alpine region (McGowan et al., 2018), which covers the last two millennia. Pollen transfer functions (e.g. Cook and van der Kaars, 2006;Fletcher and Thomas, 2010) have been developed in various Australian ecosystems; however, vegetation and temperature relationships are not always straightforward and are considered as semi-quantitative estimates.
To investigate past climate and environmental change in this vulnerable region, we sampled sediments from Club Lake: a small snow-fed alpine lake in the Australian Alps. It is well-recognised that the sediments of lakes preserve an integrated record of climate and environmental change through time and therefore act as 'sentinels of change' in the environment (Adrian et al., 2009). Here we use lipid biomarker analysis, which works on the basis that microbial communities adjust the chemical structure of their cell membranes in response to environmental temperature (Weijers et al., 2007). Glycerol dialkyl glycerol tetraethers (GDGTs), a class of microbially-derived lipids, are abundant in terrestrial, lacustrine and marine sediments. Indices derived from these components have been used in a range of palaeoenvironmental reconstructions (e.g. Baker et al., 2019;Foster et al., 2016;Powers et al., 2010;Tierney, 2012). Surface-sediment calibrations have established a strong empirical relationship between the relative distribution of GDGTs and temperature, which has enabled global and regional calibrations for reconstructing palaeotemperatures from lake sediments. In lake sediments, the methylation index of branched tetraethers from modern samples has been shown to correlate well with measured annual air temperature Tierney et al., 2010;Zink et al., 2010), while other calibrations are based on regression between GDGTs and temperature (Pearson et al., 2011). While the exact mechanisms by which biomarkers such as branched GDGTs (brGDGTs) respond to the environment are not fully understood, laboratory experiments simulating the lake environment have shown a strong response between the production of brGDGTs and temperature (Martínez-Sosa et al., 2020).
We also considered total mercury (Hg) concentration, which has recently been recognised as a useful indicator of climate change due to its sensitivity to physical transformations (volatilization and condensation) and specifically to temperature changes (Schneider et al., 2020). As an example, Schneider et al. (2020) described a 4-fold increase in atmospheric mercury deposition in lake sediments in the Venezuelan Andes, associated with a 3°C increase in temperature during the Bølling and Allerød chronozones. This work, and others (e.g. Stern et al., 2012), demonstrate that the mercury flux and deposition in sediments increases during warmer periods; however the relationship is complex, as other processes, including wind and fire, can influence the delivery and deposition of mercury (Biester et al., 2007;Martínez-Cortizas et al., 1999). Mercury deposition rates have also been shown to increase as a result of industrialisation, even in remote areas (Biester et al., 2007). Prior to industrialisation, however, mercury accumulation in natural archives represents a novel new tool for semi-quantitative palaeotemperature reconstruction. A range of other proxies, including pollen, charcoal, organic matter content, elemental geochemistry, and ancient microbial community composition, were also investigated to consider the environmental and ecosystem responses to changes in temperature. We hypothesised that these proxies would indicate links between past warming, fire activity and vegetation changes.

Study site
Located 120 m above the present-day regional timberline, Club Lake (36.413°S, 148.291°E) is a high-altitude (1955 m above sea level) alpine lake in the Snowy Mountains of southeastern Australia (Fig. 1). The lake Z.A. Thomas, S. Mooney, H. Cadd et al. Science of the Total Environment 802 (2022) 149542 is a shallow (<1.7 m deep), moraine-dammed tarn, and was formed from past cirque glaciation. The Club Lake basin is enclosed by ridges on its northern and western sides and shallow, poorly drained peat on the southern border of the basin leading to a small outflow stream. The lake is approximately 1.46 ha in area and is oligotrophic and isothermal (Timms, 1980). Mean annual precipitation at Thredbo Automatic Weather Station (AWS;1957 masl;36.40°S, 148.29°E) at the same altitude and just 1.5 km from Club Lake, is 1406 mm (BoM, 2020). Summer temperatures are mild, with a mean maximum summer temperature of 15.8°C at the Thredbo AWS, and winter temperatures regularly fall below freezing (mean winter minimum -4.6°C). Snow generally persists for around four months, though semi-permanent snow patches on the leeward side of ridges can persist throughout summer and autumn (Costin et al., 2000). As one of only 5 natural snow-fed alpine lakes on the Australian mainland, these lakes have been subject to substantial monitoring efforts (Green, 2012(Green, , 2011. Hourly water temperature was measured with a temperature logger (Tinytag plus -Gemini Data Loggers, Chichester England), suspended at a depth of 30 cm near the outlet of Club Lake from the years 2005-2016. Ice formation and ice breakup (defined as the day when ice breaks out from the central area of the lake; Green, 2011), has also been recorded since 2005.

Sampling description
Four adjacent sediment cores (CL1, CL2, CL3, CL4) were extracted from Club Lake in 2011 with a piston corer at a water depth of approximately 1 m. While care was taken to ensure that the sediment surfacewater interface remained intact in all cores, sediment was lost from the bottom of core CL3 during extraction, and was therefore discarded. The cores were frozen and split with a core splitter to avoid possible compaction or distortion by extrusion. The sediment stratigraphy of the three cores was visually consistent (Fig. S1). Core CL1 was sampled for GDGT temperature calibration, and samples for radiocarbon dating taken from CL1, CL2 and CL4. 210 Pb dating and all other proxy analyses were conducted on CL4.

Glycerol dialkyl glycerol tetraethers (GDGT) distributions
Contiguous samples of 1 cm depth were extracted from the Club Lake sediment sequence (core CL1), totalling 32 samples (plus 5 duplicates). Sediment samples were freeze-dried, homogenised and 2 g of dry sediment weighed. The GDGT fraction was extracted with DCM: MeOH (9:1) on a Dionex accelerated solvent extractor (ASE), following the method outlined in Schouten et al. (2002). This fraction was subsequently passed through aluminium oxide (Al 2 O 3 ) columns, and spiked with a deuterated C 46 standard (Huguet et al., 2006), monitoring m/z 748. A one-point calibration was used to quantify GDGTs using the internal standard. The analysis was by HPLC-MS, following Baker et al. (2019) and Schouten et al. (2002). Peak areas were integrated following the methods described in Weijers et al. (2007). It must be noted that the LC-MS method does not separate 5-methyl from 6-methyl branched GDGTs and therefore does not allow full analysis of branched GDGTs (De Jonge et al., 2014;Hopmans et al., 2016;Martínez-Sosa et al., 2021). It is, however, consistent with analysis of these compounds in previous studies of lake sediments (Loomis et al., 2012;Pearson et al., 2011;Sun et al., 2011;Tierney et al., 2010;Zink et al., 2010). Only peaks with areas >5000 were considered, as a user-selected cut-off for excluding low-concentration substances (c.f. Schouten et al., 2007). We calculated a range of different soil (Peterse et al., 2012;Weijers et al., 2007) and lake (Loomis et al., 2012;Pearson et al., 2011;Sun et al., 2011;Tierney et al., 2010;Zink et al., 2016Zink et al., , 2010 branched GDGT calibrations to investigate which method corresponded most closely with the modern instrumental temperatures from Club Lake.

Pollen and charcoal
Fourteen pollen slides were prepared from 1 cm 3 volumetric samples down the sediment column using standard palynological techniques (Faegri and Iverson, 1975) consisting of an alkali treatment (10% NaOH at 100°C for 10 mins), sieving (125 μm), acetolysis and mounting in silicone oil. A known concentration of Lycopodium spores were added to the pollen preparations to facilitate estimates of pollen and microcharcoal concentrations. Pollen, spore and charcoal particles (<125 μm: hereafter referred to as microcharcoal) were counted at ×400 magnification until a minimum of 150 terrestrial pollen grains were identified (150 represented the minimum across all available pollen slides). Pollen types were identified to family or species level where possible using a combination of print and online resources (http://apsa.anu.edu.au/). Pollen types were grouped according to life history traits and common community associations (alpine, open woodland, heath, marsh/wetland and regional) using information from herbarium records, PlantNET and Costin et al. (2000); see Table S5. Microscopic charcoal counts were converted to accumulation rates (microCHAR<125 μm: no./cm 2 /yr) using charcoal concentrations (calculated using the Lycopodium spore spike) and the deposition time (yr/cm) determined from the age-depth model. Correspondence Analysis (CA) was performed on terrestrial pollen percentage values in the 'vegan' package in R (Simpson and Oksanen, 2020;Simpson, 2007). Environmental variables (temperature, macroscopic charcoal, mercury concentration and total organic matter) were projected onto the CA ordination space to examine the correlation between pollen composition to corresponding environmental variables using the envfit function in 'vegan'.
Macroscopic charcoal was isolated from contiguous sub-samples (every 0.5 cm, n = 49) from core CL4 and quantified using the method outlined in Mooney and Tinner (2011). Volumetric (2 cm 3 ) subsamples were immersed in <4% sodium hypochlorite (bleach) for 24 h, after which they were then gently wet sieved (using a 250 μm sieve). Charcoal in the collected material was hand sorted and the concentration quantified (mm 2 /cm 3 ) using ImageJ image analysis software (Rueden et al., 2017). Macroscopic (>250 μm) charcoal accumulation (CHAR>250 μm: mm 2 /cm 2 /yr) was calculated by dividing the charcoal concentration by the deposition time estimated from the age-depth relationship.

Total organic matter
Organic content of the Club Lake sediments was estimated using loss-on-ignition (LOI). LOI was determined using methods described by Bengtsson and Enell (1986): contiguous 1 cm samples (n = 24), were oven-dried (105°C for 24 h) and then placed in a muffle furnace (550°C for 4 h). LOI was calculated as the mass lost expressed as a percentage of the oven-dried mass.

Mercury
Mercury analyses were conducted every 0.25 cm using a Milestone Direct Mercury Analyzer (DMA-80 Tri-Cell; Milestone, Bergamo, Italy). The DMA-80 determines total mercury concentration through a sequence of thermal decomposition, amalgamation, and atomic absorption spectrometry. Samples were analysed using method 7473 of USEPA (1998); two blanks and two Standard Reference Materials (SRMs) were analysed for every 40 samples. Approximately 100 mg of dry sample material was weighed into nickel boats. After every tenth sample, a replicate sample was analysed. When replicate recovery exceeded a variance of 10% compared to the original sample, a third replicate was run. Certified sediment reference materials NIST-2711a (Montana River Sediment) and SECCC WQB-1 (Lake Ontario) were analysed, and results were in agreement with published values (Table S6).

XRF core scanning and magnetic susceptibility
Non-destructive scanning geochemical analysis was undertaken using an ITRAX μXRF (micro X-ray fluorescence) core scanner (Croudace et al., 2006) at the Mark Wainwright Analytical Centre, University of New South Wales. The core section CL4 was analysed at a resolution of 200 μm with a 20s count time. X-rays were generated using a Mo-tube set at 30 kV and 55 mA. Data was quality-checked to remove potentially unreliable data points, including elements with low countrates and low signal-to-noise ratio. We therefore report the following elements: Si, K, Ca, Ti, Cr, Mn, Fe, Ni, Cu, Zn, Rb, Sr, Zr. Raw elemental data was normalised by total kilo counts per second (kcps). These data can be considered semi-quantitative and are most reliable when considering relative changes across the sequence (Croudace et al., 2006). In order to apply multivariate statistical analysis, we use a centre log ratio approach to correct for the 'closed-sum' issue of compositional datasets (Weltje and Tjallingii, 2008). Principal Component Analysis (PCA) was performed on the transformed μXRF elemental profiles from CL4 using the R 'analogue' package (Simpson and Oksanen, 2020;Simpson, 2007) to detect patterns of variability that are shared across elements. Magnetic susceptibility was measured at 5 mm intervals on the core using a Bartington Magnetic Susceptibility meter using a point sensor.

Ancient DNA analysis of microbial community composition
All laboratory work was completed in a specialized ancient DNA laboratory at the University of Adelaide. DNA was extracted using a combination of DNA extraction methods from approximately 0.25 g of each of the sediment samples (at 2, 7, 12, 17, and 21 cm depth from core CL4, corresponding to 140, 770, 1130, 1420, 1770 cal yr BP) and three controls (DNA extraction and two sampling swab controls from the test tube and bench); the first three steps of the Qiagen PowerLyzer PowerSoil kit were completed (i.e., bead beating, protein precipitation, and inhibitor removal), followed by the Rohland and Hofreiter (2007) extraction method to maximize recovery of short DNA fragments. DNA extracts were then prepared for shotgun metagenomic DNA sequencing using the Kircher et al. (2012) library preparation method; DNA libraries were sequenced on an Illumina HiSeq 10X using 150 base pair paired end DNA sequencing. An average of 50.3 M sequences per samples were obtained from sediment samples, while 2121 sequences were obtained for sampling and laboratory controls. The rates of duplicate DNA sequences were low in the sediment (averaging 1.5%), but notably higher in control samples (9.0%), as expected. Average ancient DNA length in sediment samples was consistent with samples of this age (average length of 81.2 bp in sediment vs 66.6 bp for controls). Sequences were taxonomically identified using the MALT against SILVA 132 database (Eisenhofer and Weyrich, 2019;Weyrich et al., 2017) and RMA6 files were converted using the weighted LCA algorithm in MEGAN6 CE (version 6.17.0) with the following LCA parameters: weighted LCA, lcp 80, mpi 90%, and supp 0.05 (Huson et al., 2016). Data was then imported in QIIME2 (v2020.8) (Bolyen et al., 2019) and microbes identified in controls were removed from the downstream analysis. Compositional differences were calculated using Bray-Curtis distances (rarefied to 100 sequences), and statistical differences were calculated using adonis and ANCOM (Bolyen et al., 2019) of DNA sequences < 300 bp in length to limit the contributions of modern microbial DNA.
3. 14 C and 210 Pb dating and age modeling To derive a chronological framework for the Club Lake sequence, bulk sediment was sampled for radiocarbon analysis. All radiocarbon samples were pre-treated, graphitised and measured for radiocarbon by accelerator mass spectrometry (AMS) at the Waikato Radiocarbon Laboratory, New Zealand. All samples were pretreated using an acidbase-acid (ABA) protocol (with multiple base extractions) and then combusted. The age of recent sediment accumulation at Club Lake was determined by detecting the abundance of atmosphere-derived 210 Pb. The 210 Pb analysis was conducted at the Australian Nuclear Science and Technology Organisation using alpha particle spectrometry for the determination of polonium-210 ( 210 Po) and 226 Ra activities. The method uses down-core change in 210 Pb activity to calculate an accumulation rate based on the 210 Pb half-life of 22.26 ± 0.22 yr (Appleby and Oldfield, 1992). Both the CFCS (Constant Flux: Constant Sedimentation) and the CRS (Constant Rate of Supply) models (Appleby, 2002;Appleby and Oldfield, 1978) were used to calculate calendar ages, and are in good agreement (Table S1). Here we use the CFCS model.
The 14 C and 210 Pb ages were used to develop a Bayesian age-depth model using a P_sequence deposition model in OxCal 4.4 (Bronk Ramsey, 2008;Bronk Ramsey and Lee, 2013) with the General Outlier analysis detection method (probability = 0.05) (Bronk Ramsey, 2009). The 14 C ages were calibrated against the Southern Hemisphere calibration (SHCal20) data set (Hogg et al., 2020). The model was based on 1000 iterations, with the surface (depth zero) and year of sampling (2011) as the uppermost chronological control point. The age model was based on the CL4 core, with the deeper CL2 basal age providing an age constraint for this core. Two ages from CL1 (coloured purple; Fig. 2) are not incorporated into the age model, but instead are overplotted to illustrate the consistency of the age-depth model across the three cores. Modelled ages are reported here as thousands of calendar years before present (1950 CE) or ka (Table 1 and Fig. 2). We used the mean of the modelled age solutions to estimate the age of the sediments at each sample depth.
Three ages were statistically identified as outliers, with a posterior probability of 0. High alpine environments are well known to be challenging to date due to the potential for in-wash of older carbon from the catchment. The anomalously old radiocarbon ages between 0 and 8 cm depth (identified as outliers in the age model) may have been a result of reworking in the catchment (e.g. Thomas et al., 2019), possibly by snow accumulation or increased erosion/soil trampling by cattle and sheep stock introduced after European settlement from the early-tomid nineteenth century. Following Stanley and De Deckker (2002), we included a Delta_R value of U(0,200) to take into account the possibility of older carbon being incorporated into the sediment core (e.g. allowing a ΔR of between 0 and 200 years). This value was chosen as it represents an appropriate uncertainty given the slow accumulation rate of the sediments, and effectively removes any sedimentation rate artefacts.

Club Lake temperature reconstruction: GDGTs
To derive a proxy reconstruction of the air temperature at Club Lake over the late Holocene, we calculated a range of different branched GDGT calibrations, and compared these with the modern instrumental temperatures from Club Lake (  Fig. 2. Age-depth model for the Club Lake sediments using a P_sequence deposition model in OxCal 4.4 (Bronk Ramsey, 2008;Bronk Ramsey and Lee, 2013). 1σ and 2σ age ranges (dark and light blue envelopes respectively) and probability distributions generated from the Bayesian age model are shown. The 14 C ages were calibrated against the Southern Hemisphere calibration (SHCal20) data set (Hogg et al., 2020). The black probability density functions are from core CL4, orange from CL2, and purple from CL1 (which did not contribute to the age model). Posterior and prior (set at 5) probability for each measurement being an outlier are shown in the left-hand text square brackets. Lake since 2005 (Table 2) and these suggest that the mean summer lake temperature ranges from 13.1 to 14.6°C. We also evaluated the data from the closest weather station to Club Lake (Thredbo AWS 071032; BoM, 2020) which is located at a similar elevation of 1957 m asl. The instrumental record of mean summer air temperature (MSAT) is 11.1°C (using the entire historic record from 1966 to 2020) and this has increased slightly to 11.2°C more recently (using the 30 year mean from 1991 to 2020) (BoM, 2020) and on an annual basis MSAT ranged between 9.8 and 12.7°C at the Thredbo AWS.
Although lake water temperature and air temperature are tightly coupled and often used as surrogates for each other (e.g. Edinger et al., 1968;Livingstone and Dokulil, 2001;Webb and Nobilis, 1997), it is clear that both mean annual air temperature (MAAT) and MSAT are consistently lower in the instrumental record compared to the direct measurements of Club Lake water temperature. This likely reflects the differences in how the values were calculated; the lake data logger provided hourly temperature measurements, while the instrumental temperature record used an average of the daily minimum and daily maximum temperatures. Given this situation, we considered the data from the temperature loggers as being more representative of the conditions recorded by the lake sediments.
The choice of branched GDGT calibration had a significant effect on the resulting temperature reconstructions from Club Lake (Fig. S2), ranging from negative values, to a high of 18°C (MSAT). However, all calibrations had a similar trend with generally lower temperatures at the base of the core and higher temperatures in the modern, uppermost sediments. All calibrations also showed a peak in temperature from the samples over the depth of 10-15 cm. Many of the calibrations yielded reconstructed temperatures that were unrealistic or resulted in a very large temperature range (e.g. the calibrations by Zink et al. (2010Zink et al. ( , 2016) that are considered unlikely to have occurred in the late Holocene (Fig. S2). The calibration from Pearson et al. (2011) gives modern MSAT closest to the modern instrumental temperatures measured at Club Lake, and appears to be the most appropriate to use for this study.
The Pearson et al. (2011) calibration is based on a global transect of 90 lakes from the Arctic to the Antarctic, with climatic and environmental settings of some of the lakes used in the calibration similar to that of Club Lake. Notably, the Pearson et al. (2011) GDGT calibration was also found to be the most appropriate calibration to use to calculate palaeotemperatures for Lake Mackenzie, Australia (Woltering et al., 2014). The calibration is based on summer temperature which is advantageous as some lakes are ice-covered for substantial periods over winter and so many lakes are only monitored during the summer. Where sites experience large fluctuations in temperature, MAAT is not a particularly biologically relevant parameter, for example, in terms of productivity. In the highly seasonal environment of Australia's Snowy Mountains, the MSAT temperature calibration from Pearson et al. (2011) provides estimates closest to the recent instrumental record and is used here to reconstruct temperature over the last 3500 years. Our reconstruction of temperature (MSAT) following Pearson et al. (2011) was calculated as follows:  Table 2 Mean annual lake temperature and mean summer lake temperature data from Club Lake (measured at 30 cm depth; average calculated from hourly temperature readings), and mean annual air temperature and mean summer air temperature data from Thredbo AWS (calculated from an average of daily minimum and maximum temperature). Summer temperature is calculated from December, January and February averages. NA indicates that there is missing data so a representative average cannot be calculated. Our reconstruction of MSAT varies between about 13 and 16°C, and appears to display a long-term increase in temperature averaging 0.1°C per century (Fig. 3). The experimental error of this reconstruction of temperature is largely unquantified, but is likely to be a small component of the overall error, which is approximately ±1°C at 1σ (Pearson et al., 2011). While the long-term trend is uncertain, the relative temperature changes reconstructed very likely exceed this error, providing confidence in the pattern of the temperature data. The reconstruction suggests that the period from ca. 3500 to 1500 cal yr BP experienced cooler summer mean temperatures, especially at the base of the series (e.g. from 3500 to 2600 cal yr BP). An increase towards warmer summers began just prior to 1500, with these warmer conditions persisting until ca. 1100, after which temperatures rapidly declined. The period from ca. 1000 to ca. 600 cal yr BP was cooler than the prevailing long-term trend, until another increase in temperature, lasting around 200 years. A long, slow decline in summer temperatures followed from ca. 400 cal yr BP right up into the historic period. In the last 60 years reconstructed summer temperatures reveal another steep increase: the MSAT for this period just exceeds the range of variability for the entire record, with a temperature increase of 0.03°C per year, indicating that the recent rate of increase in summer temperature exceeds anything in the entire 3500 year record.

Discussion of Club Lake temperature reconstruction
Environmental factors beyond temperature have been suggested to influence brGDGT structure, and therefore could potentially bias the temperature calibration. Here we discuss the Club Lake temperature reconstruction in the context of other palaeoenvironmental proxies from the lake, including mercury, ITRAX, organic matter content, charcoal, and pollen. Temperature calibration biases can include pH (Tyler et al., 2010), lake oxygen content (Martínez-Sosa et al., 2020;Tierney, 2012), as well as redox state and microbial community changes (Weber et al., 2018). While a strong correlation between the CBT (Cyclization of Branched Tetraethers) and MBT (Methylation of Branched Tetraethers) ratios may indicate that apparent changes in temperatures may in fact be controlled by changing pH values (Tyler et al., 2010), no correlation between the CBT and MBT ratios was observed in this study (R 2 = 0.0078). In addition, the source of GDGT input could potentially change, for example, from (allochthonous) bacteria associated with organic matter in soil, washed into the lake as a result of erosion, to autochthonous or within-lake bacteria. This is potentially problematic as the calibrations for soils and lakes are not interchangeable (Tierney et al., 2010). However, for Club Lake, other proxy analyses (Fig. 4), including the high-resolution elemental profiles from the ITRAX geochemistry and relatively stable loss-on-ignition (prior to recent European impacts) suggest that the input of clastic or detrital material from the shallow catchment soils was limited. The exception is the very top sample which may have been affected by the introduction of grazing in the catchment (as discussed in Section 4.2). The reconstructed temperature of 16.2 ± 1°C (which is higher than the contemporary Thredbo AWS mean summer temperature) must therefore be treated with caution.
Confidence in the lipid reconstruction of MSAT at Club Lake is supported by changes in mercury deposition, which is an independent proxy sensitive to temperature. This sensitivity is a result of its unique electron configuration which strongly resists removal of an electron, forming weak bonds and hence melting and volatilising at low temperatures (Schroeder et al., 1998). This is particularly the case of deep moraine lakes, at high altitudes and with small catchment area. At Club Lake, the correspondence between the reconstruction of MSAT and mercury is especially strong during the warm period identified between ca. 1500 and 1100 cal yr BP. The mercury flux does not show the same step increase at 1500 cal yr BP demonstrated in the reconstruction of MSAT, but reveals a slow increase starting from about ca. 2000 cal yr BP, though climbing fastest at 1500 cal yr BP. Notably the structure evident in the reconstruction of MSAT over the last 1000 years is absent in the flux of mercury at Club Lake, including the second period of increased temperature identified in the lipid palaeothermometer technique from ca. 600 to 400 cal yr BP.
Mercury atmospheric flux and deposition increased five-fold during the warm period between 1650 and 1000 cal years BP compared to background (Fig. 4). This increase in mercury atmospheric deposition follows closely both the temperature and charcoal reconstructions, highlighting the relationship between warmer climate and mercury atmospheric fluxes and deposition in lake sediments. No significant variations in organic matter (as determined by LOI) occurred during the warm period (ca. 1650 to 1000 cal yr BP), suggesting that the mercury Fig. 3. GDGT-inferred mean summer air temperature calibration calculated from Pearson et al. (2011). Black line shows loess smoothed fit, purple envelope shows 1σ uncertainty. Also plotted is the Thredbo AWS mean and maximum summer air temperature (averaged 2006-2015). increases did not come from the catchment, reinforcing the hypothesis that mercury deposition in Club Lake is predominantly atmospheric. Both the proportion of organic matter and the concentration of mercury rapidly increased in the Club Lake sediments during the last 300 years. Any differences between the reconstruction of MSAT using lipids and the mercury flux in the Club Lake sediments potentially relates to the altered sources, as the flux of mercury is not only sensitive to temperature but also wind (Gustin et al., 1997) and bushfire emissions (Pompeani et al., 2018). The 5-fold increase in mercury deposition in Club Lake between 1650 and 1000 years BP is clearly temporally associated with increases in charcoal (both macro and micro), although any inter-relationships between fire regimes, warmer temperature and mercury fluxes cannot be discounted. The very large increase in mercury flux in the recent period is comparable to that in Biester et al. (2007) and also likely to result from complex sources, including from coal-fired power stations, mining and other industrial activity as well as the increased temperature evident in our reconstruction and the recent instrumental record (Fisher and Nelson, 2020;Marx et al., 2010;Schneider et al., 2019, Schneider et al., 2021.

Environmental reconstructions
The reconstructed changes in temperature at Club Lake are also associated with several other palaeoenvironmental proxies, suggesting environmental responses occurred in response to the altered climatic conditions. Considerable changes in vegetation composition have occurred over the last 2000 years, with a dynamic relationship between warm/cooler MSAT and contraction and expansion of vegetation characteristic of the alpine zone. Changes in the ratio of woodland to alpine species (Table S4), suggest variability in the amount of woody vegetation occurring at high elevation sites near the Club Lake catchment. The period of elevated temperatures from 1500 to 1100 cal yr BP is coeval with the highest woodland:alpine ratio, suggesting an increase in tree species, likely Eucalyptus pauciflora, at higher elevations and a reduction in several alpine zone species (Fig. S2). Vegetation change associated with reduction of MSAT around 1000 to 800 cal yr BP was characterised by an expanding alpine zone, while an increase in the woodland:alpine ratio followed an increase in MSAT from about 800 cal yr BP. Interestingly the most recent increase in temperatures has not corresponded to an increase in the woodland:alpine ratio. In this case, the expansion of woodland species may have been reduced by other anthropogenic impacts, such as cattle grazing and fire activity. While grazing was banned in the high altitude alpine area (including the Club Lake catchment) since the mid-20th century, longer lasting impacts on vegetation communities complicates the relationship between temperature and fire. However, as Pickering et al. (2008) have noted, Australian alpine vegetation occurs in a narrow altitudinal band and any future warming has the potential to result in the extinction of alpine species, including several endemic to the area, as treelines are driven to higher elevations.
The changes in MSAT at Club Lake are also reflected in changes to the accumulation of charcoal in the sediments, both for microscopic and macroscopic charcoal, and this is particularly the case during the warmer period between ca. 1500 and 1100 cal yr BP. The increase in both charcoal size fractions tends to suggest that there was an increase in both local and regional burning during this warmer MSAT. The archaeology and use of the high alpine area by Aboriginal people are currently not well-constrained (e.g. Kamminga, 1995). Whilst there is a possibility the increase in fire activity between ca. 1500 and 1100 cal yr BP is due to an increase in anthropogenic burning, the limited evidence of human use of the Alpine region and coeval increase in temperature and biomass available to burn suggests that this increase in fire was driven by climate and environmental conditions conducive to burning.
Geochemical data generally revealed few major changes over the record, although magnetic susceptibility is elevated during the first warm period (Fig. S3). ITRAX PCA axis 1, explaining 49% of the variance (Figs. 4, S4, S5, Table S7) is strongly positively correlated with terrigenous elements (Fe, K, Mn, Rb and Ti) and negatively correlated with Cu and Ni. The surrounding geology of Club Lake is of schists, quartzite, gneiss and granulite features of Ordovician origin with distinctive chemical features, including of high concentrations of Mg, Fe and other ferro-magnesian elements, high K and concomitant Rb and Cs (Kolbe and Taylor, 1966). The effects of grazing are likely responsible for some of the recent changes (i.e. post-European settlement) recorded in the geochemical record (most notably through strong decrease in PCA Axis 1 and increase in Axis 2). Evidence of increased detrital input following the introduction of grazing is not immediately apparent from PCA Axis 1, however, a large increase in Si in the top few centimetres (contributing to PCA Axis 2) and the increase in LOI probably reflect the well documented erosion within the small catchment associated with post-European grazing. As noted, there are several recent changes in the uppermost sediments indicating the impacts of both proximal and distant (post-European) human activity: this includes the rapid increase in mercury flux at the top of the core, and recent changes in the pollen spectra, reflecting historic summer grazing and the occurrence of introduced taxa, such as Pinus pollen, observed in the top 1 cm sediment sample (Fig. S6).
Microbial community composition response to past environmental conditions and geochemical changes of the ancient sediment remains poorly understood, as sediments represent mixtures of modern and ancient microbial DNA fragments. Therefore, we assessed the microbial DNA signatures in selected core depths at different DNA fragment sizes (number of base-pairs (bp), i.e. 0-50 bp; 51-100 bp, 101-150 bp, 151-200 bp, and 201-250 bp), as ancient DNA is often much smaller in length. DNA fragment length was associated with specific community compositions (p = 0.006, F mode = 1.9778) (Fig. S7, Table S7), reflecting the mixtures of modern (living) microbial communities and signatures of more ancient species. After accounting for DNA fragment length, core depth (adonis test, p = 0.001, test = 3.57748) and predicted age of the sediment sample (adonis test, p = 0.001, F mode = 5.53421) were significant drivers of ancient microbial compositions, indicating shifts in microbial communities over time. This was largely driven by differences in the 2 cm core sample (dated at ca.1810 CE) compared to the older samples (Fig. S7). Next, we examined if microbial communities shifted according to the ecological changes identified with pollen or charcoal analyses, as microbial communities are known to be correlated to plant compositions (Thompson et al., 2017). The composition of ancient microbial communities was significantly linked to the woodland:alpine ratio (p = 0.002, F model = 3.4599) but not the charcoal flux, suggesting that these past microbial communities are linked to the overall plant ecology of the area and not burning. This was also largely driven by data from the 2 cm sample, as Neobodo designis, an uncultured Acidobacteriaceae bacterium, and an uncultured Gemmatimonas species were significantly higher in the 2 cm sample than any others (ancom; w > 120). Further, we observed more Acidobacteria in the 2 cm sample (5 taxa) than other layers (2-3 Acidobacteria taxa). Although time depth and woodland:alpine information are confounded because woodland areas decreased over time, these findings suggest that the past microbial communities in these sediments may have shifted as woodland areas retreated and temperatures warmed in the region. These findings are supported by a previous study which found a link between ancient soil bacterial microbiome and particular vegetation types (Gellie et al., 2017).

Comparison to published regional climate and environmental reconstructions
Determining the possible mechanisms behind the temperature changes and associated ecosystem responses in the Australian alpine region during the late Holocene may also provide insights for potential future vulnerability of these ecosystems. Here, we compare the Club Lake temperature and environmental reconstructions from this study with several other local and regionally important proxy records (Fig. 5). One of the most comprehensive initial pollen studies to be undertaken in the Australian highlands was at Club Lake, where Martin et al. (1986) described sharp changes in abundance of many alpine taxa after 3000 cal yr BP (and also at 1500 cal yr BP), and suggested instability induced by drought, frost, soil erosion may have been responsible. Higher summer temperatures and a longer growing season were suggested to result in an increase in the tree line over the late Holocene (Martin, 1986), which lends support to our interpretation of higher temperatures between ca. 1500 and 1100 cal yr BP at Club Lake, with an inferred treeline shift also weakly constrained by the pollen data. Stanley and DeDeckker's (2002) reconstruction of aeolian activity using the sediments from Blue Lake, located close to Club Lake in the Snowy Mountains, shows a distinct increase in the maximum size of aeolian quartz grains, interpreted to show peaks in windiness, during the periods ca. 1700 to 1400 and ca. 800 to 500 cal yr BP, which just precede (by about 200 years) our identified periods of increased temperature. Their Blue Lake age model used a reservoir correction of 900 years to correct for the presence of 'old carbon' and so any lead or lag may be an artefact of this, however. Similarly, a dust deposition record from a peat bog near the headwaters of the Snowy River also showed variations during the late Holocene, but high age model sensitivity precludes defining the exact timing of these changes (Marx et al., 2011).
On a broader regional scale, millennium-length ensemble reconstructions of regional (Australasian) temperature suggests the warmest period of the last 1000 years (prior to the industrial period), was between approximately 800-600 cal yrs. BP (Gergis et al., 2016), which directly corresponds to the slightly more muted warm period determined in the Club Lake (and other) records. While there are currently no ensemble reconstructions further back in time, a near annual resolution ẟ 13 C speleothem record from the Yarrangobilly Caves (McGowan et al., 2018), located just 80 km north of Club Lake (Fig. 1), spanning the last 2000 years, is available for comparison. Changes in ẟ 13 C were interpreted to reflect temperature and snow cover in the Australian Alps over the last two millennia. Notably, the Club Lake GDGT-derived temperature reconstruction displays a similar pattern: the clearest coherence between the two records are the warming trends between ca. 800 -400 cal yr BP, and in the centuries before 1000 cal yr BP, suggesting that both sites are recording regionally significant temperature changes.

Mechanisms of climate and environmental change at Club Lake
There are a number of climate drivers that affect the climate of south-eastern Australia, most notably the influence of the Southern Hemisphere westerly winds (SHWW), and the El Niño-Southern Oscillation (ENSO). The SHWW are one of the major drivers of climate variability across the Southern Hemisphere and exerts a strong seasonal control on regional precipitation and temperature variability (Gillett et al., 2006). Since the mid-to late-20th Century, a strengthening and poleward shift of the SHWW has been observed (Marshall, 2003;Thompson et al., 2011;Visbeck, 2009), associated with climate-carbon changes across the Southern Hemisphere, and a decline in rainfall across south-eastern Australia and Tasmania (Nicholls, 2010;Taschetto and England, 2009). Further afield, the SHWW are also linked with changes in the tropics, including the intertropical convergence zone (ITCZ) and the El Niño-Southern Oscillation (ENSO) (Cai et al., 2014;Schneider et al., 2014). Indeed, seasonal snow in the Australian Alps is known to be influenced by El Niño events (Pepler et al., 2015).
During the Holocene, changes in the seasonality of insolation have been linked to sea surface temperature changes in the tropical Pacific Ocean (Braconnot et al., 2012;Wanner et al., 2008). The subsequent southward migration of the ITCZ after the mid-Holocene (Haug, 2001) has been argued to have led to a reduced zonal sea surface temperature gradient, weakening the Walker circulation and increasing the frequency of El Niño events (Conroy et al., 2008;Koutavas et al., 2006). Rein et al.'s (2005) reconstruction of El Niño based on offshore Peruvian lithic flux, indicates maximum El Niño strength around 1500 cal. yr BP, and shows a similar pattern to the Club Lake MSAT reconstruction (Fig. 5). Intriguingly, a reconstruction of the Southern Oscillation Index (SOI) by Yan et al. (2011), which is based on the difference between moisture proxies from the western and eastern Pacific, shows a very similar pattern to the Rein et al. (2005) reconstruction over the last 2000 years, but with the opposite interpretation of a positive SOI bringing more La Nina-like conditions between 1500 and 1000 cal yrs. BP. However, the closest ENSO reconstruction from the Holocene, a rainfall record from Swallow Lagoon in northeastern Australia (Barr et al., 2019) does not show much correspondence with our temperature record.
Tropical sea surface temperatures are known to influence the intensity and location of the polar jet, with teleconnections between the SHWW and ENSO (Clem and Fogt, 2015;Ding et al., 2012;Wilson et al., 2016). Several studies have examined changing southern westerly winds over the Holocene from Southern Ocean islands. Stronger winds in the South Atlantic sector have been recorded between 1000 and 1800 cal yr BP and also between 600 and 400 cal yr BP from a range of proxy records, including exotic pollen and microcharcoal concentration from the Falkland Islands (Turney et al., 2016a), charcoal concentration from Laguna Guanaco (Moreno et al., 2014), and biogenic carbonate accumulation rate from the Palm2 marine core (Lamy et al., 2010).
Similarly, a reduction in Dracophyllum on sub-Antarctic Campbell and Auckland Islands (southwest Pacific) between 2000 and 1000 cal yr BP was interpreted to indicate stronger westerly winds (Turney et al., 2016b). The periods of increased westerlies show a remarkable correspondence with increased MSAT at Club Lake (Fig. 5), and hint at broad atmospheric shifts occurring during the late Holocene.
To consider further the relationships between increases in surface temperature in the Snowy Mountains region and dominant synoptic conditions, we used KNMI climate explorer (van Oldenborgh and Burgers, 2005). Due to the limited instrumental records in the Australian high country, we investigated the surface temperature as calculated by the ERA-Interim reanalysis over the grid box 36 to 37°S and 148 to 148.5°E  (Rein et al., 2005); E. Charcoal accumulation rate (CHAR>250) from Club Lake (this study); F. Mean Summer Air temperature (°C) from Club Lake (this study). GAM (generalised additive model) splines fitted to MSAT and CHAR data; significant positive trends are highlighted in red and significant negative trends in blue, determined from the first derivatives of the GAM splines. Dashed vertical line shows timing of the first European cattle grazing in the Australian Alps (1830s).  (corresponding to the Australian alpine region), and then correlated this with pressure and sea surface temperature fields from ERA-Interim and HadISST. We find that high summer temperatures in the Snowy Mountains are associated with high pressure systems over the broader SE Australia region and the Tasman Sea. Fig. 6a shows an alternating pattern of high and low pressure systems in the southwest Pacific region consistent with Rossby wave propagation. Correlation of summer temperatures in Snowy Mountains with HadISST shows links with sea-surface temperatures in the Tasman Sea, in the Niño 3 region of the eastern equatorial Pacific, and the equatorial Indian Ocean (Fig. 6b). Increased penetration of Rossby waves impacting Australia is thought to be related to the atmospheric circulation response to tropical sea-surface temperature forcing, with some studies linking it to increased equatorial Pacific temperatures (Ding et al., 2012;Fogt et al., 2011;Turney et al., 2017), while others emphasise the importance of the Indian Ocean (Cai et al., 2011;McIntosh and Hendon, 2018). For example, the positive phase of the Indian Ocean Dipole in conjunction with El Niño is associated with Rossby wave trains causing decreased rainfall, reduced cloud cover and increased surface temperature over Australia during summer, producing conditions favourable to bushfires and their early commencement (Wang and Cai, 2020). The warm period evident in the Club Lake MSAT between ca. 600 and 400 cal. yr BP has been linked to anomalous north-easterly winds over New Zealand and south-eastern Australia (Goodwin et al., 2014;Lorrey et al., 2008), supporting the notion of southeastwardpropagating Rossby wave trains in this time period. Previous work has suggested the influence of tropical sea surface temperatures on the midlatitudes of the Southern Hemisphere was most recently established in the 1940s (Turney and Fogwill, 2021;Turney et al., 2017). Our reconstructed temperature data from Club Lake indicates that this influence might have also been strengthened during periods of the late Holocene.

Significance and conclusions
The Club Lake MSAT record presented in this study provides valuable quantitative data on the temperature variability in Australia during the late Holocene. While the influence of tropical sea-surface temperatures is currently thought to be a permanent feature across the mid-to highlatitudes, reconstructed temperature data from teleconnected regions can help to identify changes in these relationships over longer timescales. CMIP5 model ensembles project temperature increases that range between 1 and 2°by 2050, and between 1 and 4.5°by 2100 (Fig. S8). The significance of warmer temperatures in this region in the future cannot be overstated. The combined effect of temperature increases with less precipitation falling as snow and higher rates of melting (Whetton et al., 1996) present a number of challenges. Ecologically, Australia's unique alpine fauna and flora are vulnerable to extinction and/or distribution change. While grazing was removed from the catchment in the mid-20th century, and rehabilitation and stabilisation followed (Lembit, 2002), it is possible that post-European impacts resulted in several long-lasting legacies on the landscape (Pickering et al., 2004). These impacts, the modification of disturbance regimes, and the displacement of native species, may also be amplified with climate change. Economically, the seasonal snow fields are directly linked to hydroelectric energy production, irrigation to agricultural regions through diversion of water and snow melt through the Murray and Murrumbidgee Rivers, as well as snow-related tourism. With larger rainfall anomalies and dry trends likely to be exacerbated by the projected future rise in temperatures, major impacts on seasonal snow cover and the vulnerable high-altitude flora that inhabits this region may result.
CRediT authorship contribution statement Zoë A. Thomas

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.