Spatial and temporal variation of benthic macroinvertebrate assemblages during the glacial melt season in an Italian glacier-fed stream

The biodiversity of glacier-fed streams is particularly threatened by climate change, emphasising the need of monitoring these sentinel systems. The glacier-fed Saldur stream is an International Long Term Ecological Research (ILTER) site in the Italian Central Eastern Alps. Here, we sampled benthic macroinvertebrates and measured environmental variables (discharge, suspended solids, conductivity, water temperature, and channel stability) five times at six sites (5–11 km from the glacier) during an entire glacial melt season (April–September). Our main objectives were (1) to elucidate relationships between the abiotic variables and the faunal composition, (2) to quantify and compare the spatial and temporal variability of the faunal community, and (3) to assess the composition of the benthic macroinvertebrate community in relation to conceptual models. Hosting a higher number of individuals and more diverse communities at sites with reduced glacial influence, the Saldur stream fitted well in the framework of conceptual models. Nevertheless, the spatial variability of the fauna was higher than the temporal variability. This study presents an initial characterisation of the benthic faunal assemblages in the Saldur stream, constituting a reference point for future analyses dealing with potential disruptive factors introduced by climate change and upcoming hydroelectric power production on this stream.


Introduction
High Alpine glacier-covered regions are expected to be particularly exposed to climate change effects in the next decades (Schädler & Weingartner, 2010;Jordan et al., 2016). Consequently, streams draining Alpine areas are generally identified as one of the environments most sensitive to climate change (Füreder et al., 2002;Viviroli & Weingartner, 2008; Handling editor: Marcelo S. Moretti Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10750-018-3731-8) contains supplementary material, which is available to authorized users.
Jiménez Cisneros et al., 2014). This is particularly true for glacier-fed streams and rivers. The forecasted changes in climate will eventually result in the disruption of the current sediment budgets and hydrographs Brown et al., 2018). More specifically, the discharge peak of glacier-fed rivers may shift from the summer towards the spring season, while the total annual runoff will increase in the first phase but decrease in the long term Huss et al., 2010;Huss, 2011;Huss & Hock, 2018).
Changes in the abiotic environment will have implications for aquatic biota. Locally, the decreasing percentage of glacial cover in a catchment increases the taxonomic richness of algae (Rott et al., 2006) and of benthic aquatic macroinvertebrates (Brown et al., 2007c;Milner et al., 2008;Jacobsen et al., 2010). Nevertheless, from a long-term perspective at the global level, the shrinkage of glaciers and the loss of permanent snow (Huss et al., 2017) may eventually lead to the extinction of several alpine freshwater taxa Finn et al., 2013;Giersch et al., 2015).
Given the expected changes of the cryosphere, monitoring these susceptible habitats and the organisms inhabiting them becomes highly important (Brown et al., 2006). Indeed, Khamis et al. (2014) suggested the use of invertebrates as indicators of environmental change in relation to glacial retreat, while Giersch et al. (2017) clearly linked the loss of alpine insects to the consequences of climate change, both stressing the need for further research.
With the goal of deepening the knowledge of glacier-fed stream ecosystems in the Italian Alps, the Saldur stream, a glacier-fed stream in the Italian Central Eastern Alps, is particularly interesting due to the intensification of monitoring and sampling activities in the entire catchment beginning in 2015. Indeed, between spring 2014 and autumn 2015, a ''run-of-river'' (ROR) hydropower plant designed to produce 3200 kW a year was built at * 2000 m a.s.l., and the stream and its catchment became part of the International Long Term Ecological Research (ILTER) network (site IT-25) in 2014.
Considering the 2015 melting season, the objectives of the present study were (1) to elucidate the relationships between a set of abiotic variables of the glacier-fed stream varying in time and space and the composition of macroinvertebrate assemblages, (2) to quantify and compare the effects of space (i.e. the positions of the sampling sites, 5-11 km from the glacier) and time (i.e. the time of sampling) factors on the community structure of the fauna, and (3) to assess and to characterise the composition of the benthic macroinvertebrate community and abiotic parameters in the Saldur stream before the implementation of the hydropower plant.
We hypothesise that (1) the density and diversity of macroinvertebrates is highest when and where environmental harshness is low (low glacial influence). Further, we suggest that (2) the benthic macroinvertebrate assemblage responds more to season than to site variation (month-to-month variation is higher than site-to-site variation), due to the high variability of environmental conditions during the glacial melt season. Finally, (3) as the stream flow is still unaffected in terms of quality and quantity, we assume that faunal assemblages, as well as environmental parameters, fit into the conceptual model for spatial and temporal patterns described previously Füreder, 2007;Brown et al., 2007b).

