Unravelling chironomid biodiversity response to climate change in subarctic lakes across temporal and spatial scales

We combined paleolimnological reconstructions and space-for-time substitutions to unravel chironomid biodiversity responses to climate change in subarctic mountains across temporal and spatial scales. Using sediment records, we found that long-term temporal changes in chironomid taxonomic diversity were mainly induced by the temperature tolerance/optimum of species, while little changes in functional diversity were found due to the replacement of similar functional-type taxa within the community. Overall, paleolimnological reconstructions suggested the selection of larger chironomid species by long-term climate cooling and little changes in trophic guilds. Space-for-time substitutions showed, however that low-elevation lakes with forested have more sediment-feeding taxa and larger larvae than high-elevation lakes, thus, suggesting the selection of large chironomid morphotypes with a sediment-feeding mode under warmer climate. Space-for-time substitutions and paleolimnological reconstructions, therefore, gave contrasting results for the link between climate and functional diversity of chironomid communities, likely because space-for-time substitutions failed to match the extent of both spatial and temporal climatic gradients. We suggest that future studies must address biodiversity issues across both temporal and spatial scales as an improved understanding of biodiversity responses to climate change may help us to understand how biodiversity will be affected by ongoing and future change.


Introduction
Subarctic freshwater ecosystems are exposed to accelerated warming, which potentially drives widespread ecological changes and dramatic biodiversity shifts (Smol et al., 2005), especially among aquatic macroinvertebrates that are essential for appropriate ecosystem functioning (Koltz et al., 2018). Chironomidae (Arthropoda; Diptera; Nematocera) are non-biting midges whose larvae constitute a predominant group in these ecosystems (Northington et al., 2010;Lento et al., 2019), thus, providing an ideal taxonomic group for studying biodiversity patterns in subarctic lakes. Predicting climate change effects on chironomid biodiversity remains, however, an elusive challenge for ecologists worldwide.
Abstract We combined paleolimnological reconstructions and space-for-time substitutions to unravel chironomid biodiversity responses to climate change in subarctic mountains across temporal and spatial scales. Using sediment records, we found that longterm temporal changes in chironomid taxonomic diversity were mainly induced by the temperature tolerance/optimum of species, while little changes in functional diversity were found due to the replacement of similar functional-type taxa within the community. Overall, paleolimnological reconstructions suggested the selection of larger chironomid species by long-term climate cooling and little changes in trophic guilds. Space-for-time substitutions showed, however that low-elevation lakes with forested have more sediment-feeding taxa and larger larvae than high-elevation lakes, thus, suggesting the selection of large chironomid morphotypes with a sedimentfeeding mode under warmer climate. Space-for-time substitutions and paleolimnological reconstructions, Due to the lack of long-term monitoring, spacefor-time approaches have frequently been used to identify key environmental drivers and analyse climate change effects on northern lake and river ecosystems. A space-for-time approach addresses how species are distributed across areas in landscape and relate to specific habitats and assumes that drivers of spatial gradients of species composition also drive temporal changes in diversity. However, space-fortime studies often fail to match the extent of both spatial and temporal climatic gradients, as spatial variation in climate often exceeds its temporal variation (Blois et al., 2013). Furthermore, space for time assumes stationarity in environmental conditions, but climate change is controlled by the concerted action of numerous forcing factors that operate on different timescales, thus, inducing different patterns (including warming and cooling periods) and frequencies in climate (Roberts, 2013). Climate change is, therefore, a highly non-stationary phenomenon, and temporal variability in climate and environmental conditions lead to continuous turnover of species that potentially interacts with spatial species turnover (Stegen et al., 2013).
Time also represents one of the fundamental axes that shapes ecological systems (Ryo et al., 2019), but most studies have, however, attempt to identify the mechanisms that generate biodiversity patterns over short timescales. Paleolimnological approaches are an efficient tool for reconstructing the effects of past climate change on lake ecosystems (Smol, 1992) and provide excellent opportunities to address mechanisms that generate species turnover patterns over different timescales (e.g. Castillo-Escrivà et al., 2020). Specifically, chironomid larval head capsules are well preserved in sediments which allows for the identification to morphotype level (Walker, 2001). Indeed, analysis of subfossil chironomid communities has provided valuable insights into the impacts of past environmental and climatic changes on aquatic biodiversity (Massaferro & Corley, 1998;Brodersen et al., 2001;Wohlfarth et al., 2018). Most paleolimnological studies, however, have applied the traditional taxonomy-based approach, thereby largely disregarding functional biodiversity based on species´ traits. Understanding how biodiversity has been affected by climate change in the past may also help us to understand the future effects of climate change on ecosystem functioning.
Several studies have demonstrated that spatial and temporal theories should not be treated separately, and a deeper understanding of the mechanisms underlying community assembly and biodiversity dynamics can be gained by considering combinations of these theories (Stigall, 2008;Rull, 2014). Engels et al. (2020) analysed the main drivers of chironomid biodiversity in the Northern Hemisphere by combining spatial training sets (consisting of 837 lakes located along broad climatic and environmental gradients) and numerous long-term sediment records (covering the last glacial-interglacial transition: when summer temperature changed by up to 6ºC), thus, representing the largest available metanalysis on this topic. Engels et al. (2020) confirmed the strong association between summer temperature and changes in chironomid community compositions, as frequently reported in literature (Larocque et al., 2001;Rees et al., 2008;Martel-Cea et al., 2021). However, Engels et al. (2020) highlighted more complex ecological responses for individual sites and for small amplitude change in temperature, which could likely suggest that key aspects of chironomid biodiversity were not identified and considered.
Biodiversity is composed of multiple, interrelated facets reflecting both the taxonomic composition of species, their roles in ecosystems and trophic interactions, a.k.a. taxonomic, functional and trophic diversity (Heino & Tolonen, 2017). In recent years, numerous methods have been developed to introduce these interrelated facets of biodiversity into ecology (Serra et al., 2017;Belle & Cabana, 2020). Specifically, trait-based approaches use morphological, and phenological features, with each species exhibiting a set of multiple functional traits related to their biological and physiological properties (Usseglio-Polatera et al., 2000;Poff et al., 2006). Studies of chironomid functional composition have, therefore, the potential to elucidate the underlying mechanisms of spatiotemporal patterns in biodiversity that would not be discernible from the conventional taxonomic-based approach (Stivrins et al., 2021). Few pioneering studies have investigated the relationships between both functional and taxonomic chironomid diversity and environmental changes (Luoto & Nevaleinen, 2015;Motta & Massaferro, 2019), but these studies only considered a single functional trait (i.e. feeding mode) and a unique trait category was assigned to each taxa. The issues of using binary codes for traits has already been widely discussed in the literature (see Degen et al., 2018).
In this study, we combine different spatial and temporal scales to unravel chironomid biodiversity responses to climate change in subarctic lakes. First, we quantified the role of climate-induced changes on temporal patterns in taxonomic and functional chironomid composition of a subarctic lake during the Late-Holocene. Second, we used results from Belle & Goedkoop (2020) who examined the functional diversity of chironomid communities in lakes across a 1000 m-elevation gradient, reflecting gradual spatial changes in temperature and landscape characteristics. Finally, we explored complementarity of these two approaches to study climate change effects on chironomid biodiversity in subarctic lakes.

