# Highly variable species distribution models in a subarctic stream metacommunity: patterns, mechanisms and implications

Guillermo de Mendoza, Riikka Kaivosoja, Mira Grönroos, Jan Hjort, Jari

Ilmonen, Olli-Matti Kärnä, Lauri Paasivirta, Laura Tokola, Jani Heino

► To cite this version:

Guillermo de Mendoza, Riikka Kaivosoja, Mira Grönroos, Jan Hjort, Jari Ilmonen, et al. Highly variable species distribution models in a subarctic stream metacommunity: patterns, mechanisms and implications. Freshwater Biology, 2018, 1 (63), pp.33-47. 10.1111/fwb.12993. hal-02450949

HAL Id: hal-02450949

https://univ-tlse2.hal.science/hal-02450949

Submitted on 23 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Highly variable species distribution models in a subarctic stream
metacommunity: patterns, mechanisms and implications

Guillermo de Mendoza¹,²,*, Riikka Kaivosoja³, Mira Grönroos⁴, Jan Hjort³, Jari Ilmonen⁵,

Olli-Matti Kärnä³, Lauri Paasivirta⁶, Laura Tokola³ & Jani Heino⁷

¹Centre for Advanced Studies of Blanes, Spanish National Research Council (CEAB-CSIC),

Accés a la Cala St. Francesc 14, ES-17300 Blanes, Spain.

²Laboratoire GEODE, UMR 5602 CNRS, Université Toulouse-Jean Jaurès, 5 allées Antonio

Machado, FR-31058 Toulouse, France.

³University of Oulu, Geography Research Unit, P.O. Box 3000, FI-90014 Oulu, Finland.

⁴University of Helsinki, Department of Environmental Sciences, Section of Environmental

Ecology, Niemenkatu 73, FI-15140 Lahti, Finland.

⁵Metsähallitus, Natural Heritage Services, P.O. Box 94, FI-01301 Vantaa, Finland.

⁶Ruuhikoskenkatu 17, FI-24240 Salo, Finland.

⁷Finnish Environment Institute, Natural Environment Centre, Biodiversity, Paavo Havaksen

Tie 3, FI-90570 Oulu, Finland.

* Corresponding author: guillermo.de-mendoza@univ-tlse2.fr
20 **Summary**

21 1. Metacommunity theory focuses on assembly patterns in ecological communities, 22 originally exemplified through four different, yet non-exclusive, perspectives: patch 23 dynamics, species sorting, source-sink dynamics, and neutral theory. More recently, three 24 exclusive components have been proposed to describe a different metacommunity 25 framework: habitat heterogeneity, species equivalence, and dispersal. Here, we aim at 26 evaluating the insect metacommunity of a subarctic stream network under these two 27 different frameworks.

28 2. We first modelled the presence/absence of 47 stream insects in northernmost Finland 29 using binomial generalised linear models (GLMs). The deviance explained by pure local 30 environmental (E), spatial (S), and climatic variables (C) was then analysed across 31 species using beta regression. In this comparative analysis, site occupancy, as well as 32 taxonomic and biological trait vectors obtained from principal coordinate analysis, were 33 used as predictor variables.

34 3. Single-species distributions were better explained by in-stream environmental and spatial 35 factors than by climatic forcing, but in a highly variable fashion. This variability was 36 difficult to relate to the taxonomic relatedness among species or their biological trait 37 similarity. Site occupancy, however, was related to model performance of the binomial 38 GLMs based on spatial effects: as populations are likely to be better connected for 39 common species due to their near ubiquity, spatial factors may also explain better their 40 distributions.

41 4. According to the classical four-perspective framework, the observation of both 42 environmental and spatial effects suggests a role for either mass effects or species sorting 43 constrained by dispersal limitation, or both. Taxonomic and biological traits, including 44 the different dispersal capability of species, were scarcely important, which undermines
the patch dynamics perspective, based on differences in dispersal ability between species.

The highly variable performance of models makes the reliance on an entirely neutral

framework unrealistic as well. According to the three-component framework, our results

suggest that the stream insect metacommunity is shaped by the effect of habitat

heterogeneity (supporting both species-sorting and mass effects), rather than species

equivalence or dispersal limitation.

5. While the relative importance of the source-sink dynamics perspective or the species-sorting paradigm cannot be deciphered with the data at our disposal, we can conclude that habitat heterogeneity is an important driver shaping species distributions and insect assemblages in subarctic stream metacommunities. These results exemplify that the use of the three-component metacommunity framework may be more useful than the classical four perspective paradigm in analysing metacommunities. Our findings also provide support for conservation strategies based on the preservation of heterogeneous habitats in a metacommunity context.

59

60 **Key-words**

61 beta regression, comparative analysis, insects, metacommunity theory, single-species
62 distribution models, stream macroinvertebrates, subarctic streams.

63
64 **Introduction**

65 Metacommunity theory predicts the assembly of ecological communities according to
66 different perspectives. Originally, this idea was illustrated by Leibold et al. (2004) in the form
67 of four metacommunity perspectives: (1) patch dynamics, which is based on a resource
68 competition-colonisation trade-off among species, thus taking into account species' dispersal
69 potential (Hanski, 1994); (2) species-sorting along environmental gradients, which relies on
70 differences in environmental tolerance among species (Leibold, 1995); (3) mass effects or
71 source-sink dynamics, whereby species may survive in poor-quality habitats owing to
72 constant immigration from the source populations in high quality habitats (Pulliam 1988);
73 and (4) the neutral theory, where demographic stochasticity solely explains assembly patterns
74 (Hubbell, 2001). Deciphering which of these perspectives is more suitable in the context of
75 metacommunity analysis seems difficult and may well depend on the context of analysis (e.g.
76 spatial extent, biogeographic region, ecosystem type and more; Heino et al., 2015).

77 Nevertheless, the examples of metacommunity perspectives depicted in Leibold et al.
78 (2004) are not mutually exclusive, and represent a fraction of possibilities which can be
79 expanded with the inclusion of species dispersal rates, connectivity, species interactions,
80 disturbance, priority effects, rapid local adaptation, meta-ecosystem dynamics and more
81 (Brown, Sokol, Skelton, & Tornwall, 2017; Logue, Mouquet, Peter, & Hillebrand, 2011). The
82 more recent proposal by Logue et al. (2011) claims that the metacommunity concept is better
83 generalised by three major exclusive components, which decompose the metacommunity
84 framework into (1) environmental heterogeneity, whereby habitat patches differ in
85 environmental attributes; (2) species equivalence, in terms of niche characteristics; and (3)
86 dispersal, referred to as the rate of dispersal among patches. Here, we aim at evaluating
87 species distributions in a subarctic stream insect metacommunity under these two different
frameworks (i.e., Leibold et al., 2004 versus Logue et al., 2011), specifically so as to evaluate
which of the two is more adequate for the interpretation of our observations.

Species distribution models have previously been used to predict community-level
properties such as biodiversity (Ferrier & Guisan, 2006). Their accuracy in predicting
community-level properties appears to be higher than that of community assembly models,
although at a high cost in terms of model complexity (Bonthoux, Baselga, & Balent, 2013,
Chapman & Purse, 2011). The accuracy of single-species distribution modelling, however,
may also be advantageous to test ecological theories about community assembly mechanisms.
This is because accurately modelling the distribution of single species, one at a time, provides
the opportunity to proceed with a subsequent comparative analysis across species. Using a
comparative analysis, the variation in model performance can be related, for example, to
species traits and potential phylogenetic constraints.

Stream insect species, in particular, are highly suitable to decipher community
assembly processes through the comparative analysis of single-species distribution models
(Heino & de Mendoza, 2016). This is because of the high variability among species in
tolerance of environmental conditions, as well as resource exploitation, dispersal capability,
and habit traits (Merritt & Cummins, 1996; Tachet, Richoux, Bournaud, & Usseglio-Polatera,
2010; Schmidt-Kloiber & Hering, 2015; Serra, Cobo, Graça, Dolédec, & Feio, 2016). This
variability is valuable in evaluating which community assembly mechanism dominates in
each particular context of analysis. Basically, such an analysis might shed light into the
relevance of environmental variables, spatial variables, and dispersal capability of species on
model performance. Subsequently, this information can be used as an indicator of the
preponderance of one community assembly mechanism over another (Figure 1). For example,
if many species show similar spatial patterns, and if these species share the same dispersal
potential, we can presume that the ability to disperse may be underlying the observed general
pattern for these species. This would give us hints about the adequacy to consider one
particular metacommunity theory perspective over the others. Within the classical four-
perspective framework (Leibold et al., 2004), patch dynamics would likely be suitable in this
case, as this perspective relies on the different capability of species to both disperse and
exploit resources. Within the metacommunity framework based on three exclusive
components (Logue et al., 2011), dispersal would be main driver in this case. Moreover,
stream insects are also a diverse group of species, which belong to different insect orders and
vary widely in physiological and morphological adaptations (Merritt & Cummins, 1996).
Thus, modelling the distribution of single stream insect species and subsequently proceeding
with a comparative analysis across species is also a suitable indirect practice to explore
possible evolutionary constraints on community assembly processes.