Study area
The Saldur stream drains the Matscher Valley, an eastern side valley of the Vinschgau Valley in South Tyrol, Italy. The catchment (101 km 2 ) has been an ILTER site since 2014 (IT-25), located within the Central Eastern Alps (N 46°, E 10°) and spanning between 921 m a.s.l. and 3738 m a.s.l. The geology of the area is characterised by the presence of gneisses, mica gneisses, and schists (Habler et al., 2009). The climate is relatively dry and has a total annual precipitation of approximately 500 mm measured at 1570 m a.s.l., with a clear distribution showing a minimum in winter and a maximum in summer (source: Hydrographical Office of the Autonomous Province of Bozen). The land cover is a mixture of bare soil and rocks in the upper part of the catchment, while the lower part is characterised by more forests, pastures, and grasslands. The anthropogenic activities in the catchment are mainly limited to agriculture (primarily animal husbandry, mowing, and cultivation of strawberries).
The Saldur stream originates from the Matsch Glacier, whose surface area decreased from 2.8 km 2 in 2006 (Galos & Kaser, 2014) to 2.2 km 2 in 2013. The stream has some small, intermittent, groundwater-fed tributaries. The hydrological dynamics characterising the stream are nivo-glacial, with the snow melting generally occurring between June and July and the glacial melting occurring in August. At a gauging station at 1632 m a.s.l. run by the Hydrographical Office of the Autonomous Province of Bozen, typical baseflow conditions show discharge values around 0.6 m 3 /s. During the melt season, these values span between 5 m 3 /s and 10 m 3 /s, with peaks up to 15 m 3 /s in the case of storm events.
Between spring 2014 and autumn 2015, a flowing water hydropower plant, having its weir at 2012 m a.s.l. and producing up to 3200 kW a year, was built. Hydroelectric power production started in December 2015, with the current concession (2017) allowing a production of 1600 kW a year.

Study design
For the present study, six sites along the Saldur stream were chosen. Sites 2A, 2B, and 2C were chosen adjacent to each other with a distance of 100-160 m, to be able to assess the impact of the hydropower plant in relation to the position of the weir in the near future. For the others, the denomination of the sites preserved the names assigned during previous monitoring campaigns (Fig. 1).
The investigated stretch of stream from site 1 to site 3 was 6.2 km long, and the altitude ranged between 2030 m a.s.l. at site 1 and 1645 m a.s.l. at site 3. The altitude and position of each site relative to the glacier are reported in Table 1.
On a monthly basis, a set of environmental variables and macroinvertebrate samples was collected at all sites during 2-day periods from April 2015 to September 2015. Each sample is identified by the site code coupled with the month of collection (AM April/May, JE June, JY July, AU August, SE September).