Study site for paleolimnological reconstruction
For the paleolimnological approach, we used Lake Diktar Erik (68°26′43″N, 18°36′50″E), a small lake (0.1 km 2 ) located in northern Sweden (Fig. 1). The lake is located at 375 m a.s.l. and has a maximum water depth of 16 m (Secchi depth of 6 m). Previous studies from the area have shown climate-induced changes of the surroundings of Lake Diktar Erik during the Holocene (Meyer-Jacob et al., 2017), which provides an exceptional natural experiment to study how climate and landscape changes could affect chironomid biodiversity in a subarctic lake. A 100-cm sediment core retrieved from the deepest part of the lake, and age-depth model showed that the core covered the last ca. 9500 years, and the model was already discussed in Belle et al. (2019). Briefly, the deglaciation of the surroundings of Lake Diktar Erik started at ca. 9500 cal. BP and allowed an initial phase of landscape development characterized by a rapid transition from a vegetation-free landscape to forest vegetation during early Holocene (until ca. 6500 cal. BP Meyer-Jacob et al., 2017). We, therefore, decided to focus on the upper section of the core covering approximately 6000 years (or 49 cm of sediment), when the surrounding landscape was assumed to be in a balance with prevailing climatic conditions.
Sedimentary proxies, such as organic matter concentration, carbon and nitrogen stable isotopes, carbon and nitrogen concentrations, and C/N weight ratios, were analysed to reflect both terrestrial and aquatic processes. A principal component analyses (PCA) was performed on sediment data (already published in Belle et al., 2019), and the first PCaxis scores are used as a marker of climate-induced Fig. 1 Location of the study region showing the sampling site (black circle; Lake Diktar Erik), and elevation map over the study region showing the major lakes, streams and rivers changes in organic matter dynamics, one of the key drivers of changes in chironomid communities.

