Contribution of beta diversity in shaping stream macroinvertebrate communities among hydro-ecoregions

Rivers are heterogeneous and patchy-structured systems in which regional biodiversity of aquatic communities typically varies as a function of local habitat conditions and spatial gradients. Understanding which environmental and spatial constraints shape the diversity and composition of benthic communities is therefore a pivotal challenge for basic and applied research in river ecology. In this study, benthic invertebrates were collected from 27 sites across three hydro-ecoregions with the aim of investigating patterns in α- and β diversity. We first assessed the contribution to regional biodiversity of different and nested spatial scales, ranging from micro-habitat to hydro-ecoregion. Then, we tested differences in α diversity, taxonomic composition and ecological uniqueness among hydro-ecoregions. Variance partitioning analysis was used to evaluate the mechanistic effects of environmental and spatial variables on the composition of macroinvertebrate communities. Macroinvertebrate diversity was significantly affected by all the spatial scales, with a differential contribution according to the type of metric. Sampling site was the spatial scale that mostly contributed to the total richness, while the micro-habitat level explained the largest proportion of variance in Shannon–Wiener index. We found significant differences in the taxonomic composition, with 39 invertebrate families significantly associated with one or two hydro-ecoregions. However, effects of environmental and spatial controls were context dependent, indicating that the mechanisms that promote beta diversity probably differ among hydro-ecoregions. Evidence for species sorting, due to natural areas and stream order, was observed for macroinvertebrate communities in alpine streams, while spatial and land-use variables played a weak role in other geographical contexts.