Abiotic variables
Spot measurements of water temperature (tem) with an alcohol thermometer were recorded at each site. In addition, water samples were collected in 3 polyethylene (PE) bottles of 100 ml and in 3 PE bottles of 1 l to assess conductivity (con) and suspended solids (ss), respectively. The channel stability (pfa) was estimated for each site by calculating the bottom component of the Pfankuch index (Pfankuch, 1975). Data on the hydrological regime, namely discharge (dis), were retrieved at site 3 using a stream gauging station, logging values every 10 min.
In the laboratory, electrical conductivity (EC, at 25°C) was measured using an EC metre (COND7; XS Instruments, Carpi, Italy) calibrated with a standard solution of 84 lS/cm. Suspended solids were measured after 30 min of sedimentation of each water sample in a plastic Imhoff cone (APAT & IRSA-CNR, 2003), using the volume of wet particles (ml l -1 ) as surrogate. As a final value for each site during the five different sampling events, the average of the three water samples for conductivity and suspended solids was calculated.
As additional environmental variables, we considered the following for each site: chainage of the sites (cha), defined as the distance from the glacier snout that was calculated with ArcGIS (Esri, 2014) and the Glacial Index (gla), an index of glacial influence that was calculated by combining glacier size with the distance from the glacier terminus, as defined by .

Macroinvertebrates
Benthic macroinvertebrates were sampled at each of the six stream sites on five occasions. On each sampling event, we collected a total of 12 Surber samples (0.0506 m 2 , mesh size 500 lm) at random, covering the main substrate types present at each site. Four Surber samples were combined into one sample such that each sampling from a site provided three separate replicate samples. All replicates were labelled and preserved in the field in 70% ethanol. For data treatment, the three separate replicates from one site were summed together, resulting in the sample codes mentioned before (e.g. 1-AM, 2A-AM).
In the laboratory, Surber samples were sorted, and the fauna were identified to the lowest possible taxonomic level under a stereoscopic microscope at 915 magnification, referring to Campaioli et al. (1994) and Sansoni & Ghetti (1998). Of the 27 total taxa, 9 were classified to the genus or species level and 18 (including chironomids) to the family level.