Subfossil chironomid analysis
The 49 cm core was subsampled at a 1 cm resolution. Chironomid head capsules were hand sorted from samples of wet sediment (2 cm 3 ) following Walker (2001) and identified under a microscope using Brooks et al (2007) and Rieradevall & Brooks (2001). Data are expressed as relative abundances, and only taxa occurring in at least two samples, with a maximum relative abundance of more than 2%, were included in further analysis. For further functional analyses (see details below), chironomid counts were also pooled to the genus level to fit the taxonomic resolution of the European Chironomidae Database (Serra et al., 2016). A fuzzy coding approach was used to reflect affinities of each genus to each trait category, an approach that has the advantage to acknowledge variability/diversity within each taxonomic unit and at different life stages (Chevenet et al., 1994;Bonada et al., 2007). Furthermore, the fuzzy coding approach is widely used in trait-based ecological research (Castro-Català et al., 2020) in which traits are usually coded at the genus level (Sarremejane et al., 2020).

Statistical analysis applied on paleolimnological data
Constrained hierarchical cluster analysis using a Bray-Curtis distance and a CONISS linkage method ("Rioja" package for R; Juggins, 2020) was used to identify periods of uniform chironomid communities. The number of zones that significantly differed from each other were calculated using the brokenstick model (Bennett, 1996). Detrended correspondence analysis was performed on the chironomid dataset to estimate among-sample variability and to select the appropriate ordination method, resulting in a linear approach (Legendre & Birks, 2012). A PCA was, therefore, performed on chironomid data and PC-axis significance was checked using the broken-stick model (Bennett, 1996). The functional diversity of chironomid communities was characterized using six life-history traits (i.e. 4th instar larval size, emergence season, voltinism and flight period, food preference and feeding habits), also used for the space-for-time study (see below, Table 1).
These six traits are considered responsive to climate change (Bonada et al., 2007) and represent biological features that directly influence ecosystem functions (Hulot et al., 2000;Hevia et al., 2017). Trait categories were standardized to sum 1 for each trait (depicted as the standardized trait profile following a fuzzy coding procedure, Chevenet et al., 1994), ensuring that all taxa were weighted equally in further analyses. A trait-by-sample dataset was computed as the cross-product between the relative abundance of chironomid genera and the standardized trait composition, thus, containing the relative abundance of each trait category in each sediment layer (n = 49). A Fuzzy Correspondence Analysis (FCA) was then performed to compare samples over time based on the trait composition of their chironomid communities (i.e. trait-by-samples dataset: 49 sediment samples × 30 categories of 6 traits). Cluster analysis, based on a matrix of Bray-Curtis distances calculated on lake coordinates along the first four FC-axes, was applied to identify lake clusters with similar trait composition, and samples were then grouped into successively larger groups using the CONISS method and the broken-stick model for selecting significant clusters.