In this study, we analysed the distribution of common stream insect species in the
metacommunity of a subarctic drainage basin. Species differ widely in their dispersal
capability (e.g. passive or active dispersers, aquatic or aerial adults) and tolerance of
environmental conditions such as temperature, water flow, or habitat characteristics
(Grönroos et al., 2013; Heino, 2005; Heino & Grönroos, 2014). We used environmental,
climatic and spatial variables as predictors of the distributions of single stream insect species.
Our aim was to elucidate, first, whether or not environmental and spatial factors are relevant
for explaining the distribution of stream insect species; and second, whether or not the
obtained models can be related to the different dispersal capability, site occupancy (i.e. a
gradient of rarity-commonness), and biological and taxonomic traits, of stream insect species.
Both considerations were used to evaluate which of the two different metacommunity
frameworks, either the one based on four non-exclusive perspectives (Leibold et al., 2004) or
the one based on three exclusive axes (Logue et al., 2011), is more adequate to interpret our
observations of single species distributions in stream networks (Figure 1).
141 The field work for this study (Fig. S1) was conducted in the Tenojoki drainage basin (main
142 stem length: 361 km, basin area: 16377 km², altitude of sites: from 19 to 285 m a.s.l.) in
143 northernmost Finland (70°N, 27°E). This subarctic drainage basin is close to a natural state,
144 since it is characterised by very small human populations and subsequent little impact from
145 human development. A typical feature of the area are short cool summers and long cold
146 winters (from early November to end of May). The mean annual temperature is about -2°C in
147 the continental areas of the drainage basin, and close to 0°C near the Arctic Ocean (Dankers
148 & Christensen, 2005). Annual precipitation ranges from 310 mm to 410 mm depending on
149 the location in the drainage basin (Mansikkaniemi, 1970). Most of the rainfall and snowmelt
150 enters streams and rivers, as evaporation is generally of minor importance. Vegetation is
151 dominated by mountain birch (Betula pubescens ssp. czerepanowii) woodlands at low altitude
152 and barren tundra at higher altitude, but also peatlands, heathlands and riparian meadows
153 occur commonly. Coniferous pine (Pinus sylvestris) woodlands occur only in scattered
154 locations, mostly in the southern parts of the drainage basin. Wadeable streams and rivers
155 (i.e. channel width < 25 m, water depth < 50 cm) in the area are close to a pristine state,
156 providing excellent possibilities for examining species distributions in natural environmental
157 conditions. We sampled altogether 55 tributary streams for this study (for details, see Kärnä
158 et al., 2015). All these 1st to 5th order tributaries drain into the mainstem of the River
159 Tenojoki or the River Utsjoki, and no site is located in the two mainstem rivers (Fig. S1).
We took a 3-minute kick-net sample (net mesh size: 0.3 mm) at each study site (Kärnä et al., 2015) at the same time with the environmental measurements in early and middle of June 2012 (see below). The sample for each site consisted of six 30-s subsamples that were divided between main habitats at a riffle site (ca. 50 m²) based on visual inspections of variation in depth, flow, moss cover and particle size. The six subsamples were pooled in the field to obtain a composite sample. Such a sampling method has been shown to be effective in northern streams, allowing to detect patterns in community structure (Heino, Ilmonen, & Paasivirta, 2014) and distributions of single species (Heino & de Mendoza, 2016). Samples were immediately preserved in ethanol in the field and were taken to the laboratory for further processing and identification. Animals were separated from detritus and moss fragments and identified to the lowest possible taxonomic level, mostly species (Kärnä et al., 2015).

## Species considered and species traits

We detected 107 insect taxa, of which 87 could be taxonomically determined to species or species group (Kärnä et al., 2015). Insects determined to genus level were discarded as they were considered too likely to include a few species, which is inappropriate to model single-species distributions. Then, we focused on 48 species that occurred at more than 10% of the 55 study sites, that is, that occurred in at least six sites. This is because modelling the distribution of species present in less than six sites is likely to produce spurious results and therefore the analysis of these species was considered unreliable (e.g. Pearce & Ferrier, 2000). In practice, we could model the occupancies of only 47 species because the mayfly *Baetis rhodani* occurred at all sites, so we could not use this species to model presence/absence. The 47 stream insect species considered in this study are listed in Table S1.
186

Nomenclature generally follows de Jong et al. (2014) and more specific references for the
Simuliidae (Adler & Crosskey, 2016 ; Ilmonen, 2014).

187

188

Body size class, dispersal potential, functional feeding groups and habit trait groups
were considered as species traits (Table S2). Functional feeding groups refer to exploitation
of different resources, while habit traits define modes of locomotion and attachment to
substrate (Merritt & Cummins, 1996). Body size classes and female dispersal potential
followed a previous study (Heino & de Mendoza, 2016), with additional information from
Tachet et al. (2010), Schmidt-Kloiber & Hering (2015) and Serra et al. (2016). Female
dispersal potential was characterised as being “low” or “high”. In general, all species of the
Simuliidae were considered to have high dispersal potential, owing to the fact that their
females feed as flying adults, in most cases searching for blood of vertebrates, and hence
were assumed here to generally persist much longer as active flyers than the rest of species.

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

In this regard, Baldwin, West, and Gomery (1975) often found their marked Simuliidae
females several kilometers away from their natal streams. Owing to their small size, the
Simuliidae may also be distributed long distances passively by wind (Crosskey, 1990). All
other species were considered as weak dispersers except for the caddisflies *Plectrocnemia
conspersa* and *Potamophylax cingulatus*, according to the information available for these taxa
from previous studies (Gíslason, Hannesdóttir, Munoz, & Pálsson, 2015; Hoffsten, 2004;
Müller-Peddinghaus, 2011; Müller-Peddinghaus & Hering, 2013; Schmidt-Kloiber & Hering,
2015). Although such information about dispersal abilities of stream insects is rather simple,
there is currently no better information available (Schmidt-Kloiber & Hering, 2015; Serra et
al., 2016; Tachet et al., 2010). Functional feeding and habit trait groups generally follow
Merritt & Cummins (1996).

*Local environmental, climatic, and spatial variables*
The 55 streams were surveyed during the early northern summer, between early and middle of June in 2012. We measured a set of local (i.e. proximal) environmental variables that have been found important for stream insects in northern drainage basins in previous studies (Heino et al., 2014; Kärnä et al., 2015). These comprised physical habitat and water physico-chemical variables. For physical habitat variables, we measured current velocity (m/s) and depth (cm) at 30 random spots in a riffle site. We also measured mean width of the riffle site based on five cross-channel measurements, evenly spaced across the surveyed riffle site.
Bank height and bank slope were measured at the same locations with stream width measurements. Bank height was measured as the height of the lower stream bank, i.e. the height from the water level to the edge of terrestrial vegetation. Bank slope was measured (perpendicular to the stream) as a stream bank rise (cm) over 2 m starting from the edge of terrestrial vegetation. Moss cover (%) and particle size classes (%) were visually estimated at 10 squares (1 m²) at random locations in a riffle site. We used a modified Wentworth's (1922) scale of particle size classes: sand (0.25–2 mm), gravel (2–16 mm), pebble (16–64 mm), cobble (64–256 mm) and boulder (256–1,024 mm). Based on the visual estimates for each square, we calculated mean values for each particle size class and moss cover at a site and used these mean values in species distribution modelling. We also visually estimated shading (%) by riparian vegetation and proportion of riparian deciduous trees (%). For physico-chemical properties, we measured pH, conductivity and water temperature at each site in the field using a YSI device model 556 MPS (YSI Inc., Ohio, USA) and took additional water samples during the field campaign for further analysis. Water samples were frozen at the end of the day at the Kevo Field Station situated in the northern part of the study area, and were later analysed for total nitrogen, colour, iron and manganese in the laboratory of the Finnish Environment Institute in Oulu following Finnish national standards (National Board of Waters, 1981).
236

We also included three climatic variables, including annual air temperature sum above 5°C (growing degree days), mean annual air temperature and mean July air temperature for the period 1981–2010. These variables were calculated in ArcMap 10.2 for each site from a gridded (1 x 1 km) climate data provided by the Finnish Meteorological Institute (Pirinen et al., 2012). The gridded climate data were produced using meteorological station observations and Kriging interpolation (e.g. Aalto, Pirinen, Heikkinen, & Venäläinen, 2013). The selected climatic variables are likely to be important for the distributions of insects in this subarctic area, where temperature is closely associated with insect life cycles (Danks, 2007).

237

238

239

240

241

242

243

244

Spatial variables were distance-based Moran's Eigenvector Maps (db-MEM) based on geographical distances among sites (Dray, Legendre & Peres-Neto, 2006). These spatial db-MEM variables were obtained with the function "PCNM" of the R package "PCNM" (Legendre, Borcard, Blanchet, & Dray, 2013; R Core Team, 2013). We used the largest distance in the minimum spanning tree, keeping all sites connected, as the truncation threshold. Spatial db-MEM variables represent structures of autocorrelation at all spatial scales. Only those spatial db-MEM variables showing significant positive autocorrelation were included in subsequent modelling (Borcard, Gillet & Legendre, 2011), resulting in 13 spatial variables (Figure 2). Based on eigenvalues and bubble plot maps, the spatial variables can be divided into those ranging from large-scale spatial structures (e.g. V1, V2) and those showing very small scale spatial patterns (e.g. V12, V13).