Introduction
Biodiversity, which is commonly defined as the variety of life forms (e.g. species or taxa) that live in a certain habitat (Colwell 2009), is a key definition in community ecology with direct applications for basic and applied research. Measures and indices of biodiversity are generally used to get insights on the variation of biological communities along environmental gradients (Bo et al. 2016;Falasco et al. 2018) and in response to anthropogenic pressures Doretto et al. 2018Doretto et al. , 2019.
Quantifying biodiversity is, therefore, a crucial issue in ecological studies, and different approaches have been developed, but Whittaker's idea (Whittaker 1960), and its consecutive modifications (MacArthur and Pianka 1966), of decomposing the diversity into three components: alpha (a), beta (b) and gamma (c), is still one of the most recognized and employed in community ecology (Flach et al. 2012;Hepp et al. 2012;Matsuda et al. 2015). A major reason supporting the wide use of such an approach is that when biological communities are sampled over multiple sites, disentangling the contribution of the local diversity (a) and compositional changes (b) to the regional diversity (c) provides a better understanding of how biodiversity varies within the studied area. For instance, Piano et al. (2019a) compared macroinvertebrate and diatom communities between permanent and intermittent river reaches in north-western Italy and found that flow intermittency caused a reduction in the regional diversity (c) mainly acting as an environmental filter at local scale (a), while no relevant evidence of species turnover (b) occurred between the two type of sites. Decomposing the regional diversity into its a and b components appears particularly advantageous for heterogeneous and hierarchically structured systems, like rivers, where data are usually collected at multiple and nested spatial scales, such as micro-habitats, river stretches and basins (Ward and Tockner 2001;Tornwall et al. 2015;Er} os and Lowe 2019). In addition to these natural levels of organization, there are also anthropogenic levels-transects , meso-habitats (Burgazzi et al. 2017) and hydroecoregions (Li et al. 2001)-originated for biomonitoring purposes, including the standardization of sampling procedures and comparison between observed and expected community metrics (Bo et al. 2017).
Because riverine biodiversity varies according to both local environmental conditions and spatial gradients, it is of interest to assess the contribution of each spatial scale, as done in previously published works (Schmera and Er} os 2008;Clarke et al. 2010;Ligeiro et al. 2010). For instance, Karaus et al. (2013) evaluated the contribution of lateral aquatic habitats, such as parafluvial ponds, backwaters and tributary confluences, to the diversity of macroinvertebrate communities along three river corridors. They found that b diversity among the different aquatic habitats most contributed to the total c diversity (Karaus et al. 2013).
In addition to diversity partitioning, Legendre and De Cáceres (2013) have proposed an innovative and alternative method to elucidate what are the causes of variation in b diversity at regional scale. Based on their approach, b diversity can be decomposed into two components: the species (SCBD) and local contribution to the beta diversity (LCBD). SCBD describes to which extent each species singularly affects the variation in b diversity within the community and has been tested in relation to various speciesspecific measure of ecological niche and habitat occupancy (Heino and Grönroos 2017;Siegloch et al. 2018). By contrast, LCBD is a comparative indicator of the ecological uniqueness of sites. The higher the value of a site, the greater its uniqueness in terms of community composition with respect to the other sites (Legendre and De Cáceres 2013). Hence, by scoring sampling sites based on their communities, LCBD provides direct and useful applications for river biomonitoring and management at regional scale. For instance, LCBD could be used to compare sites and as a tool to assess the achievements of restoration programmes. Recent studies, therefore, have correlated LCBD with environmental variables and spatial distances to evaluate habitat and spatial constraints on community composition (Tolonen et al. 2018;Li et al. 2020), or used LCBD as an informative indicator to identify sites to be protected to preserve riverine biodiversity (Ruhí et al. 2017).
Because rivers are among the most dynamic, varied and multifaceted environments, it is important to explore new approaches to identify the factors determining their biological diversity. Thus, linking the variation in biodiversity of aquatic community to the dendritic and hierarchically nested structure of lotic ecosystems is a current challenge for both the river conservation and application of the recent ecological frameworks (e.g. metacommunity theory). Yet, this demands for shared and robust approaches to quantify spatial and environmental factors promoting beta diversity in river networks.
In this study, benthic macroinvertebrates were collected across the three most representative hydroecoregions (HERs) in north-western Italy to examine the variation in diversity and community composition at the regional scale. Specifically, we first used diversity partitioning to assess the contribution of different and nested spatial scales, ranging from micro-habitat to hydro-ecoregion, to the regional biodiversity (aim 1). In accordance with previous studies (Ligeiro et al. 2010;Fornaroli et al. 2015;Tornwall et al. 2015), we hypothesized that when biodiversity is calculated as taxon richness, the higher spatial levels, such as site and HERs, mostly contribute to the total diversity as a consequence of the enhanced habitat heterogeneity and availability of ecological niches. On the contrary, we expected a stronger contribution of the lowest level (i.e. microhabitat) when diversity is calculated taking into account the relative abundance of macroinvertebrates because of the control imposed by local near-bed conditions on the in-stream distribution of macroinvertebrates according to their preferences. Second, we compared macroinvertebrate communities of the three HERs by testing differences in a diversity and taxonomic composition and modelling LCBD in relation to environmental and spatial variables for each HER (aim 2). Because selected hydroecoregions represent different geomorphological conditions varying from low-order mountain streams to lowland mid-order rivers, we expected significant differences in the taxonomic composition, with the contribution of environmental and spatial factors on LCBD varying according to the hydro-ecoregion.
To our knowledge, this is the first study that combines different approaches to evaluate sources and mechanisms of variation in regional biodiversity across spatial scales and hydro-ecoregions since HERs were defined for this study area, therefore providing potential information for river biomonitoring and management.

