Functional diversity of chironomid communities in subarctic lakes across gradients in temperature and catchment characteristics

Northern ecosystems are experiencing rapid and large-scale changes driven by accelerated warming, which have profound effects on the terrestrial and freshwater biodiversity. A comprehensive understanding of the distribution of aquatic biodiversity of subarctic ecosystems is therefore needed to better predict future trajectories of their unique biodiversity. In this study, we examined the functional diversity of chironomid communities in subarctic lakes across a 1000 m-elevation gradient, reflecting gradual changes in temperature and landscape characteristics. Using fuzzy correspondence analyses, we investigated spatial variability in trait composition of chironomid communities from 100 lakes in northern Sweden, and tested the hypotheses that (1) climate directly and indirectly shapes chironomid trait composition across the studied gradient, and (2) that generalist taxa with smaller body size and broader food preferences are more able to persist in cold environments. Our results showed that complex interplays between direct (e.g. temperature) and indirect climate processes (e.g. elevation-driven changes in vegetation/habitats) affect the functional diversity of chironomid communities. Specifically, traits such as larval size, food preference and feeding habits were well separated along the gradient, and this pattern revealed that low elevation lakes with forested catchments tended to have more sediment-feeding taxa and larger larvae than those above the tree line. As expected, food resource availability in lakes is strongly linked to vegetation composition/cover, and traits related to resource exploitation in chironomid communities are therefore well constrained by landscape characteristics. Furthermore, our findings suggested that short life cycles could facilitate the development of viable population in northern and high-elevation lakes where the short ice-free period is a limiting factor, thus contradicting patterns showing smaller organisms in warmer environments reported for other invertebrates. As a consequence of climate warming, the highest elevation lakes in subarctic landscapes will likely lose their typical cold-adapted chironomid taxa along with their functional attributes leading to potential impacts on the food web structure and the overall functioning of northern lake ecosystems.


Introduction
Patterns in biodiversity (i.e. taxonomic composition and functional attributes) are determined by a combination of large-scale landscape attributes and habitat properties (Townsend et al. 2003), among which climate plays a central role. Northern ecosystems are experiencing rapid and large-scale changes driven by accelerated warming (AMAP 2011), and climate change is already responsible for ecological shifts observed on freshwater biodiversity (Heino et al. 2009;Lento et al. 2019;Lau et al. 2020). Also, the recent Arctic Biodiversity Assessment by the Conservation of Arctic Flora and Fauna (CAFF) has identified climate change as the most significant threat to Arctic biodiversity and pointed out subarctic freshwaters as hotspots for ecological change (Lento et al. 2019). A comprehensive understanding of the effects of climate warming on structural and functional diversity of subarctic ecosystems is therefore needed to better predict future trajectories of their unique biodiversity.