Data treatment and statistical analysis
The benthic macroinvertebrate community was first analysed by calculating density, taxonomic richness, the Shannon-Weaver index (H 0 ), the Simpson dominance index (1-D), and the Pielou evenness index (J) for each site and month. Average agglomerative clustering through the Unweighted Pair-Group Method using Arithmetic averages (UPGMA) was applied to the faunal data after chord transformation, with the number of clusters defined by fusion levels (Legendre & Gallagher, 2001).
To compare the effects of time and space, all possible pairs of elements (months or sites, not mixed) were tested using a non-parametric Friedman test (with Bonferroni correction). We tested for differences in the faunal indices, the abiotic variables, and the taxa that scored more than 200 findings throughout the sampling period (seven in total). The same dataset was explored for Pearson correlations, again considering all possible pairs of sites and months.
Faunal metrics (ind./m 2 , taxonomic richness, % EPT, H 0 , 1-D, J) were further analysed by agglomerating samples on the basis of months and sites. Additionally, a coefficient of variation (CV) for each metric linked to each site or month was calculated (e.g. agglomerating all the samples collected throughout the season at a site or agglomerating all the samples Finally, a canonical correspondence analysis (CCA) was computed to analyse the relationship between the raw fauna and environmental data matrices. The overall analysis and the axes were tested through a Monte Carlo permutation test (n = 999). The resulting biplots were plotted using scaling 1 and scaling 2. The former approximates the position of the samples along the environmental variables (Borcard et al., 2011), and the latter shows the optimum taxa along a set of quantitative environmental factors (Borcard et al., 2011), emphasising relationships among taxa (Legendre & Legendre, 2012).

Abiotic conditions
Discharge, suspended solids, conductivity, and temperature showed clear seasonal patterns in all sites related to glacial-melting dynamics, with significant differences (Friedman tests, P \ 0.05) among the sampling events, particularly when considering the month of August (Tables 2 and 3).
The Glacial Index decreased with increasing distance from the glacier snot (cha) (r = -0.990, P \ 0.001). Values of Glacial Index overall showed a limited glacial influence in absolute terms (Table 1). The Pfankuch index showed a similar trend, increasing Cobbles, boulders Mean annual current (m 3 /s) and mean annual water level (m) were recorded only at site 3 from low channel stability at sites 1-2C and high stability at sites 2D and 3 ( The UPGMA cluster analysis based on taxon abundances revealed two distinct but intertwined patterns for grouping of the samples: both site position and month were identified as factors driving the composition of the clusters. In the structure of the six clusters highlighted in Fig. 2, two tendencies emerged: (1) in terms of position, sites tended to be  Suspended solids-ss (ml l -1 ), conductivity-con (lS cm -1 ), temperature-tem (°C), Pfankuch index-pfa. Mean ± standard deviation Abbreviations for sampling events: AM April/May, JE June, JY July, AU August, SE September clustered into two groups, the first one comprising samples from sites 1-2C and the second one grouping samples from site 2D together with the ones from site 3; (2) in terms of time, the months June and July and the months August and September were grouped together. The first sampling event, from April to May, tended to be represented in a unique cluster (Fig. 2). Density (ind./m 2 ), taxonomic richness, % EPT, H 0 , 1-D, and J showed an overall downstream increasing trend during the study period (Fig. 3). Nevertheless, neither temporal nor spatial patterns, in terms of faunal indices, were found to be significantly different among any pair of sampled months or sites.
In terms of temporal variability, the only taxa whose abundances occasionally differed between months at least at some sites (P \ 0.05) were Rhabdiopteryx sp., Protonemura sp., Chironomidae, Baetis sp., and Perlodidae (Table 4). For Rhabdiopteryx sp. and Chironomidae, the Friedman test identified significant differences in the overall model (P \ 0.05) for Protonemura sp., Baetis sp., and Perlodidae differences were found between July and September, June  Densities of four of the examined taxa (Baetis sp., Limoniidae, Limnephilidae, Perlodidae) were inversely correlated to the Glacial Index and, accordingly, with distance from the glacier (chainage). In contrast, Chironomidae was the only taxon positively correlated with the Glacial Index. The channel stability index was negatively correlated with four taxa (Baetis sp., Limoniidae, Limnephilidae, Perlodidae) and positively correlated with Rhabdiopteryx sp.; Chironomidae and Rhabdiopteryx sp. were positively correlated with discharge. Finally, a strong positive correlation was observed between the density of Rhabdiopteryx sp. and conductivity (Table 5).
Grouping the samples by month and study site, the faunal density showed higher variability in the spatial component (Fig. 4b, average CV 0.562) than in the temporal component (Fig. 4a, average CV 0.405). The same pattern was found for taxonomic richness (average CV 0.241 and 0.182, respectively). Regarding the other faunal metrics, temporal and spatial components showed similar patterns in CV (Fig. 4).
The CCA showed significance through a permutation test (F = 6.273, P = 0.001). The amount of variance explained by the collected environmental variables was 67%. The permutation test of each axis showed significance of the first four axes, explaining 61.25% of the constrained variance ( Table 6). The CCA biplots with scaling 1 reinforced the trend shown in the UPGMA clustering, with samplings tending to group based on site position and time of sampling (Fig. 5a, b, respectively). The CCA biplot using scaling 2 (Fig. 6) confirmed the results from the Friedman tests and the correlation analyses: a clear trend showing an affinity of most of the taxa for reduced Glacial Index values, greater distance from the glacier snout, and greater Pfankuch index values was evident.   The seasonal glacial melt process strongly influenced the dynamics of the benthic macroinvertebrate community, as also reported in similar European glacier-fed systems (Ward, 1994;Burgherr et al., 2002). The cyclical ordination of the sampling events shown in the CCA (Fig. 5b) (Armitage et al., 2001;Leunda et al., 2009) stemmed from a strong seasonality factor for specific groups of macroinvertebrates, whose densities were associated with the variability in the environmental factors. Moreover, the grouping of the sampling events by time was also independently identified by the classification analysis (clusters). Among the seven most abundant taxa, Rhabdiopteryx sp. and Protonemura sp. showed significant sensitivity to fluctuations of some abiotic factors (Füreder et al., 2005), with two opposite behaviours: the former preferring early-and late-season conditions (high conductivity and low discharge, channel stability, and temperature) and the latter preferring midseason conditions (high discharge and high temperature). Previous studies on the life cycle of Protonemura sp. identified strong interspecific differences in terms of thermal demands, so our taxonomic resolution cannot allow us to infer more regarding the behaviour found in our results (Haidekker & Hering, 2008). On the other hand, a previous study on Rhabdiopteryx sp. in the alpine region supports our results, clearly showing that there is larval growth during winter snow cover, which explains the prominent emergence in spring .
In terms of space, the three taxa whose abundance differed significantly among sites (Limoniidae, Limnephilidae, Perlodidae) were closely associated with sites 2D and 3, giving a first confirmation of the progressive presence of new taxa along a glacier-fed stream modelled by .

Temporal versus longitudinal variability
Disentangling the spatio-temporal variability of the benthic assemblage, the CV analysis of density and taxonomic richness revealed a greater importance of the site-to-site variation than of the month-to-month variation. In our results, there was a clear trend of diminishing CVs in taxonomic richness, with a general trend, except for site 2A, of increasing CV in density-accompanied by increasing channel stability and decreasing glacial influence of the sites (Fig. 4), as also reported in other studies (Burgherr & Ward, 2001;Lods-Crozet et al. 2001).
Despite the higher abundance and taxonomic richness of macroinvertebrates in sites 2D and 3, the high CVs exhibited by sites 1-2C suggest that the community closer to the glacier does not possess a strong core structure, as shown by the sites located more downstream, whose assemblages tended to remain more constant across the sampling period. Indeed, we observed much higher temporal dynamics in taxon turnover in the higher sites, explaining the higher CV values in comparison to the lower sites, e.g. analysing the July sampling event compared to June, at site 1, we sampled 7 new taxa, whereas at site 3, we sampled only 3 new taxa. Thus, it seems that in particular periods, although the local glacial influence may be relatively high, there is the ecological space and opportunity for new taxa to appear outside the classical ''windows of opportunity'' of spring and late summer/autumn (Uehlinger et al., 2002).

Fit to the conceptual models for glacier-fed streams
The Saldur stream fitted the habitat template for glacier-fed rivers first presented by Milner & Petts (1994), and furtherly developed by Milner et al. (2001). Indeed, in spatial terms, the model strongly relies on environmental gradients to describe the distribution of benthic macroinvertebrates along a glacier-fed stream, particularly referring to stream temperature and channel stability (Milner & Petts, 1994;Milner et al., 2001). These factors both followed a clearly defined trend during the sampling period in the Saldur stream.
As expected, water temperature generally increased with increasing distance from the glacier, mostly due to solar radiation and atmospheric heating (Maiolini & Lencioni, 2001;Kuhn et al., 2011;Khamis et al., 2015). The other main conditioning factor of the conceptual model, channel stability, behaved similarly: the Pfankuch index decreased according to the distance from the glacier (chainage), due to a more stable channel morphology associated with reduced surfaces involved in scouring and deposition processes, and according to reduced glacial influence (Milner & Petts, 1994). More specifically, a considerable shift in channel stability was identified between the four most upstream sites and the remaining two (2D and 3).
In full agreement with the conceptual model, this pattern was also observed for the descriptive metrics density and diversity of the benthic macroinvertebrate assemblages throughout the entire sampling period (Milner et al., 2010). In this framework, the whole set of calculated faunal metrics suggested that the density and the diversity of the sampled assemblages were highly influenced by the spatial patterns and ultimately by a reduced environmental harshness-also highlighted by the increased channel stability-derived by decreased glacial influence for sites 2D and 3 (Ilg & Castella, 2006;. Similarly, a more recent spatial analysis in glacierfed streams proposed by Füreder (2007) linked values of a set of faunal metrics to the percentage of glaciation in a catchment. Although our longitudinal transect covered only a limited pattern of glaciation (Table 1), its effect was visible in our results, consistent with the values of diversity, evenness, and number of taxa described by Füreder (2007), for our level of glaciation. In contrast, abundance of macroinvertebrates was consistently lower in our case. Potentially, this is due to the use of a 500-lm mesh-sized sampling net and the strong interannual variability in macroinvertebrate abundances naturally occurring in glacierfed streams (Lods-Crozet et al., 2001;Rüegg & Robinson, 2004).
The CCA graphically represented this concept in relation to the distribution and succession of taxa along the stream (Fig. 6). Indeed, quadrants I and IV of the Cartesian plane-identifiable as areas of low environmental harshness and glacial influencehosted the highest number of taxa; corroborating the affinity of most of the taxa we retrieved with low values of the Glacial Index, longer distance from the glacier snout, and higher temperature and channel stability (e.g. Castella et al., 2001;Milner et al., 2001;Füreder et al., 2005;Brown et al., 2006). Among the seven most abundant taxa, only chironomids and Rhabdiopteryx sp. showed different behaviours, both related to higher values of the Glacial Index. For Chironomidae, an affinity to high-level glacierised sites is interpreted as a consolidated behaviour (e.g. Niedrist & Füreder, 2016), whereas for Rhabdiopteryx sp., it may be interpreted as a consequence of the larval life cycle with an early emergence, after growth under the snow cover : a condition that at the time of sampling was only satisfied by the most upstream sites, associated with higher glacial influence.
In temporal terms, our results fit well with conceptual models and previous studies : 1844Brown et al., 2007b). At site 3, where discharge was measured, a clear inverse relationship between magnitude of discharge and both the abundance and number of taxa was evident for all the sampling events (see Appendix-Supplementary Material).
Finally, the divergent status of site 2A may be explained by the presence of a tiny groundwater-fed tributary close to the sampling site, whose influence represented a local discontinuity resulting in an ameliorated and more diverse habitat condition for the benthic macroinvertebrate assemblage (Fig. 3, Brown et al., 2007c;Khamis et al., 2016).
In summary, the first of the three initial hypotheses was upheld. Indeed, faunal metrics used as a surrogate to qualitatively and quantitatively describe the benthic macroinvertebrate assemblage showed their highest values at sites with reduced glacial influence.
Regarding the second hypothesis, the site-to-site variation was found to be greater than the month-tomonth variation, meaning that the spatial component had a greater effect on the benthic macroinvertebrate assemblage than did the temporal component. Based also on the fact that four out of six sites comprised a stretch whose length was * 500 m, we were expecting the opposite. Moreover, the abiotic factors related to glacier-fed streams, even within a single glacialmelting season, vary within a wide range of values. Thus, these results emphasise the fact that even at a reduced spatial scale, benthic community assemblages may be quite different in terms of both density and number of taxa constituting the community.
Regarding the third and last hypothesis, the Saldur stream showed a set of abiotic factors whose dynamics and relationships were highly consistent with what was previously recorded in other European glacier-fed rivers. Consequently, the distribution of the benthic macroinvertebrates conformed well to the spatial and temporal conceptual models of glacier-fed streams. Nevertheless, our study was limited by being focused on a longitudinal transect not including the stream reach immediately adjacent to the glacier snout. For this reason, our results are only partly comparable to studies previously conducted on other glacier-fed streams.

Conclusions
Overall, this study provides one of the few characterisations of Italian glacier-fed streams (but see Maiolini & Lencioni, 2001;Lencioni et al., 2007;Niedrist et al., 2017) and represents a good reference point for future and continuous monitoring of the Saldur stream, also in light of the potential environmental impact of the ''run-of-river'' hydropower plant and of the regular ILTER activities. Similar analyses may thus constitute the baseline from which natural variability and the variability produced by climate change and/or human activities can be discerned.
Lastly, the study identified an important issue regarding monitoring programmes: even in an extremely temporally variable ecosystem such as a glacierfed stream (where sampling year and season still constitute a great source of faunal variability), the number and position of sampling sites along a longitudinal transect might have a prominent influence when assessing the variation of density and diversity of the benthic macroinvertebrates (e.g. Brown et al., 2009). Because benthic macroinvertebrates are commonly used as powerful bioindicators, this fact must be considered when scheduling long-term monitoring programmes.