Chironomid-based temperature reconstruction for the Eemian Interglacial (MIS 5e) at Sokli, northeast Finland

The Last Interglacial (Eemian, MIS 5e) can be considered a test-bed for climate dynamics under a warmer-than-present climate. In this study we present a chironomid record from the high latitude Sokli site (N Finland), where a long continuous sediment sequence from the last interglacial has been preserved from glacial erosion. The chironomid-analysis shows a diverse fauna, with dominance of warm-water indicators and shifts in assemblage composition that can be attributed to temperature, lake depth, productivity and habitat availability. Quantitative mean July paleotemperature estimates based on the chironomid data indicate overall mean July air temperatures up to 1 °C warmer than present. Two cooling events can be discerned, the Tunturi event, dated to about 127.5 ka BP, in the lower part of the sequence, and the Värriö event, dated to about 119 ka BP, associated with the beginning of a cooling trend in the upper part of the record. Warm conditions already at the onset of the interglacial contrast with a recent chironomid-based last interglacial temperature reconstruction from Denmark, which suggests a late onset of Eemian warming. The relatively small increase in inferred temperatures compared to present day temperatures at Sokli differs from other high latitude Eemian sites, and likely reflects the influence of the Atlantic Meridional Overturning Circulation in maintaining already elevated temperatures in Fennoscandia during interglacials.


