Multitrophic biodiversity patterns and environmental descriptors of subArctic lakes in northern Europe

1. Arctic and sub-Arctic lakes in northern Europe are increasingly threatened by climate change, which can affect their biodiversity directly by shifting thermal and hydrological regimes, and indirectly by altering landscape processes and catchment vegetation. Most previous studies of northern lake biodiversity responses to environmental changes have focused on only a single organismal group. Investigations at whole-lake scales that integrate different habitats and trophic levels are cur-rently rare, but highly necessary for future lake monitoring and management. 2. We analysed spatial biodiversity patterns of 74 sub-Arctic lakes in Norway, Sweden, Finland, and the Faroe Islands with monitoring data for at least three biological focal ecosystem components (FECs)—benthic diatoms, macrophytes, phytoplankton, littoral benthic macroinvertebrates, zooplankton, and fish—that covered both pelagic and benthic habitats and multiple trophic levels.

3. We calculated the richness relative (i.e. taxon richness of a FEC in the lake divided by the total richness of that FEC in all 74 lakes) and the biodiversity metrics (i.e. taxon richness, inverse Simpson index (diversity), and taxon evenness) of individual FECs using presence-absence and abundance data, respectively. We then investigated whether the FEC richness relative and biodiversity metrics were correlated with lake abiotic and geospatial variables. We hypothesised that (1) individual FECs would be more diverse in a warmer and wetter climate (e.g. at lower latitudes and/or elevations), and in hydrobasins with greater forest cover that could enhance the supply of terrestrial organic matter and nutrients that stimulated lake productivity; and (2) patterns in FEC responses would be coupled among trophic levels.
4. Results from redundancy analyses showed that the richness relative of phytoplankton, macrophytes, and fish decreased, but those of the intermediate trophic levels (i.e. macroinvertebrates and zooplankton) increased with decreasing latitude and/ or elevation. Fish richness relative and diversity increased with increasing temporal variation in climate (temperature and/or precipitation), ambient nutrient concentrations (e.g. total nitrogen) in lakes, and woody vegetation (e.g. taiga forest) cover in hydrobasins, whereas taxon richness of macroinvertebrates and zooplankton decreased with increasing temporal variation in climate.
5. The similar patterns detected for richness relative of fish, macrophytes, and phytoplankton could be caused by similar responses to the environmental descriptors, and/or the beneficial effects of macrophytes as habitat structure. By creating habitat, macrophytes may increase fish diversity and production, which in turn may promote higher densities and probably more diverse assemblages of phytoplankton through trophic cascades. Lakes with greater fish richness relative tended to have greater average richness relative among FECs, suggesting that fish are a potential indicator for overall lake biodiversity.
6. Overall, the biodiversity patterns observed along the environmental gradients were trophic-level specific, indicating that an integrated food-web perspective may lead to a more holistic understanding of ecosystem biodiversity in future monitoring and management of high-latitude lakes. In future, monitoring should also focus on collecting more abundance data for fish and lower trophic levels in both benthic and pelagic habitats. This may require more concentrated sampling effort on fewer lakes at smaller spatial scales, while continuing to sample lakes distributed along environmental gradients.

