Meteorological validation of chironomids as a paleotemperature proxy using varved lake sediments

Study site and sampling

Nurmijärvi (61°35′N, 25°55′E; 87.7 m a.s.l.) is a boreal lake situated at the border of municipalities Sysmä and Hartola in Finland. The lake is divided into two distinct separate basins (Figure 1). Two inlets lead to the northern lake basin (max. depth 23 m) and one to the southern basin (max. depth <2 m). These inflow channels supply the lake with minerogenic material during spring floods, contributing as a typical setting for development of clastic-biogenic type of varves (Ojala et al., 2012). A major outlet stream is located in the southern basin flowing toward southeast. The catchment is characterized by boreal forest and patches of cultivated land areas. The bedrock is composed of porphyritic granites and granodiorites covered by till and fine-grained deposits. The mean July air temperature (TJul) at the study site is 16.9°C (climate normals 1981–2010).

OPEN IN VIEWER

Figure 1.

Location of the study site Nurmijärvi (61°35′N, 25°55′E) and its catchment characteristics.

The humic lake has a water color of ~58 (55–60) mg Pt L⁻¹. The lake is currently mesotrophic with summertime total phosphorus concentrations of ~20 (18–22) µg L⁻¹. The water is circumneutral with an epilimnetic pH of ~7.0 (6.7–7.3). Nurmijärvi has been introduced with pike-perch (Sander lucioperca) and northern densely rakered whitefish (Coregonus lavaretus f. pallasi).

A sediment profile was collected with a HTH-corer from the middle of the northern basin in February 2016. The core was recovered into a polyvinyl chloride (PVC) pipe that was cut lengthwise into half in the laboratory. The sediments were then subsampled at 1-cm intervals (totaling 49 samples), which usually contained several laminae, in order to gain enough chironomid head capsules for quantitative analysis. The chronology is based on repeated varve counts of four parallel cores, verified by ¹³⁷Cs determination and paleosecular characteristics of relative declination during the last centuries (Ojala et al., 2016; see age–depth model in supplementary document available online). Accordingly, counting of clastic-biogenic type of varves resulted in 250 years long chronology with chronological error estimate of ±2%.

Chironomid community analysis

Subsamples for fossil chironomid analysis were prepared applying standard methods (Brooks et al., 2007). The wet sediment was gently sieved through a 100-µm mesh, and the residue was examined using a Bogorov counting chamber under a stereomicroscope (32–40× magnification). Larval head capsules were extracted with fine forceps and mounted permanently with Euparal on microscope slides. Faunal identification was performed under a light microscope at 400× magnification. The minimum chironomid head capsule number per sample was set to 50. Identification of the chironomids was based on Brooks et al. (2007).

Principal component analysis (PCA) was used to detect the direction of the main community variance. Hence, PCA axis 1 scores were considered to reflect the chironomid response to the driving environmental variable, whereas axis 2 scores reflect secondary gradients. The PCA was run using square-root transformed relative abundances and down-weighting of rare taxa. Furthermore, constrained unweighted pair-group average (UPGMA) cluster analysis was used to group samples into local chironomid zones. Bray–Curtis similarity index was used as a measure, and the number of significant zones was assessed using the broken stick method. Temperature optima for the most common taxa in Lake Nurmijärvi were assessed using weighted averaging (WA).

Proxies of secondary gradients

Magnetic susceptibility was measured using a Bartington equipment sensor MS2C (Bartington Instruments Ltd., Witney, Oxon, UK). Loss-on-ignition (LOI) was performed from wet sediment samples, which were dried at 105°C and subsequently ignited in an oven at 550°C following Heiri et al. (2001).

Alongside chironomid analysis, fish scales and mandibles of stream insect (mayflies and caddisflies) and Chaoboridae were enumerated. These zoological macroremains are presented as a ratio to chironomids. Fish scales can be used as a proxy for food web changes (Pawłowski et al., 2016), mayflies (Ephemeroptera) and caddisflies (Trichoptera) are indicators for hydrological conditions (Luoto et al., 2013), and a Chaoborus mandible:chironomid head capsule ratio (chaob:chir) reflects state of anoxia (Quinlan and Smol, 2010). In addition to chaob:chir index, percentage of profundal chironomids, determined from taxa depth preferences in a previous Finnish dataset (Luoto, 2010), was used to indicate habitat changes in profundal areas.