245

246

247

248

249

250

251

252

253

254

255

Prior to modelling species distribution, we eliminated strongly correlated (i.e. Pearson $r$ > .7) predictor variables from the sets of local environmental and climatic variables (see Dormann et al., 2013). Hence, we removed one variable (i.e. annual temperature sum) from the climatic variables and four variables (i.e. water iron, colour, conductivity and boulders) from the stream environmental variables. The spatial variables were already not mutually correlated (Borcard et al., 2011).

256

257

258

259

260
261

262 Modelling species distributions

263 The distribution (i.e. presence/absence) of each species was modelled using binomial
264 generalised linear models (i.e. binomial GLMs with logit link function), using separately
265 local environmental, climatic and spatial variables, with the R package "Rcmdr" (Fox, 2005).
266 The deviance explained for each species was thus obtained for each binomial GLM with each
267 of these three different subsets of variables (Figure 2). The variables selected for each
268 species' model were based on forward selection and Bayesian Information Criterion (BIC),
269 separately for each variable group (i.e. environmental, climate and spatial). BIC values were
270 used because they prevented the selection of too complex models in our case, in contrast to
271 AIC (results not shown), which is often the case under large sample sizes (Burnham &
272 Anderson, 2004). Moreover, the target model under BIC selection does not depend on sample
273 size, in contrast to AIC (Burnham & Anderson, 2004). Therefore, AIC may be problematic in
274 our case as we aim at comparing model performance between species, which may differ in
275 the number of presences and absences. Also, deviating observations were removed from
276 some species' models if they had Cook's distance values > 1 and hence affected profoundly a
277 few models (Cook, 1977). For environmental variables, we registered whether the effect was
278 positive or negative on species distributions. We then used the selected variables of these
279 three subsets (i.e. local environmental, climatic and spatial) to perform variation (deviance)
280 partitioning by subtraction, similarly as performed in multivariate contexts (Legendre &
281 Legendre, 2012). Specifically, the deviance accounted for subset A, subset B, and subset A
282 and B together, was computed, so as to obtain the different fractions of variation solely
283 explained by each subset (i.e. unshared with other subsets). We eventually obtained adjusted
284 D² values (Guisan & Zimmermann, 2000; Legendre & Legendre, 2012) which could be
285 attributed to pure local environmental (E), climatic (C) or spatial effects (S), as well as to
286 total effects combining the three subsets of pure effects and their joint effects (E+C+S effects; Figure 2). Modelling methods other than GLMs could have been possible, yet species 287 probably show linear responses to the environmental predictors due to the fact that they are 288 on the edge of their geographical and ecological distributions, making GLMs adequate. 289 Adding quadratic terms to binomial models is unlikely to change results substantially in these 290 situations (e.g. Pulido, Riera, Ballesteros, Chappius, & Gacia, 2015), and increase the 291 difficulty of interpretation of the results. Also, deviance partitioning is easy to accomplish 292 when this is based on GLMs. 293

294

295 *Comparative analysis across species*

296 We performed a comparative analysis across species using beta regression (Ferrari & Cribari-297 Neto, 2004), where the adjusted *D*² values obtained with previous binomial GLMs were used 298 as the dependent variable to be explained by site occupancy, taxonomic vectors or species 299 trait vectors (Figure 2). These vectors were obtained separately from Principal Coordinate 300 Analysis (PCO). Using the taxonomic relatedness of species, a taxonomic relatedness matrix 301 was built using the function "taxa2dist" in the R package "vegan" (Oksanen et al., 2013), and 302 taxonomic vectors were handled as continuous PCO vectors with the function "pco" in the R 303 package "ecodist" (Goslee & Urban, 2007). The first four taxonomic eigenvectors were 304 selected as these had much higher eigenvalues than the rest (Fig. S2). Similarly, species trait 305 vectors were also computed using body size class, dispersal potential, functional feeding 306 groups, and habit trait groups (Table S2). Species traits were considered as regular factors, 307 except body size class which was considered as an ordered factor, to obtain a distance matrix 308 based on Gower's metric with the function "daisy" of the R package "cluster" (Maechler, 309 Rousseeuw, Struyf, Hubert, & Hornik, 2013), and eventually trait PCO vectors with the 310 function "cmdscale". The four trait eigenvectors obtained were considered for further
statistical analyses. The variation in adjusted *D*² values across species that could be attributed
to pure E, pure C, pure S, or E+C+S effects was fitted on site occupancy, the four taxonomic
and four species trait vectors selected, using beta regression with the function "betareg" of
the R package "betareg" (Cribari-Neto & Zeileis, 2010). Beta regression is adequate when the
response variable (in this case, the adjusted *D*² values) is constrained between 0 and 1.

We compared the explained variation by pure E, C, and S effects with a Kruskal-
Wallis test, with additional Mann-Whitney tests for subsequent pair-wise comparisons
between groups. Non-parametric tests were chosen since adjusted *D*² values data departed
from normality following the Shapiro-Wilk test (Zar, 1984). We also analysed the univariate
relationships between site occupancy, body size, dispersal potential, broad taxonomic insect
groups, functional feeding groups, habit trait groups, and taxonomic and trait vectors.
Depending on the continuous (e.g. site occupancy) or categorical (e.g. habit trait group)
nature of the variables involved, we followed Kruskal-Wallis tests, Mann-Whitney tests,
Fisher's exact test or Spearman correlations, as these variables were generally not normally
distributed (Zar, 1984).

326

## 327 Results

328 *Single species models*

329 Local environmental and spatial effects accounted for a higher variation in species
distributions (16.1% and 12.6% in average, respectively) than did climatic effects (5.4%; *p* <
.001, Kruskal-Wallis test), whereas the average deviance explained did not differ
330 significantly between local environmental and spatial effects (*p* = .125, Mann-Whitney test)
331 (Table S3). The local environmental factors most frequently selected in explaining species
distributions were water temperature, shading, and to a lesser extent, stream width, cobbles
332 and moss (Figure 3). The spatial variables most often selected were better represented by
333
large-scale spatial variables within the Tenojoki drainage basin (e.g. V1, V2) than by small-
spatial scale variables (e.g. V12, V13), as also shown in Figure 3. Amongst the climate
variables, mean annual temperature was significant in explaining the distribution of 32
species, and July air temperature of 20 species (not shown in Fig. 3).

The adjusted deviance explained by binomial GLMs was highly variable across species and difficult to relate to particular taxonomic groups (Table S3). For example, local environmental effects were particularly relevant for the stonefly *Siphonoperla burmeisteri* (i.e. accounting for 66.7% of adjusted *D*² values), the mayfly *Heptagenia dalecarlica* (50.2%), and the blackfly *Prosimulium hirtipes* (37.7%), whereas spatial effects were most relevant for the caddisfly *Rhyacophila nubila* (41.2%), the stonefly *Brachyptera risi* (31.5%), and the chironomid midge *Cardiocladius capucinus* (28.6%). Climate effects were also highly variable. They were generally low (see above), and accounted for more than 20% of adjusted *D*² values in only three cases: the stoneflies *Diura nanseni* and *Siphonoperla burmeisteri* (34.4% and 20.1%, respectively), and the chironomid midge *Orthocladius rivicola* (28.9%). Combining all effects, binomial GLMs explained on average 37.8% of the null deviance (Table S3).

Comparative analysis across species models

The highly variable species-local environment and species-climate relationships in binomial
GLMs were not accounted for by site occupancy, or by taxonomic and trait vectors, in the
beta regression analysis (Table 1). The deviance explained by spatial variables was, however,
significantly (i.e. *p* < .05) accounted for by site occupancy (Table 1). The influence of TAX-
PCO4 and TRA-PCO2 on the adjusted *D*² values predicted by spatial effects in binomial
GLMs was significant as well. Also, the influence of TAX-PCO3 was marginally significant
(i.e. *p* < .10), remaining like this in the binomial GLMs based on all variables combined
361 (Table 1). However, when repeating the beta regression analysis by using only the significant
362 variables selected (i.e. site occupancy, TAX-PCO3, TAX-PCO4, and TRA-PCO2), only site
363 occupancy was statistically significant ($p = .017$), but not TAX-PCO3, TAX-PCO4 or TRA-
364 PCO2 ($p = .943, p = .175$, and $p = .449$, respectively, results not shown in Table 1).

365 Analysing through beta regression the univariate relationship of these variables with the
366 adjusted $D^2$ values of binomial GLMs based on spatial effects produced a similar result (site
367 occupancy, $p = .036$, Fig. S3; TAX-PCO3, TAX-PCO4, and TRA-PCO2, $p = .760, p = .660$,
368 and $p = .524$, respectively, results not shown). This univariate relationship between the
369 adjusted $D^2$ values and site occupancy was not observed when the adjusted $D^2$ values of
370 binomial GLMs were referred to environmental or climate effects (Fig. S3). No statistical
371 significance was observed either for univariate relationships between separate species traits
372 and the adjusted $D^2$ values in binomial GLMs, with the sole exception of body size (Fig. S3).

