Temperature and spatial connectivity drive patterns in freshwater macroinvertebrate diversity across the Arctic

Lento, Jennifer; Culp, Joseph M.; Levenstein, Brianna; Aroviita, Jukka; Baturina, Maria A.; Bogan, Daniel; Brittain, John E.; Chin, Krista; Christoffersen, Kirsten S.; Docherty, Catherine; Friberg, Nikolai; Ingimarsson, Finnur; Jacobsen, Dean; Lau, Danny Chun Pong; Loskutova, Olga A.; Milner, Alexander; Mykrä, Heikki; Novichkova, Anna A.; Ólafsson, Jón S.; Schartau, Ann Kristin; Shaftel, Rebecca; Goedkoop, Willem


| INTRODUC TI ON
Rapid environmental change is occurring in the Arctic, with higher temperatures and increased precipitation leading to ecosystem alterations such as shifts in thermal regimes and terrestrial vegetation that have the potential to greatly impact freshwater biodiversity Heino et al., 2009Heino et al., , 2020Wrona et al., 2006). As mean temperature isotherms and degree-day boundaries shift northward in the Arctic, predicted responses of freshwater organisms to warming include range expansion of southern eurythermic taxa and loss of cold stenotherm taxa (Heino et al., 2009;Oswood et al., 1992;Vincent et al., 2011;Wrona et al., 2006). Such shifts in composition have previously been observed in alpine streams following warming and loss of glacial runoff (Khamis et al., 2014;Lencioni, 2018;Niedrist & Füreder, 2021). In a global meta-analysis, Parmesan and Yohe (2003) described the range shifts of 279 mainly terrestrial and marine invertebrate species in response to climate change, averaging 6.1 km/decade towards the poles. Similar processes are predicted for freshwater flora and fauna (Heino et al., 2009;Pecl et al., 2017;Shah et al., 2014), with the northward drainage of many Arctic catchments facilitating dispersal of taxa from lower latitudes as conditions at higher latitudes become more hospitable (Vincent et al., 2011).
However, the rate at which warm eurythermal taxa move north and overall diversity increases in the Arctic is dependent in part on dispersal capabilities of different organism groups and existing barriers to dispersal (Culp et al., 2012;Heino et al., 2009;Medeiros et al., 2020). For example, changes to diversity may occur at a slower rate on Arctic islands, where physical distance and a lack of geographical connectivity to the mainland acts as a barrier to dispersal and limits the northward movement of eurythermic species Culp et al., 2012;Medeiros et al., 2020).
Benthic macroinvertebrates are an ideal organism group to assess the long-term ecological effects of climate change on Arctic freshwater diversity because they provide time-integrated responses to environmental conditions, including those induced by landscape-level changes in vegetation and hydrology (Brown et al., 2018;Resh, 2008), and different taxa have varied levels of dispersal capability (Medeiros et al., 2020;Sarremejane et al., 2020).
Moreover, benthic macroinvertebrate diversity, composition, and species distribution patterns have been shown to reflect temperature gradients in the Arctic Lento et al., 2020;Medeiros et al., 2020) and alpine regions (Brighenti et al., 2019;Khamis et al., 2014;Niedrist & Füreder, 2021). Investigations of biodiversity in Arctic freshwater benthic macroinvertebrates have found strong associations between community composition and physical and chemical environmental drivers at different spatial scales, although these studies have largely focused on ecosystemspecific assessments within relatively small regions of the Arctic (e.g., richness have been found in glacier-fed streams along a latitudinal gradient in Europe , probably related to water temperature and channel stability (sensu Brown et al., 2018). Similar latitudinal diversity trends have been observed in eastern Canadian Arctic rivers, where lower temperatures and decreased terrestrial vegetation at northern latitudes appear to limit species distributions . The physiological tolerance hypothesis (Currie et al., 2004) suggests that diversity at high latitudes is limited due to cold temperatures that exceed physiological tolerance levels of most taxa. However, community compositional patterns of midges (Chironomidae) on high Arctic islands reflect both temperature and dispersal limitation caused by spatial discontinuity from the mainland (Medeiros et al., 2020). While these examples clearly document regional and local diversity patterns, they do not provide large-scale baseline information on macroinvertebrate distributions and diversity in lakes and rivers across the circumpolar region, i.e., the current state of Arctic freshwater biodiversity (Heino et al., 2009). This baseline information, or reference condition for the Arctic, is critical for the assessment of temporal, climate-driven community changes, and for predictions of future changes to species distributions and diversity Lento et al., 2019), even more so as their diversity is predicted to change rapidly with continued climate warming (Shah et al., 2014). In this context, it is also important to understand the relative importance of local and regional drivers of diversity in these systems.
Our study is part of a broader assessment of the biota of Arctic freshwaters  and includes contemporary data from more than 370 lake and 1,160 river sites from all countries of the Arctic. This research provides the first circumpolar analysis of spatial patterns in Arctic macroinvertebrate diversity, thus setting baselines of change for future assessments. Specifically, our objectives were to: (1) assess spatial patterns in diversity across the circumpolar region to identify diversity baselines for lakes and rivers; (2) compare and contrast the roles of two large-scale drivers of taxonomic composition, i.e., temperature and spatial connectivity, in influencing diversity across the Arctic; and (3) identify additional environmental correlates of macroinvertebrate community composition that will drive ongoing and future responses to climate change in Arctic lakes and rivers. While our focus is on the largest-scale assessment to date of spatial patterns in present-day freshwater macroinvertebrate diversity across the Arctic, we also address commonalities and differences among regions by evaluating biotic-abiotic relationships relevant to climate change.
Based on the results of smaller-scale studies of benthic macroinvertebrate α diversity in the Arctic, we hypothesise that (H 1 ) circumpolar diversity, measured as taxonomic richness, declines with increasing latitude and decreasing water temperature due to thermal thresholds of some taxonomic groups (e.g., physiological tolerance hypothesis; Currie et al., 2004). We predict that diversity will be lower at high latitudes and in low temperature regions, and these low-diversity systems will be dominated by cold-tolerant taxa.
Second, we hypothesise that (H 2 ) local richness is further limited on Arctic islands where distance from the mainland limits northward dispersal (e.g., Medeiros et al., 2020). We predict that diversity will be lower on Arctic islands than on the mainland at similar latitudes.
Lastly, we hypothesise that (H 3 ) limitations to local richness caused by low water temperatures and low spatial connectivity contribute to the strongest differences in composition (β diversity) between mainland and island regions and between low and high Arctic regions. We predict that island regions will contain a subset of the taxa found in mainland regions (i.e., island communities will be nested within mainland regions), and that community composition will be most strongly related to climate-influenced habitat drivers.