Temperature calibration

The chironomid-based mean July air temperature reconstruction applying the calibration-in-space approach used the expanded Fennoscandian calibration model (weighted-averaging partial least squares, WA-PLS) combining several datasets (Luoto, 2009; Luoto et al., 2014a, 2016b; Nyman et al., 2005). The temperature gradient in the training set varies from 7.9°C to 17.6°C, including Sub-Arctic, boreal, and temperate lakes. The two-component model currently includes 180 lakes and 129 taxa having a coefficient of determination (r²_jack) of 0.86, a root mean square error of prediction (RMSEP) of 0.85°C, and a maximum bias of 0.75°C.

In the initial steps of the calibration-in-time approach, different model types considered suitable for the data – such as the modern analogue technique (MAT), WA, locally weighted-WA, partial least squares (PLS), WA-PLS, and maximum likelihood (Gaussian logit model) – were tested. The cross-validated model performance and error statistics were assessed according to the r²_jack between predicted and observed values, RMSEP, and the mean and maximum biases. The observed raw temperature data of the past ~180 years from Kaisaniemi (Helsinki) was derived from the Finnish Meteorological Institute. These data were interpolated based on the geographical temperature differences over the observational period 1981–2010 in weather stations from Kaisaniemi and Jyväskylä to obtain local temperatures for our study site. The measured data were averaged based on the laminae (years) included in each of the sediment subsamples.

Temperature reconstructions

In the temperature reconstructions, sample-specific errors (eSEP) were estimated using bootstrapping cross-validation (999 iterations). Pearson correlations and associated levels of statistical significance (p ⩽ 0.05) were used to evaluate the relationship between the chironomid-inferred temperatures using both calibration approaches and the instrumental temperature data. Proxy data form natural time series used to lengthen instrumental climatic records, and may contain a significant portion of autocorrelation. Therefore, the Durbin–Watson test was used to assess serial correlation in the residuals from the linear regressions. The Durbin–Watson statistic remains always between 0 and 4, where a value of 2 means that there is no autocorrelation in the sample. Values approaching 0 indicate positive autocorrelation, and values toward 4 indicate negative autocorrelation. In addition, time series analysis was used to test whether there is lag between chironomid response (PCA axis 1) and the observed (instrumental) temperature records. A significant relationship was considered to have p ⩽ 0.01. The correlation and time series analyses were performed using the program PASW Statistics 18.

Community variance

Since the sediment stratigraphy in Nurmijärvi consisted of annually laminated layers verified by independent dating methods and enough chironomid head capsules were found for quantitative community analysis, we were able to perform temperature reconstructions that were subsequently compared with the meteorological record. However, because of variable rate of sedimentation (i.e. varve thickness) and 1-cm sampling resolution where each sample contained five varves on average (ranging between 1 and 11), the instrumental data had to be averaged accordingly.

The chironomid assemblages in the sediment core did not seem to include many true cold indicators. Only T. lugens-type and S. coracina-type, which were found most abundantly in the lower part of the core and were absent in the more recent samples (Figure 2), are found near the cold end of the gradient in the Finnish latitudinal chironomid-temperature training set (Luoto et al., 2014a). These results are logical, since the study site is located closer to the warmer, southern part of the training set. In contrast, the biostratigraphy consisted of several species with very warm temperature preference (Brooks et al., 2007), such as Endochironomus species, G. pallens-type, and T. mendax-type, which were more common in the recent parts of the profile (Figure 2). These community patterns fit well with the observed temperature trends that are evidenced by the correlation between the instrumental record and chironomid PCA axis 1 scores (Figure 4) representing the primary direction of variance.