373 The TAX-PCO3 vector showed the highest species scores for blackflies (Simuliidae)
374 and the lowest for mayflies (Ephemeroptera), and was strongly correlated ($p < .001$) to
375 dispersal potential (Table S5, Figure 4a). In contrast to this taxonomic vector, TAX-PCO4
376 showed the highest species scores for both blackflies and mayflies (Figure 4b), and was
377 strongly correlated to site occupancy ($p = .007$, Table S5). Finally, TRA-PCO2 reflects the
378 influence of functional feeding groups and body size on model performance (Figure 4c), as
379 indicated by the strong correlation of both variables (i.e. $p < .001$) with this trait vector (Table
380 S5).

381

382 **Discussion**

383 *Single species models*

384 Our results indicated that single species distributions of stream insects are highly variable in
terms of predictability, as well as the significant environmental and spatial predictors
underlying such distributions. There was no evident association between model accuracy and particular taxonomic groups (Table S3). Nevertheless, a few generalisations can be highlighted with regard to the results obtained. For example, water temperature and shading, and to a lesser extent, stream width, cobbles and moss, were more relevant as environmental predictors of species distributions than stream flow or water chemistry variables (Table S3, Figure 3). This is in line with the well-known influence of temperature and resource availability on insect life cycles at high latitudes (Danks, 2007) and indicates the influence of species sorting processes along these environmental gradients. Resource availability is represented in our case by shading, which indicates the proximity of terrestrial vegetation and hence is a surrogate of availability of allochthonous resources from terrestrial origin for aquatic insect larvae. This typically corresponds with a situation of a low-order stream which, as in our case, is influenced strongly by terrestrial material from riparian vegetation which is then taken as food resource by shredders, hence promoting their dominance (Vannote, Minshall, Cummins, Sedell, & Cushing, 1980). Shading may also be inversely related to primary productivity, but in this study, we found that the relationship of species distribution with shading was always positive (Figure 3), suggesting that rather than biofilm production, it is the external input of terrestrial material from riparian birch tree abundance what is likely driving species distributions. In our case, shading was selected as a significant variable in binomial models for some predators (*Isoperla difformis* and *Plectrocnemia conspersa*) and shredders (*Leuctra* spp.), for some collector-gatherers (*Corynoneura lobata*-type, *Eukiefferiella devonica*-group, *Orthocladius rhyacobius*-group and *Tvetenia discoloripes*), and for some collector-filterers (*Philopotamus montanus* and *Prosimulium hirtipes*) (Table S3). These latter groups perhaps benefit indirectly from the increase in potential resources that the variable “shading” represents for shredders, for example, through the enhancement of
410

nutrient re-cycling by shredding coarse plant litter (Covich, Palmer, & Crowl, 1999; Wallace

411

& Webster, 1996).

412

Spatial variables were also relevant for the distributions of some species. Specifically,

413

large-scale spatial variables were more important than small-scale variables in explaining

414

species distributions in our study (Table S3, Figure 3). At a larger spatial extent (ca. 500 km

415

latitudinal gradient), previous findings indicate a stronger relevance of environmental factors,

416

compared to spatial restrictions, on single-species distributions (Heino & de Mendoza, 2016).

417

This is perhaps not surprising because increasing the spatial extent may have a strong positive

418

effect on the relevance of niche processes through larger environmental gradients (Chase,

419

2014). However, increasing the spatial extent may also preclude species to reach

420

environmentally suitable locations owing to dispersal limitation, and thus the relative

421

contribution of both environmental and spatial constraints on species distributions does not

422

always vary predictably with spatial scale (Alahuhta & Heino, 2013).

423

424

*Comparative analysis across species*

425

Comparative analysis across the species models showed a clear relationship between model

426

performance and site occupancy. Specifically, the binomial GLMs that we built upon spatial

427

variables could be related to site occupancy, and to a lesser extent, to taxonomic and trait

428

vectors, whereas none of these variables was significantly related to model performance

429

when models were based on local environmental or climate variables (Table 1). At first

430

glance, our results also suggested both a slight influence of female dispersal potential (related

431

to the taxonomic vector TAX-PCO3), and a potential influence of functional feeding groups

432

and body size (related to the trait vector TRA-PCO2), on the performance of models based on

433

spatial variables. The taxonomic vector TAX-PCO3 perhaps relates to female dispersal

434

potential, as species scores along this vector were much higher for the blackflies than for the
rest of species, and lowest for the mayflies (Figure 4). Blackflies are possibly the best active dispersers among all the insects we considered, because females feed as flying adults and in most species they must actively search for blood meals, often several kilometers away from their natal streams (Baldwin et al., 1975). However, adult mayflies, do not feed and often have extremely short life spans (Brittain, 1990). Therefore, it seems reasonable to assume that blackflies may actively disperse better than mayflies. Site occupancy and dispersal potential were not correlated (Table S5), and both taxa were the ones with highest number of sites occupied (Fig. S4). In contrast, mayflies differed in site occupancy from non-biting midges (Chironomidae) (Fig. S4), despite species in both groups can be considered weak active dispersers, as chironomid adults are also short-lived and generally weak active fliers (Armitage, 1995). On the other hand, the trait vector TRA-PCO2 suggests an influence of feeding behaviour and body size (Figure 4, Table S5) on model performance. This is because the exploitation of food resource from terrestrial origin (i.e. shredders) would facilitate the development of more complex trophic food webs with the inclusion of predators (Figure 4). This would also contribute to the positive association of body size to TRA-PCO2 (Fig. 4), as the largest insects we found are either predators or shredders (Table S2).

Nevertheless, it is important to note that taxonomic and trait vectors had a comparatively much weaker effect on predictability by spatial variables than that of site occupancy. In fact, not only did site occupancy attain a higher statistical significance (Table 1), but it could also be partly related to the capability of the taxonomic vector TAX-PCO4 to account for the adjusted *D*² values of binomial GLMs because these two predictor variables were significantly correlated (Table S5). Moreover, when repeating the beta regression analysis by using only the significant variables selected (i.e. site occupancy, TAX-PCO3, TAX-PCO4, and TRA-PCO2), only site occupancy was statistically significant, indicating that the influence of taxonomic and trait vectors on model performance is rather weak.
Analysing through beta regression the univariate relationship of these variables with the adjusted *D*² values of binomial GLMs based on spatial effects again resulted in site occupancy as the only significant variable (see Results above). Therefore, we must conclude that any potential effect of taxonomic and trait vectors on model performance, including the effect of female dispersal potential and body size, and that of functional feeding groups, must be considered with caution: their statistical significance only appears after controlling for site occupancy and the other variables considered in the full model of beta regression. In this regard, the fact that *Baetis rhodani* is a widespread mayfly, which could not be modelled because it was present at all sites, also gives support to the idea that dispersal abilities are not so important in structuring invertebrate assemblages in high-latitude drainage basins. This is because it demonstrates that mayfly species can be widespread, despite being rather weak active dispersers. We also acknowledge that the rarest species (i.e. present in less than six sites) were not modelled because models based on such small number of presences were considered unreliable (e.g. Pierce & Ferrier, 2000). However, excluding these species does not undermine the conclusion that the distributions of most common species are better accounted for by models based on spatial variables than that of not-so-common species. In fact, we effectively modelled 47 out of the 86 taxa available at the species (most cases) or species-group (few cases) taxonomic resolution, comprising 55% of cases, which is a representative subset of species in the entire metacommunity.

*Approaching the suitability of metacommunity analysis frameworks*

With the information above about single-species distribution models and subsequent comparative analysis across species, it is possible to proceed with the evaluation of the suitability of the two different frameworks of metacommunity analysis (Figure 1) considered here: (1) the classical approach exemplified by the four different non-exclusive perspectives
485

described by Leibold et al. (2004) or (2) the three exclusive components as proposed by
Logue et al. (2011).

486

487

Among the four different metacommunity perspectives of the Leibold et al. (2004) framework, neutral theory and patch dynamics do not rely on the effect of environmental variables, in contrast to species sorting and source-sink dynamics, the latter of which also incorporates a strong influence of spatial effects (Figure 1a). In our study, single-species models often relied on the effect of environmental variables, particularly temperature and shading, while being also dependent on large-scale spatial variables (Figure 3). As environmental and spatial factors are both relevant for the distribution of species, this result suggests that either species sorting along spatially structured environmental gradients, or source-sink dynamics between populations of high-quality and low-quality habitats, are both likely as important processes driving metacommunities. Then, the comparative analysis across species showed that site occupancy is responsible for the observed differences in the relevance of spatial variables on species distributions (Table 1). This suggests that common species would be better able than rare species to maintain populations in low-quality habitats through constant immigration, favouring the source-sink dynamics perspective over species sorting.

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