Introduction
Model projections of future climate change suggest temperature increase outside the range of observational records with most pronounced warming occurring in the Arctic (above 70°N) (IPCC 2013). The Last Interglacial [roughly equivalent to Marine Isotope Stage (MIS) 5e and the Eemian Interglacial in northern Europe, 129-116 kyr ago], is the most recent interglacial period when climate was as warm as or warmer than today globally. Global mean summer temperatures during the Eemian climate optimum are estimated to have been approximately 0-2°C higher than present (IPCC 2013). Reconstructed summer temperature anomalies are larger in the circum-Arctic region, approximately 4-5°C above present for most of the Arctic, due to the various feedback mechanisms amplifying insolation-driven warming at high latitudes (CAPE-Last Interglacial Project Members 2006). Paleoclimate information from this period can contribute to our knowledge of climate dynamics and forcing under a warmer than present climate and may be used for validation of climate models. It also serves as a reference for natural climate variability during all phases of an interglacial, from the deglaciation to the thermal maximum to the instability associated with the onset of glaciation.
In order to better understand the drivers of climate dynamics in a warmer-than-present world, spatially distributed proxy data are required. Although there are numerous climate records of the Eemian from central and southern Europe, continuous records from the previously glaciated northern latitudes are scarce due to subsequent glacial erosion (CAPE-Last Interglacial Project Members 2006). The Sokli site, situated in northern Finland, however contains a thick and complete lacustrine sediment record of the Eemian (Helmens et al. 2000(Helmens et al. , 2007. The Eemian sediments are found in sequence with overlying sediments from the Weichselian and the Holocene and have been attributed to the Eemian by climate-stratigraphy and luminescence dating (Helmens et al. 2000(Helmens et al. , 2007Alexanderson et al. 2008). The Eemian sediments analyzed in this study were recovered in 2010 and a first detailed multiproxy study focusing on a cooling event (the Tunturi cooling event) in the lower half of the sequence was published by Helmens et al. (2015). Helmens et al. (2015) used high-resolution pollen analysis combined with chironomid analysis of a section of the record to reconstruct July temperatures. In the present study we increase the time-interval compared to the previous study and use chironomid analysis to encompass the complete Eemian sequence.
Chironomids (Diptera: Chironomidae), non-biting midges, are useful paleoclimate indicators that are widely distributed and abundant in lake sediments (Brooks et al. 2007). Many taxa have narrow ecological optima and their chitinized head capsules are well preserved as fossils in lake sediments (Brooks et al. 2007). Chironomids have been used to quantitatively infer paleotemperatures (mean July air temperatures) by the use of transfer functions (Walker et al. 1991) as temperature is a main driver of the broad-scale geographical distribution of chironomids (Eggermont and Heiri 2012). Chironomid-based transfer functions were initially applied to lake sediment records that covered the last glacial/Holocene transition (Walker and Cwynar 2006;Brooks 2006), but have subsequently been used to reconstruct climate conditions during the Weichselian (Engels et al. 2008(Engels et al. , 2010 and the Last Interglacial (Francis et al. 2006;Axford et al. 2011;Jones et al. 2017).
Transfer function-based reconstructions can be sensitive to the impact of confounding variables and calibration dataset properties (Brooks and Birks 2001;Birks et al. 2010;Velle et al. 2010;Juggins 2013;Salonen et al. 2013;Engels et al. 2014). Ideally, a calibration dataset should cover the entire range of environmental values (e.g. temperature) that could be represented by the fossil material, be from a similar climatic regime and sedimentary environment as the fossil material and have a good coverage of the fossil taxa . Evaluating the effect of calibration dataset properties on the resulting reconstruction may be accomplished by applying different datasets to the same fossil record (Płóciennik et al. 2011;Salonen et al. 2013;Engels et al. 2014). To identify the effects of confounding variables, the reconstructions needs to be evaluated against other paleoecological data and independent climate records . Based on multiproxy data (diatoms, non-pollen-palynomorphs, pollen and macrofossils of aquatic and terrestrial plants, geochemistry), Plikk et al. (2016) and Kylander et al. (2018) have shown that the Eemian paleolake at Sokli was a relatively nutrient-rich lake affected by infilling and associated changes in stratification, mixing and nutrient recycling over the course of lake development. The influence of these variables on the chironomid assemblage and the following reconstruction therefore needs to be thoroughly evaluated.
In this study, we present a chironomid record for the Eemian sediment sequence from the high-latitude Sokli site. We apply three independent chironomidclimate calibration datasets to the record to quantitatively reconstruct paleotemperatures during the Last Interglacial. The results are evaluated through comparison with previously published multiproxy paleolimnological data and a pollen-based temperature reconstruction from the same site . Comparisons are also made with Eemian Interglacial climate records from central Europe and the North Atlantic Ocean in order to analyze spatial variability in Eemian climate dynamics and to discuss potential drivers of past climate change.
Site description and stratigraphy Sokli (67°48 0 N, 29°18 0 E, 220 m a.s.l.) is situated in eastern Finnish Lapland (Fig. 1). The regional vegetation is northern boreal forest dominated by birch, pine and spruce. The present climate is cold boreal with a mean annual temperature of 1°C, a mean July temperature of 13°C and mean February temperature of -14°C and precipitation ranging between 500 and 550 mm/year (Drebs et al. 2002). Lakes in the area are typically ice-covered from October to the end of May. The bedrock belongs to the western-most part of the Kola alkaline province and consists of a carbonaterich magma intrusion, approximately 5 km in diameter, encircled by a hilly ring of fenites (elevations 300-350 m), called the Sokli Carbonatite Massif (Vartiainen 1980). Dispersed residual phosphorusrich deposits occur at the surface of the carbonatite (Talvitie et al. 1981).
Sediments from the entire Weichselian and the last and present interglacial occupy a hollow formed along a fracture zone in the easily weathered bedrock at the center of the massif. These sediments have likely been preserved from glacial erosion due to the location of the site within the ice-divide zone and/or to frozen-bed conditions of the ice sheet, and to the protected situation in a deep weathering depression (Helmens et al. 2007). The Sokli sequence consists of five finegrained or organic units dated to MIS 5e, MIS 5c, MIS 5a, MIS3 and the Holocene respectively, that are intercalated with till, glacio-fluvial or fluvial sediments (Helmens et al. 2000(Helmens et al. , 2007. The sequence rests on sands and till from the penultimate glacial, lying on weathered bedrock at a depth of 28 m below surface. A thick bed of diatom gyttja occurs as a marker horizon near the base of the Sokli sequence. It was first correlated with the Eemian by Ilvonen (1973) based on its pollen content. The Eemian age assignment was later supported by stratigraphic studies of the overlying Weichselian sediments and absolute dating control (Helmens et al. 2000(Helmens et al. , 2007Alexanderson et al. 2008). The Eemian bed is bracketed by optically stimulated luminescence (OSL), thermoluminescense (TL) and infrared stimulated luminescence (IRSL) dating of under and overlying sands, giving an age of the sediments of between 180-150 and 110-95 kyr (Helmens et al. 2000(Helmens et al. , 2007Alexanderson et al. 2008).
The coring site, Sokli 2010/4, is situated in the center of the massif. Coring was performed by the Geological Survey of Finland in winter 2010 by percussion drilling using a flow-through coring device. Nine meters of Eemian diatom gyttja and 3 m of underlying glaciolacustrine sands were recovered. The diatom gyttja was observed between 2500 and 1600 cm depth, and the glaciolacustrine sand between 2800 and 2500 cm depth. In the uppermost part of the diatom gyttja (1650-1600 cm) 1-2 cm thick sand layers appear, while silty gyttja is observed at two levels (2341-2324 cm and 1741-1733 cm). Thin yellowish laminations occur at intervals in the lower half of the sequence (below 1986 cm) (Fig. 2).

Timing of the last interglacial at Sokli
An age model for the Eemian lacustrine sequence based on correlation with other regional and global events has been developed by Salonen et al. (2018). The model relies on an age for the start of interglacial conditions at Sokli, as established by pollen, of 130.9 kyr BP (2517 cm) and an age for the end of the Eemian at Sokli of 117.5 kyr BP (1585 cm, in a parallel core not analyzed for chironomids) and assumes constant sedimentation rates between these points. The chronology presented by Salonen et al. (2018) and used in this study differs from the previously published chronology for the Eemian at Sokli, which was based on land-sea comparison (Helmens et al. 2015), by establishing an earlier onset of interglacial conditions at Sokli.  with an abundance of [ 5% at any one level (%), count sums (number hcs) and chironomid concentrations (number hcs/g wet sediment). Added to this diagram are also taxa (Sergentia coracina-type, Einfeldia dissidens-type, Tanypus) with lower abundances but that showed increased occurrences at defined stratigraphical levels and which have clear ecological preferences. The samples above 1631 cm yielded a low number of head capsules and individual findings are indicated with a (?). Grey bars (VE and TE) correspond to the Värriö and Tunturi cooling events, as outlined in Helmens et al. (2015), Plikk et al. (2016) and Salonen et al. (2018) The ages for the beginning and end of our record are based on the onset and end of interglacial conditions as determined from stalagmites from northern Sweden, Belgium and the Alps . A recent U/Th dating of a stalagmite from the Korallgrottan Cave in northern Sweden indicates deglaciation of southern and central Fennoscandia by 130.1 kyr BP (Personal communication: Frank, N. Institute of Environmental Physics, Heidelberg, Germany). With reference to the last deglaciation when Sokli is estimated to have been ice-free about 1 kyr before the Korallgrottan Cave, we suggest that the Sokli site was deglaciated earlier than the Korallgrottan Cave during the penultimate glaciation and use the timing of the warming recorded in the Alps at 130.9 kyr BP (Moseley et al. 2015) for the onset of the Eemian at Sokli. The date of 117.5 kyr BP for the end of interglacial conditions is based on a date from Vansteenberge et al. (2016). Based on the homogenous lithology and the stable pollen concentration, deposition during the Eemian at Sokli appears to have been generally stable. Chronological uncertainties based on age modelling are estimated by Salonen et al. (2018) to be largest in the Mid Eemian, reaching 2r values of approximately 1500 years ). Attempts to further constrain the Eemian record from Sokli, e.g. by the use of tephra layers, have so far been unsuccessful.

Chironomid analysis
Subsamples were taken from the uppermost part of the glaciolacustrine sands and the overlying 9 m of diatom gyttja (2515-1600 cm). A total of 79 wet samples (1-cm slices) were selected and analyzed for their chironomid content (sample volume ranging from approximately 0.5-2 cm 3 , weight range from 0.29 to 4.17 g). The samples were deflocculated in warm 10% KOH (20 min water bath at 85°C), then sieved over a 100 lm mesh. At least 50 head capsules per sample (Quinlan and Smol 2001) were picked with fine forceps under 259 magnification using a dissecting microscope. The head capsules were placed onto coverslips, air-dried and mounted in EuparalÓ mounting medium. The chironomids were identified mainly using Brooks et al. (2007). Samples containing low abundances of chironomids were combined at two levels (1709-1701 cm and 2235-2231 cm) and used together in all statistical analyses and the ecological interpretation. The uppermost three samples (1629-1600 cm) contained only a few chironomids and were excluded from the statistical analyses. Zonation was done with constrained hierarchical clustering using the package rioja (Juggins 2017) for R version 3.4.0 (R core Team 2017) and the significance tested using broken-stick modelling (Bennett 1996).

Temperature reconstruction
Three chironomid-temperature calibration datasets ( Fig. 1) were selected for use with the fossil Sokli data, based on (1) geographic proximity and similarity in present-day climatological regime to the Sokli site, (2) good representativity of dominant fossil taxa and (3) the relative difference in temperature gradient covered. Details of the three calibration datasets, the Norwegian (Norway and Svalbard) (Brooks andBirks 2001, 2004), northern Russian (Nazarova et al. 2008(Nazarova et al. , 2015 and Finnish (Nyman et al. 2005;Luoto 2009;Luoto et al. 2014Luoto et al. , 2016 calibration datasets, as well as the performance of the models are presented in Table 1. The reconstructions were based on 2-component weighted averaging-partial least squares (WA-PLS) regression models, and predictive power was assessed using leave-one-out cross-validation (jack-knifing; Birks 1995). Outliers, as specified in the original publications, were removed before analysis. A 1000 cycle bootstrap cross-validation was used to obtain sample specific error estimates (eSEPs) for each fossil sample. The reconstructions were made on square-root-transformed data using the program C2 version 1.5 (Juggins 2007). The variance explained by July air temperature was assessed for each calibration dataset by performing a canonical correspondence analysis (CCA) on square-root-transformed data in the program R using the package vegan (Oksanen et al. 2017) for R. The strength of the relationship between chironomid assemblages and July temperature was assessed using Detrended Canonical Correspondence Analysis (DCCA) on square root transformed data in the program CANOCO v. 4.5 (ter Braak and Š milauer 2002).
We applied a number of techniques to test the suitability of the individual calibration datasets and to test the reliability of the reconstructions. The percentage abundance of fossil taxa absent from the calibration set as well as the abundance of fossil taxa considered rare (Hill's N2 \ 5, Hill 1973) in the calibration set were calculated to assess the similarity between fossil and modern samples, using the program C2 (Juggins 2007). The fit between fossil samples and calibration samples, and the fit between fossil samples and temperature was further assessed using the Modern Analogue Technique (MAT) and the goodness-of-fit (G-O-F) statistics (Birks et al. 1990). The distance between a fossil sample and its most similar modern assemblage was compared, and the 5th percentile of all squared chord distances of the modern data was used to define the cut-off for a good analogue. G-O-F was assessed by adding the fossil samples passively into a CCA ordination with July air temperature as the only constraining variable and calculating the squared residual distances between the fossil and calibration samples and the constrained axis. Fossil samples with residual distances larger than the 90th or 95th percentile of the residual distances of all modern samples were identified as having a 'poor' or 'very poor' fit respectively. Statistical significance of the reconstructions was further tested using the method of Telford and Birks (2011). The MAT and statistical significance test were performed on percentage respectively square-root-transformed data using the programs C2 and R (Juggins 2007; R Core Team 2017), package palaeoSig (Telford and Trachsel 2015) and G-O-F statistics was accomplished with CANOCO 4.5 (ter Braak and Š milauer 2002) on square-root-transformed data with rare taxa downweighted.

Chironomid stratigraphy
A total of 80 chironomid taxa were identified and count sums ranged between 57 and 180 head capsules per sample after amalgamation. The diagram in Fig  Performance statistics (r 2 , RMSEP, max bias, average bias) are based on Jack-knifing shows the dominant taxa. No statistically significant zones were identified, which probably can be attributed to the low variation in the dominant taxa. Therefore, we divide the data into four main zones (Ch1-4) distinguished in the cluster analysis, and two subzones representing two previously reported cooling/drying events, the Tunturi and Värriö (E2) events (Helmens et al. 2015;Plikk et al. 2016;Salonen et al. 2018) (Fig. 2).

Quantitative temperature reconstructions
The reconstructions show similar trends but differ in actual temperature values. The reconstruction based on the Norwegian calibration dataset (No) shows temperatures mainly fluctuating around 12-13°C, while the Russian (Ru) and Finnish (Fi) datasets show reconstructed temperatures fluctuating mainly around 13-14 and 14-15°C respectively (Fig. 3). The Norwegian reconstruction remains below present-day mean July temperatures of 13°C. All reconstructions show slightly lower temperatures in the lower than the upper half of the record. The initial part of the Tunturi event (2345-2315 cm) is marked in all reconstructions by a drop (Norwegian 2°C, Finnish and Russian 1.5°C) in temperature. Temperatures during this phase reach down to approximately 10°C with the Norwegian, 12°C with the Russian and 13°C with the Finnish reconstruction. Temperatures then increase slightly and fluctuate around 11.5°C (No), 13°C (Ru) and 14°C (Fi) during the remainder of the event. A gradual shift to warmer temperatures occurs after 2255 cm depth (transition to zone Ch2) in all reconstructions. The warmest phase according to all records is above 1825 cm, during the initial part of Ch4, with temperatures reaching 13.6°C in the Norwegian reconstruction and 17.1°C in the Finnish and Russian reconstructions.
The Värriö cooling event (1743-1705 cm) is characterized by a two-stepped drop in temperature with 3-4°C in all reconstructions down to 10°C (No) and 13°C (Fi, Ru). The temperatures in all reconstructions then increase to levels similar to those prior to the event before they gradually decline in the uppermost part of the record to temperatures similar to, or below, those of the lowermost part of the record [i.e. 10.5°C (No), 12.5°C (Ru) and 13.5°C (Fi)].
The MAT analyses indicate that 14 (19%) of the fossil samples have no good analogues with the Russian calibration dataset samples, 38 (52%) with the Finnish dataset and 47 (64%) with the Norwegian dataset. No fossil samples had a poor or very poor fit with temperature using the Russian calibration dataset, two samples (3%) using the Norwegian dataset and 18 (24.6%) using the Finnish dataset. The maximum percentage of fossil taxa absent from the calibration datasets varies between 7.9 (Ru) and 11.3% (No) while maximum percentage of fossil taxa rare in the calibration datasets varies between 2.9 (Fi) and 7.1% (No). The Russian and Finnish reconstructions were statistically significant when compared to 999 random reconstructions (Telford and Birks 2011).

Ecological interpretation
The majority (95-99%) of the encountered specimens are lacustrine, indicating open-water conditions throughout the Eemian. With the exception of the lower-and uppermost part of the record, the major share of the encountered taxa are 'intermediate'-or 'warm-water' taxa, indicating a warm climate across most of the Eemian sequence. Most taxa of the fossil record, and especially the most common ones such as C. ambigua, C. mancus-type and T. mendax-type, frequently occur in the calibration datasets. The nearabsence of Zalutschia, a genus indicative of lower pH (Brundin 1949;Walker et al. 1985), can be seen as an indication of the interpreted high pH of this lake resulting from the particular bedrock conditions of the Sokli site (Plikk et al. 2016). The high abundance of eutrophic C. mancus-type and to a lesser extent E. albipennis-type suggests a relatively nutrient-rich lake. In addition, the low abundance of profundal taxa indicates that the lake was relatively shallow (\ 10 m) and/or that there was a low availability of welloxygenated profundal environments (Eggermont and Heiri 2012).
The onset of zone Ch1 shows decreasing abundances of profundal and cold stenothermic taxa reflecting the transition from a deep glaciolacustrine setting to a warm relatively shallow lake. The zone is dominated by C. mancus-type, T. mendax-type and C. ambigua. The former two are meso-eutrophic warm-indicating taxa, whilst the latter is an oligomesotrophic taxon associated with clear, alkaline waters with a dense cover of submerged macrophytes (Olander et al. 1999;Brodersen and Lindegaard 1999). C. ambigua is generally considered an indicator of cold conditions but has been found to dominate the assemblage also in temperate lakes (Brodersen and Lindegaard 1999). Here and in previous studies from Sokli (Plikk et al. 2016;Shala et al. 2017) it can be associated with the occurrence of Myriophyllum. The high abundance of other taxa associated with submerged macrophytes also suggests clear water conditions during zone Ch1. Slightly before and during the Tunturi cooling event an increase in oligo-mesotrophic taxa, a decrease in chironomid concentration and the presence of semi-terrestrial taxa suggest lowered productivity possibly related to a decrease in temperatures and fluctuating water levels.
The increase in E. albipennis-type, C. cylindraceustype and Pseudochironomus in zone Ch2 suggests increasing trophic levels and warming temperatures (Brooks et al. 2007). The relatively high abundance of Chironomini larvulae in this zone is intriguing. First instar larvae have a partly planktonic behavior and small size. Increasing proportions of first instar remains in the sediments are suggested to indicate anoxic conditions in the hypolimnion (Quinlan and Smol 2001). Increased anoxic conditions would also result in the elimination of bottom-dwelling larvae (Quinlan and Smol 2001), which is in accordance with the disappearance of the profundal T. lugens-type, with relatively high oxygen demands, and the decrease in profundal S. coracina-type and Chironomus anthracinus-type. Indeed, very good preservational conditions (indicating anoxia) are reflected in a high abundance of chironomid exoskeletons with an intact body in this zone. Increasing anoxia could be related to high productivity and stable stratification and thus warm temperatures.
High abundances of warm-indicating eutrophic taxa in zone Ch3, in combination with the decrease in C. ambigua and other taxa associated with submerged macrophytes suggests increasingly eutrophic and turbid conditions. Low oxygen conditions are indicated by the relatively high proportion of Chironomini larvulae and well-preserved (complete) chironomids. The sudden decrease in C. mancus-type between 1995 and 1955 cm depth and the simultaneous increase to peak levels in taxa associated with macrophytes and shallow eutrophic conditions, suggests a sudden drop in lake levels. Following this, macrophyte-associated taxa become more common, reflecting progressive infilling and overgrowth of the lake. The marked increase in C. plumosus-type, a taxon tolerant to low oxygen conditions, towards the end of the zone is interpreted as reflecting increasingly eutrophic and anoxic conditions in the littoral following the progressive infilling.
Initially high values of C. plumosus-type in Ch4, together with eutrophic and warm indicating Einfeldia dissidens-type and Tanypus, suggest a shallow eutrophic lake with poor oxygen conditions (Brooks et al. 2007). T. lugens-type has a generally northern affinity but can also be found in the profundal of deep temperate lakes (Engels and Cwynar 2011), while in the arctic and subarctic it occurs in the littoral (Brooks et al. 2007). Since there are no signs of a sudden deepening of the lake, the increasing abundance of T. lugens-type in this zone may most likely reflect the influence of colder conditions. The two distinct peaks in C. ambigua associated with the Värriö cooling event, and the decrease in C. mancus-type indicate decreasing nutrient levels and a shift to clear water conditions. The decrease in C. plumosus-type and well-preserved chironomids indicates increased oxygen levels, while the increase in oligo-mesotrophic taxa may reflect cooler conditions. After the event, there is a slight return to conditions similar to those prior to the event, before another increase in oligomesotrophic taxa indicates a return to cooler conditions. Above 1640 cm depth, sand layers start to appear in the sequence and only fragments of chironomids are found.
The influence of site-specific non-climatic variables on the temperature reconstructions The strength of July air temperature as an explanatory variable for chironomid assemblages lies in the fact that it includes the combined effects of several ecologically important variables such as length of growing season/ice cover duration, strength of stratification, catchment vegetation, productivity/respiration rate/oxygen levels, habitat/food availability (indirect effect) and chironomid metabolism/growth rate/generation rate or emergence (direct effect), to which it is related (Eggermont and Heiri 2012; Dickson and Walker 2015). One important assumption in the transfer function approach is that the relationship between the reconstructed variable and covariables remains constant through time . For the Holocene, where temperature changes have been of relatively low amplitude, it has been shown that regional reconstructions can provide distinctly varying results due to the presence of independent variation in a site-specific factor, influencing or even overriding the climatic signal (Brooks 2006;Velle et al. 2010). Careful evaluation of the fossil record is therefore important, to ensure that independent variation in variables other than the one of interest have had negligible influence on the fossil assemblage . In this section we use the results from the previously published studies of the development of the Sokli paleolake (Plikk et al. 2016;Kylander et al. 2018), to identify periods where such variations may have occurred.
Apart from temperature, primary productivity in lakes is affected by supply of nutrients which is related to local bedrock/soil conditions, catchment processes (supply of organic matter, erosion) and lake-morphometry (Nõges 2009). Plikk et al. (2016) have shown that the Sokli paleolake probably was a naturally productive lake due to the easily weathered carbonatite and phosphorus-rich bedrock, as evidenced by high occurrences of phosphorus-favored diatoms otherwise rare in northern Fennoscandian lakes. Over time, the Sokli paleolake became increasingly eutrophic due to infilling of the lake with sediments which also forced a shift in the thermal regime of the lake (Plikk et al. 2016;Kylander et al. 2018). As a result, at around 2085 cm depth, the lake changed from a stably stratified mesotrophic and dimictic lake to an increasingly shallow, unstably stratified lake where continuous resuspension events allowed nutrients to remain in circulation and to be accessible throughout the growing season. This is also evidenced in the chironomid record by the decrease in taxa favored by clear-waters and lower trophic levels and an increase in (warm) taxa associated with more turbid and nutrient rich conditions (e.g. E. albipennistype, C. mancus-type) (Brodersen et al. 2001;Langdon et al. 2010). Therefore, above 2085 cm depth, in Ch3 and especially the initial part of Ch4, the increasingly productive conditions cannot be unequivocally related to temperature but can be expected to also reflect (non-climatic) changes in water depth.
Dissolved organic carbon (DOC) has been reported as an important environmental variable affecting chironomid assemblage composition in surface sediments of humic Fennoscandian lakes (Larocque et al. 2001;Nyman et al. 2005;Luoto 2013). DOC affects the optic properties of water, pH and hypolimnetic oxygen concentrations. As a result, increased DOC can lead both to decreasing productivity in a lake through its effect on water transparency (shading) on the one hand but also to increased oxygen consumption (through increased decomposition) on the other. Increased DOC in the lake may therefore affect the midge community in ways that are difficult to separate from the impact of climate. The multiproxy record from the Sokli paleolake indicates that pH remained high throughout the main part of the lake development (Plikk et al. 2016), and typical humic indicators (e.g. Chaoborus flavicans and Zalutschia zalutschicolatype) (Luoto 2013) are lacking from the Sokli paleolake fossil record. However, in the upper part of Ch3 and especially the initial part of Ch4 increased DOC levels are suggested by the appearance of acidophilic diatoms (Plikk et al. 2016). This increase could be related to climate-related factors as suggested in Plikk et al. (2016) where it is thought to reflect an expansion of spruce in the catchment, but it may also be related to infilling and expansion of marginal wetlands. C. plumosus-type, a chironomid taxon with relatively high temperature optima in the modern data, is found to be common in high DOC lakes with low oxygen conditions (Luoto 2013). Its high abundance during Ch4 may be associated with DOC and must therefore be interpreted with caution.

Influence of calibration dataset properties on reconstructed temperatures
A key aspect in dataset selection is similarity between the modern dataset and the fossil site regarding e.g. climate and geology, minimizing the differences in ecologically important limnological and climatological variables such as pH, nutrients and continentality and increasing the likelihood of finding good modern analogues (Birks and Seppä 2004;Engels et al. 2014). The bedrock at Sokli, influencing pH and nutrient levels in the lake, differs from the generally acidic bedrock of the sites in the Finnish and Norwegian calibration datasets. The underlying geology is more varied in the Russian calibration dataset. Nevertheless, most fossil taxa are well-represented and no dominant fossil taxa have poorly defined optima in the different calibration datasets. The selected calibration datasets span a continentality gradient which has been suggested to influence the weighted average (WA) optima of the different taxa in the datasets Engels et al. 2014). Despite the differences in continentality, taxa of the Finnish calibration dataset generally have higher optima than the taxa of the more continental Russian dataset, which indicates that other factors such as length and extent of the temperature gradient and/or distribution of lakes along the gradient has an overriding influence on the WA optima of taxa (ter Braak and Looman 1986). The Finnish dataset has a considerably higher minimum temperature (7.9°C compared to 3.5 and 1.8°C) than the Norwegian and Russian datasets, while the Norwegian dataset has a lower maximum temperature (16.0 compared to 17.1 and 18.8°C) than the Finnish and Russian calibration datasets, respectively (Table 1). If the distributions of the most common fossil taxa in the three calibration datasets are plotted on the same temperature scale, it is evident that most are either warm-water taxa that occur in ranges between about 7 and 18.5°C (e.g. Cladopelma lateralis-type, C. mancus-type, E. albipennis-type, P. nubeculosum-type, C. cylindraceus-type), or taxa that have a broad distribution with respect to temperature (e.g. P. sordidellus-type, Tanytarsini undiff). High occurrences of many of these taxa in temperatures above 16°C, which is the maximum temperature of the Norwegian calibration dataset, may indicate that optimum temperatures are slightly underestimated for these taxa in all datasets, but most so in the Norwegian dataset, due to truncation at the upper end of the temperature gradient. On the other hand, the Finnish but also the Russian datasets may overestimate optima for cold-indicating taxa, due to the very/relatively low number of lakes at the cold end of the temperature gradient (ter Braak and Looman 1986).
The WA optima of the few dominant taxa (C. ambigua, C. mancus-type, T. mendax-type, Tanytarsini undiff) are likely to influence the mean temperature of the different reconstructions. C. ambigua shows clearly defined occurrences in the Norwegian calibration dataset, with no occurrences above 13°C. The narrow temperature range of this species may be questionable considering its (high) occurrence over the range of temperatures in the two other calibration datasets, and modern observations of mass-occurrence of this species in temperate lakes (Brodersen and Lindegaard 1999). This phenomenon could possibly be attributed to it being delimited by some other confining variable in the Norwegian calibration dataset. T. mendax-type on the other hand, shows well defined occurrences in a warmer interval in the Finnish calibration dataset than in the other two sets, and may likewise be limited in lakes below 12.6°C by some unknown confining variable.
The calibration dataset properties (i.e. temperature range covered and distribution of lakes along the gradient) seem to influence the mean value of the quantitative reconstructions. We do however see similar relative trends in temperature development in all three records and the presence or absence of cool and oligotrophic taxa appears to reflect cooler respectively warmer intervals as well. The Russian reconstruction has the highest number of modern analogues and all samples show a good fit to temperature. The calibration dataset covers a temperature gradient that extends further on the warm end than the other two calibration datasets, and therefore may have the best coverage of the environmental gradient for the common fossil taxa. Despite its low minimum temperature, it may not be able to fully capture the optima of cold indicating and/or deep-water taxa due to the relatively low number of lakes at lower temperatures and the dominance of shallow lakes at these temperatures (Nazarova et al. 2015). The Sokli Eemian paleolake is however relatively shallow and deep-water taxa are rare in the assemblage. The Russian dataset may thus still be more advantageous as it may possibly capture more correctly the optima of shallow water species.
All three calibration datasets have similar r 2 and RMSEP relative to the overall temperature gradient, but the Russian calibration dataset has a lower k1/k2 (0.85), suggesting that the reconstruction may respond sensitively to confounding variables (Juggins 2013). In the section above, we have tried to highlight intervals where such variables may have affected the reconstruction.
The reconstructed temperatures based on the Russian calibration dataset are between 0 and 1.5°C lower than those of a pollen-based reconstruction based on multiple calibration methods from the same site  showing a similar temperature development, which lends support to the chironomid-based reconstruction. Differences may be related to the longer temperature range, especially towards warmer temperatures, covered by the pollen calibration dataset compared to the chironomid dataset.
Eemian temperature development at Sokli: integration of multi-proxy data Based on the above discussion on the reliability of the different reconstructions and the identification of intervals where factors other than temperature might confound our results, we shall in the following section summarize the Eemian temperature development at Sokli based on the Russian reconstruction and compare to other paleoclimate data from Sokli.
Chironomid-inferred temperatures during the early interglacial (Ch1) fluctuate markedly but show a generally increasing trend from slightly below to above present day temperatures at Sokli (13°C) (Fig. 4). The rising temperature trend follows the vegetation development from an open birch-dominated forest to a mixed boreal forest reflecting climatic amelioration (Helmens et al. 2015;Plikk et al. 2016). Warming temperatures are interrupted by the Tunturi cooling event, centered at 127.5 kyr BP, when inferred temperatures drop approximately 1.5°C, down to 12°C. Multi-proxy data from the Sokli paleolake indicate that during the event, cool and dry conditions led to lower lake levels, a more open catchment and increased soil erosion (Helmens et al. 2015;Plikk et al. 2016). The event is not marked by severe changes in the chironomid record, which may be due to the dominance of taxa with relatively wide temperature tolerances; however, there is a slight increase in oligotrophic taxa likely reflecting a climate-related decrease in productivity.
Inferred temperatures again reach values higher than 13.5°C by the end of the cooling event and increase to approximately 14°C above 2165 cm depth in Ch2. The chironomid assemblage is characterized by the appearance of warm-water taxa and a decrease in clear-water oligo-mesotrophic C. ambigua. During this phase, boreal forest returned to the area and planktonic diatoms became dominant over benthic, suggesting a prolongation of the growing season and increased productivity (Plikk et al. 2016). The formation of distinct laminations indicates stable summer stratification and hypolimnetic anoxia, and there is no obvious sign of lake-level lowering at this stage. This phase may therefore represent the optimum (warm, productive conditions) of the interglacial.
The chironomid-inferred temperatures continue to be high into Ch3 before slightly decreasing and seem unaffected by the marked changes in the chironomid assemblage associated possibly with increased (biological) turbidity (Plikk et al. 2016). Pollen-inferred temperatures show a parallel trend ). However, above 1950 cm depth chironomidinferred temperatures begin to increase. This increase should be interpreted with caution since lake infilling and increasing DOC most likely influenced the chironomid record thus obscuring the climate signal.
The chironomid-inferred temperatures fluctuate but show a declining trend in the last stage of the interglacial sequence at Sokli (Ch4) (Fig. 4). The fluctuations are paralleled by changes in the diatom community (Plikk et al. 2016) and may reflect also the impact from lake-level variations, suggesting that temporary more arid conditions led to increased mixing and oxygenation of the shallow waters. Again, the initially high temperatures associated with shallow eutrophic (anoxic) conditions should be interpreted carefully. During the Värriö cooling event, about 119 kyr BP, inferred temperatures drop to approximately 13°C, and the changes in the chironomid assemblage reflect a possibly climate-related decrease in productivity. The chironomid-inferred temperatures during the event and the remaining part of the phase stayed near present-day temperatures, similar to the pollen-inferred temperatures (Helmens et al. 2015;Salonen et al. 2018), suggesting that conditions remained relatively warm towards the end of the sequence.

Wider implications
The chironomid-inferred temperatures for the Eemian at Sokli suggest that higher than present temperatures were reached already early in the interglacial, about 130 kyr BP. Inferred optimum temperatures of approximately 1°C above present, are similar to the Leveäniemi insect record (summer temperatures 1°C above present during optimum) from northern Sweden (Lundqvist 1971) and the chironomid-inferred optimum temperatures from Denmark (0.5°C above present) (Jones et al. 2017) as well as other central European insect and pollen records suggesting optimum Eemian temperatures approximately 1-2°C warmer than present (Zagwijn 1996;Aalbersberg and Litt 1998;Kühl et al. 2007;Candy et al. 2016). The Sokli paleolake record differs from other high latitude (Arctic) insect records from Baffin Island, Greenland and NE Siberia where optimum summer temperatures for the Eemian are estimated to have been between 4 and 9°C warmer than present (Bennike and Böcher 1994;Francis et al. 2006;Axford et al. 2011;Kienast et al. 2011). A temperature increase of 1°C renders a temperature at Sokli similar to the temperatures of the warmest decennium of the twentieth century (1930s) and the Holocene thermal maximum (Seppä et al. 2009), and is perhaps surprisingly low considering the difference in insolation b Fig. 4 Chironomid-inferred mean July temperatures at Sokli, northern Finland based on the Russian calibration dataset (Nazarova et al. 2015) (c) compared to marine summer seasurface temperatures (SSTs) and relative abundance of Neogloboquadrina pachyderma (sinistral) from core PS1243 from the Nordic Seas (Zhuravleva et al. 2017) and core MD03-2665 from the Eirik Drift, subpolar North Atlantic (Irvali et al. 2016) (a, b), lake-level low stands and chironomid-inferred mean July temperatures from the Hollerup site, Denmark, interglacial (Jones et al. 2017) (d) and pollen-inferred mean July temperatures from Gröbern, northern Germany (Kühl et al. 2007) (e). Summer insolation at 65°N after Berger and Loutre (1991)  The 3-kyr delay in Fennoscandian warming (early interglacial temperatures 2-7°C lower than present) suggested from the chironomid-based temperature reconstruction from Denmark (Jones et al. 2017) is not replicated in the present study (Fig. 4). Instead, as warm as or warmer than present temperatures are inferred already from the beginning of the interglacial (Fig. 4). Although this discrepancy may be related to chronological uncertainties, an influence of confounding variables on the chironomid-based temperature reconstruction at one, or both, of the sites (e.g. the influence of elevated nutrient conditions on the chironomid assemblage in the Sokli record or of deep water conditions at the Hollerup site) must be considered. The chironomid-based reconstruction at Sokli is corroborated by the pollen-based July temperature reconstruction from the same site . No pollen-based temperature reconstruction is yet presented from Hollerup, but the vegetation development follows the same trajectory as in central Europe where optimum July temperatures are recorded during the Corylus-phase. At Hollerup, this vegetation phase is coincident with low chironomidinferred temperatures and deep-water conditions, which may be a reason for concern. The interruption of warming conditions by the Tunturi cool/arid event at Sokli, centered at 127.5 kyr BP, showing similarities to a cooling event recorded in marine records from the Nordic Seas and the subpolar North Atlantic approximately 126 kyr BP (Irvali et al. 2016;Zhuravleva et al. 2017), does however support previous indications that significant meltwater impact and interruptions in the thermohaline circulation did affect the northern North Atlantic region during the early interglacial (Govin et al. 2015;Irvali et al. 2016;Jones et al. 2017;Zhuravleva et al. 2017).
The upper part of the Sokli paleolake sequence is probably truncated and the main part of glacial inception is missing from this core but is recorded in the pollen record of a parallel core (Helmens et al. 2015). However, the appearance of cold-indicating chironomid taxa in the upper part of the present record, the fluctuating temperatures and the registration of the cool and dry Värriö event around 119 kyr BP suggest that the record contains the initial part of glacial inception. Cooling events possibly connected to rapid shifts in ocean-atmospheric circulation patterns have been recorded from marine records across the North Atlantic region at the onset of glacial inception (Irvali et al. 2016). Notably, also from Hollerup dry events have been recorded in the lower respectively upper part of the record (Jones et al. 2017) (Fig. 4). Even though it is hard to tie the Värriö event to a particular marine event due to the uncertain chronology, it is an indication that the increased climatic variability of the last glacial inception can be traced also in high latitude terrestrial records.

Conclusions
In this study we present a quantitative chironomidbased July temperature reconstruction for the Eemian interglacial. Our reconstruction indicates that the interglacial in this region was characterized by temperatures as warm as and up to 1°C warmer than today already from the beginning of the interglacial. Two cooling events, which may be related to similar events in the North Atlantic region, can be distinguished during the early interglacial and the glacial inception at 127.5 respectively 119 kyr BP. Based on comparisons with qualitative data (diatom, pollen, macrofossil, non-pollen palynomorph and xrf-data) of lake development we highlighted periods where changes in variables other than temperature might have had an influence on the temperature reconstruction. The application of different calibration datasets to the same fossil record showed the importance of calibration dataset properties, especially range of temperature gradient covered and distribution of lakes along this gradient, on mean temperatures of the reconstruction, while temperature trends for the different reconstructions were similar. Pollen-inferred temperatures from the same site are approximately 1.5°C higher than the chironomid-inferred temperatures for the main part of the sequence. This may be a consequence of differences in temperature range covered by the calibration datasets, where the pollen calibration dataset covers a longer temperature range and extends further on the warm end of the gradient than the chironomid dataset. The relatively gentle warming observed in this study during the Eemian optimum compared to other high-latitude sites is related to the general warming (background) effect of the AMOC on the Fennoscandian climate.