Space-for-time substitution dataset
We re-used the published training set by Larocque et al. (2001) that describes the chironomid assemblage composition of 100 lakes in the northern boreal and Arctic/alpine ecoregions of Sweden to evaluate how spatial variation in climate affects chironomid biodiversity. Using this dataset, Belle & Goedkoop (2020) also examined the functional diversity of the chironomid communities in subarctic lakes across a 1000-m elevation gradient, reflecting gradual changes in temperature and landscape characteristics. Briefly, the functional diversity of chironomid communities was also characterized using six life-history traits (i.e. 4th instar larval size, emergence season, voltinism and flight period, food preference and feeding habits; Table 1). A FCA was also performed on community-weighted trait composition and multivariate analyses (redundancy analysis) were then performed to establish a model linking chironomid trait composition to environmental conditions.

Paleolimnological reconstruction
A total of 2786 chironomid head capsules were retrieved from the 49 sediment layers, with counts ranging from 28 to 111 HC per sample, and only four samples contained less than 50 head capsules. 13 morphotypes met the minimum frequency criteria for inclusion in further statistical analysis. Across the core, two main contrasting trends occurred between the two morphotypes with the  (Fig. 2). Furthermore, three statistically significant zones (Chi-Z1-Z3) were revealed based temporal patterns in chironomid communities (Fig. 2). The lowermost chironomid zone (Chi-Z1; from 5750 to 4720 cal. BP or from 50 to 41 cm) was predominated by Micropsectra insignilobus type (up to 60%), but its relative abundance decreased to 20% towards the upper part of the zone, while Heterotrissocladius maeaeri type taxa became more abundant (up to 40%). The next chironomid zone (Chi-Z2; from 4720 to 1990 cal. BP or from 41 to 17 cm) exhibited a steady predominance of taxa such as Micropsectra insignilobus type (up to 22%) and Heterotrissocladius maeaeri type (up to 60%), while an overall decrease in Psectrocladius septentrionalis type and Heterotrissocladius marcidus type was also found. This period also showed low relative abundances of Procladius, Microtendipes pedellus type and Paratanytarsus. The uppermost zone (Chi-Z3; from 1990 cal. BP to present or from 17 cm to the top) exhibited a predominance of Heterotrissocladius maeaeri type (Fig. 2), while the percentages of Heterotanytarsus increased and those of Micropsectra insignilobus type decreased. The first two PC-axes applied to chironomid data accounted for 17.3% and 15.8% of the total variance, respectively (Fig. 3A). PC1-axis explained the distribution of Psectrocladius septentrionalis type and Micropsectra insignilobus type and to a lesser extent that of Heterotrissocladius maeaeri type, with positive values on the PC1-axis correlating with high relative abundances of the former two taxa and low relative abundances of H. maeaeri type (Fig. 3B). The PC2-axis largely reflected the predominance of Psectrocladius sordidellus type and Heterotanytarsus, Fig. 2 Relative proportions of the most abundant chironomid taxa in the sediments of Lake Diktar Erik. Taxa labelled in grey lacked sufficient information on functional traits in the European Chironomidae Database, and were, therefore, excluded from further analyses. The dendrogram is based on the chironomid data and constructed by hierarchical cluster-ing analysis (Bray-Curtis distance, CONISS linkage method). Horizontal dashed lines indicate significant zones defined by major patterns in chironomid communities. Chironomid zones corresponding to time: Chi-Z1 = 5750-4720 cal. BP; Chi-Z2 = 4720-1990 cal. BP; Chi-Z3 = 1990 cal. BP-present with positive values representing samples with higher percentages of these taxa.
An FCA was applied to the chironomid trait compositions of the Lake Diktar Erik sediment core. Some taxa, such as Heterotanytarsus, were not included in the analysis as sufficient information was missing from the Chironomid database, and their inclusions would probably have enhanced overserved patterns during the most recent period. Overall, very low correlation ratios suggest low functional turnover in chironomid communities. The traits "4th instar larval size" (from SIZE1 to SIZE4, see Table 1 for full names of trait categories) showed the highest correlation ratios with the F1-axis (0.051), revealing that some trait categories were well separated along this axis (Fig. 4A). Along the F2-axis, some categories of the traits "Food preferences" (0.06) showed the highest, although still very low, correlation ratios (i.e. higher for the F2-axis than those for the F1-axis; Fig. 4A). Samples were mainly distributed along F1-axis, and as the trait "4th instar larval size" was separated along FCA axes, these results suggested temporal changes in size of chironomid communities towards smaller larval size, but very little changes in food preferences over time (Fig. 5). The CONISS linkage method revealed two significantly distinct groups of samples (Fig. 4B): Trai-Z1; from 5750 to 4720 cal. BP or from 50 to 41 cm; Trai-Z2; from 4720 cal. BP to present day or from 41 cm sediment depth to the top. Size trait (small larvae: SIZE1 and SIZE2) was more likely found for communities belonging of the Trai-Z1 zone, whereas intermediate size larvae (SIZE3 and SIZE4) was characteristic for the uppermost trait zone (Trai-Z2; Fig. 4A and B).