Although species-sorting processes cannot be completely discarded because of the demonstrated influence of environmental variables in many cases, deviance partitioning suggests that the pure effects of environmental and spatial factors on species distributions are stronger than their joint effects (Table S3). Also, the effect of spatial variables was better explained than that of environmental factors by our explanatory variables, particularly site occupancy, in the comparative analysis. These results slightly undermine the idea of species-sorting across spatially structured environmental gradients as the most important process shaping metacommunities. In any case, the neutral theory, which relies entirely on spatial
dynamics, is unlikely. As the dispersal potential of species has a rather weak effect on model accuracy, patch dynamics can be discarded as well as a suitable perspective of metacommunity analysis in our case. It should be acknowledged, however, that the difficulty to explain model performance with dispersal ability can also be a consequence of the coarseness of the dispersal measures currently available for freshwater invertebrates (Schmidt-Kloiber & Hering, 2015; Serra et al., 2016; Tachet et al., 2010). Moreover, the different metacommunity paradigms from Leibold et al. (2004) may always act simultaneously to a certain extent along a continuum (Figure 1a) rather than being distinct and mutually exclusive options (Brown, Sokol, et al., 2017; Gravel, Canham, Beaudet, & Messier, 2006; Logue et al., 2011).

Spatial autocorrelation may appear not only as a consequence of mass effects or species sorting along spatially structured environmental gradients when the spatial scale is not very large, but also as a consequence of dispersal limitation at very large spatial scales (Heino et al., 2015). Nevertheless, some insect species found in this study exemplify well the potential importance of the source-sink dynamics for metacommunities in subarctic streams, independently of their dispersal capability. For example, six blackfly species were examined (Table S1), of which five were present in more than 50% of sites, three of them in 75% of sites or more (Table S2). Thus, blackfly species in subarctic streams have successfully spread widely, which is advantageous to maintain metapopulations through source-sink dynamics. On the other hand, the mayflies are as widespread as the blackflies (Fig. S4), but far less capable of active dispersal. This suggests that the dispersal capability of species does not determine the metapopulation dynamics, whereas site occupancy probably does so. Spatial autocorrelation patterns have been described for the blackflies at small spatial scales, driven by strong effects of inter-specific competition for oviposition sites, and subsequent priority effects at the community level (McCreadie & Adler, 2012). The importance of priority effects
535

for the blackflies reinforces the idea of the relevance of site occupancy for community
dynamics, where rare species are in clear disadvantage for habitat recolonisation.

536

537

Alternative to the framework of Leibold et al. (2004), we can interpret our results
under the framework of Logue et al. (2011), whereby three different and mutually exclusive
components can be used to analyse metacommunities: species equivalence, habitat
heterogeneity and dispersal (Figure 1b). In our case, this alternative framework makes
interpretation of the results much easier. At the very least, we can conclude that species
equivalence is unlikely to play any role in metacommunity dynamics, similarly to discarding
neutral theory under the Leibold *et al.* (2004) framework. Dispersal can also be discarded, yet
again with caution due to the current lack of high resolution dispersal measures for freshwater
invertebrates (Schmidt-Kloiber & Hering, 2015; Serra et al., 2016; Tachet et al., 2010). Thus,
the main difference in the interpretation of the results with this alternative framework is that
we can now be certain about the role of habitat heterogeneity, while under the Leibold et al.
(2004) framework it is more difficult to discern whether species sorting or source-sink
dynamics is the dominant process. Habitat heterogeneity is indeed related to both
mechanisms. In fact, using habitat heterogeneity in space and time as the templet for
ecological strategies (Southwood, 1977) could be the framework of choice in situations
where it is difficult to discern species sorting processes from source-sink dynamics.

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

*Alternative approaches, caveats and conclusions*

555

556

557

558

559

Emergent properties at the community level are difficult to discern from field observational
data alone. In this regard, population genetics can be very useful in order to gain confidence
about the distinction between, e.g., source-sink dynamics and species sorting processes. This
is because population genetic studies could be used to estimate the relative contribution of
immigrants from nearby populations to the genetic variability of the population under study
560 (Bunn & Hughes, 1997; Hughes, Huey, & Schmidt, 2013; Hughes, Schmidt, & Finn, 2009).

561 Genetic analyses would probably provide the opportunity for a more robust interpretation of
562 our results. Genetic studies, however, are difficult to accomplish with stream insects in the
563 field when the idea is to compare many species at a time, and they are far more expensive
564 than the comparative approach of single species distributions we considered here. Therefore,
565 the comparative approach presented here can be used as a first step to explore the relative
566 contribution of environmental and spatial factors on species distributions, without using
567 expensive and time-consuming genetic analyses. In fact, by using the comparative approach
568 we can certainly conclude that the dispersal capability of species and neutral theory play little
569 role in shaping subarctic stream insect metacommunities. Rather, it is habitat heterogeneity,
570 which influences mass effects and/or species sorting processes, that matters. Subsequently,
571 the results of our study strongly recommend the preservation of habitat heterogeneity as the
572 conservation strategy to maintain biodiversity in these ecosystems.

573 Nevertheless, it should be acknowledged that one shortcoming of single-species
574 distribution modelling is that it does not consider the influence of species interactions in
575 structuring ecological communities. Stream ecology has considered that severe environmental
576 conditions may weaken the potential effects of biotic interactions in structuring communities
577 (Peckarsky, 1983). However, more recent findings pose doubts as to whether this is actually
578 true (Cadotte & Tucker, 2017; Thomson et al., 2002). In fact, biotic interactions can
579 reproduce patterns of community structure essentially identical to what it could be expected
580 from environmental filtering alone. This is because environmental changes may affect
581 population growth rates of competing species in opposite ways, and this may cause the
582 exclusion of some species that would otherwise be able to coexist (Cadotte & Tucker, 2017).
583 There exists also evidence indicating that biotic interactions limit the geographical range
584 expansion of species facing environmental changes (Pigot & Tobias, 2013; Sexton et al.,
2009). Overall, this suggests that inter-specific interactions may also play a role in our case, although the abundances of insect larvae in subarctic streams are typically low (see also Heino & Grönroos, 2017) and may thus result in weak density-dependent interactions among species (see also Morin, 2011).

Our study considered tributary streams draining into two linear sub-elements of a larger river network (Fig. S1). However, there exists growing concern about the potential role of the entire dendritic river networks in shaping biodiversity patterns, community structure and species distributions (Altermatt, 2013; Altermatt & Fronhofer, 2017; Brown, Wahl, & Swan, 2017; Jamoneau, Passy, Soininen, Lebuocher, & Tison-Rosebery, 2017; Schmera et al., 2017). For example, the consideration of whole river networks may unveil a more preeminent role for spatial factors in community assembly, undermining the role of environmental filtering. Therefore, studies conducted across whole dendritic networks could be more in line with neutral theory, as shown by Muneepeerakul et al. (2008) for fish communities, yet no environmental variable was truly considered in that study. Although we focused on tributary streams draining into the main river, the consideration of whole dendritic networks may help us to perceive more accurately the real connectivity pathways between isolated patches. This connectivity may have consequences for metacommunity stability with respect to a situation where only a linear component of this network is acting (Fagan, 2002). Also, dispersal along dendritic networks implies more variability in local richness with strong consequences also for community differentiation among patches (Carrara, Altermatt, Rodriguez-Iturbe, & Rinaldo, 2012; Seymour, Fronhofer, & Altermatt, 2015). In any case, dispersal processes in stream networks may depend on the organism group considered (Schmera et al. 2017).

In the case of stream insects, the taxa considered and the taxonomic resolution achieved prior to species-distribution modelling, may also have important consequences on
our perception of the influence of dendritic riverine networks on biodiversity patterns (Kaelin & Altermatt, 2016). Here, some taxa were discarded as it was not possible to determine the species. Provided that the influences of dendritic landscapes and biotic interactions (discussed above) may strongly affect how we understand the reality of community assemblages, it is essential to use the best taxonomic resolution possible to make accurate inferences about the mechanisms truly governing the observed patterns. In fact, the criterion of ‘best taxonomic resolution possible’ used in our modelling endeavours is a fundamental requirement to draw robust conclusions to be applied in biodiversity conservation.

Finally, for biodiversity conservation, it is essential to focus on maintaining habitat heterogeneity because it appears to determine metacommunity organization (Kärnä et al. 2015) and species distributions (Heino & de Mendoza, 2016) in streams at high latitudes. Unless habitat heterogeneity is not considered (along with potentially important effects of dendritic network structure and biotic interactions), conservation plans may fall short and not result in desired outcomes.

## Acknowledgements

We thank Sirkku Lehtinen and Marja Lindholm for help in the field or the laboratory. We also thank two anonymous reviewers and the Guest Editor Florian Altermatt for their very helpful and constructive comments on an earlier version of this manuscript. Kevo Subarctic Research Station provided facilities during the field work. This study is part of the project “Spatial scaling, metacommunity structure and patterns in stream communities” that was supported financially by a grant from the Academy of Finland. Further support was provided by grants (no: 273557, no: 267995 and no: 285040) from the Academy of Finland. The authors declare no conflict of interest.
635 **References**

636 Aalto, J., Pirinen, P., Heikkinen, J., & Venäläinen, A. (2013) Spatial interpolation of monthly
637 climate data for Finland: comparing the performance of kriging and generalized
638 additive models. *Theoretical and Applied Climatology*, **112**, 99–111.

639 Adler, P.H., & Crosskey, R.W. (2016) *World blackflies (Diptera: Simuliidae): A
640 comprehensive revision of the taxonomic and geographical inventory*. Retrieved from
641 http://www.clemson.edu/cafls/biomia/pdfs/blackflyinventory.pdf.

