First circumpolar assessment of Arctic freshwater phytoplankton and zooplankton diversity: Spatial patterns and environmental factors

1. Arctic freshwaters are facing multiple environmental pressures, including rapid cli mate change and increasing land- use activities. Freshwater plankton assemblages are expected to reflect the effects of these stressors through shifts in species distributions and changes to biodiversity. These changes may occur rapidly due to the short generation times and high dispersal capabilities of both phyto- and zooplankton. Spatial patterns and The main objectives of this study were: (1) to assess tial patterns of plankton diversity focusing pelagic (2) to assess dominant component of β diversity (turnover or nestedness); (3) to identify which environmental factors best explain diversity; and (4) to provide recommendations for future monitoring and assessment of freshwater plankton communities across the Arctic region.


| INTRODUC TI ON
The Arctic region is experiencing dramatic climate change, with both air temperature and precipitation expected to continue to increase more rapidly than in any other region of the world (Bintanja, 2018;IPCC, 2014;NOAA, 2018). This will have direct and indirect effects on freshwater ecosystems and their biodiversity. Warmer temperatures are predicted to change species distributions, allowing some species to expand their range while others will lose habitats (Heino et al., 2009;Post et al., 2009;Rolls et al., 2017). Also, human activities are growing in the Arctic regions in the form of infrastructure (roads, airports, settlements), oil and gas drilling, industrial pollution, and mining activities, which will lead to increased pressure on the landscape that may result in cumulative effects such as eutrophication and browning of freshwaters (Hayden et al., 2019;Huntington et al., 2007;Wrona et al., 2013). Latitudinal and longitudinal biodiversity studies across Arctic regions are required to understand a range of current environmental conditions and provide a baseline in order to better predict the magnitude of future changes.
A characteristic feature of the Arctic landscape is the great diversity and abundance of lakes and ponds derived from repeated glacial cycles evident in large parts of the circumpolar region, as well as in regions with discontinuous permafrost. Environmental conditions vary greatly at both larger (ecoregions) and smaller (watercourses, lakes, and ponds) spatial scales and are driven by differences in bedrocktype, geomorphological processes, and temperature regimes, among other factors (ACIA, 2005;AMAP, 1998AMAP, , 2017. Thermokarst or thaw ponds and lakes, which account for approximately 25% of the estimated total area covered by ponds and lakes in the Arctic (Grosse et al., 2013), encompass a variety of transparencies and trophic conditions, whereas most Arctic ponds and lakes are oligotrophic with low nutrient conditions (Rautio et al., 2011;Wauthy et al., 2018).
Despite this great limnological diversity, the majority of Arctic lakes and ponds are characterised by low biodiversity and relatively simple food webs Rautio et al., 2011).
These characteristics make them highly sensitive to environmental pressures that alter biodiversity and ecosystem properties Lizotte, 2008;Vadeboncoeur et al., 2003).
Due to a short growing season and typically nutrient-poor conditions, Arctic lake and pond ecosystems are expected to respond rapidly to increasing temperature and nutrient concentrations (Przytulska et al., 2017;Rouse et al., 1997). Phytoplankton and zooplankton are suitable as indicators of such environmental change because they have high taxonomic diversity, have short generation times, and respond to subtle changes in temperature and mixing regimes (Adrian et al., 2009). While species distributions and composition of phyto-and zooplankton assemblages have been well-studied in some regions of the Arctic, these biotic groups are not typically included in monitoring programmes and therefore, large-scale assessment of their circumpolar diversity is lacking (Wrona et al., 2013).
Phytoplankton species richness and biomass can vary greatly across Arctic lakes and ponds depending on the environmental conditions that regulate their essential needs for resources (primarily light, inorganic carbon and other nutrients) and on loss rates from UV radiation, viruses, and predation by zooplankton (Prowse et al., 2006;Reynolds, 2006). The number of species and overall biomass is expected to decrease with increasing latitude due to lower temperatures caused by reduced solar radiation and leading to a shorter growing season (Gaston, 2000;Stomp et al., 2011). Zooplankton communities of Arctic lakes typically include crustaceans (calanoid and cyclopoid copepods and cladocerans) and rotifers. Despite the typically oligotrophic conditions of Arctic freshwaters, zooplankton whereas isolated regions had lower taxonomic richness. Ecoregions with high α diversity generally also had high β diversity, and turnover was the most important component of β diversity in all ecoregions. 4. For both phytoplankton and zooplankton, climatic variables were the most important environmental factors influencing diversity patterns, consistent with previous studies that examined shorter temperature gradients. However, barriers to dispersal may have also played a role in limiting diversity on islands. A better understanding of how diversity patterns are determined by colonisation history, environmental variables, and biotic interactions requires more monitoring data with locations dispersed evenly across the circumpolar Arctic. Furthermore, the | 3 SCHARTAU eT Al.
abundance and biomass can be relatively high in shallow Arctic lakes because of the absence of fish and the presence of benthic algal mats that are important feeding habitats (e.g. Mariash et al., 2014;Rautio & Vincent, 2006). Many ponds and lakes in the Arctic are devoid of fish because they are shallow (<3 m) and the water column is partially or completely frozen to the bottom during the winter (e.g. . These fishless ponds and lakes are characterised by high abundances of the herbivorous zooplankton Daphnia, an important taxon for the structure and functioning of many freshwater communities and food-webs across the northern hemisphere (e.g. Jeppesen et al., 2017;O'Brien et al., 2004;Stross et al., 1980). Climate change is predicted to exert both direct and indirect effects on Arctic zooplankton assemblages. For example, warmer temperatures, a longer open-water season, and higher primary productivity will favour many cladoceran species at the expense of copepods (Rautio et al., 2008).
Climate and land-use changes in Arctic regions can jointly shift lake biological communities from a periphyton-dominated system towards food-webs driven by pelagic phytoplankton production. For zooplankton, these food-web shifts result in clear changes in diversity and composition. For example, small filter-feeding cladocerans benefit more from increasing phytoplankton biomass compared to larger cladocerans and omnivorous copepods (Hayden et al., , 2019 that rely more on their lipid storages during periods of food shortage (Grosbois et al., 2017).
Latitude, longitude, and elevation are undoubtedly important spatial factors influencing aquatic communities; however, local environmental factors mainly drive phytoplankton diversity (Henriques-Silva et al., 2016;Hessen et al., 2006;Stomp et al., 2011). For example, plankton diversity has a well-established positive relationship with lake size (Dodson, 1992;Dodson et al., 2000;Merrix-Jones et al., 2013). Global environmental change, leading to a combination of increased water temperatures and nutrient inputs to Arctic lakes, is predicted to result in increased primary productivity, an increased number of species, and a shift in community composition of plankton (Heino et al., 2009). Circumpolar plankton data from different Arctic regions would likely provide a key to understanding how climate and productivity drive biodiversity in lakes and ponds.
This study was part of the first circumpolar Arctic freshwater biodiversity assessment, led by the Freshwater Group of the Circumpolar Biodiversity Monitoring Program (CBMP), which is a part of the Arctic Council's Conservation of Arctic Flora and Fauna (CAFF) working group. The aim of this synthesis was to combine and analyse currently available data to establish a baseline for future monitoring and assessment of plankton communities in the Arctic by describing the biogeographical distribution of pelagic communities of phytoplankton and zooplankton.
Specific objectives were first to assess spatial patterns of α diversity, as taxonomic richness, and β diversity components (turnover, or change in taxa, and nestedness, or loss of taxa; Baselga, 2010) across the Arctic, and then to link regional patterns in plankton diversity to climate and other environmental variables.
We hypothesised (H1) that plankton α diversity decreases with increasing latitude, from subarctic towards high Arctic, because the harsh environmental conditions at high latitudes limit survival by some species (Henriques-Silva et al., 2016;Hillebrand, 2004;Righetti et al., 2019). Second, we hypothesised (H2) that turnover contributes more to β diversity than nestedness throughout the Arctic, due to high environmental heterogeneity and dispersal constraints caused by glaciation events in the past (Henriques-Silva et al., 2016;Righetti et al., 2019;Soininen et al., 2018). Further, we hypothesised (H3) that plankton species richness is limited by local environmental conditions such as temperature, productivity, and habitat space, and we predicted that richness would be inversely related to elevation above sea level and positively related to summer temperature (both are proxies for the length of the growing season), precipitation (proxy for nutrient flux from catchment) and lake area (more habitats) as shown across taxonomic groups (Dodson et al., 2000;Gaston, 2000;Righetti et al., 2019

| Selection of samples and description of data
Across data sources, phytoplankton were generally collected by a bottle grab for surface samples, with a volume-specific sampler for depth-specific or depth-integrated samples, or with a plankton net.
Most samples were collected by a volume-specific sampler without any subsequent sieving/filtration. If sieving was applied, a mesh size of 10 or 20 µm was used (see Supplement 1 for details and Table   S1a for the final dataset). Zooplankton samples were collected with a plankton net for depth-integrated samples, or a volume-specific sampler followed by sieving for depth-specific samples. Mesh size ranged from 20 to 156 μm, however most samples were collected by a plankton net with a 50-156 μm mesh size for crustaceans and 45-100 μm mesh size for rotifers (see Supplement 1 for details and Table   S1b,c for the final datasets). Species richness of crustaceans was not correlated with mesh size (r = −0.09), which indicated that variation in mesh size did not affect the number of crustacean species sampled. However, as sampling in most regions did not specifically target rotifers, and as they were not efficiently collected by larger mesh sizes, rotifers were excluded from zooplankton analysis. A comprehensive species list of rotifers found in the datasets is presented for information in Table S9.
Based on collection methods, all phytoplankton were considered pelagic. Zooplankton assemblages are composed of true pelagic taxa mainly found in open waters, but they can also contain benthic species, particularly in small water bodies. Since this study focused on planktonic assemblages, we excluded all samples indicated to represent benthic habitats (e.g. littoral, scrape, close to bottom). For both phytoplankton and zooplankton, depth-specific samples were combined into one sample so that all stations were represented by composite samples. The resulting datasets, which were used to give a general overview of samples and plankton taxa represented in the CBMP-Freshwater database (ABDS), included a total of 1,736 samples from 272 lakes for phytoplankton (PHYTOTAL; Figure 1a; Table   S2) and 3,105 samples from 442 lakes for zooplankton (ZOOTOTAL; Figure 1c; Table S3; see Supplement 1 for details). The lakes covered a wide geospatial gradient in latitude (55.2-83.5°N), longitude (−159.8 to −91.9°E), elevation (0-947 m a.s.l.), and lake area (<0.0001 ha-1,043 km 2 ). Further, the dataset included deep (stratified) lakes as well as shallow lakes and small ponds of different origin (tundra, bedrock, thermokarst). In this paper, we have not made any attempt to differentiate between ponds and lakes as information on lake depth was missing in many cases.
There were some differences among lakes with respect to sampling effort (e.g. single to multiple sampling events) and time of sampling. The collection time period for samples ranged from 1908 to 2015 for phytoplankton, and from 1900 to 2016 for zooplankton (see Tables S4 and S5 for (Tables S6 and S7). Temporal patterns were not assessed because time series data were too sparse. Selection of data time periods and samples is discussed further in the description of data analysis.

| Harmonisation of nomenclature
Nomenclature was harmonised according to recent taxonomic sources to adjust for any regional differences in naming conventions, and we updated nomenclature from older samples that used outdated taxonomy. For phytoplankton, we cross-referenced our taxon list with species on AlgaeBase ® (Guiry & Guiry, 2017). If a taxon in the dataset was not matched with AlgaeBase, the next higher taxonomic classification level was used (see Table S8 for the most common phytoplankton taxa across the Arctic zones). For most analyses, the number of taxa (species-level) was used, for α-and β diversity both species-level and class-level were analysed. Zooplankton nomenclature was corrected and updated based on the latest internationally accepted nomenclature (see Table S9 with references).
For taxa with unclear taxonomy, for instance due to hybridisation, we combined records of species into species groups (e.g. Daphnia longispina complex, which includes Daphnia cristata, D. cucullata, D. galeata, D. hyalina, D. longispina, and D. umbria).
Harmonisation of nomenclature resulted in a phytoplankton dataset of 1,535 taxa across 42 classes, 1,154 of which were identified to the species level. The phytoplankton diversity was characterised by 140 common taxa across 15 classes, which contributed to 85% of the phytoplankton diversity (Table S8). Many taxa were encountered only a few times and in less than 10% of the lakes.
The zooplankton dataset included a total of 363 taxa, 288 of which were identified to species (Table S9). Approximately 40% of all taxa were found in only one or two lakes. The main groups (i.e. Calanoida, Cyclopoida, and Cladocera) were represented by 29, 37, and 65 species, respectively. Rotifera were represented by a further 125 species (Table S9). Harpacticoida and Ostracoda were represented by 22 and 11 species each, whereas other groups were represented by nine taxa (Table S9).

| Harmonisation and aggregation of data prior to analysis
Assessment of spatial patterns in plankton assemblages required further data harmonisation and aggregation to minimise bias in the data resulting from differences in the timing and frequency of sampling among lakes (see Supplement 1 for details). We selected only data collected within the growing season (June-September).
To avoid confounding spatial trends with temporal trends, we omitted older data from our analysis and only selected data from 2000-2016. All data (e.g. recorded as biomass, density, or count) were converted to presence-absence due to limited data on abundance/biomass and to account for differences in sampling methodology. Prior to analysis, samples from single-taxon surveys were removed. Zooplankton samples with a high proportion of benthic species and samples with low taxonomic resolution were also removed.
For all statistical analyses, each lake appears only once in the dataset. For analysis of taxonomic richness at the lake scale, we averaged presence-absence data by lake, across temporal and spatial variation in the data (see Supplement 1 for details). The resulting datasets included 171 lakes for phytoplankton (PHYDIV; Figure 1b, Table S2) and 206 lakes for crustacean zooplankton (CRUSTDIV; Figure 1d, Table S3). For further analysis of α and β diversity, we selected one sample from each lake. This sample was selected from the last year with data (see Supplement 1 for details), resulting in a dataset of 171 lakes with 1,533 phytoplankton species belonging to 42 classes (PHYDIV) and 186 lakes with 45 crustacean taxa (species and genera; CRUSTDIV2; see Table S9).
CRUSTDIV2 is similar to CRUSTDIV except for 20 Finnish lakes, which were excluded prior to analysis (see Supplement 1). For both phytoplankton and zooplankton, samples were chosen to maximise the quantity of data for spatial analysis. We acknowledge that diversity may be underestimated in some regions if chosen samples did not reflect the full suite of taxa found across sampling years, F I G U R E 1 Maps of lakes with plankton data used in the present study. (a) Phytoplankton included in the CBMP database (PHYTOTAL); (b) phytoplankton used in analysis of phytoplankton diversity (PHYDIV); (c) zooplankton included in the CBMP database (ZOOTOTAL); (d) crustacean zooplankton used in the analysis of diversity (CRUSTDIV). Lakes with data on rotifers (ROTDIV) are also indicated. Additional maps presenting lakes with environmental data used in the analysis of Arctic plankton are shown in Figure S1. Background map downloaded from: Esri, DeLorme, GEBCO, NOAA, NGDC, and other contributors (https://www.arcgis.com/home/item.html?id=305cb 48615 d743f 89fe7 0ff3a 819e59e) but the selection of a single sample was intended to avoid introducing temporal variation into the data.

| Geographical regions
To standardise the spatial scale of the circumpolar analysis of diversity, stations were grouped based on regional conditions in biogeography, climate, and other environmental factors that could potentially affect habitat conditions within freshwater ecosystems (Olson et al., 2001). At the largest scale, stations were classified by Arctic zone (subarctic, low Arctic, and high Arctic), following the Arctic Biodiversity Assessment (CAFF, 2013), resulting in a broad south-north gradient. Data from other southern climatic zones (restricted to the Faroe Islands and Swedish alpine sites) were combined with the data from the subarctic region. This broad grouping of lakes by Arctic zone allowed us to assess broad-scale relationships between species richness and water chemistry, despite the lack of data on water chemistry from many lakes.
Lakes were further classified by ecoregion for assessment, using climate-based terrestrial ecoregions (Terrestrial Ecoregions of the World; Olson et al., 2001) rather than the larger flow-related Freshwater Ecoregions of the World (Abell et al., 2008), which do not reflect regional climate differences as strongly. To provide a standardised and automated approach to grouping lakes at smaller scales and extracting geospatial supporting variables, lakes were grouped into standardised catchments at two spatial scales by using hydrobasins from the global HydroSHEDS project (Lehner & Grill, 2013). Hydrobasins are standardly derived catchments that range in size from continental-scale basins (level 01 hydrobasins) through a hierarchy of sub-basins extracted at successively smaller scales to the smallest-sized sub-basins (level 12 hydrobasins). Midscale level 05 hydrobasins were chosen to group lakes for estimation of β diversity within each ecoregion. The smaller level 07 hydrobasins were used to extract the geospatial variables temperature and precipitation.

| Environmental data
In our analysis of spatial patterns in plankton diversity, we included several key geospatial and environmental variables: latitude, longitude, site elevation (m), lake area (km 2 ), long-term average  mean annual precipitation (mm), long-term average  mean annual air temperature (°C), long-term average  mean summer (July-August) air temperature (°C) as proxy for water temperature (which was only available for few lakes), conductivity (Cond, µS/cm), and total phosphorus (TP, mg/L). For phytoplankton we also included chlorophyll-a (Chl a, µg/L), total nitrogen (mg/L), dissolved oxygen (mg/L), pH, Secchi depth (m), calcium (Ca, mg/L), chloride (Cl, mg/L), and potassium (K, mg/L), whereas total organic carbon (TOC, mg/L) was included for zooplankton only. Some values below detection/reporting limits were reported for TP. Detection limits were often not reported and therefore we ran the analyses with the values reported (i.e. no adjustment of values).
Climate geospatial layers from WorldClim version 2 (http://www. world clim.org/) were used to calculate annual mean and monthly mean long-term average  precipitation and air temperature values for each level 07 hydrobasin (standardised catchment). Long-term average mean summer air temperatures were calculated from the average of monthly mean values for July and August. Long-term average mean summer and winter air temperatures were strongly correlated, and therefore we did not include the latter metric in our analysis. For conductivity, we used values from laboratory measurements and included only field measurements if no laboratory measurements were available. Other parameters (e.g. presence and composition of fish community, lake depth, and water chemistry) were available but only for a very limited number of sites in the CBMP database.
To understand which environmental factors potentially can explain the observed differences in taxonomic richness between Arctic zones, we compared differences in medians for all abiotic variables between Arctic zones using Kruskal-Wallis test and pairwise Wilcoxon rank sum tests with Holm-adjusted p-values for multiple comparisons. Statistical analyses were performed in R 3.4.3 (R Development Core Team, 2017).

| Analysis of spatial patterns in diversity
Our analysis focused on three aspects of freshwater biodiversity: taxonomic richness (as observed per lake); α diversity (here defined as the rarefied number of taxa found at an ecoregion scale), and β diversity (here defined as compositional dissimilarity between lakes within the ecoregion).
We assessed observed taxonomic richness (per lake) across Arctic zones (sub-, low, and high Arctic). Comparisons between zones were made with Kruskal-Wallis test and pairwise Wilcoxon rank sum tests with Holm-adjusted p-values for multiple comparisons. We used non-parametric tests since taxonomic richness data were not normally distributed and for each of the datasets, variances differed between the three Arctic zones (Bartlett's K 2 (2) = 43.21 and K 2 (2) = 20.34 for phytoplankton and zooplankton, respectively, p < 0.001 for both). Statistical analyses were performed in R 3.4.3.
We compared a standardised measure of α diversity (as rarefied taxon richness) among ecoregions to assess regional and climatic influences on diversity patterns. To account for differences in the numbers of lakes that were sampled in each ecoregion, which ranged from two to 57 and which can affect comparisons of diversity if more species are encountered as additional lakes are sampled (Colwell et al., 2012), we used rarefaction curves to estimate the average number of taxa found within each ecoregion at a chosen number of lakes (e.g. Colwell, 2013). In this procedure, sample-based rarefaction curves were generated for each ecoregion with a minimum of five sampled lakes (to ensure sufficient data to generate accurate estimates), and curves were extrapolated, if necessary, to include estimates for at least 10 lakes (see Colwell et al., 2004 for details on extrapolation). Rarefaction curves were randomised 100 times, and the average taxonomic diversity was estimated with upper and lower 95% confidence intervals. Rarefied or extrapolated α diversity estimates at 10 lakes were then extracted from the curve, with 10 lakes chosen as the extraction point to maximise the accuracy of diversity estimates for highly sampled ecoregions (e.g. those with more than 10 lakes), while avoiding excessive extrapolation of those ecoregions with 5-10 lakes. Estimation of α diversity through rarefaction curves was completed using the program EstimateS (Colwell, 2013;Colwell & Elsensohn, 2014).
We estimated β diversity and the relative contribution of each of its component parts, turnover (dissimilarity by replacement of taxa among samples) and nestedness (dissimilarity by loss of taxa among samples) following Baselga (2010). Beta diversity was estimated for each level 05 hydrobasin within each ecoregion using multiple-site dissimilarities as Sørensen dissimilarity (β SOR ) and partitioned into turnover (β SIM ) and nestedness (β NES ) components (see Baselga, 2010). A perfectly nested community means that sites with lower richness are subsets of sites with higher richness, whereas a community with high turnover can have the same richness, but different species at each site, i.e. different species composition (Almeida-Neto et al., 2012).
Hydrobasin-scale estimates of β diversity, turnover, and nestedness were averaged for each ecoregion to determine the dominant component. If neither of the components (turnover or nestedness) contributed to ≥66.7% (i.e. two-thirds, chosen as a conservative threshold) of the β diversity, they were regarded as contributing equally to β diversity. Analysis of β diversity was conducted using the betapart package version 1.5.1 in R (Baselga & Orme, 2012).

| Analysis of environmental variables
Preliminary exploration showed considerable heterogeneity in the variance of environmental variables, and substantial regional and countrywise differences in the number of taxa. Therefore, in the analysis of relationships between plankton taxonomic richness and environmental variables, we constructed mixed effects models in R 3.4.3, using the lme4 package (Bates et al., 2015). The datasets covered a limited range in temperature and other environmental variables; therefore, we are unlikely to have captured a full gradient in these variables and we chose linear models. Prior to model fitting, the data were checked for heteroscedasticity, normality, and influential outliers by examining plots of residuals (Zuur et al., 2009). For simplicity (to avoid circular statistics), we transformed latitude and longitude to the shortest geographical distance from the equator and 0/180° meridian, calculated using the WGS84 ellipsoid and the Vincenty (ellipsoid) method (Vincenty, 1975) in the R package geosphere (Hijmans, 2019). We analysed biotic-abiotic relationships for two different subsets of PHYDIV and CRUSTDIV, respectively (see Figure S1). The full mixed effects models (PHYENV1: 171 lakes; CRUSTENV1: 206 lakes) were fitted using maximum likelihood and included average number of taxa as a response variable, and distance from the equator (not in phytoplankton sets, see below), distance from 0/180° meridian, elevation, lake area, annual precipitation and average summer air temperature as fixed variables and country as a random variable with random intercepts.
In addition, we constructed models for more limited phytoplankton and zooplankton datasets where variables, indicative of the nutrient condition/productivity, were included as additional fixed variables. Total phosphorus (PHYENV2: 57 lakes) and specific conductivity (CRUSTENV2: 98 lakes) were chosen because they were best represented in the data as indicators of productivity. Our selection of predictors are anchored in existing literature: species diversity generally decreases with increase in latitude (Hillebrand, 2004), and follows decrease in annual precipitation and in temperature (Gaston, 2000;Righetti et al., 2019). Similarly, there is strong evidence that species diversity decreases with elevation (Hessen et al., 2006;Stomp et al., 2011), and some evidence that species diversity may follow longitudinal gradients (Hessen et al., 2006). In addition, planktonic species diversity is expected to correlate with lake area (Dodson, 1992) and productivity (Dodson et al., 2000). All fixed variables were centred by subtracting median from the values. Country was included as a random variable into the mixed effects models to address the (apparent) diversity of sampling methods among countries, but this variable was reduced from the full models because the number of lakes was only two or three in some of the countries and change in Akaike information criterion (∆AIC) of models with and without the random variable was <2.0. Therefore, only the full model for PHYENV1 was a mixed effect model (country as a random variable) and the full models for CRUSTENV1, CRUSTENV2 and PHYENV2 were multiple linear regression models. The full models included two-way interactions of the fixed variables only if they were informative (85% confidence intervals [CIs] of the estimates did not include zero). The full models did not include annual air temperature nor continentality (measured as temperature difference between the coldest and the warmest month) since they were strongly correlated with precipitation (Pearson's r = 0.77 and −0.70, respectively, p < 0.001). In addition, distance from equator was removed from the PHYENV1 and PHYENV2 full models because of moderate correlations with precipitation and summer temperature (Pearson's r = −0.61, p < 0.001), and because it increased substantially the variance inflation factors (>5) if included in the full models. Following the principle of parsimony, we included as few parameters as possible and simplified models according to the procedures described in Richards (2015). We identified top models, using the dredge function in the R package MuMIn (Bartoń, 2020), and ranked the models based on corrected AIC (AICc; models with ∆AICc < 6.0 from the best model were included in top models), and removed nested models. We calculated adjusted r 2 , Akaike weight, variance inflation factors and 85% CIs (Arnold, 2010;Grueber et al., 2011) for each model and each variable in top models.

| General overview of plankton taxa
The most common phytoplankton genera (PHYTOTAL) across all five classes in terms of overall number and spatial distribution were Gymnodinium (dinoflagellate), Katablepharis, Cryptomonas (cryptophyte), Chlamydomonas (chlorophyte), Cyclotella (diatom), and Dinobryon (chrysophyte; Table S8). The most common zooplankton species (ZOOTOTAL), in terms of number of lakes and spatial distribution of observations, were the cyclopoid copepod Cyclops scutifer, the calanoid copepod Leptodiaptomus angustilobus, and the cladoceran D. longispina complex. Among the rotifers, Kellicottia longispina was most common. Less than 50% (77 taxa) of all crustacean taxa (Calanoida, Cyclopoida, Cladocera) were assigned as pelagic or regularly found in either pelagic or benthic habitat; the remaining taxa may be assigned as truly benthic (Table   S9). For rotifers, 97% (131 taxa) of all taxa with known habitat preferences were regularly associated with open waters or mixed habitat (Table S9).

| Diversity patterns in different Arctic zones and ecoregions (H1)
Plankton diversity was assessed as taxonomic richness per lake compared across Arctic zones, as well as α diversity across ecoregions.
There were sparse data from the key continental regions of Arctic Russia and Canada, and the remaining coverage was geographically uneven apart from Norway.
Phytoplankton species-level taxonomic richness did not differ significantly among Arctic zones (Kruskal-Wallis test; χ 2 (2) = 2.86, p = 0.24). The median taxonomic richness was 15.0, 16.5, and 12.0 species in the subarctic, low Arctic, and high Arctic, respectively ( Figure 2a). Chrysophytes, chlorophytes, diatoms, and cyanobacteria were the most taxon-rich groups of phytoplankton across all Arctic zones (together accounting for >50% of the species richness; Figure S2a). Especially high taxonomic richness (>40 taxa/ lake) was restricted to subarctic regions of Canada, Finland, and Sweden ( Figure S3). Of all the chemical parameters, only 10 were represented by enough lakes for comparison between Arctic zones (Table S10). Among these, TP, dissolved oxygen, pH, calcium, chloride, and potassium showed significant differences between Arctic zones with highest median values in high Arctic ( Figure S4a, Table   S11). There was also a significant difference between median annual precipitation among Arctic zones, with subarctic sites being wetter than high and low Arctic (Table S11).
For comparisons among ecoregions, phytoplankton α diversity was assessed for seven ecoregions that had >5 lakes (Figure 3a, Table S14). When data were rarefied to assess taxonomic richness at 10 lakes for each ecoregion, the highest α diversity was found  Table S14).
Taxonomic richness of crustacean zooplankton differed significantly among Arctic zones (Kruskal-Wallis test; χ 2 (2) = 24.40, F I G U R E 2 Taxonomic richness (average per lake) across Arctic zones for (a) phytoplankton species-level (PHYDIV) and (b) crustacean zooplankton (CRUSTDIV). Box-plot displaying median, first and third quantiles (box), minimum and maximum (whiskers), and outliers (points). Sample size (number of lakes) is given under the zone name on the x-axis. Different letters given above boxes indicate that the zones differ significantly (p < 0.05) in taxonomic richness p < 0.001). The lowest richness was found in the high Arctic (median 3 taxa per lake) which differed significantly from subarctic (median 5 taxa per lake, Wilcoxon test, p < 0.001) and low Arctic zones (median 4 taxa per lake, Wilcoxon test, p = 0.004; Figure 2b). The taxonomic richness in the low Arctic was also significantly lower than in the subarctic (Wilcoxon test, p = 0.004). While calanoid copepods were the most species rich group in the subarctic and the low Arctic, accounting for about 50% of the total species numbers, all three groups (Calanoida, Cyclopoida, and Cladocera) were equally represented in the high Arctic ( Figure S2b). For the crustacean zooplankton dataset, geospatial and climatic variables differed significantly among Arctic zones, as did phosphorus and conductivity (Table S12). As with the phytoplankton dataset, mean annual air temperature, mean summer air temperature and mean annual precipitation were highest in the subarctic and lowest in the high Arctic, and median phosphorus concentration and conductivity were higher in the high Arctic than in other zones ( Figure S4b, Table S13). Median lake area was highest in the subarctic and lowest in the high Arctic with few lakes >1 km 2 in the latter ( Figure S4b).
For the among-ecoregion comparisons, α diversity of zooplankton was assessed for 10 ecoregions that had moderate to high numbers of sampled sites, and α diversity estimates were F I G U R E 3 Results of the circumpolar assessment of plankton diversity comparing ecoregions with >5 lakes with phytoplankton and crustacean zooplankton data. (a) Phytoplankton α diversity (number of species) with rarefaction or extrapolation to 10 lakes; (b) zooplankton α diversity (number of species) with rarefaction or extrapolation to 10 lakes; (c) dominant component of phytoplankton β diversity (turnover; nestedness, approximately equal contribution); (d) dominant component of zooplankton β diversity. The ecoregions are indicated by the following labels: AC, Arctic coastal tundra; AD, Arctic desert; BB, Brooks-British Range tundra; ES, east Siberian taiga; FI, Faroe Islands boreal grasslands; KH, Kalaallit Nunaat high arctic tundra; KL, Kalaallit Nunaat low arctic tundra; KK, Kamchatka-Kurile meadows and sparse forests; NR, north-west Russian-Novaya Zemlya tundra; SR, Scandinavian and Russian taiga; SM, Scandinavian Montane Birch forest and grasslands; TC, Taimyr-Central Siberian tundra; WI, Wrangel Island arctic desert; YG, Yamal-Gydan tundra. See Tables S14 and S15 for details extracted at 10 lakes for comparison (Figure 3b, Table S15). Among these ecoregions, α diversity of crustacean zooplankton was highest for the Taimyr-central Siberian tundra with an estimate of 20 taxa in 10 lakes followed by north-west Russian-Novaya Zemlya tundra with 16 taxa and Kamchatka-Kurile meadows and sparse forests with 12 taxa (Figure 3b, Table S15). Alpha diversity was lowest in the Faroe Islands boreal grasslands with 4 taxa, and the Arctic desert with 5 taxa, both of which had significantly lower α diversity estimates than the three most diverse ecoregions ( Figure 3b, Table S15).

| Components of β diversity (H2)
Beta diversity of phytoplankton at the class-level ranged between 0.05 to 0.71 within ecoregions (Table S14). Beta diversity was highest for the Arctic coastal tundra, followed by the East Siberian taiga, Taimyr-central Siberian tundra, and the Kalaallit Nunaat low Arctic tundra, all of which had β diversity measures of 0.50 or more. This indicates that the lakes in these regions showed the highest amonglake diversity, i.e. showed a high differentiation in phytoplankton assemblages based on class-level data. Beta diversity was lowest in the Kalaallit Nunaat high Arctic tundra (0.05) and the Scandinavian and Russian taiga (0.10), indicating that similar classes of phytoplankton were found among lakes (Table S14). Turnover was the dominant component of β diversity in only two of seven ecoregions ( Figure S5, Table   S14), accounting for at least 67.7% of the total β diversity. In contrast, species-level β diversity indicated greater differences among sites, with β diversity in nine ecoregions ranging from 0.29 to 0.88 (Table   S14). The highest β diversity was found in the Arctic coastal tundra and the East Siberian taiga (both 0.88), while the lowest β diversity was in Kalaallit Nunaat high Arctic tundra (0.29) and the Scandinavian and Russian taiga (0.49). Turnover was the dominant component of β diversity in all but the Kalaallit Nunaat high Arctic tundra (for which only 57% of β diversity was attributed to turnover; Figure 3c).
Beta diversity of crustacean zooplankton, calculated for 10 ecoregions, ranged from 0.36 to 0.89 (Table S15), indicating intermediate to low similarity among assemblages. Beta diversity exceeded 0.80 in two of the ecoregions, the Arctic coastal tundra and the north-west Russian-Novaya Zemlya tundra, indicating that these ecoregions had the highest differences in crustacean zooplankton assemblage composition among lakes. Except for two ecoregions, turnover was the most important component of β diversity in all ecoregions, accounting for 51%-97% of β diversity (Figure 3d, Table   S15). For Faroe Islands boreal grasslands and Yamal-Gydan tundra, turnover and nestedness were considered close to equal.

| Taxonomic richness and environmental variables (H3)
The three top models (∆AICc < 6) of phytoplankton taxonomic richness (PHYENV1) included elevation, distance from 0/180° meridian (a proxy for longitude), summer temperature, and lake area as fixed variables ( Table 1). The best model included elevation, lake area, and summer temperature as fixed variables (marginal r 2 = 0.10, conditional r 2 = 0.85), while the simplest of the models included only elevation and summer temperature (marginal r 2 = 0.11, conditional r 2 = 0.83). Elevation was the only predictor with strong support as it was included as a fixed variable in all three selected models. The effects of the fixed variables on the average number of taxa were weak (Table 1, Figure 4a). For example, an increase of 100 m in elevation decreased number of taxa by 1.9 (85% CI 1.3-2.5) using the best model. However, the importance of variation in climatic conditions over elevation was supported by the analysis of a smaller dataset, which included TP in addition to geospatial and climate data (PHYENV2). The average number of phytoplankton taxa was best explained using a model that included interaction between distance from 0/180° meridian and precipitation and between distance from 0/180° meridian and summer temperature (adjusted r 2 = 0.50, Table 1).
The three top models of crustacean zooplankton taxonomic richness (CRUSTENV1) included distance from equator (a proxy for latitude), summer air temperature, distance from 0/180° meridian (a proxy for longitude), and interaction of distance from equator and summer temperature as fixed variables ( Table 2). The best model included all those variables (adjusted r 2 = 0.21) while the simplest of the models included only distance from 0/180° meridian and summer temperature (adjusted r 2 = 0.18). Summer temperature was the only predictor with strong support as it was included as a fixed variable in all three selected models, but the effect was weak as an increase of 10°C in mean summer temperature increased richness by 2.5 (85% CI 1.8-3.1) using the simplest model (Table 2). In the best model, the effect of summer temperature depended on distance from the equator: the effect of summer temperature was stronger at lower latitudes and ceased at high latitudes ( Figure 4b). The support of summer temperature was weaker when considering the smaller dataset, which included conductivity in addition to geospatial and climate data (CRUSTENV2). The top models of this dataset included interaction of distance from 0/180° meridian and precipitation, summer temperature and conductivity as fixed variables. The best model included all these variables (adjusted r 2 = 0.53), but the interaction term (EW:Pr) was the only predictor with strong support as it was included in all top models (Table 2).

| D ISCUSS I ON
Taxonomic richness of pelagic zooplankton was highest in the subarctic and lowest in the high Arctic. Taxonomic diversity varied as expected across the Arctic, but the diversity relationship with latitude was weaker than expected, especially for phytoplankton. The β diversity for both phyto-and zooplankton indi-

8.40
Notes: The fixed variables are distance from 0/180° meridian (EW, unit expressed in 10 2 km), lake area (Area, 10 2 km 2 ), mean summer air temperature (Temp, °C), mean annual precipitation (Prec, mm), and elevation (Elev, 10 m above sea level); 85% confidence intervals (CI) are given for parameters in all models. All fixed variables were centered by subtracting the median, and the central value is given for each dataset (CeV). Country is used as a random variable in PHYENV1 dataset, and SD attributed to Country are given for mixed effects models, as well as standard deviation of the residual error (resSD) for all models. Marginal (mr 2 ) and conditional r 2 (cr 2 ) given for mixed effects models and adjusted r 2 (r 2 ) for multiple linear models (all r 2 values p < 0.001), and degrees of freedom (df) for all models. Akaike weights (w i ) were recalculated after removal of nested models.
a Central value for total phosphorus in PHYENV2 was 0.011. associated with elevation, lake area and summer temperature, whereas zooplankton taxonomic richness models included a proxy for longitude and the interaction between summer temperature and a proxy for latitude.

| Alpha diversity decreases towards the high Arctic (H1)
We found some evidence that plankton taxonomic diversity tended to decrease from the subarctic towards the high Arctic, which is consistent with the hypothesis that plankton richness is limited by the harsh conditions at high latitudes of the Arctic that include lower temperatures, shorter growing seasons, and longer periods of ice cover (Gaston, 2000). Broadly, phyto-and zooplankton diversity followed our hypothesised trend of less diversity at higher latitudes, which is also consistent with previous aquatic research (e.g. Henriques-Silva et al., 2016;Hillebrand, 2004;Righetti et al., 2019). Despite a shorter temperature gradient, a similar decrease in zooplankton diversity towards high latitudes has been shown for both Arctic lakes (Patalas, 1990) and ponds (e.g. Jeppesen et al., 2017;Rautio & Vincent, 2006).
Isolated lakes, those on islands, were also found to be less diverse than those on the mainland, which is in line with earlier observations (Rautio et al., 2008). A specific distribution pattern was less marked for phytoplankton than zooplankton, indicating that crustacean zooplankton may be more dispersal limited than phytoplankton. This is consistent with the prediction that dispersal rate is negatively correlated to body size among organisms that disperse passively, because abundance is the fundamental driver of dispersal and organism size and abundance are inversely related (e.g. Hillebrand, 2004;Soininen et al., 2007). These differences have implications for potential changes in species composition in response to climate change, as barriers to dispersal may slow the rate of change in species distributions (Strecker et al., 2008;Thomas et al., 2001).

| Beta diversity decreases with increasing latitude, with turnover as the most important component (H2)
Generally, β diversity was similar for zooplankton and phytoplankton (species-level data), with high β diversity for both plankton groups in coastal regions of the Canadian Arctic and in inland regions of Russia, indicating large assemblage differences among lakes in these ecoregions. Our results were partly consistent with other studies, indicating an inverse relationship between β diversity and latitude (Soininen et al., 2018).
For zooplankton and phytoplankton, species replacement (i.e. turnover) across lakes was the predominant component of  (Tables 1 and 2). Distance from 0/180° meridian is a proxy for longitude and distance from equator (dL) is a proxy for latitude. The effect of summer temperature on crustacean zooplankton depends on latitude: summer temperature has a stronger effect the lower the latitude is, and there is no effect in high latitude. The inner x-axis ticks represent the observations (lakes). Test statistics of regressions are presented in Tables 1 and 2 Notes: The fixed variables are distance from 0/180° meridian (EW, unit expressed in 10 2 km), distance from equator (dsL, 10 2 km), lake area (Area, 10 2 km 2 ), mean summer air temperature (Temp, °C), mean annual precipitation (Prec, mm), elevation (Elev, 10 m above sea level), and conductivity (Cond, 10 2 µS); 85% confidence intervals (CI) are given for parameters in all models. All fixed variables were cantered by subtracting the median, and the central value is given for each dataset (CeV). Standard deviation of the residual error (resSD) is given for all models. Adjusted r 2 (r 2 , all r 2 values p < 0.001) and degrees of freedom (df) are given for all models. Akaike weights (w i ) were recalculated after removal of nested models.
assemblage differences in all ecoregions, consistent with the results of global meta-analysis (Soininen et al., 2018 Based on calculated β diversity for Arctic lakes, a large part of global variability is indeed found in Arctic lakes (Soininen et al., 2018). This may indicate that both historical events, such as repeated glaciations (Hewitt, 2000), as well as continuing harsh environmental conditions contribute to heterogeneity in β diversity (Soininen et al., 2018;Whittaker et al., 2001). A number of lake features could also contribute to these patterns, including local selective pressures such as the presence, composition, and abundances of fish populations, which may have a strong effect on zooplankton assemblage composition (Brooks & Dodson, 1965;Hessen et al., 2006). For example, large-bodied Daphnia species are often the dominant species in many fishless ponds and lakes in the Arctic (Edmondson, 1955;Jeppesen et al., 2017;van Geest et al., 2007), whereas daphnids are usually absent in lakes with fish (O'Brien et al., 2004). Our results also showed this pattern, as D. pulex was found in nearly 50% of our studied lakes in Svalbard, Greenland, and Alaska, whereas they were not found in any of the larger lakes with fish in Fennoscandia. However, a better understanding of how diversity patterns are determined by colonisation history, environmental variables, and biotic interactions, requires more monitoring data, especially on fish assemblages and other local stressors (Lau et al., 2020).

| Climatic variables are the most important factors affecting taxonomic richness (H3)
Taxonomic richness of both phytoplankton and crustacean zooplankton was positively related to mean summer air temperature, closely following the results of other analyses (Gaston, 2000;Hessen et al., 2007;Righetti et al., 2019). This supports other findings that taxonomic richness is positively related to the length of the growing season , which is directly linked to the mean summer air temperature. The negative effect of elevation on phytoplankton richness also points in the same direction. However, the positive effect of summer temperature on zooplankton richness decreased with increasing latitude. This may be attributed to dispersal limitation exceeding that of temperature in the highest latitude sites, such as lakes at Svalbard and Canadian Arctic archipelago, which are islands, preventing higher richness. Furthermore, the effects of mean summer air temperature and other environmental variables were weak, which may have been due to interacting effects of parameters, since temperature will also affect ice cover, incoming dissolved organic carbon, and nutrients. Temperature when paired with higher dissolved organic carbon was found to increase phytoplankton abundance and uptake of TP (Weidman et al., 2014). Zooplankton at higher latitudes may be more limited by food availability rather than temperature (de Senerpont Domis et al., 2013). Evidence of an east-west gradient (measured as distance from 0/180° meridian) on taxonomic richness was weak, which warrants further study of data that are more evenly distributed circum-polarly. The random effect of the country was evident in the phytoplankton dataset indicating substantial country-dependent differences in the sampling and/or identification methods.
We did not find any strong evidence for a positive relationship between taxonomic richness and lake area for either phytoplankton or zooplankton. This result was in contrast to numerous other studies (Dodson, 1992;Dodson et al., 2000;O'Brien et al., 2004;Stomp et al., 2011), but consistent with the findings of Strecker et al. (2008) who found a negative relationship between richness and lake size in high Arctic lakes on Ellesmere Island in Canada.
Other studies indicated that lake size is a poor predictor of taxonomic richness, and that other environmental variables such as lake productivity, plankton community composition, fish community structure, and differences in immigration and extinction are more important for zooplankton species richness (Hessen et al., 2006;Merrix-Jones et al., 2013). Differences between studies could also be explained by differences in the range of lake area data and species richness, with smaller range and lower species numbers in the Arctic lakes (O'Brien et al., 2004).
Furthermore, large lakes were rare in our datasets and this may have contributed to the weak relationship between richness and lake area in our plankton datasets. Environmental data (especially lake water chemistry) were generally sparse, signifying a clear need for a large spatial study from the Arctic including climatic, lake chemistry, lake morphometry, and biotic variables to test our results.

| CON CLUS I ON S AND RECOMMENDATIONS
Our findings confirmed that there is considerable variation in plankton community composition across the Arctic. Arctic lakes are identified as sentinels of climate change and collecting detailed data across environmental gradients and across regions, from the subarctic to the high Arctic, will improve our understanding of the changes to biodiversity across the northern hemisphere. Currently, gaps in plankton monitoring coverage include continental areas of Arctic Canada and Siberia. However, monitoring in large areas of the Arctic is very challenging, as these locations are only accessible by air, and remote sensing methods (e.g. high-resolution satellite data) should therefore be considered in combination with more conventional monitoring methodology . For more accessible sites, remote sensing data should be accompanied by a set of basic water chemistry measurements.
Future monitoring efforts should also endeavour to improve consistency in sample processing methods and taxonomic resolution using the most updated nomenclature. Further, it is important to include both impacted and pristine sites in different ecoregions and zones for long-term monitoring, in order to differentiate between climate-induced and multi-stressor changes. We acknowledge that the current dataset used to draw these conclusions does not cover all areas of the Arctic, nor does it include all data that have been collected to date. To facilitate more conclusive circumpolar datasets in future biodiversity assessments, CAFF has created a data upload facility with data policy statement https:// abds.is/index.php/contr ibute -data. We hope that this initial effort would promote future monitoring efforts and provision critical data to the next circumpolar biodiversity assessment with higher regional resolution.

ACK N OWLED G M ENTS
This study is a result of the first circumpolar assessment of Arctic fresh- Brianna Levenstein (University of New Brunswick, Canada) for assistance with data formatting, nomenclature and data analysis. We wish to thank Sarah Laske (USGS Alaska Science Center) and three anonymous reviewers for commenting on a previous version of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The metadata and (where allowed by data contributors) raw data used in this paper will be made available on the CAFF ABDS; http:// abds.is. The ABDS is the Arctic GBIF (Global Biodiversity Information Facility).