Although the primary long-term community variance in Nurmijärvi is clearly linked with the temperature oscillations, the latter part of the chironomid stratigraphy is surprisingly patchy (Figures 2 and 4), considering that the sample intervals only represent a time span between a year and a decade. This is also well reflected by the cluster analysis, which shows that all the significant zone borders are located in the upper sediment section representing the past 20 years (Figure 2). Previous studies have shown that the chironomid survival and distribution can be regulated by the availability of hypolimnetic oxygen (Brodersen and Quinlan, 2006). In this study, the sediment samples were derived from the hypolimnetic part of the lake, which most likely suffers from summertime oxygen deficiency allowing the preservation of laminated sediments. Therefore, the alternating presence and absence of the taxa in the upper half of the profile probably has something to do with oxygen availability, but the potential oxygen changes are difficult to track based on the chironomid assemblage composition since all the occurring taxa can be found from low oxygen lakes. There also seems to be no oxybiontic taxa present in the upper half of the record. In contrast, in the initial part of the record, the abundant and consistent occurrence of T. lugens-type suggests milder hypoxia (Millet et al., 2010). The chaob:chir values, which reflect hypolimnetic oxygen depletion (Quinlan and Smol, 2010), suggest that mild anoxia may have occurred through the record in Lake Nurmijärvi (Figure 5). During the past 20 years, the inferred anoxia appears to change concurrently with percentage of profundal chironomids (Figure 5) that also indicate that the profundal areas may have been temporarily inhabitable. It is noteworthy though that since development of visible seasonal laminae, such as in the entire examined Nurmijärvi sediment profile, requires anoxic hypolimnetic conditions, the deepwater areas have remained more or less stable throughout the record. In any case, changes in oxygen conditions and subsequent reflectance on chironomid communities do not mean that the chironomid-temperature signal would necessarily be compromised. A previous study from southern Finland (Luoto and Salonen, 2010) showed that the oxygen conditions in two boreal lakes with different limnological conditions (eutrophic vs oligotrophic) oscillated according to changes in air temperatures. Consequently, higher hypolimnetic oxygen values were found during cold climate periods, such as the ‘Little Ice Age’, whereas lower oxygen availability was associated with warm episodes, including the Medieval Climate Anomaly and the present warming. These results are similar to the current chironomid assemblage succession, suggesting that hypoxia has increased in the lake alongside climate warming during the past century.

In comparison to chironomid-temperature optima between the CiT and CiS datasets, the largest offset was found in Thienemannimyia, Procladius, T. lugens-type, and S. coracina-type that had >3°C lower optima in the CiS dataset (Figure 3). Thienemannimyia and Procladius are free-living and can tolerate wide temperature ranges at the genus level that probably explains the offset, but the differences in estimated temperature optima of T. lugens-type and S. coracina-type may be related to lake depth. Both T. lugens-type and S. coracina-type live in the profundal and are typical cold-indicating taxa (Brooks et al., 2007). Since Lake Nurmijärvi is deeper than most of the lakes in the CiS training set, its profundal area provides a colder habitat for these taxa compared with the shallower lakes of the CiS dataset at the same air temperature regime.

When setting the results into a context of direct human-induced catchment changes, there seems to be no clear response in specific chironomid taxa (Figure 2) to the known human activities, such as the artificial lowering of the lake level related to saw-mill activity in the 1920s and 1930s (Markkanen, 1968). Nevertheless, there is a concurrent change during the 1930s in PCA axis 1 scores and magnetic susceptibility, which reflects catchment originated changes possibly related to operation of the saw-mill. In the later part of the core, there is a synchronous increase in sediment organic matter and changes in chironomid PCA axis 1 and 2 scores following the 1990s (Figure 5). There is a commonly observed eutrophication-like response of aquatic communities from temperate (Manca and DeMott, 2009; Visconti et al., 2008) to arctic (Luoto et al., 2014c, 2015, 2016a) lakes superimposed on climate change. Based on the chironomid biostratigraphy and simultaneous increases in the share of organic matter and measured air temperatures, this phenomenon appears to be operating also in Nurmijärvi. In addition, recent fish introductions, reflected by the increase in fish scales (Figure 5), may have caused changes in the chironomid communities as suggested by PCA axis 1 scores. Since the number of stream insects (mayflies and caddisflies) remained very low during the entire examined time period (Figure 5), it appears that hydrological changes (stream conditions) close to the sampling site have remained small. Considering all the changes that have occurred in potential secondary environmental factors, it is apparent that most significant confounding influences on chironomid-based temperature reconstructions have occurred in the most recent sediment samples. Therefore, the inferences of the past ~20 years should be considered with caution.

