Spatial and Temporal Patterns of Macroinvertebrate Assemblages in the River Po Catchment (Northern Italy)

Abstract

In the last decade, large scale biomonitoring programs have been implemented to obtain a robust understanding of freshwater in the name of helping to inform and develop effective restoration and management plans. A comprehensive biomonitoring dataset on the macroinvertebrate assemblages inhabiting the rivers of the Po Valley (northern Italy), comprised a total of 6762 sampling events (period 2007–2018), was analyzed in this study in order to examine coarse spatial and temporal trends displayed by biotic communities. Our results showed that macroinvertebrate compositions and derived structural and functional metrics were controlled by multiple environmental drivers, including altitude and climate (large scale), as well as habitat characteristics (local scale). Altitude proved to be the primary geographic driver, likely due to its association with thermal and precipitation regimes, thus explaining its overriding influence on macroinvertebrate assemblages. Significant temporal variations were observed across the study period, but notably in 2017, the overall taxonomic richness and diversity increased at the expense of Ephemeroptera, Plectoptera and Trichoptera taxa during an unprecedented heatwave that occurred across southern Europe. The detail of this study dataset allowed for important environmental attributes (e.g., altitude, habitat characteristics) shaping biotic communities to be identified, along with ecologically vulnerable regions and time periods (e.g., extreme climatic events). Such research is required globally to help inform large-scale management and restoration efforts that are sustainable over long-term periods.

1. Introduction

Assessing the ecological status of riverine ecosystems is fundamental to informing the conservation and management of freshwater ecosystems [1]. A large amount of data has been generated through biomonitoring efforts globally, including in the USA, Australia, South Africa and various nations within Europe [2]. Such data often span decadal time periods and cover broad spatial scales, making them ideal to uncover spatio-temporal effects of environmental drivers on biological communities. Examining data generated from biomonitoring programs is thus crucial to guide river management strategies and achieve ecological improvements through identifying primary environmental drivers shaping riverine ecosystems, as well as ecologically important regions and time periods. Moreover, such information can be useful in informing water managers and environmental regulators on how to enhance efficiency of the existing network of biomonitoring sites (e.g., biologically important locations).

Rivers serve many societal functions and are one of the most intensively human influenced ecosystems worldwide [3]. The benefits of water provision for various human purposes, including agriculture, energy production, drinking water supply, and the displacement of pollutant loads often impair riverine ecosystems and biodiversity [4], with potentially serious societal costs. Anthropogenic activities threaten riverine ecosystems through habitat loss and degradation [5] such as physical modification of riverine environments, deforestation of pristine wildernesses, pollution and introduction of non-native species [6,7,8]. Riverine environments, due to their societal and ecological importance, are integral to environmental policies such as the Water Framework Directive (WFD—[9]), and benthic invertebrate fauna (macroinvertebrates) are one of the predominant Biological Quality Elements [9] used for their ecological assessment. Macroinvertebrates are intermediately positioned in freshwater food webs, consuming basal resources like primary producers, organic matter and detritus, and serving as a food source for amphibians, birds, reptiles, fish and humans.

Both biodiversity assessments and biomonitoring programs utilizing macroinvertebrate assemblages routinely process and examine single metrics targeting different facets of community assemblages, such as taxonomic richness (e.g., [10]); taxonomic diversity (e.g., [11,12]); ecological response guild summaries (e.g., [13,14,15]); community abundance (e.g., [16]) and functional diversity (e.g., [17]). However, metrics characterizing different community properties respond differently to environmental stressors. Thus, single-metric approaches can miss important information about the assemblage structure and its response to environmental gradients [18], highlighting the importance of adopting multi-metric approaches to increase the likelihood of detecting ecological responses to different environmental conditions and stressors. In summary, variations in environmental conditions driven by natural and anthropogenic processes can mediate its influence on different structural and functional metrics. As such, other multivariate approaches (e.g., Non-Metric Multidimensional Scaling—NMDS) may thus be used along with metrics to better represent assemblage structure variations.

Despite the large amount of data collected on Italian river macroinvertebrates since the implementation of WFD programs, only recently has such data become available for researchers, policy makers and water managers (e.g., [19,20,21]). Fornaroli and coauthors [21] collected and homogenized data on the presence, distribution and abundances of macroinvertebrate taxa inhabiting the River Po catchment (northern Italy) in the last decade, providing the first checklist of macroinvertebrates occurring in this area. This data source represents a comprehensive spatial and temporal data set on the macroinvertebrate assemblages of the rivers of the Po Valley to help overcome the lack of analyses in the region and to help guide effective restoration and management efforts.

The objective of this study is to examine spatial and temporal trends of the macroinvertebrate assemblages inhabiting the River Po catchment. We tested the following hypotheses: (H1) macroinvertebrate assemblages show a geographic pattern that is primarily driven by altitude; (H2) there are differences among the assemblages inhabiting the pristine and anthropogenically altered study sites; (H3) there are families associated with the human alteration status; (H4) macroinvertebrate response metrics display significant temporal trends. 2. Materials and Methods 2.1. Study Area and Samples Selection

For this study we selected suitable data from a recently released dataset [21] collating the data collected by the Environmental Agencies using the Italian national standardized method [22] for the implementation of WFD activities. The underpinning data were collected between 2007 and 2018 and comprised a total of 6762 sampling events. All study rivers lie within the River Po catchment (North Italy), which covers 71,000 km² and crosses different Administrative Regions (Aosta Valley, Piedmont, Lombardy, Liguria, Veneto, Trentino, and Emilia-Romagna) including mostly the Subalpine area and the Po Plain, (Biogeographic region: Alpine, Continental and Mediterranean, [23]). Following the national standardised monitoring protocol [22], three samples were collected each year with a Surber-net enclosing an area of 0.05 or 0.10 m², with a mesh size of 0.5 mm. Each sample consisted of ten replicates collected proportionally to reflect predominance of the most representative microhabitats (mix of substratum type, current velocity, presence/absence of macrophytes, etc.). The different sampling efforts (0.5 or 1 m²—i.e., all replicates of each Surber-net summed) were standardized by reporting all the abundances to 1 m². The protocol includes sampling each site in at least two different seasons (the most related to the local life cycle of insects) [11]; consequently, the winter season is underrepresented in this dataset and excluded from the analyses.