Space-for-time substitution
The RDA model that linked trait categories to environmental variables explained 16.9% of the total variation in chironomid trait compositions of lakes. Most of the variance contained in both the first and second RDA axes is described by elevation and TOC (e.g. low elevation and humic sites inducing positive RDA1 and negative RDA2 scores; Fig. 6), suggesting a strong influence of the elevation/humic gradient on chironomid trait compositions. Overall, the traits "4 th instar larval size", "Food preference" and "Feeding habits" were well separated along RDA axes, and chironomid communities were mainly composed of large-bodied larvae belonging to the shredder and sediment-feeding trophic guilds at low-elevation sites, whereas chironomid larvae in higher elevation lakes tended to be smaller and commonly had predatory feeding mode (Fig. 6).

Discussion
Taxonomic changes in chironomid communities of the Lake Diktar Erik sediment core appeared to be mainly induced by temperature tolerance/optimum of species (Larocque et al., 2001;Heiri et al., 2011), with decreasing trends of warm-water taxa Psectrocladius septentrionalis type, Micropsectra insignilobus type and increasing proportion of the intermediate temperature-taxon Heterotrissocladius maeaeri type. These findings confirm the long-term cooling previously reported from the area (Fig. 5; Barnekow & Sandgren, 2001;Seppä et al., 2009) and also show good agreement with the reconstructions by Larocque & Hall (2004) who produced chironomid-based temperature inferences from lakes in the area. However, the most recent large taxonomic shift towards a predominance of Heterotanytarsus (observed in PC2 instead of PC1; Fig. 5) and concurrent changes in sedimentary proxies (see "PC1 Sediments" in Fig. 5) have not previously been reported in the area (Larocque & Hall, 2003) and could suggest the influence of processes that occurred after the sample collection of earlier studies (i.e. late 1990s). Unlike other taxa, Heterotanytarsus is typically found in boreal lakes with high DOC concentrations associated with terrestrial organic matter inputs (Medeiros & Quinlan, 2011), and the observed, intriguing taxonomic change in the recent years could suggest unprecedented ecological changes that needs to be  (Bigler et al., 2002). Symbols represent significant zones previously reported using hierarchical clustering analysis confirmed in the future. Additional environmental conditions (such as brownification and tree-line dynamics) should, therefore, contribute to shape chironomid communities, thus, potentially confounding the well-established association between chironomid composition and climate. Our results, therefore, stress ecological processes that drive biodiversity dynamics are likely multidimensional and involve several factors operating in concert with climate change. Species trait frameworks have already been used to identify complex driving factors of chironomid communities (Serra et al. 2017;Stivrins et al. 2021). Due to the scarcity of reliable ecological information at the morphotype level, we had to aggregate some morphotypes (see "Subfossil chironomid analysis" section) and use fuzzy coding approach to incorporate intra-taxonomic variability and account for trait affinity differences within a genus (Bonada et al., 2007). Using fuzzy correspondence analyses, we investigated temporal variability in trait composition of chironomid communities of Lake Diktar Erik during the last 5750 years. The fact that not all taxa were included in the trait database, especially the absence of Heterotanytarsus for which feeding preference information was missing (Serra et al., 2016), could have induced significant bias in ecological interpretations of our observed reported pattern during the most recent period. Overall, we observed little temporal changes in functional diversity of chironomid communities (see temporal patterns of F2 in Fig. 5), which suggests a large degree of functional redundancy (Mayfield et al., 2021). The first FA-axis mainly represented temporal changes in larval size, whereas the second FA-axis mainly represented their feeding habits (Fig. 4). Overall, warmer conditions favoured the selection of smaller larvae, and no specific trend in feeding habit was observed throughout the record (Fig. 5). Our results, therefore, suggest that climate cooling (and associated changes) observed after the Holocene Thermal Maximum tended to favour the selection of larger chironomid taxa but had little effect on their feeding habits and food preference in Lake Diktar Erik. This finding of a selective pressure towards larger chironomid species upon climate cooling with little changes in trophic guilds conflicts with previous findings from other studies that report that small larval size and a predatory feeding mode are favoured in cold environments ( Fig. 6; Saether, 1979;Nolte & Hoffmann, 1992;Belle & Goedkoop, 2020). As previously reported by Blois et al. (2013), space-for-time substitutions often cover too large ecological gradients (in this case encompassing the above/below tree limit, up to 7.7 °C range temperature, etc. Belle & Goedkoop, 2020) which might not be relevant for inferring small changes occurring at finer temporal scales.