642 Alahuhta, J., & Heino, J. (2013) Spatial extent, regional specificity and metacommunity
643 structuring in lake macrophytes. *Journal of Biogeography*, **40**, 1572–1582.

644 Altermatt, F. (2013) Diversity in riverine metacommunities: a network perspective. *Aquatic
645 Ecology*, **47**, 365–377.

646 Altermatt, F., & Fronhofer, E. A. (2017) Dispersal in dendritic networks: Ecological
647 consequences on the spatial distribution of population densities. *Freshwater Biology,
648 63*, 22–32. https://doi.org/10.1111/fwb.12951

649 Armitage, P. D. (1995) Behaviour and ecology of adults. In P. D. Armitage, P. S. Cranston, &
650 L. C. V. Pinder (Eds.), *The Chironomidae: Biology and Ecology of Non-Biting
651 Midges* (pp. 194–224). London: Chapman & Hall.

652 Baldwin, W. F., West, A. S., & Gomery, J. (1975) Dispersal pattern of black flies (Diptera: Simuliidae) tagged with ³²P. *The Canadian Entomologist*, **107**, 113–118.

654 Bonthoux, S., Baselga, A., & Balent, G. (2013) Assessing community-level and single-species models predictions of species distributions and assemblage composition after 25 years of land cover change. *PLoS ONE*, **8**, e54179.

657 Borcard, D., Gillet, F., & Legendre, P. (2011) *Numerical Ecology with R*. New York:
658 Springer.
659 Brittain, J. E. (1990) Life history strategies in Ephemeroptera and Plecoptera. In I. C.

660 Campbell (Ed.), Mayflies and Stoneflies (pp. 1–12). Dordrecht: Kluwer Academic

661 Publishers.

662 Brown, B. L., Sokol, E. R., Skelton, J., & Tornwall, B. (2017) Making sense of
663 metacommunities: dispelling the mythology of a metacommunity typology.

664 *Oecologia*, **183**, 643–652.

665 Brown, B. L., Wahl, C., & Swan, C. M. (2017). Experimentally disentangling the influence
666 of dispersal and habitat filtering on benthic invertebrate community structure.

667 *Freshwater Biology*, **63**, 48–61. https://doi.org/10.1111/fwb.12995

668 Bunn, S. E., & Hughes, J. M. (1997) Dispersal and recruitment in streams: evidence from
669 genetic studies. *Journal of the North American Benthological Society*, **16**, 338–346.

670 Burnham, K. P., & Anderson, D. R. (2004) Multimodel inference: understanding AIC and
671 BIC in model selection. *Sociological Methods & Research*, **33**, 261–304.

672 Cadotte, M. W., & Tucker, C. M. (2017) Should environmental filtering be abandoned?
673 *Trends in Ecology & Evolution*, **32**, 429–437.

674 Carrara, F., Altermatt, F., Rodriguez-Iturbe, I., & Rinaldo, A. (2012) Dendritic connectivity
675 controls biodiversity patterns in experimental metacommunities. *Proceedings of the
676 National Academy of Sciences of the United States of America*, **109**, 5761–5766.

677 Chapman, D. S., & Purse, B. V. (2011) Community versus single-species distribution models
678 for British plants. *Journal of Biogeography*, **38**, 1524–1535.

679 Chase, J. M. (2014) Spatial scale resolves the niche versus neutral theory debate. *Journal of
680 Vegetation Science*, **25**, 319–322.

681 Cook, R. D. (1977) Detection of influential observation in linear regression. *Technometrics*,
682 **19**, 15–18.
Covich, A. P., Palmer, M. A., & Crowl, T. A. (1999) The role of benthic invertebrate species in freshwater ecosystems: Zoobenthic species influence energy flows and nutrient cycling. *BioScience, **49***, 119-127.

Cribari-Neto, F., & Zeileis, A. (2010) Beta regression in R. *Journal of Statistical Software*, **34**, 1-24.

Crosskey, R. W. (1990) *The Natural History of Blackflies*. John Wiley & Sons, Chichester.

Dankers, R., & Christensen, O.B. (2005) Climate change impacts on snow coverage, evaporation and river discharge in the sub-arctic Tana basin, Northern Fennoscandia. *Climatic Change*, **69**, 367-392.

Danks, H. V. (2007) How aquatic insects live in cold climates. *The Canadian Entomologist*, **139**, 443-471.

de Jong, Y., Verbeek, M., Michelsen, V., Bjørn, P. P., Los, W., Steeman, F., Bailly, N., Basire, C., Chylarecki, P., Stloukal, E., Hagedorn, G., Wetzel, F. T., Glöckler, F., Kroupa, A., Korb, G., Hoffmann, A., Häuser, C., Kohlbecker, A., Müller, A., Güntsch, A., Stoev, P., & Penev, L. (2014) Fauna Europaea: All animal species on the web. *Biodiversity Data Journal*, **2**, e4034.

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., García Marquéz, J. R., Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013) Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. *Ecography*, **36**, 27-46.

Dray, S., Legendre, P., & Peres-Neto P. R. (2006) Spatial modelling: A comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). *Ecological Modelling*, **196**, 483-493.
707 Fagan, W. F. (2002) Connectivity, fragmentation, and extinction risk in dendritic
708 metapopulations. *Ecology*, **83**, 3243–3249.

709 Ferrari, S. L. P., & Cribari-Neto, F. (2004). Beta regression for modelling rates and
710 proportions. *Journal of Applied Statistics*, **31**, 799–815.

711 Ferrier, S., & Guisan, A. (2006) Spatial modelling of biodiversity at the community level.
712 *Journal of Applied Ecology*, **43**, 393–404.

713 Fox, J. (2005) The R Commander: A Basic Statistics Graphical User Interface to R. *Journal
714 of Statistical Software*, **14(9)**, 1–42.

715 Gíslason, G. M., Hannesdóttir, E. R., Munoz, S. S., & Pálsson, S. (2015) Origin and dispersal
716 of *Potamophylax cingulatus* (Trichoptera: Limnophilidae) in Iceland. *Freshwater
717 Biology*, **60**, 387–394.

718 Goslee, S.C. & Urban, D.L. (2007) The ecodist package for dissimilarity-based analysis of
719 ecological data. *Journal of Statistical Software*, **22(7)**, 1–19.

720 Gravel, D., Canham, C. D., Beaudet, M., & Messier, C. (2006) Reconciling niche and
721 neutrality: The continuum hypothesis. *Ecology Letters*, **9**, 399–409.

722 Grönroos, M., Heino, J., Siqueira, T., Landeiro, V. L., Kotanen, J., & Bini, L. M. (2013)
723 Metacommunity structuring in stream networks: Roles of dispersal mode, distance
724 type, and regional environmental context. *Ecology and Evolution*, **3**, 4473–4487.

725 Guisan, A., & Zimmermann, N.E. (2000) Predictive habitat distribution models in ecology.
726 *Ecological Modelling*, **135**, 147–186.

727 Hanski, I. (1994) A practical model of metapopulation dynamics. *Journal of Animal Ecology*,
728 **63**, 151–162.

729 Heino, J. (2005) Positive relationship between regional distribution and local abundance in
730 stream insects: a consequence of niche breadth or niche position? *Ecography*, **28**,
731 345–354.
732 Heino, J., & de Mendoza, G. (2016) Predictability of stream insect distributions is dependent
733 on niche position, but not on biological traits or taxonomic relatedness of species.
734 *Ecography*, **39**, 1216–1226.

735 Heino, J., & Grönroos, M. (2014) Untangling the relationships among regional occupancy,
736 species traits and niche characteristics in stream invertebrates. *Ecology and Evolution,*
737 **4**, 1931–1942.

738 Heino, J., & Grönroos, M. (2017) Exploring species and site contributions to beta diversity in
739 stream insect assemblages. *Oecologia*, **183**, 151–160.

740 Heino, J., Ilmonen, J., & Paasivirta, L. (2014) Continuous variation of macroinvertebrate
741 communities along environmental gradients in northern streams. *Boreal
742 Environment Research*, **19**, 21–38.

743 Heino, J., Melo, A. S., Siqueira, T., Soininen, J., Valanko, S., & Bini, L. M. (2015)
744 Metacommunity organisation, spatial extent and dispersal in aquatic systems:
745 Patterns, processes and prospects. *Freshwater Biology*, **60**, 845–869.

746 Hoffsten, P.-O. (2004) Site-occupancy in relation to flight-morphology in caddisflies.
747 *Freshwater Biology*, **49**, 810–817.

748 Hubbell, S. P. (2001) *The Unified Neutral Theory of Biodiversity and Biogeography*.
749 Princeton, NJ: Princeton University Press.

750 Hughes, J. M., Huey, J. A. & Schmidt, D. J. (2013) Is realised connectivity among
751 populations of aquatic fauna predictable from potential connectivity? *Freshwater
752 Biology*, **58**, 951–966.

753 Hughes, J. M., Schmidt, D. J., & Finn, D.S. (2009) Genes in streams: Using DNA to
754 understand the movement of freshwater fauna and their riverine habitat. *BioScience,*
755 **59**, 573–583.
756

Ilmonen, J. (2014) Checklist of the family Simuliidae (Diptera) of Finland. *ZooKeys*, **441**,