To allow for the evaluation of temporal trends, only sites sampled for at least three years and with at least two sampled seasons annually were retained. Among the selected sites, none belong to the Emilia-Romagna administrative region because they do not fulfill the minimum requirements of replicates that we imposed for temporal trend analyses, as samples collected before 2013 are not present in the original database. The 270 remaining sites (Figure 1) had an average (median) of 11 samples, collected over an average period of 7 years. The considered rivers are natural, artificial (channel) or partially modified by different forms of infrastructure (e.g., dams, weirs) constructed along the rivers. For the selected samples sites, we extracted percentage representation of sampled microhabitats, classified in nine mineral substrate classes (e.g., cobble and gravel) and eight biological substrates (e.g., emerging macrophytes and coarse particulate organic matter) as included in the national standardized method [22], the presence of human impacts (i.e., each site was a priori classified as pristine or altered), the WGS-84 coordinates (as decimal degrees) and altitude (as m a.s.l.). The altitude of the sampling sites ranges from 29 to 2280 m a.s.l.

Identification of organisms was performed at different taxonomic rank levels by operators of the Environmental Agencies, thus, we homogenized the taxonomic level to the least common denominator (mostly family level) using the “biomonitoR” package [24] within R software [25].

2.2. GIS Analyses

We extracted several climatic variables from the WorldClim database (version 1.4) at the highest resolution (30 arc-s) using the function getData from the “raster” package [26] for each study site. These variables represented mean monthly and annual values of temperature and precipitation. For subsequent analyses, we considered 18 bioclimatic variables derived from the monthly temperature and rainfall values in order to generate more biologically meaningful variables. The bioclimatic variables represent annual trends (e.g., mean annual temperature, annual precipitation) seasonality (e.g., annual range in temperature and precipitation) and extreme or limiting environmental factors (e.g., temperature of the coldest and warmest month, and precipitation of the wet and dry quarters).

2.3. Macroinvertebrate Metrics