Conclusions
In this study, we used a combination of paleolimnological reconstructions and space-for-time substitutions to unravel chironomid biodiversity responses to climate change in subarctic lakes across temporal and spatial scales. The analysis of the Lake Diktar Erik sediment core suggested the selection of larger chironomid taxa by long-term climate cooling with only marginal changes in trophic guilds. This contradicts findings from previous space-for-time investigations that instead reported larger larvae and more sediment-feeding taxa in low-elevation lakes with forested catchments than those above the tree line. Our space-for-time substitutions and paleolimnological reconstructions gave contrasting results for the link between climate and functional diversity of chironomid communities, likely because natural variation in climate cannot be easily captured by space-for-time substitutions. Furthermore, the observed, ongoing climate-induced changes induced a shift towards a new chironomid assemblage composition that have not been seen in the recent past.
Our study suggests that none of the two approaches alone can provide reliable insights into biodiversity responses to climate change. This implies that synergies between existing theories of biodiversity are needed to better include the full spectrum of ecological processes that drive continuous and multidimensional biodiversity patterns. Combining landscape paleolimnology and space-for-time substitutions will, therefore, not only strengthen our understanding of how biodiversity is/was affected by climate change, but also help us to understand how it would be affected by future climate change.
Author contributions SB, TV and WG designed the study. SB and FK analysed samples. SB wrote the paper with substantial contribution from all co-authors. Data will be available upon reasonable request.
Funding Open access funding provided by Swedish University of Agricultural Sciences.

Conflict of interest
The authors declare that they have no competing interests. Data will be available upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.