757

91–95.

758

Jamoneau, A., Passy, S. I., Soininen, J., Lebuocher, T., & Tison-Rosebery, J. (2017) Beta

759

diversity of diatom species and ecological guilds: Response to environmental and

760

spatial mechanisms along the stream watercourse. *Freshwater Biology*, **63**, 62–73.

761

https://doi.org/10.1111/fwb.12980

762

Kaelin, K., & Altermatt, F. (2016) Landscape-level predictions of diversity in river networks

763

reveal opposing patterns for different groups of macroinvertebrates. *Aquatic*

764

*Ecology*, **50**, 283–295.

765

Kärnä, O.-M., Grönroos, M., Antikainen, H., Hjort, J., Ilmonen, J., Paasivirta, L., & Heino, J.

766

(2015) Inferring the effects of potential dispersal routes on the metacommunity

767

structure of stream insects: as the crow flies, as the fish swims or as the fox runs?

768

*Journal of Animal Ecology*, **84**, 1342–1353.

769

Legendre, P., & Legendre L. (2012) *Numerical Ecology (3rd English edition)*. Elsevier,

770

Amsterdam.

771

Legendre, P., Borcard, D., Blanchet, F.G., & Dray, S. (2013). *PCNM: MEM spatial*

772

*eigenfunction and principal coordinate analyses*. R package version 2.1-2/r109.

773

Retrieved from http://R-Forge.R-project.org/projects/sedar/

774

Leibold, M. A. (1995) The niche concept revisited: mechanistic models and community

775

context. *Ecology*, **76**, 1371–1382.

776

Leibold, M. A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J. M., Hoopes, M. F.,

777

Holt, R. D., Shurin, J. B., Law, R., Tilman, D., Loreau, M., & Gonzalez, A. (2004)

778

The metacommunity concept: a framework for multi-scale community ecology.

779

*Ecology Letters*, **7**, 601–613.
780 Logue, J. B., Mouquet, N., Peter, H., & Hillebrand, H. (2011) Empirical approaches to metacommunities: a review and comparison with theory. *Trends in Ecology and Evolution*, **26**, 482–491.

783 Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., & Hornik, K. (2013) *cluster: Cluster analysis basics and extensions*. R package version 1.14.4. Retrieved from http://CRAN.R-project.org/package=cluster.

786 Mansikkaniemi, H. (1970) Deposits of sorted material in the Inarijoki-Tana river valley in Lapland. *Reports of Kevo Subarctic Research Station*, **6**, 1–63.

788 McCreadie, J. W., & Adler, P. H. (2012) The roles of abiotic factors, dispersal, and species interactions in structuring stream assemblages of black flies (Diptera: Simuliidae). *Aquatic Biosystems*, **8**, 14.

791 Merritt, R. W., & Cummins, K. W. (1996) *An Introduction to the Aquatic Insects of North America*, 3rd ed. Dubuque: Kendall/Hunt Publishing.

793 Morin, P. J. (2011) *Community Ecology*, 2nd ed. Oxford: Wiley-Blackwell.

794 Müller-Peddinghaus, E. (2011) *Flight-morphology of Central European caddisflies (Insecta: Trichoptera) in relation to their ecological preferences*. PhD Thesis, University of Duisburg-Essen, Essen.

797 Müller-Peddinghaus, E., & Hering, D. (2013) The wing morphology of limnephilid caddisflies in relation to their habitat preferences. *Freshwater Biology*, **58**, 1138–1148.

800 Muneepeerakul, R., Bertuzzo, E., Lynch, H. J., Fagan, W. F., Rinaldo, A., & Rodriguez-Iturbe, I. (2008) Neutral metacommunity models predict fish diversity patterns in Mississippi-Missouri basin. *Nature*, **453**, 220–222.
803 National Board of Waters (1981) *Vesihallinnon analyysimenetelmät*, Tiedotus 213. Helsinki:
804 Vesihallitus.

805 Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O'Hara, R. B.,
806 Simpson, G. L., Solymos, P., Stevens, M. H. H., & Wagner, H. (2013) *vegan:*
807 Community Ecology Package. R package version 2.0-9. Retrieved from
808 http://CRAN.R-project.org/package=vegan.

809 Pearce, J., & Ferrier, S. (2000) An evaluation of alternative algorithms for fitting species
distribution models using logistic regression. *Ecological Modelling*, **128**, 127–147.

811 Peckarsky, B. L. (1983) Biotic interactions or abiotic limitations? A model of lotic
812 community structure. In T. D. Fontaine III, & S. M. Bartell (Eds.), *Dynamics of Lotic
813 Ecosystems* (pp. 303–323). Ann Arbor: Ann Arbor Science Publishers.

814 Pigot, A. L., & Tobias, J. A. (2013) Species interactions constrain geographic range
815 expansion over evolutionary time. *Ecology Letters*, **16**, 330–338.

816 Pirinen, P., Simola, H., Aalto, J., Kaukoranta, J. P., Karlsson, P., & Ruuhela, R. (2012)
817 Climatological statistics of Finland 1981–2010. *Finnish Meteorological Institute
818 Reports*, **1**, 1–96.

819 Pulido, C., Riera, J. L., Ballesteros, E., Chappuis, E., & Gacia, E. (2015) Predicting aquatic
820 macrophyte occurrence in soft-water oligotrophic lakes (Pyrenees mountain range).
821 *Journal of Limnology*, **74**, 143–154.

822 Pulliam, H. R. (1988) Sources, sinks, and population regulation. *American Naturalist*, **132**,
823 652–661.

824 R Core Team (2013) *R: A language and environment for statistical computing*. R Foundation
825 for Statistical Computing, Vienna, Austria. Retrieved from http://www.R-project.org/.

826 Schmera, D., Árva, D., Boda, P., Bódis, E., Bolgovics, Á., Borics, G., Csercsa, A., Deák, C.,
827 Krasznai, E. Á., Lukács, B. A., Mauchart, P., Móra, P., Sály, P., Specziár, A., Süveges
828

K., Szivák, I., Takács, P., Tóth, M., Várbíró, G., Vojtkó, A. E., & Erős, T. (2017).

829

Does isolation influence the relative role of environmental and dispersal-related

830

processes in stream networks? An empirical test of the network position hypothesis

831

using multiple taxa. *Freshwater Biology*, **63**, 74–85.

832

https://doi.org/10.1111/fwb.12973

833

Schmidt-Kloiber, A., & Hering, D. (2015) www.freshwaterecology.info – An online tool that

834

unifies, standardises and codifies more than 20,000 European freshwater organisms

835

and their ecological preferences. *Ecological Indicators*, **53**, 271–282.

836

Serra, S. R. Q., Cobo, F., Graça, M. A. S., Dolédec, S., & Feio M. J. (2016) Synthesising the

837

trait information of European Chironomidae (Insecta: Diptera): Towards a new

838

database. *Ecological Indicators*, **61**, 282–292.

839

Sexton, J. P., McIntyre, P. J., Angert, A. L., & Rice, K. J. (2009) Evolution and ecology of

840

species range limits. *Annual Review of Ecology, Evolution, and Systematics*, **40**, 415–

841

436.

842

Seymour, M., Fronhofer, E. A., & Altermatt, F. (2015) Dendritic network structure and

843

dispersal affect temporal dynamics of diversity and species persistence. *Oikos*, **124**,

844

908–916.

845

Southwood, T. R. E. (1977) Habitat, the templet for ecological strategies? *Journal of Animal*

846

*Ecology*, **46**, 337–365.

847

Tachet, H., Richoux, P., Bournaud, M., & Usseglio-Polatera, P. (2010) *Invertébrés d'eau*

848

*douce: systématique, biologie, écologie (nouvelle édition revue et augmentée).*

849

Paris: CNRS.

850

Thomson, J. R., Lake, P. S., & Downes, B. J. (2002) The effect of hydrological disturbance

851

on the impact of a benthic invertebrate predator. *Ecology*, **83**, 628–642.
852 Vannote, R. L., Minshall, G. W., Cummins, K. W., Sedell, J. R., & Cushing, C.E. (1980) The river continuum concept. *Canadian Journal of Fisheries and Aquatic Sciences*, **37**, 130–137.

855 Wallace, J. B., & Webster, J. R. (1996) The role of macroinvertebrates in stream ecosystem function. *Annual Review of Entomology*, **41**, 115–139.

857 Wentworth, C. K. (1922) A scale of grade and class terms for clastic sediments. *Journal of Geology*, **30**, 377–392.

859 Zar, J. H. (1984) *Biostatistical Analysis, Second Edition*. Prentice Hall, Englewood Cliffs.
861 Supporting Information

862 Additional Supporting Information may be found in the online version of this article.

863 Table S1. Insect species considered.

864 Table S2. Species traits considered and site occupancy.

865 Table S3. Results of binomial GLMs.

866 Table S4. Taxonomic and trait vectors from Principal Coordinate Analysis (PCO), with
corresponding scores for each species.

868 Table S5. Statistical significance of the correlations among site occupancy, species traits, trait vectors, and taxonomic vectors.

870 Figure S1. A map of the study area located in the Tenojoki drainage basin.