Area of study
We selected 27 sampling sites (Table 1) belonging to three different hydro-ecoregions (HERs): Western Alps (HER1), Southern Alps (HER4) and Padana Plain (HER6). These HERs represent relatively homogeneous areas in terms of climatic, geological and topographic conditions according to the European Water Framework Directive (WFD 2000/60/EC). Such categorization is a fundamental requisite introduced by this normative act to make river biomonitoring type-specific within the European Union (Heiskanen et al. 2004;Moog et al. 2004). The hydro-ecoregions Western Alps (HER1) and Southern Alps (HER4) account for low-order mountain streams (Fig. 1): The average elevation of the HER1 and HER4 sites was 1065 ± 443 m a.s.l. (mean ± SD) and 1004 ± 500 m a.s.l. (mean ± SD), respectively. The stream Strahler order (Strahler 1952(Strahler , 1957 of the HER1 and HER4 sites ranged from 3 to 5 and from 2 to 5, respectively. The sites of these hydro-ecoregions had similar land uses, represented mainly by natural areas (range: 82-87%) and followed by agricultural (range: 10-17%) and urbanized areas (range: 1-3%). By contrast, HER6 includes lowland rivers ( Fig. 1): the average elevation of the selected sites was 233 ± 100 m a.s.l. (mean ± SD), while the stream Strahler order varied from 4 to 7. On average, agricultural areas were the main land-use category (67%), followed by natural (21%) and urbanized areas (12%). The three HERs here considered are representative in terms of spatial extent of the variety of hydro-ecoregions present in north-western Italy (Fig. 1).

Macroinvertebrate sampling
From 2014 to 2016, each site was sampled once in two different seasons: autumn or late summer. However, a preliminary analysis showed that the taxonomic composition of macroinvertebrate communities did not vary seasonally (PERMANOVA: F 1,26 = 1.125; P = 0.302); hence, the factor season was not included in the statistical analyses. At each site, ten Surber samples were collected according to the official sampling protocol adopted by the Italian legislation for wadeable rivers (Bo et al. 2017), for a total of 270 samples (10 samples 9 9 sites 9 3 HERs). The application of the Italian methodology requires a multi-habitat proportional sampling whereby individual samples are collected proportionally to the microhabitat (i.e. substrate type) occurrence. A Surber net (0.05 m 2 , 250 lm mesh-size) was placed flush with the river bottom, and benthic macroinvertebrates were collected by perturbing and washing the surface of substrate elements inside the metal frame of the net (Doretto et al. 2020a). Benthic invertebrates were preserved in the field in 90% ethanol labelled jars, transported to the laboratory where samples were sorted and all individuals were counted and systematically identified to family level under a microscope (NIKON SMZ 1500 light microscope; 60-100 X), following the taxonomic keys for the Italian macroinvertebrate fauna (Campaioli et al. 1994(Campaioli et al. , 1999. Such taxonomic resolution is the level officially required for river biomonitoring in Italy (Bo et al. 2017).

Statistical analyses
We used the Additive Diversity Partitioning (Veech et al. 2002;Crist et al. 2003) to assess the contribution of each spatial scale to the total diversity of macroinvertebrate communities. This technique is explicitly designed for nested datasets and decomposes the regional diversity (c) measured at the higher scale into the sum of a and b components measured at the lower spatial scales. In this study, c represented the total diversity obtained by pooling together all the samples, while a was the diversity within the lowest spatial scale (i.e. Surber samples), b1 was the diversity among Surber samples, b2 was the diversity among sites and b3 was the diversity among HERs. The null hypothesis was that the observed values for alpha and beta components at each spatial scale did not differ significantly from those predicted by a random distribution of taxa across samples at all spatial scales. An individual-based randomization process was used to randomly assign the number of individuals and taxa at all spatial scales, keeping their initial abundance equal (Crist et al. 2003). This procedure generated 999 null distributions that were used to calculate the null alpha and beta values at each spatial scale. The statistical significance was calculated as the proportion of null values that were, respectively, higher or lower than the observed values. This analysis was performed separately for two community metrics: the total taxon richness and Shannon-Wiener index.
Significant differences in the community metrics and taxonomic composition among HERs were tested with different statistical analyses. The Kruskal-Wallis test was applied to taxon richness and macroinvertebrate density, while ANOVA and Tukey test were used for the Shannon-Wiener index. To avoid pseudoreplication, Surber samples (N = 10) from each site were pooled together prior to calculating these community metrics.
Significant changes in the taxonomic composition among HERs were visually and statistically examined by means of non-metric multidimensional scaling (NMDS) and permutational analysis of variance (PERMANOVA), respectively. Surber samples were pooled together for each site, and the Bray-Curtis dissimilarity index was applied on macroinvertebrate abundances. Moreover, indicator species analysis (Dufrêne and Legendre 1997) was performed to determine whether individual macroinvertebrate families were indicators of HERs.
Finally, differences in the local contribution to beta diversity (LCBD) were calculated following the approach of Legendre and De Cáceres (2013). This statistical technique allowed us to quantify the proportion by which each site-specific macroinvertebrate community contributed to the total diversity in the area of study. A site with a high LCBD value hosts a very unique macroinvertebrate community compared to the other sites and contributes more than others to the regional biodiversity. For each site, we pooled together all the Surber samples to obtain an overall community, and macroinvertebrate abundance data were Hellinger-transformed prior to calculating LCBD, for which the Bray-Curtis dissimilarity index was used.
Generalized linear models (GLM) were run separately for each hydro-ecoregion to evaluate the effects of environmental and spatial variables on LCBD. For each sampling site, the former set of explanatory variables included: elevation (m a.s.l.), stream Strahler order and the percentages of natural, agricultural and urbanized areas (1 km radius). Geographical coordinates were used to produce a new set of six spatial variables by means of the principal coordinates of neighbour matrices (PCNM) methodology (Borcard and Legendre 2002). This technique belongs to a broader category of spatial analyses (Moran's eigenvector maps-MEM- Dray et al. 2006) and applies Euclidean distance to geographical coordinates of sampling sites (i.e. longitude and latitude) to generate a new set of orthogonal variables (i.e. eigenvectors) describing spatial patterns among sites. Model selection was performed according to the AIC (Akaike's information criterion) to keep only the most relevant variables, and variance partitioning analysis was then performed to separate the pure and joint fractions of variance in LCBD explained by the environmental and spatial variables, respectively.
All the statistical analyses (significance threshold = 0.05) were performed with the software R (R Core Team 2019) using basic functions and the following packages: vegan (Oksanen et al. 2015

Results
A total of 30,243 macroinvertebrates belonging to 69 families were collected: 44 families were found in the HER1, while the total number of macroinvertebrate families in each of the HER4 and HER6 was 58. Sample-based taxa accumulation curves indicated that representative taxonomic communities were collected for each HER (Supplementary materials Fig. S1). The average number of taxa per sample (± SD) was 12 (± 3.80), 14 (± 5.44) and 11 (± 3.73) in the HER1, HER4 and HER6, respectively. The mean number of individuals per sample was 101 (± 76.54) in the HER1, 96 (± 75.31) in the HER4 and 139 (± 121.54) in the HER6. In general, the most abundant families were Chironomidae (18%), Baetidae (16%), Gammaridae (10%) and Heptageniidae (7%). These families together represented 51% of the community in this study.
Contribution of spatial scales to the macroinvertebrate diversity Results of the additive diversity partitioning showed that the b2 component (i.e. diversity among sites) alone accounted for 41% of the total variance in the taxon richness (Fig. 2a) and was significantly higher than predicted (Table 2). On the contrary, the other three components contributed almost equally to c diversity: a (i.e. diversity within Surber samples) and b1 (i.e. diversity among Surber samples) explained 17.5% and 19.9%, respectively, and their percentages of explained variance were significantly lower than predicted ( Table 2). The hydro-ecoregion (b3) accounted for 21.6% of the total diversity, and its contribution was significantly higher than expected (Table 2).
When diversity was calculated using the Shannon-Wiener index, we found differences in the contribution of the components (Fig. 2b). Although the observed variance explained by the a component was significantly lower than predicted, it alone accounted for 61.8% of the total diversity (Fig. 2b, Table 2). By contrast, the observed percentages of variance attributable to b1, b2 and b3 were significantly higher than predicted, but these components generally contributed for a minor extent to the c diversity. b1 and b2 accounted for 12% and 18%, respectively, while the hydro-ecoregion (b3) explained 8.2% of the variance associated with the Shannon-Wiener index (Fig. 2b).

Comparison of macroinvertebrate communities among HERs
No significant differences were observed in the total taxon richness and macroinvertebrate density among hydro-ecoregions (Fig. 3a, b). By contrast, Shannon-Wiener index significantly varied among HERs, but pairwise comparisons showed significant differences only between HER4 and HER6 (Fig. 3c).
Results of the NMDS ordination (Fig. 4) and PERMANOVA showed significant differences in the taxonomic composition of benthic communities among HERs (PERMANOVA: F 2,267 = 12.8, P \ 0.001). Moreover, indicator species analysis identified 39 macroinvertebrate families as indicator taxa of one hydro-ecoregion or combinations of them (Table 3).
HER1 had two exclusive indicator taxa, namely Chloroperlidae and Blephariceridae, while the number of macroinvertebrate families significantly and exclusively associated with HER4 and HER6 was 6 and 12, respectively. In HER4, the majority of these were EPT taxa, such as Ephemeridae, Ephemerellidae, Goeridae and Oligoneuriidae, whereas HER6 had a higher   (Table 3). When looking at combinations of HERs, the highest number indicator taxa was found for HER1-HER4 group, with ten EPT families out of 14 macroinvertebrate families (Table 3). On the contrary, two (i.e. Naididae, Empididae) and three (i.e. Ceratopogonidae, Goeridae and Gyrinidae) families were statistically associated with the HER1-HER6 and HER4-HER6 groups, respectively (Table 3). LCBD varied from on site to another one (Supplementary Materials Fig. S2), but, on average, no significant differences were found among HERs (Fig. 5a). Differences were observed, instead, in the relative effects of environmental and spatial variables on LCBD. Based on the model selection, four  (Table 4). Together these variables explained up to 88% of the total variance in LCBD for this hydro-ecoregion, with the pure percentage of environmental (i.e. stream Strahler order and percentage of natural areas) and spatial variables (i.e. PCNM4 and PCNM2) equal to 57% and 21%, respectively. Their combined percentage was 10%, while 12% of the variance was unexplained (Fig. 5b). On the contrary, percentage of urbanized areas was the environmental variable selected within HER4 along with four spatial variables: PCNM3, PCNM5, PCNM6 and PCNM4, but none of these affected significantly the LCBD ( Table 4). The overall percentage of variation explained by the selected predictors (47%) was almost entirely attributable to the spatial descriptors, which accounted for 46% of the total variance in LCBD in this hydro-ecoregion, while the percentage of urbanized areas explained 1%. The combined proportion of these variables was negligible and 53% of the variance was unexplained (Fig. 5c).
Two explanatory variables were included in the statistical analysis for HER6 after the model selection: the percentage of agricultural areas and PCNM6, but their effects on LCBD were not significant (Table 4). These two variables accounted for 36% of the total variance in LCBD, while a relatively high percentage (64%) remained unexplained (Fig. 5d). The amount of variation explained by the percentage of agricultural areas (24%) was twofold higher than PCNM6 (12%), while their combined effect was negligible (Fig. 5d).

Discussion
Lotic ecosystems are highly heterogeneous and patchy-structured habitats in which regional biodiversity of benthic communities is expected to vary according to both local environmental conditions and spatial gradients (de Mendoza et al. 2018;Alther et al. 2019;Harvey and Altermatt 2019). In this study, variation in the a and b diversity of stream macroinvertebrate communities across hydro-ecoregions was examined and our results support previous studies.
We first assessed how biodiversity varied over multiple and nested scales (aim 1) and found that diversity metrics responded differently across spatial scales. Variation among sites (i.e. b2 in the diversity partitioning analysis) and hydro-ecoregions (i.e. b3) were the components that accounted for the largest amount of explained variance in taxon richness. This finding suggests that site and hydro-ecoregion are the spatial scales that mostly contribute to the total richness, probably through their habitat heterogeneity which, in turn, provides different ecological niches   (Tornwall et al. 2015;Bo et al. 2016;Piano et al. 2019b). In a multi-scale study, Ligeiro et al. (2010) evaluated patterns in macroinvertebrate family richness among four nested spatial scales, ranging from Surber sample to river segment. Similar to our study, authors found that larger spatial scales explained the greater proportions of variance in the total diversity (Ligeiro et al. 2010).
When the relative abundance of benthic invertebrates was taken into consideration by using the Shannon-Wiener index, we found an opposite trend. Around three-quarters of the variation in c diversity were explained by the a and b1 components, which corresponded in our analysis to the micro-habitat scale (i.e. Surber samples). In a field study in Brazilian headwater streams, Siegloch et al. (2018) found that sampling site, followed by ecoregion, was the spatial scales that mostly contributed to the total family richness in aquatic insects metacommunities. These patterns are similar to our findings for taxon richness. However, the contribution of a component increased when the diversity partitioning was performed using the Shannon-Wiener index (Siegloch et al. 2018), as we observed in our study. Unlike taxon richness, these findings indicate that the relative abundance of stream macroinvertebrates responds more than taxon richness to small-scale factors. According to previous studies, near-bed conditions (i.e. hydraulic forces and substrate characteristics) act as strong predictors of the microhabitat distribution of macroinvertebrates because of their filter effects on the species-specific traits (Beisel et al. 1998;Lamouroux et al. 2004;Sagnes et al. 2008;Fornaroli et al. 2015).The availability of food resources in streams is often patchy distributed, depending mostly on the substrate size and stability (Minshall 1984;Fenoglio et al. 2005Fenoglio et al. , 2006Doretto et al. 2016Doretto et al. , 2017. Similarly, water velocity and depth affect macroinvertebrate density by acting on some life-history habits, including respiration (Brooks & Haeusler 2016), aquatic dispersal (Gibbins et al. 2016) and feeding strategy (Mérigoux and Dolédec 2004). In a field experiment, Burgazzi et al. (2018) collected macroinvertebrates from 50 sampling plots located in five longitudinal 10-m transects. The authors found that macroinvertebrate abundance was spatially autocorrelated and negatively affected by water depth (Burgazzi et al. 2018).
When macroinvertebrate communities were compared among HERs (aim 2), the multivariate analysis showed clear changes in the taxonomic composition with 39 benthic invertebrate families significantly associated with one or two hydro-ecoregions. Although family level is often considered as a coarse taxonomic resolution compared to genus or species level, it was successful in capturing compositional changes among hydro-ecoregions in this study. Moreover, these results confirm that hydro-ecoregions contribute to the regional diversity by hosting some distinctive taxa. Both HER1 (Western Alps) and HER4 (Southern Alps) include low-order mountain streams with analogous topography and slight geological and climatic differences. Such a comparability in abiotic conditions likely explains the high number of shared EPT taxa, which are a key faunal component of alpine streams (Doretto et al. , 2020b. By contrast, HER6 (Padana Plain) encompasses lowland and large rivers, including the downstream sections of many HER1 and HER4 sites here considered (Fig. 1). Nevertheless, we found distinct macroinvertebrate communities with a greater occurrence of non-EPT and non-insect taxa.
In accordance with previous research on the longitudinal distribution of benthic invertebrates (Vannote et al. 1980;Usseglio-Polatera et al. 2000;Tachet et al. 2010;Manfrin et al. 2013), clear up-to downstream shifts from EPT-to non-EPT-dominated macroinvertebrate communities were observed in our watercourses. In the light of this, our findings highlight the importance of a reference-based river biomonitoring in which the ecological status is calculated as the ratio between observed and reference communities (Heiskanen et al. 2004;Moog et al. 2004;Verdonschot 2006). In addition, the indicator families here identified may serve as target taxa or groups for these hydroecoregions in relation to specific river management, monitoring and restoration goals. For instance, specific indices or metrics could be developed based on these indicator taxa and, consequently, make river biomonitoring even more type specific. Similarly, their habitat preferences may be used as primary criteria for river restoration programmes because our results proved that these indicator taxa are key faunal components of the selected hydro-ecoregions.
The analysis of the local contribution to beta diversity (LCBD-sensu Legendre and De Cáceres 2013) enabled us to quantify the ecological uniqueness of sites. Each HER included sites with high and low LCDB values, but results of the variance partitioning showed the differential effects of environmental and spatial constraints on LCBD according to HER, showing that the mechanisms that promote b diversity probably differed among hydro-ecoregions. The percentage of natural areas and stream Strahler order played a significant role for HER1 sites, followed by spatial variables. These findings suggest that, in this hydro-ecoregion, benthic invertebrate communities are primary shaped by the local habitat conditions, especially the land use and stream size, while spatial processes have a minor role (i.e. species sorting scenario- Leibold et al. 2004;Brown et al. 2011;Heino et al. 2015). On the contrary, a greater proportion of unexplained variance in the LCBD was observed within HER4 and HER6 sites, with minor and not significant contributions of the spatial variables and percentage of agricultural areas, respectively.
Our findings corroborate previous works focused on macroinvertebrate metacommunities and alpine streams (Kitto et al. 2015;Scotti et al. 2020). For instance, in a large-scale study encompassing 50 Scandinavian streams, Tolonen et al. (2018) partitioned the environmental and spatial controls on the LCBD; with the original hypothesis that a greater proportion of explained variance due to local factors should indicate the effect of environmental filters (i.e. species sorting), while whether LCBD is mainly accounted for by spatial variables then stochastic or dispersal-related processes (e.g. mass effect) are the main drivers of the metacommunity dynamics. They found weak correlations between LCBD and spatial descriptors, while environmental variables at site scale explained generally a greater proportion of the total variance in LCBD (Tolonen et al. 2018). Tonkin et al. (2016) used the catchment area, elevation and land use as proxy variables of the river network position and tested their relationships with LCBD in 124 streams belonging to five hydro-ecoregions. The authors found highly context-dependent results, with correlations varying according to the hydro-ecoregion (Tonkin et al. 2016).
To conclude, this work sheds light on the patterns and mechanisms underlying the variation in diversity and community composition of macroinvertebrates among the three most-represented hydro-ecoregions in north-western Italy. Although further investigations are necessary to implement the results here obtained, this study enhances our basic knowledge on what factors shape benthic communities and also provides useful information of practical interest. First, our results show that a multi-scale approach is a good choice for investigating biodiversity of benthic invertebrate assemblages. Site is a key spatial scale, but also lower or higher hierarchical levels (i.e. micro-habitat or hydro-ecoregions) could have a relevant role depending on the type of diversity metric. Moreover, these findings stress the importance of using multihabitat sampling protocols to make data collection as much representative as possible (Böhmer et al. 2004;Barbour et al. 2006;Mondy et al. 2012). Second, environmental and spatial controls on macroinvertebrate communities are context dependent and vary from one hydro-ecoregion to another one. Under a reference-based river biomonitoring context, this finding may help in defining the best strategies for the conservation and management of lotic ecosystems and their biota, preserving either the habitat conditions or the connectivity and species dispersal pathways, depending on the relative role of environmental or spatial controls on the benthic communities.
Availability of data and material Data that support the findings of this study are available from the corresponding author upon reasonable request.

Compliance with ethical standards
Conflict of interest No potential conflict of interest was reported by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.