The results also indicated that there is a slight statistically significant lag in the chironomid response to air temperature change. This may be related to the length of the chironomid lifecycle, since though multivoltine, bivoltine, and univoltine lifehistories are most common, many northern species require up to 7 years to complete their development (Pinder, 1986). The observed lagtime of ~4–8 years is nonetheless much shorter than with pollen, for example, which may lag behind several decades because of soil forming processes. It should also be noted that the lag appears not to be visually consistent through the record (Figure 4) and there is a possibility that the lag is partly related to uncertainties in the chronology (chronological error estimate ±2%).

Temperature calibration and reconstructions

The close relationship between chironomid stratigraphy and temperature development allowed the construction of a CiT model. The CiT model is based on less sites compared with the latitudinal CiS model and it also has shorter temperature gradient and contains significantly less taxa. The temperature optima of the most common chironomid taxa in the Lake Nurmijärvi record showed similar characteristics between the CiT and CiS datasets, though there was a systematic trend where the CiT optima were generally ~2°C warmer (Figure 3). This is related to the temperature gradients of the datasets, since in the CiT, the lowest measured temperatures reach only 13.9°C compared with 7.9°C in the CiS dataset. When considering this offset in a downcore temperature reconstruction, it is probable that the dataset would infer similar temperate trends but the CiT would provide lower values than the CiS approach. Based on the descriptive statistics of the models (Table 1), the CiT model has lower RMSEP, whereas the CiS model outperforms in its r²_jack and maximum bias. These results describe a rather typical pattern, since the prediction power of a good model performance is usually related to the gradient length and the number of included samples and taxa, while lower modeling errors can originate from the opposite characteristics (Olander et al., 1999). Also in a previous comparison between chironomid-based CiS and CiT models from Switzerland (Larocque-Tobler et al., 2011), the prediction power of the spatial calibration was higher but the temporal calibration had improved prediction error, indicating that these differences are inherent. In the current data, when the RMSEPs of the CiS and CiT models are compared relative to the training set temperature gradients, they come out very similar (Table 1), suggesting that the gradient is indeed in a major role determining the modeling errors.

In the reconstructions of air temperatures using the CiT and CiS, both approaches showed correlation with the instrumental record and also with the chironomid PCA axis 1 scores. In case of the CiT, these relationships are of course intrinsic. When comparing the models’ output, the CiT model appears to be more accurate approach to reconstruct paleotemperatures from the Nurmijärvi shortcore, since it follows more closely the instrumental temperature trend over the past ~180 years (Figure 4). In the lower part of the profile, before 1930 CE, the CiS approach underestimates the temperature values, and in the samples from 1940s onwards, the values show more variability than expected based on the measured data (Figure 4). These features are not present in the CiT reconstruction and may be related to issues dealing with the modeling technique. The tendency of WA-PLS models to yield more extreme reconstructions in comparison with simple WA is known from previous studies (Toivonen et al., 2001; Walker et al., 1997). This phenomenon appears to be more common in low-diversity samples (Luoto et al., 2010). Although the ability of WA-PLS models to extrapolate the results can be beneficial when working with samples close to or beyond the training set environmental gradient, such as often during the early-Holocene or other extreme climate periods, it may also cause unrealistic variability to a reconstruction, as demonstrated by the present results (Figure 4). Likewise, the applicability of the CiT model, based on simple WA, can be compromised if the fossil chironomid communities drastically change earlier during the Holocene compared with the communities prevailing during the calibration period. In such case, the CiS model having larger taxa coverage and potentially better modern analogues together with the ability to extrapolate beyond the calibration set temperature gradient probably outcompetes the CiT model in its reliability.

Consistent with previous findings (Larocque-Tobler et al., 2011), the present results suggest that the CiT approach is probably most suitable for shorter time intervals, whereas the regional CiS training sets can be more applicable for Holocene scale reconstructions. For periods dating beyond the Holocene, continental CiS training sets (Heiri et al., 2011; Lotter et al., 1999) and training sets that take the influence of continentality into account (Engels et al., 2014; Self et al., 2011) are probably the most useful. Based on the findings, it can justified to use both CiT and CiS approaches, though long-term and large-scale changes should be considered with caution in CiT, whereas the CiS most likely require regression-based smoothing, such as LOESS.