871 Figure S2. Eigenvalues from taxonomic Principal Coordinate Analysis (PCO) based on
872 taxonomic distances between species.

873 Figure S3. Results of binomial GLMs in relation to site occupancy and species traits.

874 Figure S4. Comparison of site occupancy values between different insect groups.

875
<table><caption>**Table 1.** Results of beta regression showing the effects of site occupancy, biological trait vectors and taxonomic vectors on different fractions of variation (adjusted deviance, Adj. *D*²) explained by binomial GLMs: local environment (E) effects, climate (C) effects, spatial (S) effects, and combined (E+C+S) effects. Significant values (*p* < .05) are shown in boldface; marginally significant values (*p* < .10) in italics.</caption><thead><tr><th></th><th>Estimate</th><th>SE</th><th>z</th><th>p</th><th>Log-likelihood</th><th>Pseudo R<sup>2</sup></th></tr></thead><tbody><tr><td colspan="7"><strong>Adj. D<sup>2</sup> of E effects</strong></td></tr><tr><td>(Intercept)</td><td>-1.6002</td><td>0.2220</td><td>-7.209</td><td><strong>&lt;0.001</strong></td><td>49.06</td><td>0.2500</td></tr><tr><td>Site occupancy</td><td>-0.0042</td><td>0.0091</td><td>-0.469</td><td>0.646</td><td></td><td></td></tr><tr><td>TAX-PCO1</td><td>-378.2072</td><td>323.7056</td><td>-1.168</td><td>0.243</td><td></td><td></td></tr><tr><td>TAX-PCO2</td><td>91.4160</td><td>161.7362</td><td>0.565</td><td>0.572</td><td></td><td></td></tr><tr><td>TAX-PCO3</td><td>225.8850</td><td>218.0393</td><td>1.036</td><td>0.300</td><td></td><td></td></tr><tr><td>TAX-PCO4</td><td>-112.1618</td><td>101.9824</td><td>-1.100</td><td>0.271</td><td></td><td></td></tr><tr><td>TRA-PCO1</td><td>-0.5557</td><td>0.8628</td><td>-0.644</td><td>0.520</td><td></td><td></td></tr><tr><td>TRA-PCO2</td><td>1.3525</td><td>1.0712</td><td>1.263</td><td>0.207</td><td></td><td></td></tr><tr><td>TRA-PCO3</td><td>0.3959</td><td>0.9475</td><td>0.418</td><td>0.676</td><td></td><td></td></tr><tr><td>TRA-PCO4</td><td>-1.3110</td><td>1.3459</td><td>-0.974</td><td>0.330</td><td></td><td></td></tr><tr><td colspan="7"><strong>Adj. D<sup>2</sup> of C effects</strong></td></tr><tr><td>(Intercept)</td><td>-2.9806</td><td>0.2682</td><td>-11.115</td><td><strong>&lt;0.001</strong></td><td>94.66</td><td>0.2309</td></tr><tr><td>Site occupancy</td><td>0.0029</td><td>0.0103</td><td>0.278</td><td>0.781</td><td></td><td></td></tr><tr><td>TAX-PCO1</td><td>-579.9412</td><td>375.5543</td><td>-1.544</td><td>0.123</td><td></td><td></td></tr><tr><td>TAX-PCO2</td><td>83.1028</td><td>185.9703</td><td>0.447</td><td>0.655</td><td></td><td></td></tr><tr><td>TAX-PCO3</td><td>352.6843</td><td>252.8857</td><td>1.395</td><td>0.163</td><td></td><td></td></tr><tr><td>TAX-PCO4</td><td>-1.9501</td><td>120.4650</td><td>-0.016</td><td>0.987</td><td></td><td></td></tr><tr><td>TRA-PCO1</td><td>-1.6362</td><td>1.0045</td><td>-1.629</td><td>0.103</td><td></td><td></td></tr><tr><td>TRA-PCO2</td><td>0.4327</td><td>1.2476</td><td>0.347</td><td>0.729</td><td></td><td></td></tr><tr><td>TRA-PCO3</td><td>0.6454</td><td>1.1129</td><td>0.580</td><td>0.562</td><td></td><td></td></tr><tr><td>TRA-PCO4</td><td>-2.1768</td><td>1.5519</td><td>-1.403</td><td>0.161</td><td></td><td></td></tr><tr><td colspan="7"><strong>Adj. D<sup>2</sup> of S effects</strong></td></tr><tr><td>(Intercept)</td><td>-2.4436</td><td>0.2222</td><td>-10.997</td><td><strong>&lt;0.001</strong></td><td>60.19</td><td>0.2137</td></tr><tr><td>Site occupancy</td><td>0.0208</td><td>0.0083</td><td>2.498</td><td><strong>0.012</strong></td><td></td><td></td></tr><tr><td>TAX-PCO1</td><td>309.6453</td><td>360.2543</td><td>0.860</td><td>0.390</td><td></td><td></td></tr><tr><td>TAX-PCO2</td><td>-4.0452</td><td>162.3764</td><td>-0.025</td><td>0.980</td><td></td><td></td></tr><tr><td>TAX-PCO3</td><td>365.5554</td><td>205.9633</td><td>1.775</td><td>0.076</td><td></td><td></td></tr><tr><td>TAX-PCO4</td><td>217.5296</td><td>102.0912</td><td>2.131</td><td><strong>0.033</strong></td><td></td><td></td></tr><tr><td>TRA-PCO1</td><td>-1.0051</td><td>0.8333</td><td>-1.206</td><td>0.228</td><td></td><td></td></tr></tbody></table>
<table><tr><td>TRA-PCO2</td><td>2.5228</td><td>1.1999</td><td>2.102</td><td><strong>0.036</strong></td></tr><tr><td>TRA-PCO3</td><td>1.4151</td><td>0.9765</td><td>1.449</td><td>0.147</td></tr><tr><td>TRA-PCO4</td><td>-0.7984</td><td>1.3342</td><td>-0.598</td><td>0.550</td></tr></table>

Adj. D² of E+C+S effects

<table><thead><tr><th></th><th>Estimate</th><th>SE</th><th>z</th><th>p</th><th>Log-likelihood</th><th>Pseudo R<sup>2</sup></th></tr></thead><tbody><tr><td>(Intercept)</td><td>-0.7133</td><td>0.2331</td><td>-3.060</td><td><strong>0.002</strong></td><td>17.76</td><td>0.1791</td></tr><tr><td>Site occupancy</td><td>0.0111</td><td>0.0095</td><td>1.164</td><td>0.244</td><td></td><td></td></tr><tr><td>TAX-PCO1</td><td>-140.2227</td><td>352.1883</td><td>-0.398</td><td>0.691</td><td></td><td></td></tr><tr><td>TAX-PCO2</td><td>32.5012</td><td>178.6572</td><td>0.182</td><td>0.856</td><td></td><td></td></tr><tr><td>TAX-PCO3</td><td>407.2860</td><td>228.4095</td><td>1.783</td><td>0.075</td><td></td><td></td></tr><tr><td>TAX-PCO4</td><td>81.0016</td><td>110.2921</td><td>0.734</td><td>0.463</td><td></td><td></td></tr><tr><td>TRA-PCO1</td><td>-1.3052</td><td>0.9142</td><td>-1.428</td><td>0.153</td><td></td><td></td></tr><tr><td>TRA-PCO2</td><td>1.6446</td><td>1.1935</td><td>1.378</td><td>0.168</td><td></td><td></td></tr><tr><td>TRA-PCO3</td><td>0.9742</td><td>1.0367</td><td>0.940</td><td>0.347</td><td></td><td></td></tr><tr><td>TRA-PCO4</td><td>-1.0409</td><td>1.4575</td><td>-0.714</td><td>0.475</td><td></td><td></td></tr></tbody></table>
882 **Figure legends**

883 **Figure 1** Conceptual representation of (a) the four non-exclusive classical approaches in metacommunity studies (Leibold et al. 2004), and (b) the more recent framework of metacommunity analysis based on three exclusive components (Logue et al. 2011); according to the relative relevance of the variables used in this study (axes): spatial variables (x-axis), environmental variables (y-axis) and the different dispersal capability of species (z-axis).
884
885
886
887
888 Circles represent the theoretical location where the emphasis of each approach is situated across the three axes.
889

890 **Figure 2** Flow chart of the statistical analyses performed in this study.

891 **Figure 3** Frequency of local environmental variables (top) and spatial variables (bottom) selected as significant in explaining species distributions through binomial GLMs. Spatial variables are arranged from small-scale (i.e. V13) to large-scale extent (i.e. V1), and include V11 which was never selected. Climate variables are not shown (mean annual temperature was selected 32 times, and July air temperature 20 times). Species-environment relationships are shown in black when positive and grey when negative. Specific information for each species can be found in Table S2.
892
893
894
895
896
897

898 **Figure 4** Species scores on taxonomic vectors TAX-PCO3 (a) and TAX-PCO4 (b), and on trait vector TRA-PCO2 (c), arranged from lowest to highest values. For the trait vector, the location of the different functional feeding (FFG) and habit trait groups (HTG), is indicated, as well as the four different body size classes considered (BS, represented by columns of four different sizes), and those insects considered as of high female dispersal potential (H).
899
900
901
902
903
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /></figure>