K E Y W O R D S
climate change, fish, freshwater, macroinvertebrates, macrophytes, monitoring baseline, phytoplankton, zooplankton patterns) of northern lakes O'Reilly et al., 2015).
Large-scale climate-and land-use-induced changes in landscape biogeochemical processes and catchment vegetation (Elmendorf et al., 2012;Wrona et al., 2016) will also affect nutrient and carbon transport (Creed et al., 2018;Hayden et al., 2019;Larsen, Andersen, & Hessen, 2011). For example, many Swedish lakes have experienced dramatic declines in total phosphorus concentrations since the mid-1990s, and these declines could be attributed to the combined effects of greening, climate-driven changes in soil properties and terrestrial organic matter input to lakes, and catchment recovery from acidification (Huser, Futter, Wang, & Fölster, 2018). In contrast, increasing temperature and nutrient concentrations in sub-Arctic Finland are evident in lakes of forested areas (Hayden et al., 2019;Hayden, Myllykangas, Rolls, & Kahilainen, 2017). Changes in chemical and physical habitat characteristics will ultimately affect the biological assemblages of these lakes and the ecosystem services they supply.
Changes in thermal and hydrological regimes are recognised as the major stressors of northern Fennoscandian lakes based on a recent assessment (Lento et al., 2019). These stressors, both individually and collectively, are expected to affect the structure and function of lake biological communities. For instance, shortened ice-cover period can alter lake internal nutrient dynamics (e.g. reduce nitrification in winter; Powers et al., 2017) and phytoplankton biomass and composition Weyhenmeyer, Peter, & Willen, 2013).
Warming will also reduce suitable habitats for cold-adapted taxa particularly those with a narrow thermal tolerance (i.e. cold stenotherms; Wrona et al., 2006). While species with a wider thermal tolerance (i.e. eurytherms) may have a greater capacity to accommodate to climate change, they will still be influenced by warming as their growth and vital physiological processes are temperature dependent (Culp et al., 2012). Range expansion of southern eurythermic species and losses of cold-stenothermic species are also expected (Culp et al., 2012;Vincent et al., 2011;Wrona et al., 2006). The concerted action of these changes will impact on all trophic levels of the lake food webs. The environmental stressors will thus be likely to cause shifts to the productivity, species distribution patterns and food-web structure of Arctic freshwaters Reist et al., 2006;Wrona et al., 2016).
Current approaches for investigating climate-change impacts on northern freshwater communities mainly include field studies along climate (e.g. temperature, precipitation) and environmental gradients (e.g. water colour, ambient nutrient concentrations ;Hayden et al., 2017;Scott, Barton, Evans, & Keating, 2011;Vadeboncoeur et al., 2003;Weyhenmeyer et al., 2013), experimental studies (Hansson et al., 2013;Petchey, Timon McPhearson, Casey, & Morin, 1999), and sediment records (e.g. diatoms and chironomids; Rosén, Cunningham, Vonk, & Karlsson, 2009;Smol et al., 2005). While each of these approaches has its own merits, the sampling of sites along natural gradients can provide specific insights into biological responses at different spatial scales ). This approach is particularly important for assessing biodiversity patterns related to ecosystem productivity, as these are known to depend on spatial scale (Chase & Leibold, 2002). In lakes at lower latitudes and/ or elevations that represent a warmer and wetter climate and allow higher ecosystem productivity, higher taxon richness has been observed in phytoplankton (Weyhenmeyer et al., 2013), macrophytes (Heino & Toivonen, 2008), zooplankton (Hessen, Faafeng, Smith, Bakkestuen, & Walseng, 2006;Shurin et al., 2007), benthic macroinvertebrates (Heino, 2009;Johnson & Goedkoop, 2002), and fish . Studies focusing on a single organismal group may provide important insights into the environmental variables that shape assemblages of this specific group, but they disregard the fact that responses of a single group are not shaped only by environmental descriptors but also by the biotic interactions with taxa at other trophic levels (Seibold, Cadotte, MacIvor, Thorn, & Müller, 2018) and their habitat (Hayden et al., 2017;Johnson & Goedkoop, 2002). Evaluations of Arctic and sub-Arctic freshwater biodiversity responses using multiple organismal groups at different trophic levels and habitats are rare (but see Hayden et al., 2017Hayden et al., , 2019, but essential for our understanding of how biodiversity at the whole-ecosystem level will respond to a rapidly changing environment, and how monitoring programmes for Arctic and sub-Arctic lakes can be further developed and optimised. We studied high-latitude lakes in Norway, Sweden, Finland, and the Faroe Islands with data on at least three biological focal ecosystem components (FECs)-benthic diatoms, macrophytes, phytoplankton, littoral benthic macroinvertebrates, zooplankton, and fish-that covered both pelagic and benthic habitats and multiple trophic levels (Table S1). We analysed biodiversity patterns among these lakes which are distributed across large environmental and spatial gradients, and investigated whether their spatial biodiversity patterns were correlated with a set of abiotic and geospatial variables. This research differs from previous studies in which diversity patterns of individual FECs were investigated separately from different systems (for most recent evaluation see Kahlert et al., 2022;Laske et al., 2022;Lento et al., 2022;Schartau et al., 2022). We hypothesised that biodiversity of individual FECs among these lakes would be higher in warmer and wetter regions (e.g. at lower latitudes and elevations), and in hydrobasins with greater relative forest cover that supplies greater inputs of allochthonous organic matter (and thus provides larger habitat heterogeneity) and of nutrients that stimulate lake productivity (Finstad, Helland, Ugedal, Hesthagen, & Hessen, 2014;Hayden et al., 2019;Tanentzap et al., 2014Tanentzap et al., , 2017. We also hypothesised that biodiversity responses of individual FECs and trophic levels, particularly those in similar habitats, were tightly coupled. In addition, we identified hotspots and coldspots of lake biodiversity, and provided contemporary baselines (e.g. taxon richness) for future comparisons.

| ME THODS
Northern Europe was completely glaciated 10,000-12,000 YBP and the icecap started to melt from east to west, where the last part melted in the Swedish mountains (Svendsen et al., 2004). The eastern part of the sub-Arctic in northern Europe is lowland, where retreating glaciers formed many ice-dammed lakes and temporarily shifted the direction of flow in rivers due to glacial inflow. This caused major alterations to dispersal routes, with many fish species recolonising western regions from eastern refuges (Nesbø, Fossheim, Vøllestad, & Jakobsen, 1999;Østbye, Bernatchez, Naesje, Himberg, & Hindar, 2005). For example, fish species diversity in north-eastern Norway is significantly higher than in central Norway, where immigration occurred mainly via marine routes used by salt-tolerant species (Lehtonen, Rask, Pakkasmaa, & Hesthagen, 2008;Tammi et al., 2003). In the current dataset, the north-eastern watercourses (i.e. Paatsjoki-Pasvik, Alta-Kautokeino) discharging into the Barents Sea were colonised from eastern refuges (Østbye et al., 2005;Siwertsson et al., 2010), whereas the large Tornio-Muoniojoki watercourse draining north to south between Finland and Sweden had strong connections to southern colonisation routes from the Baltic Sea, both after glacial melting and even today (Hayden et al., 2017). Lakes in Swedish mountains drain towards the Baltic Sea and are influenced by both deglaciation and contemporary recolonisation processes. This complexity of different colonisation routes makes longitude potentially as important as latitude in determining species distribution patterns.
Lakes in northern Europe are distributed across various ecoregions and/or landscape gradients (Olson et al., 2001), such as from forested lowlands to mountainous areas characterised by shrub vegetation and bare rock, and from maritime to inland areas. Lakes occur frequently in these landscapes, yet most (c. 90%) are relatively small (<10 km 2 ) and shallow (Lehner & Döll, 2004), and sub-Arctic lakes in northern Europe are particularly sensitive to climate change and associated species dispersal northward (Lento et al., 2019).

| Study design
In this study we quantified biodiversity patterns across multiple organismal groups and trophic levels, among lakes distributed across large environmental and spatial gradients. We selected sub-Arctic lakes in Norway, Sweden, Finland, and the Faroe Islands containing data on at least three FECs (i.e. benthic diatoms, macrophytes, phytoplankton, littoral benthic macroinvertebrates, zooplankton and fish). The dispersal abilities of these FECs differ, as the smallest organisms (phytoplankton) may disperse directly via wind, macrophytes and invertebrates via waterbirds or flight (insects), but fish dispersed mainly via active colonisation after deglaciation. Thus, fish are much more dependent on historical and contemporary watercourse connectivity than other FECs. Based on presenceabsence data, we calculated relative taxon richness for each FEC (hereafter richness relative ) by dividing the taxon richness of a FEC in a lake by the total taxon richness of this FEC in all selected lakes (see below). We also calculated the among-FEC average relative richness (average richness relative ) for each lake. This is the sum of richness relative of individual FECs divided by the number of FEC, and it indicates overall biodiversity in these lakes. This approach avoids potential bias due to differences in the size of the species pool among FECs, i.e. a more diverse FEC group has stronger influences on wholeecosystem biodiversity than the less diverse FEC groups. This also allowed us to identify hotspots and coldspots of overall lake biodiversity, i.e. not for a single FEC, but including all FECs. Then, we investigated whether the observed spatial biodiversity patterns were correlated with a set of abiotic and geospatial variables in order to identify key environmental descriptors of biodiversity.

| Data collection
This study used data collected by the Freshwater Group of the Circumpolar Biodiversity Monitoring Program (CBMP), part of the Conservation of Arctic Flora and Fauna (CAFF) which is the biodiversity working group of the Arctic Council. Data were identified and acquired for the CBMP-Freshwater database (stored at the Arctic Biodiversity Data Service; www.abds.is) which contains data for lakes and rivers generally north of the CAFF-border. These data originate from national monitoring databases, academic research programmes, and published records, although data collection was not exhaustive for all countries due to time and funding restrictions. This study used data of benthic diatoms, phytoplankton, macrophytes, zooplankton, macroinvertebrates, and fish from lakes in Norway, Sweden, Finland, and the Faroe Islands.
To facilitate the analysis of data collected from multiple sources, it was necessary to ensure sampling methods were comparable among datasets. The Scandinavian countries and Finland have long histories of monitoring freshwater biotic and abiotic ecosystem components with standardised sampling and analytical methods (e.g. Appelberg et al., 1995;Friberg et al., 2006). More recently, intercalibration of assessments of lake ecological status has been conducted within the European Water Framework Directive (Kelly et al., 2014;Poikane, Kelly, & Cantonati, 2016).
Sample collection methods are briefly described here, but are provided in detail elsewhere (Kahlert et al., 2022;Laske et al., 2022;Lento et al., 2022;Schartau et al., 2022; Arctic Biodiversity Data Service, www.abds.is). Benthic diatoms were collected by scraping multiple rocks from the littoral zone to create a composite sample.
Phytoplankton were generally collected with a volume water sampler (e.g. Ruttner type). Macrophyte data were collected through either visual surveys (photographs or aquascope and rake) or sampling along transects and within quadrats. Zooplankton were sampled by using volume samplers without filtering or with plankton nets of mesh size ranging 20-90 µm (i.e. the common size range of zooplankton used in monitoring). Sampling methods for littoral benthic macroinvertebrates included kick nets, Ekman grabs, and stone scrubs, with mesh size ranging 250-500 µm. Fish were collected by net sampling (gill nets, seine nets, trawls) for abundance quantification, and occasionally supplemented with minnow traps (Faroe Islands) or lures and long lines (Finland) for detecting presence-absence of fish species. Methods employed were judged most suitable for the sampled habitats.
To avoid differences in taxonomic identity due to region-specific sample processing and/or age of samples, nomenclature for each FEC was harmonised to standardise regional naming conventions and any outdated nomenclature was corrected. Where there was ambiguous taxonomy or potential for errors in identification to species, taxonomic complexes were used to group species (e.g. see Kahlert et al., 2022;Lento et al., 2022;Schartau et al., 2022).
For example, many of the lakes sampled in the northeast of the study area (Paatsjoki watercourse) contain European whitefish (Coregonus lavaretus) morphs (Østbye et al., 2005;Siwertsson et al., 2010;Thomas et al., 2016), but in this paper all morphs were pooled as a single species. Taxonomic identification was dictated by freshwater monitoring programmes in the Nordic countries and was at the species or genus level for benthic diatoms, phytoplankton, macrophytes, zooplankton, and fish. For macroinvertebrates, we generally used genus level for insects and family or order for non-insects. Supporting data included a combination of site-specific environmental variables from the field and geospatial variables (Table S2). Few site-specific supporting variables (e.g. water chemistry, physical habitat descriptors) were measured in all sampled lakes, limiting the number of variables that could be used for analysis. Variables that were available for most or all lakes included lake area, latitude, longitude, elevation, and concentrations of total nitrogen (TN) and total phosphorus (TP) in lake water.
We regarded latitude as an important variable as it summarises climate, vegetation, proximity to the ocean, and other associated environmental characteristics. For TN and TP, we used interannual averages (number of years ranged 1-53) based on annual means from measurements in 1-12 months in each year. Other variables, such as lake depth, conductivity, pH, concentrations of dissolved or total organic carbon and calcium, and spectral absorbance, were available for <70% of the selected lakes, and their availability was not consistent in individual lakes (e.g. a lake had the depth data but not pH). Hence, these variables were excluded from our analyses, though they may explain diversity patterns of specific FECs.
Geospatial variables were calculated for all lakes using ArcGIS (version 10.3, ESRI). Because catchments could not be delineated for all sites in the CBMP-Freshwater database (due to the large number of sites), geospatial variables were summarised across hydrobasins, which are standardly derived flow-based catchments (Lehner & Grill, 2013). Hydrobasins have global coverage and are available as geospatial layers at different scales, ranging from level 01 hydrobasins (continent-scale catchments) to level 12 hydrobasins (smallest scale). For this study, geospatial variables were summarised across level 07 hydrobasins (mid-level spatial scale) for all lakes <150 km 2 in area and across level 05 hydrobasins for lakes >150 km 2 in area (i.e. Inarijärvi in Finland), to ensure the hydrobasins also included the catchment and did not just describe the lake itself.
The geospatial variables summarised for each hydrobasin included area, as well as climate, land cover, bedrock geology, and ecoregion variables (Table S2). Climate variables included longterm average  temperature and precipitation from WorldClim Version 2 (Fick & Hijmans, 2017; http://world clim. org/version2). Long-term average mean values of temperature and precipitation were averaged across all cells in a hydrobasin to calculate monthly means for each hydrobasin. Monthly means and standard deviations were then used to calculate annual mean and annual standard deviation for both temperature and precipitation, as well as an annual coefficient of variation (CV) for precipitation. Land-cover data (percent of each land-cover class for each hydrobasin) were extracted from the CAFF Land Cover Change Index (http://caff.is), and focused on land-cover data from 2010 classified using the coding of the International Geosphere Biosphere Programme with multiple land-cover classes. Bedrock geology data were derived from the Geological Map of the Arctic (Harrison et al., 2011). For each hydrobasin, the percent coverage of each geologic setting type (i.e. lithology) was calculated.

| Lake selection and data description
A total of 74 lakes from the Faroe Islands (five), Finland (33), Norway (19), and Sweden (17) were selected as they had data available for three or more FECs (Table S1). The majority were sub-Arctic lakes, except for Lake Solvatnet (Norway) which was in the Arctic (78.9°N).
Total phosphorus and TN data were absent for six lakes. We generally lacked data for maximum and mean lake depths, precluding calculation of habitat coverages such as proportion of lake surface covered by the littoral zone. Because not all selected lakes contained data for the full set of FECs and environmental and geospatial variables (Table S1; Figure 1), our inferences on the interactions between FECs (e.g. the influences of fish composition and diversity on macroinvertebrate and zooplankton assemblages) are necessarily conservative.
Both presence-absence and relative abundance (biovolume for phytoplankton) data were available for analysis of the FECs for many lakes (details below) with the exception of macrophytes for which only presence-absence data were available. Among the selected lakes, relative abundance data of benthic diatoms and phytoplankton were collected in 6 and 38 lakes, respectively (Table S3).
Similarly, littoral macroinvertebrates, zooplankton, and fish were sampled from 53, 47, and 21 lakes for relative abundance, respectively. Presence-absence data of these three FECs were available in 11, 7, and 38 additional lakes, respectively. However, we included a few presence-absence data for macrophytes and macroinvertebrates that were collected before 1950 (i.e. <4% of the total data of these two FECs; Table S4).
Richness relative calculated using these data did not appear as outliers in the macrophyte and macroinvertebrate richness relative distributions for the whole time period, i.e. from before 1950 to 2015.
Sampling frequencies differed among lakes as the data were obtained from different monitoring programmes (Table S3). For F I G U R E 1 (a) Number of focal ecosystem components (FECs) sampled in each of the 74 selected Fennoscandian lakes and (b) lake frequencies for individual FECs in each region. The border defined by the Conservation of Arctic Flora and Fauna is indicated in (a). Monitoring efforts for the FECs differed both within and among countries. BMI, benthic macroinvertebrates phytoplankton, most lakes (76%) were sampled in multiple years.
However, for other FECs the majority of lakes were sampled in a single year (i.e. 59% [macroinvertebrates] to 71% [fish] in the presence-absence data set, and 52% [fish] to 77% [zooplankton] in the relative abundance data set), although there might be multiple sampling stations and occasions within a single year. The FECs were mostly sampled in summer, especially those that were sampled only once in individual lakes in a year. Because of low sampling frequencies for most FECs, we were unable to rarefy the species richness data. This might have resulted in underestimation of richness relative for individual FECs in lakes with low sampling frequencies compared to those in lakes with higher sampling frequencies. Also, as the time and length of individual seasons vary along the geospatial gradients (e.g. latitude and elevation), and because sampling was conducted only once a year for some FECs and lakes, the effects of season on the FECs were not analysed in our study.
Relative abundance data were used to calculate the biodiversity metrics, i.e. taxon richness, inverse Simpson index (diversity), and taxon evenness (diversity/richness). Calculations of richness relative (based on presence-absence data) and biodiversity metrics were first conducted for each sample (or subsample) of individual FECs in the lakes. Had the samples been collected at multiple stations and time periods (e.g. months and years), the values were averaged by stations, months, and years. This allowed us to obtain interannual averages for each lake, averaging across temporal variation in the data. The use of interannual averages was necessary to minimise bias in the data resulting from differences in the timing and frequency of sampling among lakes and FECs. However, relatively few lakes were sampled in multiple years, and most sampling (i.e. ≥75% occasions for individual FECs) was conducted from 1990-2015 (Tables S3 and S4). Our approach does mean that the spatial data represent different, unknown amounts of temporal variation, and the influence this has had on the spatial patterns reported here cannot be quantified. Also, the biodiversity of larger lakes may have been underestimated in instances where organism distribution was patchy. For example, two stations within a lake might differ in community composition but have similar richness relative and biodiversity metrics for a specific FEC.

| Data analysis
Principal component analysis (PCA) was used to investigate among-lake variability in richness relative (all FECs; based on the presence-absence data) and biodiversity metrics (excluding macrophytes; based on the relative abundance data) as well as correlations between FECs. For most lakes, we lacked information on the full set of FECs. Therefore, we used restricted maximum likelihood to estimate variances and covariances in PCA. Restricted maximum likelihood produces less biased estimates than normal maximum likelihood method when there are missing values (Kenward & Roger, 1997). This approach allowed us to include all 74 lakes in the PCA.
Redundancy analysis (RDA) was used to examine relationships between biodiversity patterns (i.e. richness relative and biodiversity metrics) and a set of abiotic and geospatial variables (Table S2; ter Braak & Prentice, 1988;Jongman, Braak, & Tongeren, 1987). Since complete datasets of both response and explanatory variables are required for RDA, and only a few lakes contained data on benthic diatoms (Figure 1), we excluded benthic diatoms in order to maximise the number of lakes used in RDA.
We conducted two RDAs using the richness relative dataset: (1)  For the biodiversity metrics, the exclusion of benthic diatoms and macrophyte data still resulted in a small lake sample size for RDAs (i.e. eight lakes; Table S2), since fewer lakes contained abundance data (e.g. for macrophytes and fish) than those that had presence-absence data (Table S3), and fewer lakes were sampled for multiple FECs. Thus, we limited RDAs to biodiversity metrics of (1) macroinvertebrates and zooplankton (37 lakes), and (2) macroinvertebrates and fish (14 lakes; Table S2), to determine whether patterns and correlations with environmental variables were similar to those observed in the RDA of richness relative data. Lakes for these two RDAs were situated at 62.1-70.5°N and 62.1-69.0°N, respectively.
In all RDAs, the environmental predictor variables were selected using the model-building procedure in the vegan package (Oksanen et al., 2017) and the step function in R (version 3.3.3; R Core Team, 2017). The step function uses Akaike's information criterion as the selection criterion. The procedure began with an unconstrained model (i.e. without any environmental variables) and used stepwise forward selection to choose the best model. At each step when an environmental variable (i.e. constraint) was added, ANOVA-like permutation tests were implemented to assess if the constraints were significant (p < 0.05). We also checked the variance inflation factors (VIFs) of the selected variables. Only significant environmental constraints with low VIFs (i.e. ≤10; Oksanen et al., 2017) were included in the final model, and a permutation test was used to determine if the final model was statistically significant (p < 0.05). We also included country as a dummy variable, but this inflated the variance substantially because the vegetation reflected by land cover and ecoregion data already showed differences between countries (i.e. country and vegetation were redundant). Consequently, country was not added in our RDAs. Based on our lake selection criteria we did not include the single lake on Svalbard (Solvatnet, 78.9°N) in our RDAs. This lake had low biodiversity (Tables 1 and S2), probably because it was shallow (<4 m) and froze completely during winter. Thus, our inferences about latitudinal influences should be applied to sub-Arctic lakes situated at or below 71.0°N (i.e. continental area only). TA B L E 1 Top five biodiversity hotspots and coldspots identified based on relative taxon richness (%Richness) of individual focal ecosystem components (FECs) or all FECs combined, and average %Richness (Avg%Richness) among FECs

TA B L E 1 (Continued)
The RDAs were run without limiting the number of predictor variables that could be included in the final analysis, but we recognised that this could contribute to overfitting the models. We also conducted RDAs using a reduced initial number of environmental variables (i.e. 9; with no or weak correlations) to investigate if this would change the biodiversity patterns and predictors ( Figure S1). Results and major predictors of the subset RDAs were similar to those in RDAs that did not limit the number of environmental variables. Overall, the use of entire set or reduced number of initial environmental variables in RDAs did not lead to different conclusions on the biodiversity patterns and their predictors. We proceeded with the RDAs that did not and the Breusch-Godfrey tests were conducted using the vegan and the lmtest packages (Zeileis & Hothorn, 2002) in R, respectively.

| RE SULTS
The first two principal component (PC) axes captured most of the variance in both richness relative (80%; Figure 2a,b) and biodiversity metrics (64%; Figure 2c,d) of FECs among the lakes. Principal component analysis results of the richness relative dataset showed that lake ordinations mainly followed two nearly orthogonal gradients.
The first represented a gradient in average richness relative and the richness relative of macrophytes, phytoplankton, and fish, and was strongly correlated with the first PC axis that explained >50% of total variance. The second gradient covered the richness relative of organism groups at intermediate trophic levels, i.e. macroinvertebrates and zooplankton, and was more associated with the second PC-axis explaining 29% of total variance. The second PC axis primarily separated some Swedish lakes from the major lake group (Figure 2a,b).
The PCA of biodiversity metrics showed that diversities of all FECs were positively correlated with the first PC axis, suggesting that they together explained a majority of variance among the lakes (Figure 2c,d).
Taxon evenness of the intermediate trophic levels, i.e. macroinvertebrates and zooplankton, was negatively correlated with the second PC axis (Figure 2d). Taxon evenness of most FECs generally decreased with increasing diversity or richness, as they showed almost opposite directions in the multivariate space. However, evenness for benthic diatoms increased with their diversity. A number of lakes had similar diversity and showed ordination scores from −5 to −2 on the first PC-axis and from −2 to 1 on the second axis (Figure 2c). Some Swedish and Finnish lakes were outside these ranges and separated along both axes.
In all RDAs, 52-88% of total variance in the FEC biodiversity data (i.e. total inertia) was constrained by the environmental predictors, with large proportions (i.e. 73%-96%) of this explained by the first two axes. Based on the RDA results of richness relative among five FECs, the Finnish and Norwegian lakes were separated from lakes in Sweden and Faroe Islands on both axes (Figure 3a). Although fewer lakes were included in this RDA, the results showed two orthogonal gradients in richness relative (Figure 3b) similar to those observed in the PCA (Figure 2a,b). We therefore interpreted the apparent country-grouping as a result of intrinsic differences in lake biotic and environmental attributes (e.g. latitude and vegetation) rather than of artefacts due to differences in sampling methods among countries.
Both gradients in the RDA were strongly associated with climate descriptors and hydrobasin vegetation that could be related to climate and/or human influences. For instance, the average richness relative and the richness relative of macrophytes, phytoplankton, and fish generally increased with increasing latitude and annual variability (i.e. coefficient of variation) in precipitation, and with decreasing annual mean air temperature and precipitation. These FEC and average richness relative were also higher in lakes when the percentage cover of woody savannas (i.e. woody vegetation in grasslands) increased or that of grasslands decreased in hydrobasins. Macroinvertebrate and zooplankton richness relative increased with increasing elevation and decreasing cover of evergreen needle-leaf forests in hydrobasins.
A separate RDA that included only consumers (macroinvertebrates, zooplankton, and fish) showed a similar ordination pattern that was explained by a primary gradient in macroinvertebrate and zooplankton richness relative , and a secondary gradient in fish richness relative and the average richness relative (Figure 3c,d). These gradients were strongly associated with the first and second axes respectively, which captured similar proportions of variance (i.e. 50 and 46%). The first gradient (or axis) corresponded to variability among all Fennoscandian lakes, primarily separating Swedish and Faroe lakes that were positively associated with macroinvertebrates and zooplankton from those in Norway and Finland. The second gradient mainly separated some Finnish lakes (i.e. ≥0.5 on the second axis) from the Fennoscandian lake group. Latitude instead of elevation appeared to be the main descriptor for macroinvertebrate and zooplankton richness relative , such that these were higher at lower latitudes. The RDA results also provided further insights into the local environmental descriptors for the FEC richness relative within climate regions. For instance, fish richness relative and the average richness relative increased with increasing lake TN and percentage cover of Scandinavian and Russian taiga in the hydrobasins, and the richness relative of the invertebrate groups were positively correlated with lake area and proportions of extrusive igneous rock in hydrobasins (Figure 3c,d), although all these FEC attributes also generally changed with latitude and/or elevation (Figure 3a,b). The tight association between fish richness relative and the average richness relative observed in both RDAs suggested that fish could be a structuring FEC for overall biodiversity in Fennoscandian lakes.
The RDA of macroinvertebrate and zooplankton biodiversity was predominated by data from Finnish lakes, but lakes from other regions also had ordination patterns similar to the Finnish lakes Percentage land cover of Scandinavian and Russian taiga was a strong predictor for fish diversity and richness (Figure 4d) as for fish richness relative (Figure 3d). However, fish diversity and richness also increased with increasing hydrobasin area, annual variability in precipitation, and intrusive igneous bedrock proportions in hydrobasin, while those of macroinvertebrates decreased with increasing lake area.

| D ISCUSS I ON
We found that the biodiversity of FECs from different trophic levels in northern European sub-Arctic lakes were strongly F I G U R E 3 Redundancy analysis of richness relative (%Richness) among five focal ecosystem components (phytoplankton, macrophytes, zooplankton, benthic macroinvertebrates, and fish) covering three trophic levels in 13 Fennoscandian lakes situated at 62.1-69.3°N (a and b) and among three focal ecosystem components comprised of only consumers in 39 Fennoscandian lakes at 62.1-71.0°N (c and d). Eigenvalues of the constrained inertia were 5.756 and 2.985 (i.e. 96% and 75% of total inertia) respectively. Adjusted R 2 of the models were 0.88 and 0.67 respectively. The upper panels show lake ordinations, while the bottom panels show explanatory environmental variables (red arrows) indicated by permutation tests (p < 0.05). Percentages of variance explained by the canonical axes are given in parentheses. AnnualCV_Precip and AnnualMeanPrecip, annual coefficient of variation in precipitation and annual mean precipitation; AnnualMeanTemp, annual mean temperature; %EvergreenNLF, %Grasslands, %OpenShrublands, %ScandRusTaiga, %Water, and %WoodySavannas, percentage cover of evergreen needle-leaf forests, grasslands, open shrublands, Scandinavian and Russian Taiga, water bodies, and woody savannas, respectively; %IgneousExtrusive, percentage of extrusive igneous rock in geological setting in hydrobasin. See Figure 2 for other abbreviations affected by climatic variables (i.e. temperature and precipitation) indicated by geographical position (e.g. latitude) and hydrobasin vegetation. At lower latitudes representing a warmer and wetter climate, the richness relative of phytoplankton, macrophytes and fish-key players in bottom-up and top-down controls in lake food webs-decreased, but that of zooplankton and macroinvertebrates-intermediate trophic levels in pelagic and benthic food chains-increased. Among the Finnish lakes, fish and the overall average richness relative also increased with TN, which was probably linked to increased forest cover in the hydrobasin. Lake productivity and biodiversity could be enhanced by a greater supply of organic matter and nutrients originating from terrestrial vegetation and forestry practices at least for oligotrophic and mesotrophic lakes (Hayden et al., 2019;Tanentzap et al., 2014Tanentzap et al., , 2017, while larger inputs of terrestrial organic matter can inhibit lake productivity (e.g. Creed et al., 2018;Finstad et al., 2014).
Biodiversity responses of individual FECs possibly also depended on the interactions between FECs, though they were co-determined by trophic level and the environmental variables. Our findings suggest that, under climate change, the cold-adapted species may be replaced by more southern species that better tolerate warmer water, and conclusions of the climate impacts on lake biodiversity may differ depending on the organismal groups investigated.

F I G U R E 4
Redundancy analysis of biodiversity metrics of benthic macroinvertebrates and zooplankton (a and b; 37 lakes situated at 62.1-70.5°N) or fish (c and d; 14 lakes at 62.1-69.0°N). Eigenvalues of the constrained inertia were 4.729 and 5.157 (i.e. 79% and 86% of total inertia) respectively. Adjusted R 2 of the models were 0.52 and 0.77 respectively. The upper panels show lake ordinations, while the bottom panels show explanatory environmental variables (red arrows) indicated by permutation tests (p < 0.05). Percentages of variance explained by the canonical axes are given in parentheses. TN and TP, lake total nitrogen and total phosphorus concentrations; %ClosedShrublands, percentage cover of closed shrublands; %IgneousIntrusive and %Supracrustal, percentages of intrusive igneous and supracrustal rocks in geological setting, respectively; %KolaPeninsulaTundra and %ScandMontBirchForestGrasslands, percentage of Kola Peninsula Tundra and Scandinavian montane birch forest and grasslands in hydrobasin, respectively. See Figures 2 and 3 for other abbreviations

| Responses of individual trophic levels
Based on the biodiversity trends along the decreasing latitudinal and/or elevation gradients of our sub-Arctic lakes, we infer that, as future climate becomes warmer and wetter, the richness relative of fish, macrophytes and phytoplankton will decrease, but that of the intermediate trophic levels will increase, due to concurrent climateinduced changes in the lakes' physicochemical environment, catchment vegetation, species distributions, and biological interactions.
The richness relative of fish, macrophytes and phytoplankton contrasts to the observed biodiversity changes along larger latitudinal gradients (Heino & Toivonen, 2008;Reist et al., 2006;Weyhenmeyer et al., 2013), and could be driven by the low richness relative of these FECs in the lakes in Faroe Islands and Sweden. The climate effects on macroinvertebrates and zooplankton are expected to be stronger in high-latitude lowland lakes where the richness relative was lower than in lakes at lower latitudes and/or in mountain areas. Thermal preferences and cold tolerance vary among invertebrate taxa (Danks, 1992;Wrona et al., 2013). At high latitudes, taxa intolerant to extreme cold conditions are particularly rare and some common taxa there are generally not determined to genus-or species-level (e.g. chironomids), leading to reduced overall macroinvertebrate diversity estimates (Johnson & Goedkoop, 2002;Scott et al., 2011). In a warmer climate, cold-tolerant species are predicted to be gradually outcompeted by eurytherms and their southernmost distribution limit may be driven northwards. Mountain lakes (especially those in geographically separated mountain ranges), however, can have greater habitat heterogeneity that supports more species-rich and diverse macroinvertebrate assemblages than do their lowland counterparts (Finn & Poff, 2005;Heino, 2009). However, due to data deficiency, our analysis did not include local habitat variables (e.g. substratum type, riparian vegetation), which could be the major descriptors of macroinvertebrate assemblages (Johnson & Goedkoop, 2002). Local habitat data should be collected in future monitoring programmes to assess the changes in habitat characteristics on the geographic gradients (e.g. latitude and elevation) and how these changes will affect the macroinvertebrate diversity patterns.
The latitudinal and elevation effects on macroinvertebrates and zooplankton in this study are generally consistent with those previously reported for northern lakes. For instance, taxon richness and diversity of most macroinvertebrate functional feeding groups (i.e. grazers, shredders, collector-gatherers, and predators) increased from Swedish Arctic-alpine to boreonemoral or nemoral ecoregions that correspond to a decreasing latitudinal gradient (Johnson & Goedkoop, 2002;Johnson, Goedkoop, & Sandin, 2004). Latitudinal declines in zooplankton species richness are similarly evident in Norwegian and North American lakes (Hessen et al., 2006;Schartau et al., 2020;Shurin et al., 2007). We also found that intra-annual fluctuations in temperature and precipitation were important descriptors for taxon richness of zooplankton (decreased), macroinvertebrates (decreased), and fish (increased). These fluctuations may support species coexistence by enhancing temporal niche partitioning and preventing competitive exclusion (Shurin et al., 2010).
However, large temporal environmental variation also can reduce diversity as species may have lower fitness and higher stochastic extinction risks (Adler & Drake, 2008 [Linnaeus, 1758]) while those at higher latitudes (67-71°N), but lower elevations, generally had higher fish richness relative including a range of cold-cool and warm-water species (e.g. Salmonidae, Cyprinidae, Percidae, Gasterosteidae). These patterns are evidently related to the postglacial colonisation history of the area (Laske et al., 2022;Tammi et al., 2003), where the northern-most Fennoscandian inland was colonised from eastern refuges by many cold-adapted species, while southern Finnish Lapland was colonised from southern refuges (Nesbø et al., 1999;Østbye et al., 2005;Siwertsson et al., 2010). We cannot, however, exclude that historical fish introductions may have influenced the observed patterns (Hammar, 1989;Lehtonen et al., 2008). Our finding suggests that most cold-water fish in the lakes were stenothermal, i.e. with narrow thermal tolerances, and more prevalent in northern Fennoscandia . In contrast, warm-and cool-water fish are more restricted to the southern sub-Arctic as their distributions also depend on local climate and environmental conditions related to catchment vegetation or other factors (e.g. adjacent brackish waters, direction of river flows, river connectivity; Hayden et al., 2017;Heino et al., 2009;Reist et al., 2006). This probably explains why fish richness relative and diversity were positively related to lake TN and woody savannas or Scandinavian and Russian taiga in the hydrobasin. In our study, the shrub and forest cover contained broadleaf deciduous species (e.g. birch and aspen) that are more nutrient-rich and have faster litter decomposition rates than conifer species. In addition to postglacial colonisation, we postulate that terrestrial organic matter and nutrient inputs from taiga forest had enhanced allochthonous subsidies and in lake primary productivity, and consequently supported more diverse fish communities particularly in mesotrophic lakes (Dodson, Arnott, & Cottingham, 2000;Finstad et al., 2014;Tammi et al., 2003).

| Potential interactions between trophic levels
The richness relative of macrophytes and phytoplankton were positively correlated with fish but not with macroinvertebrates or zooplankton. Macrophytes can provide habitats to support aquatic consumer production and enhance habitat complexity for epiphyte growth and protection of animals from predators (Burks et al., 2006;Carpenter & Lodge, 1986), but they may primarily favour the small macroinvertebrate taxa without increasing their overall richness (McAbendroth, Ramsay, Foggo, Rundle, & Bilton, 2005). The diversity of invertebrate detritivores may be promoted through the macrophyte decomposition process, although many macroinvertebrate taxa lack the ability to digest and assimilate detrital cellulose (Carpenter & Lodge, 1986). Our finding that increasing macrophyte diversity did not necessarily enhance the overall macroinvertebrate diversity corroborated the findings by McAbendroth et al. (2005), probably because macroinvertebrate taxa responded differently to climatic changes along latitudinal and/or elevation gradients, which were more strongly associated with patterns in their biodiversity.
The positive relationships between richness relative of fish, macrophytes, and phytoplankton could be attributed to similar associations between FECs and environmental variables, but also to the habitat structure provided by macrophytes that supports fish diversity and production (through provision of food and predation refuges), which in turn promotes higher densities (and probably more diverse assemblages) of phytoplankton and benthic diatoms through trophic cascades (Burks et al., 2006;Carpenter & Lodge, 1986). Besides, it is possible that the relatively equal proportions of littoral and pelagic habitats in many of the studied mesotrophic lakes had contributed to these patterns (Hayden et al., 2017). A similar correlation between fish and phytoplankton diversity is also evident in oligotrophic boreal lakes (Lau, Vrede, & Goedkoop, 2017).
The positive macrophyte-fish-phytoplankton interactions are likely to be stronger in nutrient-poor and clear-water lakes (like most of our study lakes), which are more abundant and have supported more macrophyte taxa in northern Fennoscandia. These interactions also may partly underlie our observation that fish richness relative increased with latitudes.
An alternative explanation for the positive correlations between richness relative of fish, macrophytes and phytoplankton is their co-dependence on lake productivity (Dodson et al., 2000). Hayden et al. (2017) reported changes in community and size structures of macroinvertebrates, zooplankton, and fish along a temperature and productivity gradient in a long south-flowing sub-Arctic watercourse between Finland and Sweden. Higher fish density and biomass were found in warmer and more nutrient-rich lakes, but macroinvertebrate and zooplankton responses were nonlinear partly due to habitat-specific top-down control by fish. If this productivity-diversity relationship applies to these lakes, then our results suggest that productivity was controlled by climate (reflected by latitude) at regional scales (i.e. among all Fennoscandian lakes) and by ambient nutrients (i.e. lake TN and/or TP) and hydrobasin vegetation (e.g. woody savannas, taiga) at local scales (i.e. among Finnish lakes).
The strong, positive correlations between the richness relative of macrophytes, phytoplankton, and fish with the average richness relative suggest that these FECs had greater influences on the overall biodiversity than did the intermediate trophic levels and benthic diatoms, although data for the latter were limited. A higher average richnessrelative could be attributed to increased habitat diversity and complexity. In our study, lakes with greater fish richness relative also tended to have greater average richness relative (i.e. biodiversity hotspots;  Thomas et al., 2016;Vander Zanden, Shuter, Lester, & Rasmussen, 1999;Vander Zanden & Vadeboncoeur, 2002). Top-down control of lower trophic levels by fish also affects ecological interactions that can alter community structure and pelagic-benthic reliance along environmental gradients in northern lakes (Hayden et al., 2017(Hayden et al., , 2019. However, Arctic catchments and lakes are highly diverse at both small and large spatial scales due to historical and contemporary factors, which complicates predictions of changing food-web structures and energy sources along latitudinal or climatic variables.

| CONCLUSIONS
Our study provided strong evidence that the biodiversity trends of northern European lakes observed along gradients of climate, geographic locations and hydrobasin vegetation were trophic-level specific: patterns for intermediate trophic levels differed from those for fish and the primary producers. Climate, including temperature and precipitation (related to latitude and elevation), was the variable most strongly associated with spatial patterns. Catchment vegetation, which integrates climatic influences as well as human impacts (e.g. land use), was also important, probably through regulating supplies of nutrients and allochthonous organic matter subsidies to lakes, thereby ultimately affecting their productivity. While patterns were observed in the current lake and catchment data, we strongly recommend consistent recording and reporting of supporting variables in future monitoring efforts, such as basic lake morphometry (e.g. mean and maximum depths, relative littoral area) and physicochemical variables (e.g. Secchi-depth, light, pH, conductivity, concentration of dissolved organic carbon), so that variables that are known to affect lake biodiversity elsewhere can be investigated for their role in sub-Arctic lake biodiversity. The correlations between certain trophic levels (e.g. richness relative of phytoplankton, macrophytes and fish, and diversity of most FECs) suggest that biodiversity responses of one trophic level affect those of the others. Hence, a food-web perspective that integrates multiple trophic levels should be an ideal approach for the monitoring and management of Arctic and sub-Arctic lake biodiversity (see also Duffy et al., 2007;Seibold et al., 2018). Thus far, the availability of abundance data was strongly unbalanced between pelagic (i.e. phytoplankton, zooplankton) and benthic food chains (i.e. benthic diatoms, macrophytes, macroinvertebrates), and few data were available for fish, benthic diatoms, and macrophytes. Future biodiversity monitoring of Arctic and sub-Arctic lakes should also focus on fish (probably by non-destructive sampling of environmental DNA) and lower trophic levels in both benthic and pelagic habitats. This requires more concentrated sampling effort on fewer lakes at smaller spatial scales, while continuing to sample lakes distributed along environmental gradients.
Monitoring programmes should also be adapted to better cover the biodiversity of benthic habitats in northern lakes by the inclusion of molecular methods (e.g. genetic barcoding) that potentially increase the taxonomic resolution of key groups with a highly complex taxonomy, such as the chironomids. Such an approach will contribute to better estimates of biodiversity and ultimately lead to better assessment criteria. In the case of destructive sampling of lake food webs, additional analyses beyond diversity and abundance such as stable isotopes and pollutant concentrations would allow more detailed comparisons among circumpolar monitoring programmes (e.g. CAFF and the Arctic Monitoring and Assessment Programme). We advocate the use of multiple trophic levels for holistic assessment of lake biodiversity and community responses to environmental changes in future monitoring programmes. We also recommend that future investigations address whether diversity trends of a FEC are mainly determined by the colonisation history, environmental variables, or biotic interactions with other trophic levels that concurrently change with the environment. In addition, we advocate a more consistent temporal sampling regime, so that temporal variations in biodiversity patterns and environmental variables could be captured to quantify the effects of climate change.

ACK N OWLED G M ENTS
We thank Richard K. Johnson and the anonymous reviewers for comments to improve an earlier version of the manuscript. We also thank Stefan Hellgren for help in preparing the maps. This study was partly supported by grants from the Academy of Finland (projects 1140903 and 1268566) to Kimmo K. Kahilainen, and from the Norwegian Agency of Environment and the Norwegian Institute for Nature Research to Ann Kristin Schartau.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used in this study are stored and available at the Arctic Biodiversity Data Service (www.abds.is).