We extracted the relative density (i.e., log₁₀(individual per m² + 1) of three macroinvertebrate families (Plecoptera: Nemouridae, Ephemeroptera: Heptageniidae and Trichoptera: Limnephilidae) that were selected based on these families being: (i) identified as driving assemblage differentiation along altitudinal gradient; (ii) widely used as good environmental indicators for lotic ecosystems (e.g., [13,14]); and (iii) being the more abundant family within the respective order during the study period. We also evaluated three assemblage metrics commonly used in the bioassessment of riverine environments: (i) the richness of all families; (ii) the percentage representation of Ephemeroptera, Plecoptera and Trichoptera (%EPT) individuals; and (iii) the Shannon diversity index (H’) [11]. Such metrics and indices represent different parameters of the benthic assemblages such as taxa abundance, sensitivity and diversity.

Based on the same list of macroinvertebrates, each family in the database was linked to its auto-ecological preferences downloaded from the freshwaterecology.info database [27,28]. Freshwaterecology.info is an online tool that unifies, standardizes and codifies more than 20,000 European freshwater organisms and their biological properties and ecological preferences: through this database, we obtained traits of more than 140 macroinvertebrate families. The nomenclature of functional traits is reported herein by their “grouping features” and “traits” (see [29]). Grouping features represent a functional trait category (e.g., “substrate preference” and “food”), while traits signify modalities residing within these (e.g., substrate preference—“gravel”, “macrophytes”; food—“fine detritus”,“ dead animal”). A total of 21 grouping features comprising 113 traits were utilized from the functional traits database in subsequent analyses (Table S1). We derived functional traits from Tachet and coauthors [30] and, subsequently, aggregated them in two functional diversity indices.

Firstly, the Functional Richness (FRic) was derived, which represents the amount of functional space filled by the community [31]. FRic is defined by the trait extremes and thus reflects the potential maximum functional dissimilarity. FRic is calculated as the hypervolume enclosing the functional space filled by the community, the resulting FRic variable is standardized by its maximum, ranging from 0 to 1. Secondly, Rao’s quadratic entropy [32] was used to estimate Functional Diversity (FDiv) because it has been considered more appropriate than other indices [33]. This index is estimated using abundance data and standardized by the maximum value to constrain the values within the range 0–1.

2.4. Data Analyses

We used non-Metric Multidimensional Scaling (NMDS), a gradient analysis approach based on a distance or dissimilarity matrix, to visualize differences in taxonomic structure of the macroinvertebrate assemblage [34]. NMDS is an iterative procedure and was performed using the function metaMDS from the “vegan” package [35]. For this, we used the Bray-Curtis dissimilarity distance, a non-Euclidean distance used to quantify the compositional dissimilarity between two different samples, and log₁₀(x + 1) transformed abundance in each sample. To examine the gradient of effect of the geographic variables, we fitted smooth surface splines using the ordisurf function from the “vegan” package. This procedure uses generalized additive models (GAMs) to overlay a smoothed response surface, which allows a more detailed interpretation than a simple linear vector [36].

To assess the correlation between macroinvertebrate assemblages and environmental variables, we ran variance partitioning [37,38] for the macroinvertebrate assemblage matrix and the metrics, using the varpart function in “vegan” while the significance of fractions was tested using the dbrda function. This method partitions the variation between the pure effects of each variable, or group of variables and the shared variance explained. First, we considered the three groups of variables (i.e., geographic, climatic or habitat) together, and afterwards we ran separate analyses for each group. The pure effect of each geographic variable (i.e., altitude, latitude or longitude) was considered, while for habitats the variables were divided in two groups (distribution in classes of mineral and biological substrates abundance) as well as for climatic variables (temperature and precipitation related variables).

Differences in assemblage composition among alteration classes were quantitatively explored, as well as seasonal controls, with the additive effects of “sampling season” and “alteration status” being tested within a permutational multivariate analysis of variance (PERMANOVA) via the adonis function within the “vegan” package. Indicator species analysis (IndVal.g; [39]) was employed to search for significant indicator species discriminating the pristine versus altered study sites. Raw species relative density data were used within the function multipatt function from the “indicspecies” package [39] that incorporated a correction for unequal group sizes. IndVal combines the information on the concentration of species abundances in a particular group (A) and the degree of occurrence in that particular group (B). Thus, ideal indicator species are those that are always present at sites in a given group and never occur in other groups [40]. The indicator values range from 0–1, the value presented in this study corresponding to the square root of the product between A and B. The significance of the indicator values was tested by 999 Monte Carlo permutations where the observed indicator value was tested against those derived from randomized data, alpha was set to 0.001.

Temporal trends of different macroinvertebrate metrics (three taxonomic, two functional and three selected taxa) were explored separately for each season as their differing taxonomic compositions (driven by organism’s life-cycles) can respond differently to long-term changes in environmental conditions (e.g., [41]). Invertebrate metrics were modeled as a smoothed response over time through GAMMs via the gamm function in the “mgcv” package [42]. For this, cubic splines were implemented and a Gaussian distribution was modeled. Each study site was input as a random effect to account for communities from each sampling location potentially being correlated over time [8] and a first-order autoregressive process was adopted to account for temporal autocorrelation (correlated errors within years were accounted for); these criteria yielded the optimal model compared to respective models (GAMs or GAMMs) without a random effect or a first-order autoregressive process in all instances (indicated by an Akaike Information Criterion value of ≤−2). Then, residual diagnostics were inspected for each optimal GAMM to ensure the assumptions of homogenous variances and normal distributions were met; where these criteria were satisfied, the significance of the time smoothing parameter was obtained. In addition, periods where invertebrate responses significantly increased or decreased were explored by the first derivatives of each GAMM using the method of finite differences derived from 200 equally spaced time points over the study period [43]. When residual checks revealed the assumptions of GAMMs were not satisfied, temporal changes were instead visualised through a “LOcally Estimated Scatterplot Smoothing” (LOESS) function.

We performed all statistical analyses using R project software [25].

3. Results 3.1. Similarity among Macroinvertebrate Assemblages

The NMDS plot showed that macroinvertebrate assemblages inhabiting the River Po catchment displayed a strong geographic pattern (Figure 2). The main influence of longitude and altitude were observed along the NMDS axis 1 (Figure 2a,b) while the effect of latitude was more important along NMDS axis 2 (Figure 2c). Sites positively loaded on axis 1 tended to be characteristic of lowland sites (<400 m a.s.l.), while sites negatively loaded were typical of high-altitude sites (>1800 m a.s.l). We identified a clear altitudinal pattern in the position of the macroinvertebrate families in the ordination space (from right to left and from bottom to top Figure 2c). Assemblages inhabiting high-altitude sites were characterized mostly by families (with average NMDS1 score < −0.5) belonging to the Plecoptera order (i.e., Taeniopterygidae, Chloroperlidae, Perlodidae, Perlidae and Nemouridae), but also some Trichoptera (i.e., Glossosomatidae, Philopotamidae and Limnephilidae), Coleoptera (i.e., Scirtidae and Hydraenidae), Diptera (i.e., Blephariceridae, Thaumaleidae, Pediciidae and Psychodidae), Oligochaeta (i.e., Haplotaxidae) and Platyhelminthes (i.e., Planariidae). Macroinvertebrates inhabiting lowland sites were characterized mostly by non-insect families (with average NMDS1 score > 1.2) (e.g., Unionidae, Acroloxidae, Salifidae, Haemopidae, Corbiculidae and Hydridae), although some exceptions to this were observed with certain insect families (e.g., Sciomyzidae, Coenagrionidae and Hydrometridae). The effects of longitude and latitude on assemblage structure are less representative but coherent with the geographic characteristics of the studied catchment and partial covariates with the observed pattern of altitude.

In the variance partitioning framework, when examining climatic, habitat and geographical variables simultaneously, climatic variables alone explained a much higher proportion of the variability than habitat or geographic variables. The total variance explained was 32.2%, with the pure effect of the climate accounting for 9.4% and, combined with the shared influence of habitat and geographic variables, this increased to 28.0% (Figure 3a). The pure effect of the habitat variables was 3.1% and, combined with the shared influence of climatic and geographic variables, this increased to 15.1%. Similarly, the pure effect of geographic variables was 1.0% and, combined with the shared influence of habitat and climatic variables, this increased to 15.1%. When examining variable types (climate, habitat and geographic) separately, the variance explained by climatic variables was mostly due to the shared influence among temperature and precipitation related variables (16.4% Figure 3b). The total variance explained by habitat was 15.1%, with the mineral substrate yielding the greatest statistical influence (pure effect 11.2%, Figure 3c). Similarly, geographic variables explained 15.1% of the statistical variation whereby altitude proved to be the most influential variable (pure effect = 8.1%, pure plus shared effect = 14.0%, Figure 3d).

3.2. Presence of Human Impacts

PERMANOVA highlighted significant differences among the assemblages inhabiting pristine sites versus those subjected to some form of human alteration (F = 167.97, R² = 0.05, p < 0.001). The two clusters cover most of the altitudinal and latitudinal ranges of the River Po catchment, while the pristine sites were better represented in its western part (Figure S1).

Indicator species analysis highlighted a range of macroinvertebrate families that differed significantly between the two alteration classes, with many insect families associated with pristine sites and non-insect families associated with altered sites (Table 1). Twenty-three macroinvertebrate families, spanning across several taxonomic orders, displayed greater affinities for pristine sites, while 12 families showed greater affinity for the sites impacted by anthropogenic environmental alterations (Table 1). For all the selected families, the A term was bigger than the B term, suggesting these taxa are good indicators of the altered group, but are present only occasionally.

3.3. Metrics and Indices

The density of macroinvertebrates ranged between 12 and 36,820 individuals *m⁻² (1673 ± 2260, mean ± standard deviation) in the whole study area, whereas family richness varies between 3 and 44 (19 ± 7). The macroinvertebrate families selected on the basis of their wide geographic distribution and relevance to bioassessment (i.e., Nemouridae, Heptageniidae and Limnephilidae) were often present in the considered samples (48%, 77% and 42%, respectively) and their density ranged between 0–830, 0–1660 and 0–6760 (respectively) per square meter. The percentage of Ephemeroptera, Plecoptera and Trichoptera individuals (%EPT) within a sample was 51 ± 25%, while the Shannon diversity index (H’) ranged between 0.01 to 2.96. The Pearson’s correlation between metrics and indices was significant in the majority of instances (Figure 4), but the only two pairs of highly correlated indices (r > 0.70) were taxonomic richness versus functional richness (r = 0.76, p < 0.001) and H’ versus functional diversity (r = 0.80, p < 0.001).

In the variance partitioning framework, the relative density of the three selected macroinvertebrate families (Nemouridae, Heptageniidae and Limnephilidae) and functional diversity were less influenced by the geographic position, climate and habitat than the multivariate assemblage structure and other univariate metrics (Table 2). The total variance explained for structural indices and functional richness ranged between 27 and 31% and the influence of climatic variables was always the greatest (24–28%), as observed also for the multivariate assemblage structure. The variance explained by climatic variables was mostly due to the shared influence with habitat and geographic position than to its pure effect for all metrics (Figure S2).

3.4. Temporal Trends of Metrics and Indices

Using GAMMs, it was possible to assess the significance of the temporal trends for four out of eight tested metrics and all proved to be significant (p < 0.05) for all the seasons with coherent patterns. Over the first few years of monitoring (c. 2008–2011), family richness, Shannon and functional diversity declined significantly (Figure 5a–c,e). The taxonomic richness and Shannon notably increased between c. 2011–2017, decreased significantly thereafter (Figure 5a–c,g–i); these patterns were broadly mirrored by the functional diversity response, although fewer significant changes were detected for this metric (Figure 5j–l). The %EPT broadly contradicted the temporal trends observed for the other metrics, although they largely displayed a lack of significant temporal variations. One exception to this was a significant increase in values in the fall months after 2017, whereby an abrupt decline (c. 10–15%) occurred across the two years prior and then completely recovered by the end of 2018.

4. Discussion

Findings from this study highlighted that the structure of macroinvertebrate assemblages in the River Po catchment (Northern Italy) were controlled by multiple drivers acting at both coarse (e.g., altitude and climate) [44] and local spatial scales (e.g., habitat characteristics) [45]. The variance partitioning analyses highlighted that the variation in macroinvertebrate assemblages was mainly explained by the shared influence of the three groups of environmental variables and indicated that climate defined most of the variance explained, followed by habitat and then geographic position. Altitude has a well-documented control of flow and thermal regimes, which are widely recognized as primary determinants shaping macroinvertebrate assemblage structures (e.g., [46,47,48]). Mineral substrate composition showed greater association with macroinvertebrate than biological substrates, which is likely due to the latter being almost completely absent in sites above 500 m a.s.l. [21], while coarse substrates are commonly associated with steeper channel slope, and hence occur more widely at higher altitudes. Our findings highlighted that altitude was the main geographic driver within the study area, as reported in many other studies (e.g., [49,50]), thus confirming our H1 hypothesis.

The results of NMDS analyses showed that the macroinvertebrate assemblages of the present study were strongly driven by geographic variables and particularly by altitude, confirming again our H₁ hypothesis. High altitude assemblages, as well as those inhabiting pristine sites, were dominated by Plecoptera families. The association between Plecoptera and altitude is well known (e.g., [51,52]), as is their vulnerability to climatic changes through their restricted altitudinal niche, which has been previously highlighted [53] and represents a new threat for this order in the study region (alongside pollution) [54]. Families belonging to other orders such as Trichoptera, Coleoptera, Diptera, Oligochaeta and Platyhelminthes showed similar altitudinal gradients (i.e., strong association with high-altitude sites). These results are due, at least partially, to their biogeographic distribution and not the effect of human alterations. Among the families representative of lowland sites, some are commonly associated with the potamal stretches (longitudinal distribution according to [30]) of rivers and streams (e.g., Unionidae, Corbiculidae), or with lentic waters (e.g., Hydrometridae, Acroloxidae), whereas others are more habitat generalist (e.g., Salifidae, Haemopidae). While some of these are considered indicators of good environmental status (e.g., Unionidae and Hydrometridae), Corbiculidae is considered an invasive family in Italy and was associated with altered sites in this work. In summary, high altitude assemblages were dominated by insects, while lowland assemblages were dominated by non-insects, and were more susceptible to biological invasions [55].

According to PERMANOVA, significant differences occurred among the assemblages inhabiting the pristine versus altered sites, confirming our H₂ hypothesis. Moreover, IndVal analysis identified various families associated with the different alteration status, supporting our H₃ hypothesis. The macroinvertebrate families associated with pristine sites were typically indicators of good water quality, as highlighted within Biological Monitoring Working Party (BMWP—[13]). Conversely, most families associated with altered sites are generally classified indicators of impaired water quality. Such findings agree with the extensive literature available about the tolerance of macroinvertebrate families to water quality impairment (e.g., [56,57,58]); nonetheless, many families commonly classified as sensitive to nutrient enrichment (e.g., Leuctridae, Heptageniidae and Leptoceridae) did not display significant differences between their abundance in pristine and altered sites. In addition, many of the families associated with pristine sites are classified as rheophilic within the Lotic-invertebrate Index for Flow Evaluation metric (LIFE—[14]), which is widely used to define hydroecological relationships. Conversely, many of the families associated with altered sites possess LIFE scores that indicate they are reflective of slower flow conditions. Such ecological trends are critical within biomonitoring assessments, and as such further explorations of the relationships among selected pressures and the responses of biological indices required to help guide specific management actions [59].

The results of the variance partitioning analyses applied to structural indices, functional richness and those of multivariate assemblage structure were broadly comparable; while the absolute abundance of Nemouridae, Heptageniidae and Limnephilidae and the functional diversity were less responsive to the considered environmental drivers. We processed a limited number of metrics which responded differently to different stressors. Given that individually, these metrics only characterize a small part of the assemblage, they may lack sensitivity to changes in certain environmental conditions or when multiple anthropogenic stressors are present [60], which may have resulted in incongruent responses and certain ecological trends not being detected [61]. Future studies that focus on the effect of specific anthropogenic modifications in this catchment (as well as studies using biomonitoring data elsewhere) could build upon this study and focus on macroinvertebrate indices and metrics developed to assess a species stressor (e.g., BMWP for nutrient enrichment and LIFE for hydrological alterations). Notwithstanding, univariate indicators used in this study consider complete macroinvertebrate assemblages which reliably characterize assemblage structure and are capable of capturing ecological responses to environmental gradients. The relative density of selected families are likely controlled both by stochasticity and by natural and anthropogenic drivers that act at smaller temporal and spatial scales than considered in this research, such as short-term climate variations and/or local pollution events [62].

In accordance with our last hypothesis (H₄) we identified significant temporal trends for some of the studied metrics. Inter-annual variations in macroinvertebrate assemblages can be due to environmental changes such as climatic variations driving river flow and thermal regimes [63], changes in anthropogenic pressures (e.g., river restoration), to the establishment of invasive species [8], or to the interactive effects of those drivers [41]. The first few years of the study period witnessed evident declines in various response metrics; namely, taxonomic richness and diversity, during a time when discharges were increasing following a major drought in 2007 [64]. Comparably, significant increases in the same metrics were observed in 2017, when an unprecedent heatwave occurred across southern Europe that triggered a severe drought in many areas, persisting into fall [65]. These significant temporal trends indicate that macroinvertebrate assemblages inhabiting the Po catchment thrive during low-flow periods, likely due to many taxa within the regions displaying strong tolerance mechanisms to drought in a region known for its harsh, dry summers. Despite this, Ephemeroptera, Plecoptera and Trichoptera taxa displayed the opposite trend, most notably during autumn months. Such findings are unsurprising, as at the order level these taxa are negatively associated with low-flow periods [66], while autumn months depicts a critical time for the larval recruitment development. As such, significant temporal variations displayed in the Po catchment are most likely attributed to inter-annual flow and thermal regime variations. It is unlikely that temporal trends within this large catchment were driven by modifications of the anthropogenic pressures because those would have varied spatially across the Po catchment and generally do not present abrupt changes like the ones that we observed for the macroinvertebrate assemblages [67].

5. Conclusions Final Remarks Our results provide evidence of the link between local climate and macroinvertebrate assemblages in the catchment of River Po (North Italy). The detail of this study dataset allowed for important environmental attributes (e.g., altitude, habitat characteristics) shaping biotic communities to be identified, along with ecologically important and vulnerable regions and time periods (e.g., extreme climatic events). Such research is required globally to help inform large-scale management and restoration efforts that are sustainable over long-term periods. While this study represents the first of its kind across the region, further studies need to focus both on the effects of inter-annual climate variability on macroinvertebrate assemblages to further explore temporal trends described here, as well as to cultivate a better understanding of primary stressors shaping biotic assemblages through the use of specifically developed metrics (e.g., nutrient enrichment—BMWP; flow—LIFE,). To achieve a better understanding of macroinvertebrate assemblage dynamics within this catchment, the biomonitoring network can be improved by focusing on “reference” sites (largely free from human interference), although this study highlighted the difficulty finding such conditions in lowland environments. Further analyses should focus on the response of the families which we have identified representative of pristine sites in this study; in medium-high altitude sites in particular, such research could unveil differences in ecological dynamics between pristine and impacted sites. Moreover, from a management perspective, there remains a global need to incorporate measures of different anthropogenic impacts within (or alongside) ecological databases, as this would allow key drivers and stressors to be identified and integrated into effective management strategies.

Figure 1. Map of the study area with the location of sampling sites.

Figure 2. Nonmetric multidimensional scaling (NMDS) ordination plot for aquatic macroinvertebrate assemblages where the a-priori identified alteration status are colored. Dark gray denotes the assemblages that belong to pristine sites (n = 534) while light gray denotes those assemblages that belong to sites where an impact is present (n = 2563). Shaded ellipses represent the 95% confidence interval surrounding the centroid of each impact class in the ordination space. Overlaid smooth surface display blue splines using generalized additive models (GAMs) for Altitude (a), Longitude (b) and Latitude (c). Numbers on the splines represent respectively meter above sea level and WGS84 coordinates. Macroinvertebrate families are positioned in the ordination space with uppercase labels, as weighted averages. 3D stress = 0.16.

Figure 3. Results of variance partitioning on macroinvertebrate assemblage matrix. Panel (a) shows the results for the three considered groups of variables together while other panels show the detailed results for each group. Panel (b) refers to climatic variables, panel (c) to habitat variables and panel (d) to geographic ones. Values displayed are the adjusted R2 and negative values are not shown. The unexplained portion is shown in the bottom right of panels (Residuals). Significance of the pure effects are show as asterisks. *** p < 0.001.

Figure 4. Correlation among macroinvertebrate communities metrics and indices. The diagonal represents the density function of each metric or index, below the diagonal the pairs plot are presented and above the diagonal the Pearson's correlation between each pair of metrics and indices is reported along with the significance levels. Levels of statistical significance are as follows: ***, p < 0.001; **, p < 0.01; ns, p > 0.05.

Figure 5. Temporal trajectories of (a-c) family richness, (d-f) %EPT, (g-i) Shannon diversity index (H') and (j-l) functional diversity as a function of the generalized additive mixed-effect models (GAMMs) outputs for the different seasons. The significance of the time smoothing parameters is reported within each panel. Significant increases or decreases in slope are represented respectively in blue and red colors. Corresponding graphs, only to visualize temporal trends via LOESS (Locally Estimated Scatterplot Smoothing) for the other metrics and indices are shown in Figure S3.

Taxa Group	Family	A	B	IndVal
Pristine Sites
Oligochaeta	Haplotaxidae	0.76	0.90	0.82
Plecoptera	Chloroperlidae	0.80	0.70	0.75
Trichoptera	Polycentropodidae	0.76	0.71	0.74
Diptera	Dixidae	0.71	0.73	0.72
Ephemeroptera	Leptophlebiidae	0.80	0.59	0.69
Plecoptera	Taeniopterygidae	0.76	0.52	0.63
Oligochaeta	Enchytraeidae	0.78	0.49	0.62
Diptera	Psychodidae	0.76	0.46	0.59
Diptera	Rhagionidae	0.70	0.48	0.58
Plecoptera	Perlodidae	0.86	0.32	0.52
Trichoptera	Philopotamidae	0.82	0.31	0.51
Trichoptera	Glossosomatidae	0.75	0.32	0.49
Trichoptera	Beraeidae	0.71	0.28	0.44
Coleoptera	Scirtidae	0.80	0.20	0.40
Trichoptera	Odontoceridae	0.79	0.20	0.40
Diptera	Thaumaleidae	0.69	0.19	0.36
Plecoptera	Perlidae	0.83	0.16	0.36
Diptera	Blephariceridae	0.81	0.16	0.36
Platyhelminthes	Planariidae	0.83	0.15	0.35
Diptera	Pediciidae	0.72	0.16	0.34
Diptera	Athericidae	0.84	0.09	0.28
Coleoptera	Hydraenidae	0.90	0.04	0.19
Plecoptera	Nemouridae	0.96	0.02	0.15
Altered sites
Oligochaeta	Tubificidae	0.80	0.70	0.59
Hirudinea	Erpobdellidae	0.76	0.71	0.55
Crustacea	Gammaridae	0.71	0.73	0.51
Crustacea	Asellidae	0.80	0.59	0.41
Gastropoda	Physidae	0.76	0.52	0.37
Odonata	Calopterygidae	0.78	0.49	0.30
Gastropoda	Lymnaeidae	0.76	0.46	0.30
Hirudinea	Glossiphoniidae	0.70	0.48	0.27
Bivalvia	Corbiculidae	0.86	0.32	0.27
Heteroptera	Aphelocheiridae	0.82	0.31	0.24
Gastropoda	Hydrobiidae	0.75	0.32	0.23
Gastropoda	Planorbidae	0.71	0.28	0.20

Metric/Index	Total Explained Variance	Geographic	Climatic	Habitat
Assemblage distance matrix	32.2	15.1	28.0	15.5
Nemouridae density	16.2	13.2	14.1	5.4
Heptageniidae density	10.3	4.9	8.4	3.6
Limnephilidae density	6.0	2.8	5.9	0.1
Family richness	31.3	8.6	27.6	11.0
EPT proportion	30.2	13.8	25.7	19.8
Shannon	27.6	7.6	24.6	11.3
Functional Richness	26.7	13.9	24.6	6.6
Functional Diversity	11.3	4.4	9.4	2.5

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/12/9/2452/s1, Table S1: Macroinvertebrate functional traits examined within this study, with biological traits in non-italicized text and ecological traits being italicized, Figure S1: Distribution of pristine and altered sites along the geographic gradients, Figure S2: Results of variance partitioning on macroinvertebrate metrics and indices. Values displayed are the adjusted R2 and negative values are not shown. The unexplained portion is shown in the bottom right of panels (Residuals), Figure S3: Temporal trajectories of (a-c) log10 transformed Nemouridae absolute abundance, (d-f) log10 transformed Heptageniidae absolute abundance, (g-i) log10 transformed Limnephilidae absolute abundance and (j-l) Functional Richness visualized using LOESS (Locally Estimated Scatterplot Smoothing) for the different seasons.

Author Contributions

Conceptualization, R.F., A.B. and A.L.; Methodology, R.F. and J.C.W.; Software, R.F., J.C.W. and A.L.; Formal Analysis, R.F. and J.C.W.; Data Curation, R.F.; Writing-Original Draft Preparation, R.F., A.B., J.C.W. and A.L.; Visualization, R.F. and J.C.W.; Supervision, A.B. and A.L.; Writing-review and editing, R.F., A.B., J.C.W. and A.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The Authors would like to thank all technicians and operators of the environmental agencies for the sharing of data collected during water quality monitoring programs over the past decade, without which this work would not have been possible. We also thank the two anonymous reviewers for their constructive comments on an earlier version of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Vugteveen, P.; Leuven, R.S.E.W.; Huijbregts, M.A.J.; Lenders, H.J.R. Redefinition and Elaboration of River Ecosystem Health: Perspective for River Management. Hydrobiologia 2006, 565, 289-308.

2. Schmutz, S. Riverine Ecosystem Management; Springer International Publishing: New York, NY, USA, 2018; ISBN 9783319732497.

3. Wohl, E.; Lane, S.; Wilcox, A. The Science and Practice of River Restoration. Water Resour. Res. 2015, 51, 5974-5997.

4. Vörösmarty, C.J.; McIntyre, P.B.; Gessner, M.O.; Dudgeon, D.; Prusevich, A.; Green, P.; Glidden, S.; Bunn, S.E.; Sullivan, C.A.; Liermann, C.R.; et al. Global threats to human water security and river biodiversity. Nature 2010, 467, 555-561.

5. Allan, J.D.; Flecker, A.S. Biodiversity Conservation in Running Waters. Bioscience 1993, 43, 32-43.

6. Lewin, I.; Jusik, S.; Szoszkiewicz, K.; Czerniawska-Kusza, I.; Ławniczak, A.E. Application of the new multimetric MMI_PL index for biological water quality assessment in reference and human-impacted streams (Poland, the Slovak Republic). Limnologica 2014, 49, 42-51.

7. Nilsson, C.; Berggren, K. Alterations of riparian ecosystemes caused by river regulation. Bioscience 2000, 50, 783-792.

8. Mathers, K.L.; White, J.C.; Fornaroli, R.; Chadd, R. Flow regimes control the establishment of invasive crayfish and alter their effects on lotic macroinvertebrate communities. J. Appl. Ecol. 2020, 57, 886-902.

9. EC European Community. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 Establishing a Framework for Community Action in the Field of Water Policy; European Community: Brussels, Belgium, 2000; Volume L327/1.

10. Ofenböck, T.; Moog, O.; Gerritsen, J.; Barbour, M. A stressor specific multimetric approach for monitoring running waters in Austria using benthic macro-invertebrates. Hydrobiologia 2004, 516, 251-268.

11. Hering, D.; Moog, O.; Sandin, L.; Verdonschot, P.F.M. Overview and application of the AQEM assessment system. Hydrobiologia 2004, 516, 1-20.

12. Pinto, P.; Rosado, J.; Morais, M.; Antunes, I. Assessment methodology for southern siliceous basins in Portugal. Hydrobiologia 2004, 516, 191-214.

13. Armitage, P.D.; Moss, D.; Wright, J.F.; Furse, M.T. The performance of a new biological water quality score system based on macroinvertebrates over a wide range of unpolluted running-water. Water Res. 1983, 17, 333-347.

14. Extence, C.; Balbi, D.M.; Chadd, R.P. River flow indexing using British benthic macroinvertebrates: A framework for setting hydroecological objectives. Regul. Rivers Res. Manag. 1999, 15, 543-574.

15. Extence, C.; Chadd, R.P.; England, J.; Dunbar, M.J.; Wood, P.J.; Taylor, E.D. The assessment of fine sediment accuulation in rivers using macroinvertebrate community response. River Res. Appl. 2013, 29, 17-55.

16. Suren, A.M.; Jowett, I.G. Effects of floods versus low flows on invertebrates in a New Zealand gravel-bed river. Freshw. Biol. 2006, 51, 2207-2227.

17. Merritt, R.W.; Wallace, J.; Higgins, M.; Alexander, M.; Berg, M.; Morgan, W.; Cummins, K.; Vandeneeden, B. Procedures for the functional analysis of invertebrate communities of the kissimmee River-floodplain ecosystem. Fla. Sci. 1996, 59, 216-274.

18. Reynoldson, T.B.; Norris, R.H.; Resh, V.H.; Day, K.E.; Rosenberg, D.M. The reference condition: A comparison of multimetric and multivariate approaches to assess water-quality impairment using benthic macroinvertebrates. J. N. Am. Benthol. Soc. 1997, 16, 833-852.

19. Erba, S.; Cazzola, M.; Belfiore, C.; Buffagni, A. Macroinvertebrate metrics responses to morphological alteration in Italian rivers. Hydrobiologia 2020, 847, 2169-2191.

20. Calabrese, S.; Mezzanotte, V.; Marazzi, F.; Canobbio, S.; Fornaroli, R. The influence of multiple stressors on macroinvertebrate communities and ecosystem attributes in Northern Italy pre-Alpine rivers and streams. Ecol. Indic. 2020, 115, 106408.

21. Fornaroli, R.; Agostini, A.; Arnaud, E.; Berselli, A.; Bettoni, E.; Boggero, A.; Brolandelli, C.; Cadrobbi, G.; Cason, M.; Castelli, L.; et al. A ten-year geographic data set on the occurrence and abundance of macroinvertebrates in the River Po basin (Italy). Biogeographia 2020, (in press).

22. ISPRA Istituto Superiore per la Protezione e Ricerca Ambientale. Manuali e Linee Guida 107/2014. Report in Italian; ISPRA Linee guida per la valutazione della componente macrobentonica fluviale ai sensi del DM 260/2010; ISPRA: Roma, Italy, 2014; p. 88.

23. EEA European Environment Agency. Biogeographical Regions in Europe. Internet Database. Available online: https://www.eea.europa.eu/data-and-maps/figures/biogeographical-regions-in-europe-2 (accessed on 31 August 2020).

24. Laini, A.; Bolpagni, R.; Daniel, B.; Burgazzi, G.; Guareschi, S.; Cédric, M.; Tommaso, C. biomonitoR: An R package for calculating biomonitoring indices of running waters. XXVIII Congr. Ital. Ecol. Soc. 2018. Available online: https://www.biomonitor.it/ (accessed on 31 August 2020).

25. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019.

26. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling. R package Version 3.0-12. Available online: https://CRAN.R-project.org/package=raster (accessed on 31 August 2020).

27. Schmidt-Kloiber, A.; Bremerich, V.; De Wever, A.; Jähnig, S.C.; Martens, K.; Strackbein, J.; Tockner, K.; Hering, D. The Freshwater Information Platform: A global online network providing data, tools and resources for science and policy support. Hydrobiologia 2019, 838, 1-11.

28. Schmidt-Kloiber, A.; Hering, D. www.freshwaterecology.info-An online tool that unifies, standardises and codifies more than 20,000 European freshwater organisms and their ecological preferences. Ecol. Indic. 2015, 53, 271-282.

29. Schmera, D.; Podani, J.; Heino, J.; Erös, T.; Poff, N.L.R. A proposed unified terminology of species traits in stream ecology. Freshw. Sci. 2015, 34, 823-830.

30. Tachet, H.; Richoux, P.; Bournaud, M.; Usseglio-Polatera, P. Invertébrés d'eau Douce, 2nd ed.; CNRS Editions: Paris, France, 2010; ISBN 2-271-05745-0.

31. Villéger, S.; Mason, N.W.H.; Mouillot, D. New multidimensional functional diversity indices for a multifaceted framework in functional ecology. Ecology 2008, 89, 2290-2301.

32. Rao, R.C. Diversity and Dissimilarity Coefficients: A Unified Approach. Theor. Popul. Biol. 1982, 21, 24-43.

33. Botta-Dukát, Z. Rao's quadratic entropy as a measure of functional diversity based on multiple traits. J. Veg. Sci. 2005, 16, 533-540.

34. Clarke, K.R.; Ainsworth, M. A method of linking multivariate community structure to environmental variables. Mar. Ecol. Prog. Ser. 1993, 92, 205-219.

35. Oksanen, J.; Blanchet, F.G.; Friendly, M.; Kindt, R.; Legendre, P.; Mcglinn, D.; Minchin, P.R.; O'Hara, R.B.; Simpson, G.L.; Solymos, P.; et al. Vegan: Community Ecology Package; R package Version 2.5-6. Available online: https://CRAN.R-project.org/package=vegan (accessed on 31 August 2020).

36. Tonkin, J.D. Drivers of macroinvertebrate community structure in unmodified streams. PeerJ 2014, 2, e465.

37. Borcard, D.; Legendre, P.; Drapeau, P. Partialling out the Spatial Component of Ecological Variation. Ecology 1992, 73, 1045-1055.

38. Peres-Neto, P.R.; Legendre, P.; Dray, S.; Borcard, D. Variance partitioning of species data matrices: Estimation and comparison of fractions. Ecology 2006, 87, 2614-2625.

39. De Cáceres, M.; Legendre, P. Associations between species and groups of sites: Indices and statistical inference. Ecology 2009, 90, 3566-3574.

40. Dufrêne, M.; Legendre, P. Species assemblages and indicator species: The need for a flexible asymmetrical approach. Ecol. Monogr. 1997, 67, 345-366.

41. Durance, I.; Ormerod, S.J. Trends in water quality and discharge confound long-term warming effects on river macroinvertebrates. Freshw. Biol. 2009, 54, 388-405.

42. Wood, S.N. Package 'mgcv'. R package Version 1.8-33. Available online: https://CRAN.R-project.org/package=vegan (accessed on 31 August 2020).

43. Monteith, D.T.; Evans, C.D.; Henrys, P.A.; Simpson, G.L.; Malcolm, I.A. Trends in the hydrochemistry of acid-sensitive surface waters in the UK 1988-2008. Ecol. Indic. 2014, 37, 287-303.

44. Vaughan, I.P.; Ormerod, S.J. Large-scale, long-term trends in British river macroinvertebrates. Glob. Chang. Biol. 2012, 18, 2184-2194.

45. Lamouroux, N.; Dolédec, S.; Gayraud, S. Biological traits of stream macroinvertebrate communities: Effects of microhabitat, reach, and basin filters. J. N. Am. Benthol. Soc. 2004, 23, 449-466.

46. Jackson, H.M.; Gibbins, C.N.; Soulsby, C. Role of discharge and temperature variation in determining invertebrate community structure in a regulated river. River Res. Appl. 2007, 23, 651-669.

47. Chinnayakanahalli, K.J.; Hawkins, C.P.; Tarboton, D.G.; Hill, R.A. Natural flow regime, temperature and the composition and richness of invertebrate assemblages in streams of the western United States. Freshw. Biol. 2011, 56, 1248-1265.

48. White, J.C.; Hannah, D.M.; House, A.; Beatson, S.J.V.; Martin, A.; Wood, P.J. Macroinvertebrate responses to flow and stream temperature variability across regulated and non-regulated rivers. Ecohydrology 2017, 10, e1773.

49. Lewin, I.; Szoszkiewicz, K.; Jusik, S.; Ławniczak, A.E. Influence of selected environmental factors on macroinvertebrates in mountain streams. Open Life Sci. 2015, 10, 99-111.

50. Suren, A.M. Macroinvertebrate communities of streams in western Nepal: Effects of altitude and land use. Freshw. Biol. 1994, 32, 323-336.

51. Kamler, E. Distribution of Plecoptera and Ephemeroptera in relation to altitude above mean sea level and current speed in mountain waters. Pol. Arch. Hidrobiol. 1967, 2, 29-42.

52. Hynes, H.B.N. Biology of Plecoptera. Annu. Rev. Entomol. 1976, 21, 135-153.

53. De Figueroa, J.M.T.; López-Rodríguez, M.J.; Lorenz, A.; Graf, W.; Schmidt-Kloiber, A.; Hering, D. Vulnerable taxa of European Plecoptera (Insecta) in the context of climate change. Biodivers. Conserv. 2010, 19, 1269-1277.

54. Fochetti, R. Endemism in the Italian stonefly-fauna (Plecoptera). Zootaxa 2020, 4722, 381-388.

55. Fenoglio, S.; Bonada, N.; Guareschi, S.; López-Rodríguez, M.J.; Millán, A.; Tierno De Figueroa, J.M. Freshwater ecosystems and aquatic insects: A paradox in biological invasions. Biol. Lett. 2016, 12, 20151075.

56. Paisley, M.F.; Trigg, D.J.; Walley, W.J. Revision of the Biological Monitoring Working Party (BMWP) score system: Derivation of present-only and abundance-related scores from field data. River Res. Appl. 2014, 30, 887-904.

57. Wilding, N.A.; White, J.C.; Chadd, R.P.; House, A.; Wood, P.J. The influence of flow permanence and drying pattern on macroinvertebrate biomonitoring tools used in the assessment of riverine ecosystems. Ecol. Indic. 2018, 85, 548-555.

58. Guareschi, S.; Laini, A.; Sánchez-Montoya, M.M. How do low-abundance taxa affect river biomonitoring? Exploring the response of different macroinvertebrate-based indices. J. Limnol. 2017, 76, 9-20.

59. Clews, E.; Ormerod, S.J. Improving bio-diagnostic monitoring using simple combinations of standard biotic indices. River Res. Appl. 2009, 25, 348-361.

60. Laini, A.; Bolpagni, R.; Cancellario, T.; Guareschi, S.; Racchetti, E.; Viaroli, P. Testing the response of macroinvertebrate communities and biomonitoring indices under multiple stressors in a lowland regulated river. Ecol. Indic. 2018, 90, 47-53.

61. Yuan, L.L.; Norton, S.B. Comparing responses of macroinvertebrate metrics to increasing stress. J. N. Am. Benthol. Soc. 2003, 22, 308-322.

62. Robinson, C.T.; Minshall, G.W.; Royer, T.V. Inter-annual patterns in macroinvertebrate communities of wilderness streams in Idaho, U.S.A. Hydrobiologia 2000, 421, 187-198.

63. Monk, W.A.; Wood, P.J.; Hannah, D.M.; Wilson, D.A. Macroinvertebarte community response to inter-annual and regional river flow regime dynamics. River Res. Appl. 2008, 245, 988-1001.

64. Montanari, A. Hydrology of the Po River: Looking for changing patterns in river discharge. Hydrol. Earth Syst. Sci. 2012, 16, 3739-3747.

65. Kew, S.F.; Philip, S.Y.; Van Oldenborgh, G.J.; Otto, F.E.L.; Vautard, R.; Van Der Schrier, G. The exceptional summer heat wave in southern Europe 2017. Bull. Am. Meteorol. Soc. 2019, 100, S49-S53.

66. Datry, T.; Larned, S.T.; Fritz, K.M.; Bogan, M.T.; Wood, P.J.; Meyer, E.I.; Santos, A.N. Broad-scale patterns of invertebrate richness and community composition in temporary rivers: Effects of flow intermittence. Ecography 2014, 37, 94-104.

67. Randin, C.F.; Ashcroft, M.B.; Bolliger, J.; Cavender-Bares, J.; Coops, N.C.; Dullinger, S.; Dirnböck, T.; Eckert, S.; Ellis, E.; Fernández, N.; et al. Monitoring biodiversity in the Anthropocene using remote sensing in species distribution models. Remote Sens. Environ. 2020, 239, 111626.