3
Functional traits reflect how species interact with their environment and with co-occurring species (Dıáz and Cabido 2001) and can be used to estimate the functional diversity of communities. The effects of environmental changes on ecosystem functions are partly set by the functional diversity in a community (Hulot et al. 2000;Reynolds et al. 2002;Hooper et al. 2005), and trait composition of communities therefore helps to define ecological resilience (McLean et al. 2019). Although functional diversity responses to changing environmental conditions have been well documented for temperate lakes (see also Schmera et al. 2017), we are still lacking information about their effects on ecosystem properties of subarctic lakes.
Reponses of species distribution to environmental changes are mainly controlled by the ability of species to disperse and to persist in suitable environments (Lavergne et al. 2010). For example, short generation times (represented by traits such as multivoltinism and small body size) may facilitate dispersal, while the ability to exploit available food resources (e.g. feeding preference, food type) influences the establishment and persistence of viable populations (Urban et al. 2012). Accordingly, generalist species with smaller body size and broader food preferences should be more likely to find suitable habitats in extreme environments (Angert et al. 2011). In freshwater ecosystems, invertebrate trait composition responses to climate change are still debated, and most studies have considered only a few traits (see also Gardner et al. 2011). According to Bergmann's rule, smaller invertebrates should be found in warmer environments (see also Daufresne et al. 2009;Hayden et al. 2017), but this rule might not be valid for all taxonomic groups of invertebrates and/or ecosystem types (see also MacLean and Beissinger 2017). Furthermore, traits related to resource exploitation have been largely neglected in this framework but should be considered as they could supplement our knowledge of community trait responses to climate change.
Chironomidae (Arthropoda; Diptera; Nematocera) are non-biting midges whose larvae constitute a predominant group in northern freshwater ecosystems (Northington et al. 2010;Lento et al. 2019), thus providing an ideal model for studying biodiversity patterns in subarctic environments. Chironomid species are well studied and have been good indicators of environmental change, as they have well defined ecological requirements across environmental gradients, such as temperature, nutrient and oxygen concentrations (Brodersen and Quinlan 2006;Verbruggen et al. 2011;Heiri et al. 2011). Furthermore, chironomid species exhibit distinct functional traits related to their biological and physiological properties (Usseglio-Polatera et al. 2000;Serra et al. 2016). As a single taxon may have several traits, minor changes in taxonomic composition of chironomid assemblage can shift in assemblage trait composition and ecosystem functioning (c.f. van Kleef et al. 2015;Serra et al. 2017;Desrosiers et al. 2019). Hence, novel functional (rather than taxonomic) approaches must be developed to explore functional diversity changes of chironomid communities across environmental gradients (Motta and Massaferro 2019).
Due to the lack of long-term monitoring in most ecosystems, a comprehensive understanding of the effects of climate change on the functional diversity of chironomid communities of northern ecosystems is lacking, but spacefor-time substitutions could help to overcome this problem. Relationships between the functional attributes of chironomid communities and environmental conditions can be studied using calibration training data sets that describe the spatial distribution and abundance of taxa across a temperature gradient (Eggermont and Heiri 2012). Training sets are produced by collecting sediment samples from a large number of lakes across latitudinal or elevation gradients and by identifying the chironomid exoskeleton remains (i.e. head capsules) in these sediments, and provide insight into the composition of the contemporary chironomid communities of these lakes . For this purpose, surface sediments are usually retrieved from the largest depth of lakes where chironomid remains from the entire lake accumulate (Heiri 2004;van Hardenbroek et al. 2010). Training sets could thus give valuable insights into biodiversity patterns in subarctic environments and their relationships with environmental conditions. This study used samples collected across an elevation gradient (169-1183 m.a.s.l.) of subarctic lakes in Sweden (Larocque et al. 2001) to evaluate the effects of temperature and catchment characteristics on chironomid trait diversity. We tested the hypotheses (1) that climate directly and indirectly shapes trait composition of chironomid communities across the studied 1000 m-elevation gradient, and consequently (2) that generalist taxa with smaller body size and large ranges in food preferences are more able to persist in cold environments than specialist taxa.

Study sites and chironomid training set
We used the previously published training set (Larocque et al. 2001) describing chironomid assemblage composition of 100 lakes in the northern boreal and Arctic/alpine ecoregions of Sweden (see Table S1 for limnological and geographic characteristics). Due to influences of elevation (ranging 169-1183 m a.s.l.), the area is characterised by a strong gradient in temperature (estimated mean July air temperatures ranging 7.0-14.7 °C, Larocque et al. 2001), and the ice-free season ranges from ca. 6 months at low elevation lakes to < 2 months at the highest elevations (Table S1). Moreover, catchment vegetation ranges from coniferdominated (Pinus sylvestris, Picea abies) boreal forests at low elevations, through mountain-birch woodland (Betula pubescens ssp. tortuosa) and alpine meadows, to areas with sparse alpine-tundra vegetation (lichens and mosses) or barren ground at the highest elevations. The lakes are generally small (< 20 ha), shallow (average maximum depth of 6 m) and oligotrophic (total phosphorus concentrations ranging 8.3 ± 7 µg L −1 ).
Surface sediment samples (0-1 cm) were retrieved from the deepest part of each lake using a modified Kajak-gravity corer. In the lab, chironomid head capsules were handsorted from samples of wet sediment (3-5 mL) following recommendations from Walker (2001) and identified under a microscope using Wiederholm (1983). Chironomid counts were then pooled to the genus level to fit the taxonomic resolution of the European Chironomidae Database (Serra et al. 2016). Data are expressed as relative abundance, and only taxa occurring in at least two lakes were used.