| Data collection
This study is focused on lake and river benthic macroinvertebrate communities across the Arctic, including Alaska, Canada, Greenland, Iceland, Faroe Islands, Norway, Sweden, Finland, and Russia ( Figure 1). This work is part of an international effort by the Freshwater Group of the Circumpolar Biodiversity Monitoring Program (CBMP) of the Conservation of Arctic Flora and Fauna (CAFF), the biodiversity working group of the Arctic Council Lento et al., 2019). Freshwater biodiversity was assessed within the CAFF geographic area (Figure 1) for all Arctic countries, and included data from mainland regions of the Arctic as well as islands such as Greenland, Iceland, Faroe Islands, Svalbard (Norway), Wrangel Island (Russia), and islands in the Canadian Arctic Archipelago. Description of the environmental regimes across the circumpolar region can be found in Culp et al., (2022).
Benthic macroinvertebrate data and supporting variables for lakes and rivers were compiled as part of the extensive circumpolar CBMP-Freshwater Database (housed at CAFF's open data portal, the Arctic Biodiversity Data Service; abds.is). In Fennoscandia (Finland, Norway, and Sweden), data were primarily acquired from national monitoring programmes of lakes and rivers. In Canada, river macroinvertebrate data were available from the Canadian Aquatic Biomonitoring Network (CABIN) of Environment and Climate Change Canada; however, few lake monitoring data were available through this programme. A subset of data from Alaska originated from the National Aquatic Resource Surveys funded by the U.S. Environmental Protection Agency. Remaining data from Canada, U.S.A., and Fennoscandia, as well as data from other countries and regions (Faroe Islands, Iceland, Greenland, Russia, and Svalbard), primarily originated from academic research and monitoring related to industry. See Appendix S1 for more details.

| Supporting variables
Water quality and physical habitat data were available for approximately half the sites, and only conductivity and total phosphorus were consistently collected, making it necessary to supplement the analysis with geospatial variables. We used climatic and land cover geospatial datasets that covered the entire circumpolar region (see Table S1 and Appendix S1 for more details). Climate data were summarised as the long-term average (LTA;1970-2000 maximum August air temperature, as a measure of summer high temperatures, and the LTA annual precipitation coefficient of variation (CV), as a standardised measure of variability in precipitation (from WorldClim Version 2; http://world clim.org/version2; Fick & Hijmans, 2017). Air temperature was used as a proxy for water temperature as recommended by Yang et al., (2021), as water temperature data were generally not available. Permafrost extent was summarised from the Circum-Arctic Map of Permafrost and Ground-Ice Conditions (Version 2; https://nsidc.org/data/ggd318; Brown et al., 2002). From the same geospatial layer, the relative area of glacier cover was calculated (Table S1). Bedrock geology data were available for most countries from the circumpolar Geological Map of the Arctic (Harrison et al., 2011) that extends south to 60°N, and data for North American river sites below that latitude were obtained from a North American bedrock geology layer (Garrity & Soller, 2009). We ensured similar data were obtained for sites within the boundary of both data layers (Alaska, northern Canada, Greenland, and Iceland), and only calculated the relative area of geology classes that were common to both layers for the analysis (see Table S1 for details).
Geospatial variables were summarised across catchments from existing global hydrological basin layers (Lehner & Grill, 2013). These

| Sample selection
Subsets of lake or river sites were selected for analysis based on sampling methods and sample timing (see Figure S1 and Appendix S1 for details). Despite the variety of data sources, there was a reasonable degree of consistency in sampling methods across datasets.
In lakes, macroinvertebrates were generally collected from littoral rocky habitats with similar sampling approaches (dip net or kick net, Surber sampler, or stone scrubs/scrapes; see Table S2 for a summary of methods used in each country). Canadian lake data were omitted because littoral sampling was rare and was focused on soft sediment sampling using grab samplers. Mesh size for lake samples ranged from 230 to 500 μm among data sources (Table S2). River macroinvertebrates were collected with a kick net using a time-limited approach in most countries (although samples in Greenland, Iceland, and Russia were collected with Surber samplers or stone scrubs/ scrapes; see Table S2). Similar to lakes, mesh size ranged from 250 to 500 μm (Table S2). There were some regional differences with respect to measurement of sampling effort (e.g., whether samples were area-limited or time-limited) and sample processing methods (e.g., subsample size), which commonly differ among countries (Buss et al., 2015). The effects of this regional variation on diversity estimates were minimised by: (1) restricting data to those collected by particular sample devices to ensure similar habitats and macroinvertebrate groups were targeted (as described above); and (2) converting data to presence/absence for all analyses rather than using total/ relative counts or densities.
To minimise temporal influence on our spatial analysis, only sites with data from 2000 or later were retained. The exceptions were lake sites in Greenland and Iceland, where the majority of sites were sampled in the 1980s and 1990s, respectively, and river sites from the low Arctic in Greenland, which were sampled in the 1980s.
These data were included to ensure that all countries were represented, but were interpreted with caution, as they represented older baseline data. Time of year for sample collection varied across datasets, largely due to differences in timing of ice-off and emergence across latitudinal and longitudinal gradients. While this resulted in sampling times across the circumpolar region ranging from June to early October, samples were collected within a relatively short period (1 month) within each specific geographic region (i.e., sites at similar latitudes/longitudes). Moreover, repeated sampling at a site (within or among years) was rare for both lakes (8.8% of sites) and rivers (9.2% of sites). To ensure there was no temporal influence on our spatial analysis, only data from the most recent sample date (or most recent sample date with supporting water chemistry data) were included in the analyses.

| Taxonomic harmonisation
Taxonomic nomenclature was updated to correct outdated naming conventions, remove taxa that were not consistently identified and enumerated in all countries (i.e., Acari, Collembola, Hydra, Nematoda, Turbellaria), combine or remove data with mixed-level taxonomy, and harmonise differences in nomenclature across datasets, ensuring data comparability across the circumpolar region (see Appendix S1 for more details). Only samples with insects identified to genus or family were retained, and genus-level data were combined at the family level for insects (family or lower taxonomic resolution for non-insects) prior to analysis. Although genus-level analysis would have provided a more accurate assessment of diversity, exclusion of data that were at family level would have severely limited the spatial coverage of the dataset and resulted in the exclusion of entire countries. Preliminary assessment of a subset of sites with genus-level data indicated that although diversity estimates were higher than at the family level, spatial patterns in diversity were similar. Family-level data retain ecological response patterns (sensu Bowman & Bailey, 1997); however, such data must be interpreted conservatively as higher taxonomic levels may not represent species level-responses to environmental gradients (Heino, 2014).

| Rarefied α diversity
To compare α diversity spatially despite differences in sample density (i.e., samples per region) across the Arctic, we assessed rarefied α diversity, which we defined as the taxa richness of lakes or rivers within a catchment (hydrobasin), rarefied to the same number of sites in all catchments. We conducted spatial analysis of rarefied α diversity among hydrobasins as a function of latitude and temperature (see Figure S1 and description below).

| Grouping and rarefying data
Rarefied α diversity was calculated for a set number of sites within each hydrobasin. Rarefaction is a prerequisite for sound largescale comparisons of diversity (Colwell et al., 2004;Gotelli & Colwell, 2001), as taxa richness for any given area increases with an increasing number of samples (until an asymptote is reached; Gotelli & Colwell, 2001). Sites were grouped by level 05 hydrobasins and rarefaction curves were created in EstimateS (Colwell, 2013;Colwell & Elsensohn, 2014) by estimating diversity at increasing numbers of sites through a randomisation procedure with 100 iterations. Rarefaction curves were extrapolated (see details in Colwell et al., 2004) if necessary to ensure they exceeded a minimum of 10 nodes (i.e. a diversity estimate at 10 sites), and only hydrobasins with four or more sites were retained in the analysis to avoid extrapolation from a small number of sites. In total, 346 lake and 1,079 river sites were retained in rarefied α diversity analysis ( Figure 1).

| Spatial analysis
We assessed rarefied α diversity in relation to geographic location, isolation, and temperature to assess latitudinal and thermal gradients in diversity while accounting for the influence of dispersal barriers. The mean rarefied α diversity estimated to be found at 10 sites (±SE) was initially plotted for each hydrobasin as a function of the average latitude of all sites in the hydrobasin to explore latitudinal patterns. To evaluate whether latitudinal trends in rarefied α diversity were driven by connectivity to the mainland, we tested lakes and rivers separately with the Analysis of Covariance (ANCOVA) model: diversity = latitude + isolation factor + latitude * isolation factor, where the isolation factor was a categorical term indicating whether a hydrobasin was on an island or on the mainland. The interaction term (latitude * isolation factor) tested whether the slope of the relationship between rarefied α diversity and latitude differed between island and mainland hydrobasins. If the term was not significant (at α = 0.05), a reduced model was run to focus on differences in mean diversity among groups (categorical term). If neither the interaction nor the categorical term was significant, a least-squares linear regression model (diversity = latitude) was run with data pooled for mainland and island hydrobasins. The relationship between rarefied α diversity and long-term average maximum August air temperature was also tested for lake and river mainland and island hydrobasins.
Because mainland and island hydrobasins covered different temperature ranges with little overlap, separate least-squares linear regressions of the diversity-temperature relationship were run for these two groups (model: diversity = temperature). Prior to analyses, residuals were examined to ensure assumptions of ANCOVA and regression were met. Statistical analysis of rarefied α diversity was run in Systat 12 (San Jose, CA, U.S.A.) and plots were created using Systat 12 and the ggplot2 package (Wickham, 2016) in R Studio version 3.6.0 (R Development Core Team, 2015).

| Community composition
Assessment of compositional differences included regional assessments of dissimilarity using β diversity analysis and partitioning of β diversity into its component parts (nestedness and turnover), as well as a site-scale evaluation of taxonomic composition and its response to environmental drivers (see steps in Figure S1 and description below). Community composition analysis was conducted using 372 lake and 1,164 river sites ( Figure 1).

| Classifying data
We assessed large-scale patterns in community composition across the Arctic by calculating β diversity and its component parts at a regional scale, i.e., through pairwise dissimilarity of geographic regions. 2.5.2 | Regional analysis of β diversity and its components Pairwise dissimilarity between geographic regions (listed above) was estimated using β diversity analysis methods. Data for sites in a geographic region were combined, resulting in a single set of presence/ absence data for each region. Because differences in sampling effort among regions may have affected β diversity estimates (e.g., if one region included more taxa due to a greater number of sites being sampled), subsets of 20 sites were selected at random and combined for analysis for any regions with data from more than 20 sites. Thus, the combined data used for β diversity analysis represented the full composition of taxa across 20 or fewer sites in each region. Beta diversity was estimated as Sørenson dissimilarity (β SOR ) through pairwise comparisons of geographic regions (following Baselga, 2010;. Regional β diversity was also divided into its component parts of nestedness and turnover to evaluate which contributed the most to observed pairwise differences among regions. Nestedness reflects regions that contain a subset of the taxa found within other regions (i.e., a reduced taxa pool with no new taxa beyond those found in the more taxa-rich region), while turnover describes compositional differences among regions that result from finding unique taxa in each region. Following Baselga (2010) and , we calculated the share of β diversity that was due to nestedness and turnover as the nestedness-resultant component of Sørenson dissimilarity (β NES ) and as Simpson dissimilarity (β SIM ), respectively. The proportion of β diversity that was due to nestedness and turnover for each pairwise comparison was used to create heat maps that visualised differences among regions for lakes and rivers. Proportions were calculated to compare the relative importance of each component while controlling for regional differences in total β diversity. Beta diversity and its components were calculated in R Studio version 3.6.0 (R Development Core Team, 2015) using the beta.pair function in the package betapart ).

| Site-level multivariate analysis of composition and drivers
Dissimilarity among sites was assessed within and among regions through a site-level, multivariate assessment of spatial compositional patterns and associated environmental drivers using linear methods (see steps in Figure S1), which are appropriate when there is low turnover across samples (Legendre & Legendre, 2012). As both lake and river data contained large numbers of zeroes, data underwent a Hellinger transformation, which removes bias in ordinations of community data with Euclidean distance (Legendre & Gallagher, 2001).
We then used principal components analysis (PCA) to separately analyse spatial patterns in community composition for lakes and rivers, to identify major gradients in composition and dominant taxa across the circumpolar region.
We used redundancy analysis (RDA) to relate the Hellingertransformed community data to environmental variables and identify correlates of community composition patterns identified in the PCA. Analysis was run with: (1) the full set of lake and river sites (372 and 1,164 sites, respectively) and only geospatial variables (see Table S1 for list of chosen variables); and (2) the subset of lake and river sites (91 and 440, respectively) that had both geospatial variables and site-specific water chemistry variables (conductivity and total phosphorus). The RDA used centring and standardisation to account for differences in scale of environmental variables, and all variables were retained in the model. Monte Carlo permutation tests were used to determine the significance of each RDA axis and of individual environmental variables. For tests of variable significance, marginal effects (the individual effects of each variable with all others held constant) were used to determine the order of inclusion in the model and significance was tested on conditional effects of environmental variables. All ordinations were run in Canoco Version 4.55 (ter Braak & Šmilauer, 2002).

| Rarefied α diversity
Rarefied α diversity of lake communities ranged from 7 to 46 taxa (estimate of the diversity found in 10 sites), and generally declined with increasing latitude above 65°N (Figure 2a). However, diversity patterns in lakes differed markedly between mainland hydrobasins and islands (ANCOVA, interaction p < 0.05, Table 1 Table 2A).
Rarefied α diversity of river communities ranged from 2 to 55 taxa, and declined consistently with increasing latitude (Figure 2c).
River diversity was low in most hydrobasins on islands, with a more gradual decline above 55°N (Figure 2c). The slope of the relationship between rarefied α diversity and latitude did not differ among mainland and island river hydrobasins (ANCOVA interaction p = 0.61, Table 1), but mean diversity was significantly higher for mainland hydrobasins (intercept p < 0.001; Table 1). A combined regression of diversity as a function of latitude that included all mainland and island river hydrobasins (to test for common slope) revealed a significant decline in diversity with increasing latitude that was equivalent F I G U R E 2 Average (±SE) rarefied α diversity at 10 sites for all hydrobasins with four or more sites plotted as a function of average latitude (°N) in (a) lakes and (c) rivers, and as a function of long-term average (LTA) maximum August temperature (°C) in (b) lakes and (d) rivers. Colours and shapes of hydrobasin symbols indicate geographic location, as described in the legends. Latitude and temperature range along the x-axis differs between plots to allow for better visual assessment of lake trends. The temperature axis is plotted from high to low temperatures for consistency with latitudinal plot to a loss of approximately six taxa with an increase of 5° latitude (slope = −1.29; t = −5.73, p < 0.001).
Rarefied α diversity for both lakes and rivers was more strongly related to air temperature than to latitude, with clear declines related to lower long-term average maximum August air temperature (Figure 2b,d). For lakes, island hydrobasins had lower temperatures than most mainland hydrobasins (Figure 2b), which contributed to the strong relationship of diversity with temperature. Diversity declined significantly with decreasing temperature in mainland hydrobasins (regression slope = −2.99; Table 2B, Figure 3a), equivalent to a loss of 15 taxa with a decline of 5°C. Diversity in island hydrobasins was lower at the coldest temperatures (Figure 3c), but the slope of this relationship (−0.43) was not significant (Table 2B). For rivers, rarefied α diversity strongly declined with decreasing air temperature across the full temperature gradient covered by mainland and island hydrobasins (20°C to 1°C; Figure 2d). For mainland river hydrobasins, which covered a temperature gradient from 20°C to 10°C, the slope of the relationship of diversity with temperature was −2.48, equivalent to a loss of approximately 12 taxa with a decline of 5°C (Table 2; Figure 3b). Island hydrobasins, which covered a temperature range from 10°C to 1°C, had a shallower slope of −1.17 (Table 2; Figure 3d), equivalent to a loss of approximately six taxa with a decrease of 5°C.

| Regional β diversity patterns and partitioning
Beta diversity (β SOR ) for lakes ranged from 0.22 to 0.74 among geographic regions, indicating moderately low to high dissimilarity in composition. Dissimilarity between regions was highest on average in comparisons between islands (i.e., Iceland, Greenland, Faroe Islands, or Wrangel Island) and mainland (i.e., Fennoscandia, northwest Russia, and U.S.A.; mean ± SD β SOR = 0.61 ± 0.07), and lowest between mainland regions within Fennoscandia (mean ± SD β SOR = 0.38 ± 0.14), where geographic connectivity between regions is high.
Beta diversity between regions was more variable for rivers than for lakes, ranging from extreme similarity (β SOR = 0.09) to extreme dissimilarity (β SOR = 0.93). The strongest dissimilarity was between mainland and island regions (mean ± SD β SOR = 0.72 ± 0.15), particularly between Svalbard and mainland Fennoscandia or North America (β SOR ranging from 0.86 to 0.93), whereas strong similarity was found in comparisons between any two mainland regions (mean ± SD β SOR = 0.30 ± 0.12).
Beta diversity partitioning of lake data revealed that most regions had a moderate to high percent contribution of nestedness in comparisons with Finland, Sweden, and northwest Russia (Table 3A), which indicated that island regions largely contained a subset of the taxa found in these mainland regions. Lakes in Alaska shared a high degree of nestedness with the islands in closest proximity to North America (Greenland and Wrangel Island; Table 3A). Comparisons among all island regions were dominated by the turnover component of β diversity (mean percent turnover in island comparison = 83%), indicating that compositional differences among islands were primarily due to replacement of taxa.
Beta diversity in rivers was predominantly attributed to nestedness in most comparisons of Arctic islands with mainland or island regions (83% of β diversity on average), particularly comparisons involving Svalbard, Greenland, Ellesmere Island, and Wrangel Island (Table 3B), which suggested that islands generally contained a subset of the taxa found in mainland regions. Turnover contributed more to β diversity in comparisons between mainland regions. TA B L E 1 Results of lake and river ANCOVA models comparing the relationship between α diversity and latitude among mainland and island hydrobasins (island factor). The degrees of freedom, F-ratio, and p-value are reported for the full model interaction term, which tests whether slopes differ between mainland and island hydrobasins, and for the reduced model (without interaction) categorical term, which tests whether intercepts differ between mainland and island hydrobasins. Reduced model results are only presented when the interaction term in the full model was not significant (at α = 0.05), as the reduced model cannot be calculated when there is a significant difference in slopes

| Community compositional variation in relation to environmental drivers
Lake ordinations of community composition revealed marked geographic differences and notable separation of three groups of regions: Analysis of biotic-abiotic associations for lakes with both geospatial and water chemistry data indicated that neither TP nor conductivity contributed significantly to distinguishing among sites ( Figure S2b). As gradients in TP (range 0.001-0.074 mg/L, with one outlier) and conductivity (range 7.6-316 μS/cm, with six outliers) were relatively short and many sites lacked water chemistry data, we focused on assessing the response to geospatial correlates using all lake sites. Air temperature dominated the first RDA axis (which The PCA of all river sites separated more diverse sites in Fennoscandia, Canada east mainland, and Canada west mainland from the less diverse island regions (Canada east and west island, Greenland, Iceland, Svalbard, Wrangel Island) and Alaska along the first axis (which explained 13.7% of composition variation; Figure S3a). Most taxa were associated with these lower-latitude mainland sites at the positive end of the first axis, whereas the higher-latitude island sites were primarily associated with coldtolerant taxa such as Chironomidae, Oligochaeta, and Tipulidae.

F I G U R E 3
Linear regressions of average rarefied α diversity at 10 sites as a function of long-term average (LTA) maximum August temperature (°C) for hydrobasins in (a) mainland lakes, (b) mainland rivers, (c) island lakes, and (d) island rivers. Solid lines are least-squares linear regression lines and dashed lines are 95% confidence intervals, the latter of which is presented only for slopes that differ from zero. Temperature axis is plotted from high to low temperatures for consistency with latitudinal plot TA B L E 3 Comparison of (A) lake and (B) river benthic macroinvertebrate community composition among geographic regions, indicating the proportion of pairwise β diversity (differences in community composition) that is attributable to turnover (upper portion of table, shaded from low turnover in light green to high turnover in dark green) and nestedness (lower portion of table, shaded from low nestedness in light blue to high nestedness in dark blue). A high proportion of turnover (approaching 1) indicates that differences between regions are due more to replacement of taxa, whereas a high proportion of nestedness indicates that differences between regions are due more to loss of taxa. Proportions were used to control for regional differences in total β diversity. Regions in (A) are arranged from east to west, and regions in (B) are arranged as mainland regions east to west followed by island regions east to west The second PCA axis (explaining 9.1% of variance), showed a separation of many eastern and western Canada mainland sites, with the latter associated with families of stoneflies including Nemouridae, Capniidae, Chloroperlidae, and Perlodidae, as well as the mayfly family Heptageniidae ( Figure S3a). The PCA of the subset of sites that had water chemistry data indicated similar patterns, with a strong separation of Fennoscandia and Canada mainland sites from island and Alaskan sites along the first axis (explaining 17% of the variation; Figure 5a). Taxon associations along the first axis were similar, with the additional association of Simuliidae with higherlatitude sites. The second PCA axis (explaining 9.4% of the variation) provided a similar separation of eastern and western Canada mainland sites (Figure 5a). In both the ordination of all sites and the reduced subset with water chemistry data, Fennoscandian sites were spread across much of the second axis gradient, indicating similarity to both mainland Canada east and Canada west.
Water chemistry variables played a larger role in differentiating river communities, and both TP (range 0.001-0.740 mg/L, with one outlier) and conductivity (range 7.6-884 μS/cm, with one outlier) were significant in the RDA model. Therefore, assessment of biotic-abiotic relationships for rivers focused on the subset of 440 sites with water quality variables. This RDA showed that temperature was the most important variable (conditional effect: F = 11.09, p = 0.002), and was highly correlated with the first axis, which explained 47.2% of the constrained and 3.5% of the unconstrained variance. High air temperature was positively associated F I G U R E 4 (a) Principal components analysis ordination of lake benthic macroinvertebrate communities and (b) redundancy analysis of lake benthic macroinvertebrate communities and geospatial variables, with sample symbols coloured by geographic region. A selection of taxon codes (for ease of interpretation) is located on each ordination near where taxon points plotted (see Table S3 for codes). Axis labels indicate (a) % variance in communities and (b) % constrained variance in communities explained by axes. Geospatial variables: Glaciers, % glacier; IgneousExt, % extrusive igneous bedrock; P-Continuous, % continuous permafrost; P-DisCont, % discontinuous permafrost; P-Isolated, % isolated permafrost; PrecipCV, long-term average precipitation CV; Sediment, % sedimentary bedrock; Supracrust, % supracrustal bedrock; TempMaxAug, long-term average maximum August temperature with most taxa, and separated relatively summer-warm regions in Fennoscandia and Canada east and west mainland from cooler regions on higher-latitude islands and in Alaska (Figure 5b). Highlatitude sites were also positively associated with variability in precipitation, continuous permafrost, specific conductivity, and relative area of glaciers, all of which had significant conditional effects in the model (F = 1.73, 3.31, 5.26, and 3.2, respectively, p ≤ 0.004 for all; Figure 5b). Along the second axis, which explained 17.7% of the constrained and 1.3% of the unconstrained variance, sites were primarily separated based on the relative area of isolated permafrost, which was negatively associated with some Ephemeroptera, Plecoptera, and Trichoptera taxa (Figure 5b). Patterns in the RDA of all sites were similar, but without conductivity and TP in the ordination, permafrost coverage and geology played a larger role in differentiating among sites ( Figure S3b).

| D ISCUSS I ON
Our spatial analysis of over 1,500 Arctic lake and river sites revealed that air temperature, as a proxy for water temperature, was the most important correlate of macroinvertebrate diversity, but also pointed to the importance of biogeographical patterns. Rarefied α diversity of benthic macroinvertebrates in both lakes and rivers declined with increasing latitude, with this pattern strongly linked to temperature.
Rates of latitudinal decline in diversity changed at approximately 70°N, which coincides with the northernmost border of mainland regions in the Arctic. We therefore conjecture that diversity of remote, high latitude Arctic islands is limited by dispersal processes. Similar to trends in rarefied α diversity, compositional differences among both lakes and rivers related most strongly to summer air temperature, with colder regions predominated by cold-tolerant taxa. Regional dissimilarity was F I G U R E 5 (a) Principal components analysis ordination of river benthic macroinvertebrate communities and (b) redundancy analysis of river benthic macroinvertebrate communities and environmental variables (water chemistry and geospatial), with sample symbols coloured by geographic region. A selection of taxon codes (for ease of interpretation) is located on each ordination near where taxon points plotted (see Table S3 for taxon codes). Note that orders are presented where a large number of family points plotted, for ease of interpretation.
Axis labels indicate (a) % variance in communities and (b) % constrained variance in communities explained by axes. Environmental variables: Glaciers, % glacier; P-Continuous, % continuous permafrost; P-DisCont, % discontinuous permafrost; P-Isolated, % isolated permafrost; PrecipCV, long-term average precipitation CV; Sedimentary, % sedimentary bedrock; SpCond, specific conductivity; Supracrust, % supracrustal bedrock; TempMaxAug, long-term average maximum August temperature; TP, total phosphorus consistently highest when mainland regions were compared with island regions, which suggests that poor connectivity also contributed to compositional differences. Together, these results have implications for the type and rate of change that might be expected in Arctic freshwaters with continued warming.

| Alpha diversity patterns are mainly driven by temperature (H 1 )
Rarefied α diversity of lakes and rivers exhibited a strong decline with increasing latitude, most evident above 65-70°N, consistent with previous studies of latitudinal trends in macroinvertebrate richness (e.g., Culp et al., 2019;Oswood, 1989;Scott et al., 2011). Latitudinal declines in diversity have been described across the globe (Hillebrand, 2004), and mechanisms proposed to explain these widespread patterns include hypotheses of physiological tolerance (Currie et al., 2004). In our study, the major driver associated with rarefied α diversity was temperature, a finding that is similar to earlier studies across wider gradients from equatorial to temperate latitudes (Jacobsen et al., 1997) and that supports hypotheses of diversity limitation due to physiological tolerance constraints (e.g., inability of most taxa to survive in the harsh environmental conditions at high latitudes). The observed lower diversity at high latitudes and in colder regions of the Arctic supports our first hypothesis (H 1 ) and underscores well-known physiological and ecological relationships of ectothermic invertebrates, e.g., that the basal metabolic rate of ectotherms scales with temperature (Woodward et al., 2010) and that this relationship differs among taxa .
Temperature has been shown to play a key role in structuring macroinvertebrate communities in Arctic and alpine regions (Lencioni, 2018;Milner et al., 2001;Niedrist & Füreder, 2021), and where temperatures are below the lower physiological thresholds of most taxa, diversity may be severely limited. Variation in temperature along the latitudinal gradient in our study resulted in a strong spatial pattern in macroinvertebrate diversity, but variability in the diversity-temperature relationship at particular latitudes appeared to reflect longitudinal differences in thermal regimes, i.e., climatic variation among geographic regions at similar latitudes. For example, some coastal areas were warmer in summer than inland areas at similar latitudes.
Furthermore, the diversity-latitude relationship reflected historical differences in warming across the Arctic, as diversity at similar latitudes was higher in western Canada than in eastern Canada, where temperatures have historically remained more stable and colder .
The limited taxa found at the highest latitudes and coldest regions in our study were cold-tolerant taxa, consistent with our prediction for H 1 .
For example, the highest latitudes were dominated by extremely coldtolerant Chironomidae and Oligochaeta (c.f. Milner et al., 2001;Niedrist & Füreder, 2021), both of which were associated with the lowest temperatures and highest percent glacier cover in the hydrobasin. Moreover, Chironomidae made up more than 85% of the macroinvertebrate populations in rivers on Ellesmere Island (highest latitude in the Canadian archipelago), Svalbard (high latitude islands in Norway), and Greenland.
Several studies of glacier-fed systems Lencioni, 2018;Milner et al., 2001;Niedrist & Füreder, 2021), as well as high-latitude lakes and rivers in Russia (Baturina et al., 2020;Baturina & Loskutova, 2010;Kondratjeva et al., 2014), have noted the high tolerance to extreme cold conditions by Chironomidae (especially subfamilies Diamesinae and Orthocladiinae), Oligochaeta, Tipulidae, and Simuliidae, all of which were also associated with colder, high-latitude rivers in our study. Furthermore, taxa of Ephemeroptera, Plecoptera, and Trichoptera were correlated with higher maximum August temperatures across the circumpolar region, consistent with patterns observed in regional studies Johnson & Goedkoop, 2002;Scott et al., 2011). The observed shifts in both rarefied α diversity and community composition across thermal gradients indicate the strength of temperature as a structuring force in the Arctic.

| Spatial discontinuity limits α diversity (H 2 )
At the northern end of the latitudinal gradient, observed declines in rarefied α diversity were likely to be related to both temperature and dispersal limitation to islands (i.e., marine barriers), because the poleward expansion of invertebrate species distributions is restricted by spatial discontinuities (Heino et al., 2009;Hickling et al., 2006;Medeiros et al., 2020;Thomas et al., 2001). Such barriers are especially prevalent in North America, where the Canadian Arctic archipelago creates dispersal limitations; however, these barriers also exist among islands off the Eurasian north coast where the remoteness of the few islands in the Arctic Ocean is even larger. In our study, the observed decoupling of the macroinvertebrate diversity relationship with latitude was related to geographic location, in support of our second hypothesis (H 2 ). We hypothesise that the lower levels of lake and river macroinvertebrate diversity observed on islands reflect limitations imposed by dispersal. Dispersal capability varies widely among macroinvertebrate taxa, and is dependent upon factors such as body size, reproduction, locomotion, and dispersal strategies of different life stages (active and/or passive as juveniles or adults; Sarremejane et al., 2020). Dispersal capability and spatial distance between suitable regions and habitats limit the climate-change-induced rate at which species may move northward (Heikkinen et al., 2006;Thomas et al., 2001), although the extent of this limitation in different taxonomic groups is difficult to measure and remains largely unclear (Heino et al., 2009;Hickling et al., 2006;Sarremejane et al., 2020).
Such biogeographically-driven dispersal limitation will delay or impede the response of species distributions to climate warming in island regions compared to regions with high connectivity.

| Community composition reflects temperature differences and geographic separation (H 3 )
Spatial differences in composition of macroinvertebrates appeared to reflect regional proximity and isolation, and the highest pairwise dissimilarity among regions for both lakes and rivers was found between mainland and island regions, consistent with our third hypothesis (H 3 ). We observed a high contribution of nestedness to dissimilarity between island and mainland regions across the Arctic, consistent with our prediction for H 3 . In part, taxa that were predominant on islands at high latitudes (e.g., cold-tolerant groups) were also ubiquitous across the entire circumpolar region (albeit with low taxonomic resolution), contributing to the observed high degree of nestedness. Dissimilarity was lower when estimated between two mainland regions, which suggested that the higher degree of spatial connectivity and regional proximity contributed to stronger similarity in taxa and compositional patterns. The loss of taxa at high latitudes and on islands, coupled with the predominance of cold-tolerant taxa in these depauperate regions supports the predicted dual roles of temperature and geographic connectivity as drivers of Arctic macroinvertebrate diversity.

| Implications of temperature and spatial discontinuity relationships
The strong relationship between temperature and macroinvertebrate diversity in both lakes and rivers, along with the predicted continued warming of Arctic freshwaters, emphasise the need to monitor diversity and ecological change, including the early establishment of invasive species and loss of endemic species   Shah et al., 2014). Furthermore, the longer growing season and potential for increased nutrient levels at high latitudes (e.g., Huser et al., 2020) may lead to higher primary production within Arctic freshwater systems, further bolstering the food webs despite different responses among trophic levels (Lau et al., 2020).
We conjecture that early examples of such change (e.g., species migrations, vegetation change) can first be seen in the sub-Arctic zone (e.g., Johnson & Goedkoop, 2002;Lau et al., 2020).
Predicted scenarios of biodiversity change in response to warming indicate that changes in taxonomic richness in the Arctic will depend on differences between the rates of northward range expansion of taxa from warmer southern regions (Shah et al., 2014) and the extirpation of cold stenotherms (Bálint et al., 2011;Culp et al., 2012;Lento et al., 2019). Given the importance of spatial discontinuity in our study, the rate of loss of cold stenotherms may depend upon the magnitude of warming and the dispersal capability of the species, particularly for those on islands that may be unable to move to more suitable habitats. Several cold stenothermic species may become extinct with continued warming of Arctic regions, but our family-level analysis showed that compositional differences between low and high latitude river communities at coarser taxonomic resolutions were largely driven by nestedness, suggesting that high Arctic macroinvertebrate families were also found at lower latitudes. Thus, while these families may have adapted to and predominate in extreme cold conditions, they may not require those conditions to survive (Muhlfeld et al., 2020). With continued warming, macroinvertebrate communities may be expected to experience a net increase in family-level richness as taxa from warmer regions move north (Shah et al., 2014), although taxonomic richness at the species level may more strongly reflect losses of cold stenothermic taxa. The northward movement of taxa may also contribute to lower dissimilarity among regions, as northern community composition begins to more closely resemble that at lower latitudes. The arrival of taxa from lower latitudes will probably contribute to the loss of the community composition that is unique to the Arctic, including communities predominated by cold-tolerant species that are poor competitors in a warming environment and that are easily preyed upon (Flory & Milner, 1999). As cold-tolerant species are disfavored by warmer conditions, the expected increase in diversity is not necessarily positive, as unique species will go locally extinct.
Within our family-level analysis, we identified clear patterns in diversity and community composition across the Arctic. Although it was not unexpected to find such compositional differences by geographic region, the extent to which these regions differed at the family level was greater than might be anticipated, given the relatively coarse taxonomic level. The strength of our results indicates that spatial variation is visible even without high taxonomic precision, and suggests that a response to climate change will be evident even as broad shifts in macroinvertebrate communities. Our results emphasise the strength of these regional patterns as well as the utility of family-level assessments. Although assessment of a subset of data indicated that diversity patterns with genus-level data were similar to those found with family data (albeit with higher diversity; results not shown), an advantage of genus-level analysis would be the opportunity to detect specific cold-stenotherm taxa that may be at risk from warming, particularly species of Chironomidae, which is a diverse and ecologically important group in the Arctic (e.g., Medeiros & Quinlan, 2011;Oliver & Dillon, 1997). In the future, the application of DNA barcoding techniques could improve the taxonomic resolution (and diversity estimates) in Chironomidae, as well as in other complex groups including the Tipulidae, Oligochaeta, and Simuliidae.
Our study identified temperature and spatial connectivity as major drivers of spatial patterns in benthic macroinvertebrate diversity and community composition. Our study is the most extensive assessment to date of spatial diversity patterns in Arctic freshwater benthic macroinvertebrates, and provides a baseline to support continued assessments of climate change effects. Given the observed strong relationship of both diversity and community structure with temperature, we conclude that continued warming in the Arctic will enforce community shifts toward less dominance of cold-tolerant taxa and higher α diversity. However, our observations of strong differences between mainland and island regions suggest that barriers to dispersal will limit the rate at which these changes occur. Our circumpolar spatial analysis suggests that continued warming and its associated impacts (e.g. land use change, exploitation) will alter the unique biodiversity of Arctic freshwaters.

ACK N OWLED G EM ENTS
J.L. initiated and designed the study, coordinated data collection, contributed data, harmonised data, conducted analysis, created tables/figures and drafted the paper. J.M.C. and W.G. initiated and designed the study, contributed data, and drafted the paper. B.L.
contributed data, harmonised data, conducted data analysis, and contributed to writing/editing paper. All other authors contributed data and contributed to writing/editing paper. The authors wish to thank the CAFF Secretariat and the co-leads of the CBMP for their support. Thank you to Jani Heino for reviewing and improving the manuscript through constructive comments and editing. The authors also thank three anonymous reviewers and the Associate Editor for their suggestions that greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
The authors assert that there are 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 Arctic