Functional trait database
The functional diversity of chironomid communities was characterized using four life-history traits (i.e. 4th instar larval size, emergence season, voltinism and flight period), as well as traits describing food preference, and feeding habits obtained from the European Chironomidae Database published by Serra et al. (2016; 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.  . 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 traitby-lake 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 lake.

Data analysis
An overview of the data analysis workflow adopted in this study is provided in Fig. S1. A principal component analysis (PCA) was performed on lake characteristics to provide an overview of environmental gradients in the training set (step 1 in Fig. S1). Functional diversity of chironomid communities was analysed following a procedure similar to that by Desrosiers et al. (2019). In this procedure, correlation between functional traits (Table 1) and chironomid taxa found in the training set was analysed using a fuzzy correspondence analysis (FCA) applied to the taxa-by-trait category dataset (18 taxa × 33 categories of 6 traits; (step 2.1 in Fig. S1). The contribution of each trait to total variation was quantified by correlation ratios representing the amount of variance explained by each trait along the first two FCA axes (Picazo et al. 2012). Cluster analysis, based on a matrix of Euclidean distances among the taxa calculated on taxon coordinates along the first four FCA axes, was applied to define functional groups of taxa with similar trait combinations (step 2.2 in Fig. S1). Chironomid genera were then grouped into successively larger groups using the Ward's method (Ward 1963) and the significance of the functional groups was assessed using the broken-stick model (Bennett 1996). An FCA was also performed to compare lakes along the elevation gradient based on the trait composition of their chironomid communities (i.e. trait-by-lake dataset: 100 lakes × 33 categories of 6 traits; (step 3.1 in Fig. S1).
Cluster analysis, based on a matrix of Euclidean distances calculated on lake coordinates along the first four FCA axes, was applied to identify lake clusters with similar trait composition, and lakes were then grouped into successively larger groups using the Ward's method and the broken-stick model for selecting significant clusters (step 3.2 in Fig. S1). In addition, the "IndVal" method (Dufrêne and Legendre 1997) was applied to identify indicator categories of traits associated with each lake cluster. Only significant categories (p < 0.05) were selected as indicator trait category in the given group of lakes.
Statistical comparisons of functional group relative abundances and site elevation among lake clusters were also made by Kruskal-Wallis tests followed by multiple comparison tests (step 4.1 in Fig. S1). Finally, a redundancy analysis (RDA) using Monte Carlo unrestricted 999 permutation tests (Legendre and Legendre 2012) was performed to establish a model associating chironomid trait compositions (100 lakes × 33 categories of 6 traits) to environmental characteristics of lakes (Table S1; step 4.2 in Fig. S1). All statistical tests and graphical displays were done using R 3.5.2 statistical software (R Core Team 2018).

Distribution of subfossil chironomid along the elevation gradient
The first two PCA-axes of lake characteristics accounted for 34.5% and 21.2% of the total variance, respectively (Fig. 1). PCA revealed the existence of two distinct gradients: the first axis explained elevation, total organic carbon concentration in water (TOC), organic matter concentration in sediments (as loss on ignition; LOI) and climate parameters (mean air temperature of January and July: TJan and TJul, respectively), whereas the second axis was strongly correlated with specific conductance, pH, and concentrations of major  Fig. 1). The first gradient was interpreted as an elevation/humic pattern, where high elevation and clear water induced negative PCA1 scores and positive PCA2 scores, along which most lakes were spread (Fig. 1). The second gradient consisted of a conductivity/alkalinity pattern, where low-alkaline water and low conductivity induced negative PCA1 and PCA2 scores.
As expected, chironomid communities showed large changes across the elevation gradient (Fig. 2a), similar to that in the original taxonomic resolution applied by Larocque et al. (2001). Several genera including Pseudodiamesa, Abiskomyia and Corynoneura were found only in high-elevation lakes while others such as Heterotanytarsus, Dicrotendipes, Microtendipes, Cladotanytarsus and Polypedilum were associated with warmer lakes at low-elevation (Fig. 2a). A third group of genera including Tanytarsus, Paratanytarsus, Micropsectra, Corynocera, Sergentia and Procladius were instead evenly distributed across the elevation gradient (Fig. 2a). In total, 18 chironomid genera had sufficient information regarding their ecological requirements and biological traits in the European Chironomidae Database (plotted in black; Fig. 2a), and were kept for the second part of the study. For the remaining 11 taxa (plotted in grey), the trait database contained very little information, and these taxa were therefore excluded from further analyses. The exclusion of these taxa should, however, not have biased subsequent analyses as they consisted of rare taxa (accounting for less than ca. 20% of the total chironomid counts on average, Fig. 2b), and were evenly distributed along the elevation gradient (Fig. 2b).

Relationships among traits and trait-based functional groups
The first two axes of the FCA (F1-F2) applied to chironomid taxa and their functional traits (taxa-by-trait category dataset: 18 taxa × 33 categories of 6 traits) explained 15.4% and 12.8%, respectively, and accounted for 28.2% of the total variance (Fig. 3a). The traits "4 th instar larval size" (SIZE1 and SIZE2 vs. SIZE4 and SIZE5, see Table 1 for full names of trait categories), and to a lesser extent "Food preference" (primary producers vs. other food resources) and "Feeding habits" showed the highest correlation ratios along the F1-axis (0.502, 0.152 and 0.134, respectively). These results revealed that some categories of these traits were well separated along this axis (Fig. 3a). Along the F2-axis, some categories of the traits "Feeding habits" (0.372, predators vs. other trophic guilds) and "Food preference" showed the highest correlation ratios (0.281; animals vs. other food resources; Fig. 3a), i.e. higher for the F2-axis than those for the F1-axis (Fig. 3a).

Trait-based site clustering and relationships with environmental gradients
The first four axes of the FCA applied to the chironomid trait compositions of lakes (trait-by-lake dataset: 100 lakes × 33 categories of 6 traits) explained 22.5% of the total variability among lakes (Fig. 4). The Ward's linkage method revealed four significantly distinct groups of lakes (i.e. cluster 1-4 in Fig. 4). Cluster 1 and 4 consisted of 32 and 11 lakes, respectively, and their locations were opposed along the second FCA axis (Fig. 4). Cluster 2 was composed of 54 lakes and located near the centre of F1-F2 biplot. Cluster 3 was  composed of only 3 lakes and was opposed to other clusters along the first FCA axis (Fig. 4). The "IndVal" method highlighted significant indicator trait categories within each lake cluster, categories of traits characterising each cluster were listed in Fig. 4.
Relative abundances of genera from all the functional groups (FG1-4) significantly varied among all four lake clusters ( Fig. 5a and Table 2). Chironomid larvae from FG3 were the most abundant in all lake clusters. Chironomid communities in clusters 1 and 2 were strongly predominated  by taxa of FG3, while other functional groups were poorly represented (FG1, FG2 and FG4). The genera within FG1 were mostly represented in lakes of cluster 4 (Fig. 5a). Chironomid communities in cluster 3 were instead predominated by genera belonging to FG2 and FG3, while FG1 was poorly represented and FG4 not found ( Fig. 5a and Table 2). Elevation also significantly varied among lake clusters, being higher for clusters 1 and 3 than for lake clusters 2 and 4 (Fig. 5b).
Overall, the RDA model linking 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 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 elevation/humic pattern in chironomid trait compositions among lakes. The traits "4 th instar larval size", "Food preference" and "Feeding habits" were separated along RDA axes, and chironomid trait compositions of lakes were well distributed across the elevation gradient. Some indicator traits previously highlighted by the "IndVal" approach, such as large larvae (SIZE5), detritivore feeding habits (DEFEE, SHR) and food preference (SEDMIC and WOOD; Fig. 6), were typical for the chironomid communities of the low-elevation lakes (Clusters 2 and 4; Fig. 6). In contrast, several trait categories, including small larvae (SIZE1), predator and herbivore feeding habits (PRED, FFEEDT and SCR) and associated food preference (MICPHY, MACINVERT and MICINV), were characteristic for lakes high elevations (clusters 1 and 3). The remaining trait categories were judged indifferent, i.e. not clearly associated with any environmental descriptor (Fig. 6).

Discussion
We examined how chironomid biodiversity of subarctic lakes is distributed across an elevation gradient in temperature and landscape characteristics. Based on spatial patterns of functional diversity of chironomids in subarctic lakes, we show that some categories of the traits larval size, food preference, and feeding habits distinctly separate along the studied gradients. Our study therefore suggests complex interplays between direct (e.g. temperature, ice-free period) and indirect climate-related processes (e.g. elevation-driven changes in vegetation/habitats) on the functional diversity of chironomid communities. Resource availability in lakes is strongly linked to vegetation composition/cover, and predominant traits related to resource exploitation in chironomid communities are therefore strongly constrained by landscape characteristics. Finally, our results also show that short life cycles could facilitate the development of viable population in northern and high-elevation lakes where the ice-free period is a limiting factor, thus contradicting patterns showing smaller organisms in warmer environments reported for other invertebrates.

Trait-based functional groups of chironomid taxa
In our training-set, four different groups of chironomid genera were defined by three functional traits, i.e. larval size, feeding habits, and food preferences. These results support previous findings that traits such as food types and larval size are the main contributors to trait variability among aquatic invertebrates, (see also Usseglio-Polatera et al. 2000;Poff et al. 2006;Desrosiers et al. 2019), and also for chironomid larvae (Cummins 1973;Berg 1995;Luoto and Nevalainen 2015). As expected, these trait-based groupings corresponded poorly with known phylogenetic links among main chironomid subfamilies (i.e. Diamesinae, Orthocladiinae, Tanypodinae and Chironominae; Serra et al. 2016). Finally, these traits represent biological features that directly influence ecosystem functions (Hulot et al. 2000), such as energy flow within food webs, making these functional groups excellent indicators of the effects of climate change on the structure and function of subarctic lakes.

Direct and indirect climate processes shape assemblage structure
Overall, our results suggest complex interplays between direct (e.g. temperature) and indirect climate processes (e.g. elevation-driven changes in vegetation/habitats) on the functional diversity of chironomid communities. The observed strong effects of elevation in structuring trait composition of chironomid communities demonstrate that both temperature and large-scale landscape attributes (e.g. land use and vegetation composition/cover) strongly influence the functional diversity of northern aquatic invertebrate assemblages, confirming similar findings for temperate lakes (Poff 1997;Richards et al. 1997;Statzner et al. 2004;Townsend et al. 2003;Motta and Massaferro 2019). However, our results also show a strong influence on trait composition of habitat characteristics, such as water depth, suggesting that some functional features of chironomid communities also showed strong relationships with ecosystem types and local environmental conditions (c.f. Richards et al. 1996, Motta andMassaferro 2019). This is corroborated by Johnson and Goedkoop (2002), who found that local-scale habitat descriptors had a strong unique effect on the functional composition of lake's macroinvertebrate communities.

Effects of elevation gradient on trait composition
Specifically, our study shows that some categories of the traits larval size, food preference, and feeding habits are significantly separated along an elevation/humic pattern, and trait composition of chironomid communities were therefore distributed across the elevation gradient in temperature and landscape characteristics. At low elevation sites, communities were mainly composed of large bodied larvae belonging to the shredder and sediment feeding trophic guilds, whereas larvae in higher elevation lakes tended to be smaller and commonly had predatory feeding mode. Our results thereby support studies that have shown that temperature might play both a direct and/or indirect role in selection for smaller chironomid taxa in high-elevation and cold-water lakes (Hamerlík et al. 2017;Motta and Massaferro 2019). Furthermore, low-elevation lakes with forested catchments typically have more allochthonous organic material that largely covaries with the forest cover ( Fig. 1), and therefore tend to have more sediment-feeding taxa than lakes above the tree line (Luoto and Nevalainen 2015). As expected, resource availability in lakes is strongly linked to catchment characteristics, and traits related to resource exploitation in chironomid communities are therefore well constrained by landscape characteristics. Our results, however, conflict with patterns previously reported for other aquatic invertebrate groups showing that organisms tend to be smaller and better dispersers in warmer environments (c.f. Chown and Gaston 2010;Hayden et al. 2017). Indeed, these studies suggested that increasing body size in colder environments is the result of prolonged lifecycles (slower growth) in response to low temperatures (Angilletta et al. 2004). However, short life cycles could also facilitate the development of viable populations in subarctic and high-elevation lakes where the ice-free period is a dispersallimiting factor (Nolte and Hoffmann 1992). In addition, the ability to feed on available food resources, indirectly linked to catchment characteristics (Heino 2008;Kivilä et al. 2019) might also influence the persistence of taxa. Interestingly, our analysis also show that lake water depth contributes to shape structure and trait composition of chironomid communities, but additional studies are needed to better understand how lake type (including deep vs. shallow lakes) could affect the trait composition of chironomid communities and their responses to climate change.

Implications for climate change effects, and future challenges
Predicted climate change and concurrent climate-induced changes in catchment characteristics will affect Arctic and subarctic lakes and their unique biodiversity (Heino 2008;Lento et al. 2019). Using a space-for-time approach, our study can help to predict future trajectories of freshwater biodiversity in these subarctic lakes. As a consequence of climate warming, migration of chironomid species towards higher elevation sites will be induced by both temperature tolerance of taxa and changes in catchment characteristics (i.e. vegetation composition/cover and associated patterns in organic matter run-off). Hence, the highest elevation lakes in subarctic landscapes might lose their typical cold-adapted chironomid taxa along with their functional attributes, when these are dispersal-limited and replaced by better adapted species from lower elevations or southern latitudes (see Lento et al. 2019).
In addition to space-for-time substitutions, paleolimnological studies are powerful approaches to reconstruct the effects of past environmental changes archived in lake sediments (Smol 1992) and provide excellent opportunities to explore the responses of functional diversity to past climate change. In fact, these paleo-ecological studies cover far more than the variation in climate that we have seen so far since the 1990s, and analysis of subfossil chironomid communities has provided valuable insights on past environmental and climatic impacts on lake's biodiversity (see also Heiri et al. 2003;Luoto and Nevalainen 2018). Most of these studies have, however, applied a taxonomy-based approach, and impacts on the trait composition have been so far largely neglected. Our study, combining trait-based and subfossil communities, takes this one step further by quantifying effects on trait compositions, and contributes to the growing hierarchical understanding of the roles of environmental and climatic changes in structuring functional diversity of chironomid communities.
Trait plasticity is well documented for insects (Cummins 1973;Nylin and Gotthard 1998) leading to differences between actual and theoretical functional traits of species (Bennett et al. 2015) and a potential limitation to the use of chironomid traits as indicators of environmental change. Moreover, most aquatic insects are assumed omnivorous in their early instars (Cummins and Klug 1979), and classification based on feeding strategies of later instars may therefore mismatch with those of all larval stages. Serra et al. (2017) also reported large differences in chironomid traits between North America and Europe and suggested alternative adaptation to similar ecological environments. Functional plasticity implies that many combinations of traits appear for a single species depending on environmental conditions (Jourdan et al. 2019). A better understanding of functional trait plasticity and how it is affected by environmental changes could strengthen the reliability of ecological inferences of traitbased chironomid properties.

Conclusion
Our study shows that complex interplays between direct (e.g. temperature) and indirect climate processes (e.g. elevationdriven changes in vegetation/habitats) affect the functional diversity of chironomid communities. As expected, traits related to resource exploitation (i.e. food preference and feeding habits) in chironomid communities are indirectly constrained by landscape characteristics as these strongly impact on resource availability. Furthermore, our findings suggest that short generation times could facilitate the development of viable populations in subarctic and high elevation lakes where the short ice-free period is a limiting factor, thus contradicting patterns reported for other invertebrates. These results also highlighted how species traits can contribute to a mechanistic understanding of functional changes along climate gradients and provide new knowledge about functional changes in subarctic lakes under climate-warming scenarios. We envision that functional approaches in paleolimnological studies have the inherent strength to forecast future trajectories of freshwater biodiversity in high-latitude lakes by learning from the past changes archived in their sediments.