Role of the Hyporheic Zone in Increasing the Resilience of Mountain Streams Facing Intermittency

Abstract

We investigated the impact of intermittence in previously-perennial Alpine stream reaches, targeting the role of the hyporheic zone in increasing the resilience of these aquatic systems. We selected a perennial and an intermittent site in a reach of the Po River (North-Western Italy). We installed piezometers reaching −1 m (permanent and intermittent site), and −3 m (intermittent site) and monitored three supraseasonal droughts over a period of three years. We classified the hyporheic fauna into three categories of increasing affinity to life in the hyporheic (stygoxene, stygophile, stygobite), and used communities composition, abundance, beta-diversity and functional groups: (1) to compare assemblages at the same depth but with different hydrological characteristics, as well as assemblages from two depths at the intermittent site, and (2) to assess how the connection with surface water and the direction of the vertical aquifer flow determined the faunistic assemblages. Different taxonomic groups responded differently to intermittence, the hyporheic zone acted as a refuge increasing the resilience of the system, but resilience decreased with increasing degree of affinity to hyporheic life. Disentangling the effects of intermittence on the different faunistic component in the hyporheic zone can help guiding effective protection and restoration measures of river systems with temporary reaches.

1. Introduction

The climate of Earth is changing rapidly [1], posing challenges for species and habitat conservation. Temperature increases and changes in precipitation amounts, patterns and seasonality are leading to habitat loss and fragmentation, changes in species phenology and enhanced rates of biodiversity loss [2,3]. The main effects of climate change on lotic ecosystems are the increased frequency and magnitude of hydrological extremes, with more frequent and extended droughts predicted for the mid-latitudes [4,5]. The Alpine area is strongly impacted by climate change, and it is also under the increasing pressure of water abstraction; as a result, Alpine and perialpine streams are extremely sensitive to the effects of droughts [6]. Riverbed desiccation is one of the most pressing environmental issues related to climate change [7]. Drought is a ‘ramp’ disturbance that disrupts lateral, longitudinal and vertical hydrological connectivity [8], representing a major threat for stream invertebrates [9,10,11,12].

Recently, numerous studies have investigated the mechanisms of resistance and resilience used by benthic invertebrates to confront droughts (e.g., [13,14]). By definition, resistance is the capacity of a taxon, a community or an ecosystem to remain unchanged when being subjected to a disturbance, and resilience is the ability to recover from a disturbance and return to the pre-disturbed state [15]. Resistance and/or resilience are conferred by traits related to the physiological, morphological and life-history features of the organisms [16,17,18]. Stream biota generally exhibits low resistance and variable resilience to supraseasonal droughts [19]. Benthic communities in Alpine streams are generally considered more resilient than resistant [20,21], especially if compared with the aquatic biota of geographical regions where the drying phase is a natural part of the annual flow regime, such as the Mediterranean area [13,16,22]. In the case of intermittent streams, resistance involves a range of physiological adaptations allowing an organism to survive within dry riverbed sediments or remnant pools. For lotic invertebrates, traits that allow in situ desiccation-resistance include: diapause, desiccation-resistant eggs, cocoons or cells, body armoring and aerial respiration [16,23]. Resilient responses to flow intermittence are related to fast reproduction/growth rates or high dispersal ability, and include small body-size, asexual reproduction, active aerial dispersion and invertebrate drift [16,24]. Resilience, therefore, requires the use of specific habitats which retain free water or high humidity as refuges, where the impacts are reduced and survival is enhanced [25,26,27,28]. In intermittent rivers, the accessibility of refuges for aquatic invertebrates during dry periods, and dispersal pathways from these refuges following rewetting may differ due to the variable spatial arrangement of temporary and perennial reaches [29]. Vertical connectivity, and the use of the hypoheic zone as a refuge, becomes a dominant resilience mechanism: benthic invertebrates actively enter refuges as the stream shrinks and dries and the streams are recolonized by invertebrates that survived within the hyporheic zone during stream drying. Both or part of these mechanisms contribute to the resilience of the benthic community [30].

The hyporheic zone may retain water after streambed drying [31] and is well known to act as a temporary habitat for benthic invertebrates [28,30,32,33,34], which use the hyporheic zone as a nursery zone, for the deposition and incubation of eggs and the growth of young instars [35]. The hyporheic zone is used as well as a refuge against droughts [36,37], high superficial temperatures [31], strong sheer stress during high-discharge events [38] and catastrophic floods (e.g., [39,40]). The faunistic contingent seeking refuge in the hyporheic is predominantly composed by insect larvae. These are ecologically classified [41] as stygoxenes, organisms that have no affinities for groundwater systems where they occur only accidentally. Two other ecological groups, which use the hyporheic zone as a non-refugial habitat, can be present: stygophiles and stygobites. The former are those species that actively exploit the resources of the groundwater environment for part of their life cycle, and can be further divided into three categories: (1) the occasional hyporheos consists mainly of benthic insect larvae, the early instars of which reside in the hyporheic zone, but which can also spend all their life in the surface environment; (2) amphibite species complete part of their life cycle (typically the nymphal stage) in the sediment, they include a taxonomically variable group of stoneflies; (3) the permanent hyporheos consists of many organisms of meiofaunal size (<1 mm), which can spend all their life cycle either in subsurface or in surface water, and are represented mainly by crustaceans [41]. Stygobites are specialized subterranean forms that complete their whole life cycle exclusively in subsurface water (almost exclusively crustaceans). Colonization of the hyporheic zone by stream benthos is probably a mixture of active immigration and passive transport [42,43,44]: stygoxenes tend to move downwards, penetrating into the interstitial hyporheic zone during increased disturbance intensity, and move upwards emerging from the sediment after suitable superficial conditions are re-established [45]. The passive transport into and from the hyporheic depends on the direction and strength of hydrologic exchange, and contrasting communities are known to characterize upwelling and downwelling zones ([28], and references therein). Stygobitic taxa typically dominate upwelling groundwater zones, whereas downwelling surface water facilitates stygoxenes and occasional hyporheos (see review in [30]). Body size is another factor which affects the movement and colonization mechanism: meiobenthic invertebrates (organisms passing through a sieve of 500-µm mesh size but retained on a 44-µm mesh [46]) can actively move through the streambed and some taxa are able to swim; however, their dispersal along the surface of the riverbed is mainly passively due to drift ([47], and references therein). Ref. [44] investigated the movement of meiofauna and macrobenthos in the shallow hyporheic of two headwater gravel streams, and recorded predominantly vertical movements for temporary meiofauna (i.e., insect larvae).

We examined the use of the hyporheic by benthic (i.e., stygoxene) and hyporheic (stygophile and stygobite) taxa in an intermittent reach affected by supraseasonal droughts in previously-perennial stream reaches. We hypothesized that the persistence of communities in temporary reaches is primarily associated with resilience mechanisms, and particularly so in stygoxene and stygophile taxa, and that the direction of surface/groundwater exchange would be an important driver of changes in taxonomical and functional characteristics of the biological assemblages. To do so, we investigated community measures of diversity and resilience/resistance functional traits to evaluate taxonomical, compositional and functional changes of benthic and hyporheic communities in relation to the aquifer variation during a three-year period, which encompassed three long supra-seasonal droughts. We tested the following hypotheses:

(1)

different faunistic groups respond differently to intermittence, due to their different degree of specialization to life in the hyporheic; in particular, stygobites would be less affected by drought than stygophile and stygoxenes;

(2)

the connection with surface water and the direction of the vertical aquifer flow determine the faunistic composition at different depths.

(3)

the hyporheic acts as a refuge increasing the resilience of the benthic communities to intermittence, but resilience decreases with increasing degree of specialization to life in the hyporheic.

2. Materials and Methods

2.1. Study Area

The selected area was a low order reach of the Po River in the South-Western Alps hydroecoregion (HER 4, Piemonte, North-Western Italy [48]). The hydrological regime of streams feeding the Po River in this area is nivo-pluvial; hence, several stretches are facing seasonal hydrological alterations in summer due to reduction in precipitation and the subsequent increase in water abstraction. High variability in precipitation, and an increase in mean annual temperature, were recorded in 2017–2019 (data from Piemonte Environmental Agency, https://www.arpa.piemonte.it/rischinaturali/tematismi/clima/rapporti-di-analisi/annuale.html). 2017 was the 3rd warmest year, and the 4th driest year of the last 60 years; 2018 was the 2nd warmest year, and the 5th wettest year of the last 61 years (Figure S1). Finally, 2019 was the 5th warmest year, and the 9th wettest year of the last 62 years (Figure S1). During the extreme summer drought of 2017, the riverbed in piedmont sections of the Po River, including the intermittent sites here considered, completely dried from July [49] to January 2018. In 2018, the riverbed was dry from July to the end of October, and from the end of November to the end of the sampling period in June 2019. Sampling was conducted from 27 July 2017 to 4 June 2019. The only significant rainfall event of the entire sampling period occurred on 27 October–7 November 2018.

We selected a weakly braided reach of about 5 km length in the upper basin of the Po River, where the river runs on alluvial, fluvioglacial and megafan deposits from the Middle–Upper Pleistocene. The surrounding hills and mountains are formed from Lower Triassic siliciclastic units [50]. We selected two sampling stations at 4.9 km distance from each other, with different levels of hydrological permanence. The first one (near Sanfront village, 44°39′16” N, 7°19′27” E, 490 m a.s.l., named hereafter “permanent” station) was in the upstream section, in the main channel, where the river width always ranged around 30 m and the flow was permanent for the whole year. The second sampling station (near Martiniana Po village, 44°39′16” N, 7°19′27” E, 351 m a.s.l, named hereafter “intermittent” station) was in a downstream section, in braided section which experienced non-flow periods, and corresponded to a lateral branch of about 10 m width. The selected reach is one in a set of 13 reaches, which were recently investigated to assess the responses of benthic communities to recurrent dewatering events [51,52,53]; it was selected for accessibility and because it was located on the main river rather than on tributaries.

In order to investigate when, how and to which extent invertebrates used the hyporheic zone by invertebrates during supraseasonal droughts, on 11 July 2017, metal piezometers (internal diameter 12 cm) were inserted into the riverbed using a mechanical drill rig. One piezometer (with holes of 1.4 cm diameter in rows in the last 20 cm of the pipe) was installed in the permanent station in the middle of the channel, reaching 1 m depth below the surface. At the intermittent station, we installed two piezometers at a distance of 10 m from each other. One piezometer, built as the one of the permanent site, was installed in the channel, reaching −1 m depth; a second piezometer reaching −3 m depth, with rows of holes (1.4 cm diameter) from 50 cm below the surface to the bottom of the piezometer was installed at the margin of the river bank. This second piezometer allowed for the collection of integrated samples of all fauna present from the depth of about −1 m from the riverbed surface to the maximum depth of −3 m.

2.2. Hydrological and Physical-Chemical Data

The −3 m piezometer was instrumented with a HOBO U20L-1 water level and temperature datalogger, plus a second datalogger to compensate barometric variations in atmospheric pressure. Data were recorded at 4 h time intervals from 7 November 2017 to the end of the sampling period. The aquifer oscillation rate was calculated from water level data recorded at the −3 m piezometer at the intermittent site. In order to assess temporal trends (i.e., the time intervals required to smooth peaks in recession rates and over which time spans recession rates varied with a correlated trend), we calculated time series of recession rates averaged over different time windows (2, 5, 10, 15 and 20 days) previous to the sampling date; we computed a correlation matrix for the 5 oscillation rates including all data. Positive rates indicated that the aquifer level was raising, negative values that it was lowering. Rainfall (total daily values, measured at Saluzzo monitoring station, 3.5 km downstream of the study reach) and hydrometric levels (daily averages, measured at Villafranca monitoring station, 22 km downstream) were downloaded from the Piemonte Environmental Agency website (https://www.arpa.piemonte.it/). Surface and hyporheic water temperature, conductivity and oxygen concentration were measured before sampling with a multiparametric probe (Hydrolab model Quanta). Due to the instrument malfunctioning, oxygen measurements taken in the last two sampling dates were not considered reliable, and were omitted from the analysis.

2.3. Invertebrate Sampling

Hyporheic fauna was collected with a flexible plastic hose, connected to an electric pump and an electric generator, which was inserted into the piezometer to reach the bottom. The amount of water collected in each sampling occasion varied according to the hydraulic gradient and water level, which affected the pumping power of the sampling device and water availability. The collected water was filtered with a 100 micron mesh plankton net; samples were fixed in the field with 90% ethanol and carried to the laboratory for further identification. Samples were collected at key periods during the hydrological cycle (i.e., at the onset/end of a drought period, during droughts, low and high flows) when the aquifer level allowed for the operation of the pump (Table S1); however, samples were not collected on 20 April 2018 from the shallow intermittent piezometer, and on 8 June 2018 from the permanent and shallow intermittent piezometer because the surface water level was above the piezometer; and on 3 July 2018 from the deep intermittent site for technical problems.

Copepoda, Amphipoda, Isopoda and Cladocera were classified to the species level for adults, and genus for juvenile stages, following [54,55,56,57,58,59,60,61]. All other taxa were identified to the lowest possible level (genus or family) following [62,63,64,65]. Specimens of each taxon were measured with the use of a graduated dish with a scale interval of 1 mm and divided into five length classes: 0–1 mm, 1–2 mm, 2–3 mm, 3–4 mm, >4 mm.

2.4. Data Analysis

All identified invertebrate taxa were classified as macrofauna, permanent and temporary meiofauna based on their size (i.e., ≤1 mm for meiofauna, >1 mm for macrofauna), and as stygoxene, stygophile, stygobite, based on the definition of [41], but considering also [66]. The different filtered volume was controlled by expressing abundances as the ratio of counts to volume (i.e., abundances were expressed as N. ind. L⁻¹).

We investigated in detail the effects of intermittence by comparing the communities collected: (i) at same depth but in distant sites (thus excluding faunistic exchanges between the two sites) with different hydrological characteristics (i.e., intermittent vs. permanent site); (ii) at different depths (shallow and deep, i.e., −1 and −3 m from the surface) of the same intermittent site, where the proximity of the two piezometers allowed to detect the downwards movements of invertebrates when the available habitat contracted (i.e., the aquifer level was lowering); (iii) over time, i.e., assessing the effects of the temporal variations of the aquifer depth. Hence, we categorized the sampling dates/sites (Table S1) according to the following factors:

(1)

Factor “station”: perm = permanent station; int-1 = intermittent station reaching −1 m depth; int-3 = intermittent station reaching −3 m depth;

(2)

Factor “aquifer phase” (for the two intermittent stations only): recession without downwelling, recession with downwelling, rewetting with downwelling, rewetting without downwelling. These phases were based on the aquifer recession rate calculated over relevant time intervals preceding the sampling date (see results: hydrology, for explanation of the time-interval selection), and the presence/absence of surface water. In detail: downwelling flow occurred when surface water was present and recharged the aquifer, and no downwelling when surface water was not present; recession when the water table level was lowering (negative recession rate); rewetting when the water table level was rising (positive recession rate).

(3)

Factor “hydrological phase” based on the condition of surface water at the site: drought (i.e., no surface water), low flow (i.e., surface water level below the threshold of 15 cm at perm and int-1, and 10 cm at int-3), high flow (i.e., surface water level above the 15/10 cm thresholds).

2.4.1. Community Metrics

We used Bray–Curtis dissimilarity matrices based on log(x + 1) transformed data of abundances with each taxon divided, when applicable, into size classes (≤1 mm: meiofauna; 1–2 mm, 2–3 mm, 3–4 mm, >4 mm: macrofauna). Each taxon/size class was further categorized based on its affinity to the hyporheic (i.e., stygoxene, stygophile, stygobite). We calculated the Bray–Curtis similarity matrices between pairs of sites for each of the three affinity classes (named faunistic groups hereafter), adding a dummy variable to correct for denuded samples [67]. We performed multifactorial Permutational Analysis of Variance (PERMANOVA) analyses to test for differences and estimate components of variation due to each of the three factors, and SIMPER analysis to identify the species which most contributed to the Bray–Curtis dissimilarity among samples from different groups [68]. We tested for differences in abundances among and between groups using Kruskall–Wallis and Mann–Whitney tests for the same factors as for the PERMANOVAs. We used a Principal Coordinates Analysis (PCoA) to look for patterns of similarities among groups for the PERMANOVAs with significant results. To assess if the communities converged during periods of low flows (and low aquifer levels) we ran PERMANOVAs to test for differences in community composition between the three sites during low flow and drought, and during high flow and baseflow.

To better understand which processes drive changes in the community composition between the two shallow (permanent and intermittent) sites and the two (shallow and deep) intermittent sites throughout the study, total beta diversity was decomposed into the nestedness and turnover components. Nestedness represents a condition in which differences in the taxonomic composition between sites (or samples, as in this case) are explained by the gain or loss of taxa [69]. For this reason, nestedness generally indicates effects of environmental filters and/or gradients. By contrast, turnover represents a condition in which two sites differ in their composition owing to species replacement, which is usually an indicator of high habitat heterogeneity between sites. For each piezometer, we used presence–absence data of invertebrate taxa in all the available sampling occasions to calculate the nestedness and turnover components, applying the Jaccard coefficient of dissimilarity. We calculate pairwise dissimilarity between the permanent and shallow intermittent sites, and between the deep intermittent and the shallow intermittent sites, for each of the components of beta diversity (total, nestedness, turnover). We calculated the mean of each of these components, and expressed it as percentage of total diversity. All beta diversity analyses were carried out for the whole community and separately for the three main faunal groups (i.e., stygoxene, stygophile, stygobite).

2.4.2. Hydrology/Community Metrics Relationships

We calculated the Pearson’s correlation of the main metrics (total abundance, total number of taxa/size classes, Shannon Diversity Index) for each ecological group (styoxenes, stygophiles, stygobites) with the variation of the aquifer recession rate and average aquifer level in the two intermittent sites to assess which variations of the aquifer determined changes in the communities characteristics. We also focused on periods in the aquifer cycle relevant to explain the possible vertical movements of invertebrates in the intermittent sites and the recolonization timeframe during rewetting; i.e., the periods when water resumed in the shallow piezometer after a drought period (the corresponding dates were 18 January 2018 and 4 December 2018), and analyzed the relation with the aquifer phase.

2.4.3. Species Traits

Species traits were selected and categorized as traits for resilience and resistance [14,24]. We selected a set of traits which can be applied to invertebrates living in gravel-bed streams, and which are related to the ability to enter the hyporheic habitat. Resilient traits are related to dispersal capacity, which governs the rate at which taxa enter hyporheic refugia from the surface, and return to a rewetted river channel from the hyporheic refugia. Resistance traits are the physiological adaptations allowing organisms to survive within dry riverbed sediments or the hyporheic zone. The list of traits is provided in Table S2. The trait values for each taxon were queried from the freshwaterecology.info database and from relevant literature [24,70], complementing when necessary (i.e., for crustaceans) from [66] and from expert knowledge of each taxon (Rossetti G., Stoch F., Cottarelli V. com pers. and Bruno M.C., unpublished). We used a matrix with the total abundance for each taxon not divided into size classes, because we considered it unfeasible to assign traits to the length sub-groups within each taxon. We classified taxa as benthic (if the benthos is their primary habitat, i.e., all insects, Nematoda, Oligochaeta, Ostracoda, Hydrachnidia, Tardigrada) or hyporheic (all crustaceans, roughly corresponding to stygobites) because, without using size-classes, it was not possible to separate stygoxenes from stygophiles.

We followed the approach proposed by [24] to describe the communities in terms of resistance-resilience traits. Presence of trait states was coded as present (=1) or absent (=0) when assigning each trait to each taxon. There were no taxa which had no resistance, nor resilience traits. For each sample, we calculated the total number of resilient and resistant traits overall, and for benthic and hyporheic assemblages separately. We calculated the average number of taxa for each trait category for each site, and for each hydrological phase of the intermittent sites. These calculations were performed for the whole assemblages, and for the benthic and hyporheic assemblages separately. Overall and pairwise statistical differences were tested with non parametric tests (Kruskall–Wallis H and Wilcoxon T, respectively).

All analyses were performed using the following software: PRIMER 6 version 6.1.18 and PERMANOVA + version 1.0.8 [71], STATISTICA 64 [72] and the package BAT [73] in R [74]. Significant thresholds were always set at p = 0.05.

3. Results

3.1. Hydrology and Physical-Chemical Variables

Groundwater (GW) level (at the deep intermittent site) was significantly and positively correlated with total daily rainfall (p < 0.001) and surface water level (p < 0.001 for both monitoring stations), and negatively to GW temperature (p < 0.001); total rainfall was also positively and significantly correlated with surface water level (p < 0.001) and negatively with GW temperature (p = 0.03). These results and the times series in Figure 1A indicate that rainfall events recharged the aquifer and led to quick increases in GW level. These short-term events were superimposed on slower trends in GW level. During prolonged drought periods (e.g., summer 2017 to winter 2018), relatively scarce precipitation resulted in a continuous drop of the GW table that fell below the bottom of the −3 m piezometer after December 2017. During this period, no surface water was observed at the intermittent site, and isolated rainfall events were only able to temporarily reverse this trend. After January 2018, repeated rainfall events—and the return of surface water in the intermittent channels—gradually replenished the aquifer, and the GW level increased approximately 2 m. The GW table then alternatively switched between stationary phases during wet periods with sustained precipitation (e.g., spring 2018) and recession phases when rainfall events were less frequent and intense (e.g., winter 2019). GW temperature exhibited a smoother behavior than GW level, with evident seasonal variations (temperature amplitude around 15 °C) and smaller (1–2 °C) high-frequency fluctuations. It should be noted that heat transfer through the metallic case of the piezometer may have slightly increased the recorded temperature on sunny days. Moreover, when the GW level was below the piezometer bottom (<−3 m; autumn-winter 2017) the recorded temperatures may not fully represent the temperature of the aquifer, due to the absence of water in the piezometer, even though the thermal capacity of the porous medium should have limited temperature fluctuations recorded by the datalogger.

The time series of recession rates calculated over different time windows (Figure 1B) portrayed the highly dynamic behavior of the GW table, which is characterized by an alternation of positive and negative recession rates. Positive recession rates show higher peak values compared to negative rates due to the rapid recharge and GW table rise after precipitation events. The magnitude of these peaks progressively smoothed out with increasing time windows. Depending on the chosen time window, median recession rate ranged between 2.4 cm d⁻¹ (2 days averaging) and 1.5 cm d⁻¹ (20 days averaging). These values exemplify typical values of recession rates at the intermittent site.

Figure 1C summarizes the permanence of GW at different depths during the monitoring period, expressed as the fraction of monitoring days during which the GW table was recorded above a given depth. The piezometer depth allowed to monitor GW levels for the majority (95%) of the period, except for the prolonged drought during Autumn and Winter 2017. The median GW depth was 1.62 m below the ground. Shallow levels (>−0.5 m) of the GW table occurred less than 5% of the time, showing that shallow hyporheic sediments were often in unsaturated conditions. The change in slope in Figure 1C at −2.0 m indicates that the GW level was relatively less dynamic below this depth compared to the overlying sediments. The −1 m piezometer was dry 14% of the recorded time, and the −3 piezometer 5%.

Recession rates were significantly and positively correlated (Figure 1D) only over short time spans, i.e., when correlating 2 and 5 days prior to sampling, 5 and 10 days, 10 and 15 days, 15 and 20 days. Over longer time spans, recession rates varied with a non-correlated trend (i.e., when correlating 2 and 10–20 days prior to sampling, 5 and 15–20 days, 10 and 20 days). Therefore, because short time intervals described similar aquifer fluctuations, we retained 2 and 5 days prior to the sampling date as timeframe to analyze the community metrics. We therefore assumed that the community collected in a certain day was composed by the individuals/taxa living constantly at the collection depth, plus those which arrived in the water column following a period time (of 2–5 days) of constant aquifer movement (raising or lowering of the aquifer level).

Temperature was not significantly different between surface and hyporheic water at the intermittent piezometer, and significantly higher on the surface at the permanent piezometer; conductivity never differed; oxygen concentration and % saturation were significantly higher in surface water for the permanent and deep intermittent piezometers (Figure 2A). Surface and groundwater never significantly differed in temperature, conductivity, oxygen concentration and % saturation among the three stations (Friedman ANOVA by ranks) (Figure 2B–D).

3.2. Hp1: Responses to Intermittence of the Three Faunistic Groups Differ: Community Metrics.

Mean abundances (

Table 1. Community metrics for each factor, divided into the three faunistic groups. Stx = stygoxenes; stp = stygophiles; stb = stygobites; Aquifer phase: rec-no-dwn = recession without downwelling; rew-no-dwn = rewetting without downwelling; rec-dwn = recession with downwelling; rew-dwn = rewetting with downwelling. 2d, 5d: time intervals used for calculation of the hydrological phases (days). Station codes as in Figure 2.

	Stb—Mean Abundance	Stp—Mean Abundance	Stx—Mean Abundance	Stb—Mean Number of Taxa	Stp—Mean Number of Taxa	Syx—Mean Number of Taxa
Perm	0.7	0.32	0.29	13	18	20
Int-1	0.84	1.22	0.65	16	13	18
Int-3	0.23	0.21	0.14	10	13	15
Rec-no-dwn 2d	0.17	0.45	0.21	11	11	18
Rew-no-dwn 2d	0.60	0.42	0.26	5	8	9
Rec-dwn 2d	0.84	0.78	0.53	11	8	13
Rew-dwn 2d	0.10	0.14	0.08	6	5	6
Rec-no-dwn 5d	0.16	0.45	0.23	11	11	19
Rew-no-dwn 5d	0.63	0.42	0.24	5	7	6
Rec-dwn 5d	0.18	0.43	0.14	10	7	7
Rew-dwn 5 d	1.39	0.68	0.83	7	6	12

) were never significantly different among sites for each faunistic group (

Table 2. Results of PERMANOVAs (for factor aquifer phase, only the relevant comparisons are shown, i.e., same phase with and without downwelling, and different phase with/without downwelling). Factor aquifer tested for intermittent sites. 2 days: time interval used for calculation of the hydrological phases. Only significant results are shown. K-W: Kruskall–Wallis test; M-W: Mann–Whitney test; ns: not significant. Station codes as in Figure 2, aquifer codes as in Table 1.

	Overall	Stygobites
	PERMANOVA	PERMANOVA	K-W, M-W
Factor: Station	ns	ns	ns
	ns	ns	int-3 vs. int-1, p = 0.019
	ns	ns	ns
Factor: Aquifer	p = 0.021 (2 days)	p = 0.016 (2 days)	ns
	rec-no-dwn vs. rew-no-dwn, p = 0.016 (2 days)	rec-no-dwn vs. rew-no-dwn, p = 0.014 (2 days)	ns
	rec-no-dwn vs. rec-dwn, p = 0.013 (2 days)	rec-no-dwn vs. rec-dwn, p = 0.014 (2 days)	ns
	Stygophiles	Stygoxenes
	PERMANOVA	K-W, M-W	PERMANOVA
Factor: Station	p =0.01	ns	p = 0.001
	perm vs. int-1, p= 0.012	perm vs. int-1, p = 0.006,	perm vs. int-1, p = 0.002
	int-3 vs. int-1, p = 0.008	int-3 vs. int-1, p = 0.006	int-3 vs. int-1, p = 0.009
Factor: Aquifer	ns	ns	ns
	ns	ns	ns
	ns	ns	ns

); mean abundances were always lower, and number of taxa higher at all three sites for stygoxenes. For all faunistic groups, mean abundances of taxa were highest at the shallow intermittent site, mean abundance and number of taxa were lowest at the deep intermittent site (Table 1). Number of taxa was highest at the permanent site for stygophiles and stygoxenes, and at the shallow intermittent one for stygobites (Table 1). Community composition differed among sites for stygoxenes and stygophiles; the SIMPER analysis indicated the species of these two groups which most contributed to diversity (90% total contribution) were: Bryocamptus (R.) cuspidatus, Ostracoda, Eucyclops serrulatus, Chironomidae >4 mm for the permanent site; Nematoda 1–2, 2–3 and >4 mm, Chironomidae 1–2 mm, Limoniidae 1–2 and 2–3 mm, Ceratopogonidae 1–2 mm, Oligochaeta 0–1, 1–2, 2–3, 3–4 mm; Naididae 0–1 mm, Epactophanes richardi for int-1; Baetis sp. 1–2 mm for int-3 (Table S3). Bray–Curtis similarity was higher when comparing the two intermittent sites than when comparing the two (intermittent and permanent) shallow sites (Table S4) for all faunistic groups; i.e., the three communities were more similar at different depths of the same (intermittent) site than at the same depth but in distant sites with different hydrological characteristics (intermittent vs. permanent). For all groups, the shallow intermittent site was more variable in composition (i.e., lowest Bray–Curtis similarity, Table S4), and the deep intermittent site, less variable (highest Bray–Curtis similarity, Table S4). The comparisons of community metrics over space (shallow intermittent and permanent sites) and over depth (shallow vs. deep intermittent sites) are described in detail below.

3.2.1. Community Metrics: Effect of Intermittence in The Shallow Hyporheic (−1 m)

We compared the metrics calculated for the communities collected at −1 m depth in the permanent and intermittent sites (called perm-1 and int-1 therein) to assess the effects of intermittence on communities living at the same depth. Mean abundances were always lower, and number of taxa higher in both sites for stygoxenes than for the other two faunistic groups; stygophiles were the most abundant faunistic group at int-1, stygobites at perm-1 (Table 1). Mean abundances were higher at int-1 than at perm-1 for all groups, the number of taxa was highest at perm-1 for stygophiles and stygoxenes, at int-1 for stygobites (Table 1 and Table 2).

Bray–Curtis similarity between the two sites was slightly higher for stygoxenes (Table S4), nonetheless, community composition differed significantly for stygoxenes and stygophiles (Table 2), as also shown by the PCoA (68.8% and 71.7% variance explained, respectively, Figure 3A,B). Differently from the previous two faunistic groups, differences in assemblages were not significant for stygobites. For all groups, the shallow intermittent site was more variable in composition (i.e., lowest Bray–Curtis similarity, Table S4) than the permanent site.

3.2.2. Community Metrics: Effect of Intermittence Over Depth

We compared the metrics calculated for the communities collected at −1 m and at −3-m depth at the intermittent sites (called int-1 and int-3 therein) to assess the effects of intermittence on communities living at different depths. In fact, the int-3 piezometer collected a depth-integrated sample of fauna living from about −1 to −3 m, i.e., all the animals which usually live at −3 m depth, plus those which migrated deeper into the aquifer when its level was lowering. At the deep intermittent site, total abundances within the hyporheic zone increased during the latter stages of the drought phases; during the first supraseasonal drought (July 2017–January 2018), the proportion of stygoxenes + stygophiles increased while the proportion of stygobites remained similar; a similar phenomenon occurred in the first of the two following and shorter droughts (August–September 2018; December–June 2019, the two separated by a flood), while in the second abundances increased but were due to an increase of stygobites (2 species of Cyclopoida, 2 of Niphargidae) (Figure 4). Mean abundances were lowest and mean number of taxa highest at both sites for stygoxenes; stygophiles were the most abundant faunistic group at int-1, stygobites at int-3 (Table 1). Mean abundances and number of taxa were higher at int-1 for all three faunistic groups, and mean abundances were significantly different between sites for stygobites and stygophiles (Table 2). Community composition differed significantly for stygoxenes and stygophiles when comparing the two intermittent sites (Table 2; PCoA 68.8 % and 71.7 % variance explained, respectively, Figure 3A,B). Differently from the previous two faunistic groups, the two sites were not significant different in composition for stygobites.

Only during four sampling occasions were paired samples collected (i.e., water was present at both piezometers): 18 January, 1 February, 16 March, 4 December 2018. In these occasions, three stygoxenes taxa were exclusively collected at −3 m (specifically, Chironomidae 2–3 mm, Leptophlebidae 3–4 mm, Bryocamptus (Rheocamptus) zschokkei). The remaining stygoxene taxa were either shared between the two depths (Oligochaeta, Chironomidae and Baetis sp. 1–2 mm) or present only at −1 m (the remaining 11 taxa). Stygophiles and stygobites taxa were present either at both depths on each date (Oligochaeta, Naididae, Chironomidae, Ceratopogonidae, unidentified Ephemeroptera 0–1 mm and Acanthocyclops cf magistridussarti, Speocyclops cf franciscoloi, Niphargus microcerberus, respectively), or only in the more superficial samples (5 and 11 taxa, respectively); there were not any stygophiles and stygobites taxa exclusively collected from deep samples.

3.3. Hp 2: Hydrological Connectivity and Aquifer Flow Dynamics Drive Faunal Patterns: Hydrology/Community Metrics Relationships

The communities did not converge in composition during periods of low flow and drought (PERMANOVA p = 0.003), but converged during high flows and baseflow (PERMANOVA p = 0.187). The responses of the hyporheic communities to the aquifer variation rate were assessed only for the intermittent sites. The correlations of the main metrics (total abundance, total number of taxa/size classes, Shannon Diversity Index) of each ecological group with the variation of the aquifer oscillation rates calculated over the 2 and 5 preceding days were never significant. Conversely, the aquifer level at the deep intermittent site did correlate with several metrics (Table S5): the total abundances of stygoxenes and stygophiles were negatively correlated with the level averaged over 2 and 5 days, the diversity of stygobites was positively correlated with the level calculated over 5 days intervals.

Mean abundances were not significantly different among hydrological phases (Table 2). Stygobites were the only faunistic group for which community composition differed for aquifer phases calculated over 2 days; in particular, recession without downwelling (very low mean densities, high number of taxa) differed in composition from recession with downwelling (highest mean densities, same high number of taxa), and from rewetting without downwelling (which had higher average densities, but less taxa) (Figure 5, first axis PCoA, Table 2). Composition was most similar for rewetting with downwelling, i.e., when surface water was present but the aquifer lowering. Stygobitic taxa characterizing these different aquifer phases (SIMPER analysis, 90% contribution) were Acanthocyclops cf magistridussarti, Diacyclops zschokkei, Speocyclops sp. 2, Niphargus microcerberus for the recession without downwelling phase; Niphargus transitivus, Speocyclops cf franciscoloi, Diacyclops antrincola for rewetting without downwelling; Speocyclops cf franciscoloi, Stammericaris sp. 2, Phreatalona protzi for recession with downwelling.

Lastly, to get a better grasp of the recolonization timeframe during rewetting, we analyzed in detail the effects of rewetting after the first drought (from July 2017 to January 2018), comparing the last drought sample (16 November 2017) collected only at the permanent and intermittent -3 site, with the sample collected on January 10, after 3 days of intense rainfall (total precipitation = 2.44 cm), which raised the aquifer of about 1.8 m, from −3 to −1.2 ma below the surface. On that date, the shallow intermittent piezometer was still dry although surface water was present (at the deep intermittent site, the level rose about 1.9 m in three days), the rewetting rate was 24.22 cm/day, with rates calculated over a two-day time interval (Table S1); the aquifer level rose further in the following days, even if the rain ceased, and on January 18 the shallow piezometer was again with water, after about 7 months of drought. The recovery was fast, as abundance of stygobites, stygophiles, stygoxenes were, on the 18 of January, 0.4, 1.02, 0.02 ind L⁻¹, respectively.

3.4. Hp 3: the Hyporheic Use Increases the Resilience of the System, with Faunistic-Specific Resilience Responses

3.4.1. Analysis of Nestedness and Turnover

To assess the effects of intermittence in the shallow hyporheic, we analyzed the percentage contribution of the taxa nestedness and turnover to the total diversity, comparing perm and int-1 stations. For the whole community and for stygobites the measured total diversity between the two sites of 0.85 and 0.86 (respectively) was due almost equally to loss/gain of taxa (nestedness, 50.8 and 51.1%, respectively) and to species replacement (turnover, 49.2 and 48.9%). Conversely, the diversity of 0.75 and 0.86 measured respectively for stygophiles and stygoxenes was due for 66.9 and 65%, respectively, to nestedness, and species turnover contributed for only 33.1 and 35% (Figure 6A). The number of taxa of these two faunistic groups was higher at the permanent station (Table 1), suggesting that the assemblages at int-1 are predominantly a subset of those at the permanent site, although about 25% amount of taxa replacement did occur.

The analysis of taxa replacement and turnover between int-1 and int-3 sites allowed assessing the effects of intermittence along a gradient of aquifer depth. For the whole community the measured total diversity between the two sites scored 0.87 and was due equally to loss/gain of taxa (nestedness) and to species replacement (turnover); the diversity of 0.90, 0.82 and 0.87 measured respectively for stygobites, stygophiles and stygoxenes was due 65.9, 63.7, 73% to nestedness (Figure 6B). The number of taxa of these faunistic groups was higher at the shallow intermittent station (Table 1), suggesting that the assemblages at int-3 are predominantly a subset of those recorded at int-1, although there was also a small amount of taxa replacement.

3.4.2. Species Traits

Table 3. Results of Kruskal–Wallis ANOVA by ranks and Wilcoxon matched pair tests. Comparisons for hydrological phases for intermittent sites only. Resil = resilience; resist = resistance; BT = benthos; HR = hyporheos. Station codes as in Figure 3.

	Among Sites			Among Hydrological Phases
	Valid N	Kruskal–Wallis H	p-Value	Valid N	Kruskal–Wallis H	p-Value
Resilience	36	5.167	0.076	20	3.237	0.198
Resistance	36	6.949	0.031	20	3.787	0.151
Resilience BT	36	0.578	0.749	20	1.074	0.585
Resilience HR	36	2.743	0.254	20	0.648	0.723
Resistance BT	36	0.329	0.848	20	0.290	0.865
Resistance HR	36	4.299	0.117	20	2.856	0.240
	Over all samples			Over all hydrological phases
	Valid N	Wilcoxon T	p-value	Valid N	Wilcoxon T	p-value
Resist vs. resil	35	36.5	<0.001	19	1.5	<0.001
Resil BT vs. resil HR	35	291.5	0.700	19	61.5	0.178
Resist BT vs. resist HR	34	161.5	0.020	18	25.5	0.009
Resist BT vs. resil BT	29	12.0	<0.001	17	7	<0.001
Resist HR vs. resil HR	32	0.0	<0.001	17	0	<0.001
	Int-1			Int-3
	Valid N	Wilcoxon T	p-value	Valid N	Wilcoxon T	p-value
Resist vs. resil	5	0.0	0.043	14	1.5	0.001
Resil BT vs. resil HR	4	2.0	0.273	15	42.5	0.320
Resist BT vs. resist HR	4	0.0	0.068	14	20.0	0.041
	Perm			High flows
	Valid N	Wilcoxon T	p-value	Valid N	Wilcoxon T	p-value
Resist vs. resil	16	16.0	0.007	4	0.0	0.068
Resil BT vs. resil HR	16	32.5	0.066	3	2.0	0.593
Resist BT vs. resist HR	16	61.0	0.717	3	2.5	0.789
	Low flows			Drought
	Valid N	Wilcoxon T	p-value	Valid N	Wilcoxon T	p-value
Resist vs. resil	3	0.0	0.109	12	1.5	0.003
Resil BT vs. resil HR	3	2.0	0.593	13	22.5	0.108
Resist BT vs. resist HR	3	0.0	0.109	12	8.5	0.017

shows the results of the statistical comparisons of number of resilience/resistance taxa over all samples, for each station and hydrological phases, for the whole assemblages and for the hyporheic and benthic assemblages separately. The box-plots of resilience and resistance traits are shown in Figure 7. Resilience traits were always more abundant on average than resistance ones (1.4 times overall; 1.2 times for benthos, 1.7 times for hyporheos), and significantly so for the whole community, benthos and hyporheos (p < 0.001). The number of resilient taxa was overall higher (1.1 times) and the number of resistant taxa significantly lower (0.8 times, p = 0.020) in the hyporheos than in benthos.

There were more resilience than resistance traits at each station for benthos and hyporheos, and significantly so for the whole community. Resistance traits were the only traits category which differed significantly in number among stations for the whole community (p = 0.031), and were significantly lower in the hyporheos than in benthos at the deep intermittent site (0.5 times, p = 0.041). Both types of traits were more abundant at the shallow intermittent site and less abundant at the deep intermittent site for both benthos and hyporheos.

There were significantly more resilience than resistance traits for each hydrological phase and significantly so for benthos and hyporheos (p < 0.001), and for the whole community overall (p < 0.001) and during droughts (p = 0.003). Number of resilient and resistant taxa were always higher in benthos than in hyporheos during each of the hydrological phases, and significantly so for resistance traits (p = 0.009). Resilience and resistance traits were as average more abundant in benthos and hyporheos during low flows, and less abundant during drought. During high flows and drought, the number of resilient taxa was 1.3 times higher for the hyporheos compared to benthos, whereas the number was lower (0.9 times) during low flows. The number of resistant taxa was significantly lower (0.5 times, p = 0.017) in the hyporheos than in benthos during drought.

4. Discussion

In this study, repeated sampling along different intermittence and aquifer depth conditions allowed for the assessment of the relationships between invertebrate assemblages’ diversity and composition, and water permanence. Our main results are that flow intermittence is a primary driver of aquatic communities in these rivers, and community responses are largely due to resilience rather than resistance mechanisms, supporting the findings of [24]. We targeted different faunistic groups collected only in the hyporheic zone (sensu [74]), and their responses to disturbances differed based on their affinity with the hyporheic habitat, i.e., if they were already present in the interstitial space before the disturbances occurred [30], or if they infiltrated into the HZ during disturbance.

4.1. Hp1: Responses to Intermittence of the Three Faunistic Groups Differ

The first hypothesis we tested stated that different taxonomic groups would respond differently to intermittence, with stygobites being less affected by drought than stygophile and stygoxenes. Indeed, the communities of the three faunistic groups were more different in distant sites (which, however, are hydrologically-connected during periods of surface flow) with different hydrological characteristics (intermittent vs. permanent) than at different depths of the same (intermittent) site. This supports the hyporheic refuge theory according to which taxa migrate into the hyporheic zone to avoid unfavorable conditions when intermittence occurs [19,27,75]. Gradients in invertebrate communities, shifting from stygoxenes via stygophiles to stygobites occur with increasing distance from surface water ([76], and references therein). Our results suggest that when intermittence occurred, surface taxa and the temporary hyporheos (stygoxenes + stygophiles) moved into the shallow intermittent site (they were more abundant there compared with the permanent site). The presence of intermittent or permanent surface water had a strong effect on the invertebrates collected in this shallow hyporheic zone, highlighting its use as a refuge from flow intermittence, as reported in numerous other studies (e.g., [28,75,77,78]). In fact, we found more individuals, and more temporal variability, at the shallow intermittent sites than at the shallow permanent site for the more surface-related faunistic groups. For these taxa the use of the hyporheic zone represents a survival strategy, because the intermittent site experienced several consecutive months of complete drought. However, a reduced number of taxa of these categories were able to survive there, because there was often no connection with surface water. In our study, copepods and numerous insect larvae of different size classes characterized these assemblages. [79] investigated the hyporheic assemblage responses to variation in flow permanence and surface–subsurface exchange along a 52-km long flow-permanence gradient, and reported that hyporheic abundances varied significantly between temporary and perennial sites in direct relation with surface flow permanence. Recent studies [51,52] on intermittence on surface habitats in the Po River watershed, which included the area of this study, showed that during the supraseasonal drought, surface refugia (pools, wet woody debris, etc.) in the vicinity of the piezometer were not present, and surface benthic invertebrates communities were significantly altered in intermittent sections, with lower diversity and density than the upstream permanent reaches.

We can assume that, as stygoxenes and most of stygophiles colonize the hyporheic from the surface, the impacts on supraseasonal drought on surface water habitats and communities are transferred to the hyporheic zone. We, in fact, recorded that the assemblages of stygoxenes and stygophiles collected at the shallow intermittent site were predominantly a subset of those living at the permanent site. This underlines the impacts of drought on the potential colonizers of the hyporheic habitat and, hence, the role of hyporheic interstitial space in benthic invertebrate recovery. In fact, both benthic invertebrates which actively enter refuges as the stream dries, and those early larval stages that already are present in the hyporheic during stream drying, can contribute to the resilience of the benthic community [30]. Differently for the two previous groups, stygobite assemblages at the same shallow depth were quite similar even if the two sites were distant, and structured by the same amount of taxa replacement and nestedness. In fact, stygobitic taxa are adapted to groundwater and their communities are not structured by colonization from the surface.

We also assessed the effects of intermittence at different depths of the aquifer at the intermittent site. The deep piezometer, being perforated along its entire length, allowed collecting the animals living from the top of the aquifer to a depth of 3 m from the surface, such as those that followed the recession of the water table and were able to survive at lower depth. In fact, evidence indicates that the hyporheos migrates into deeper, saturated sediments during drying of the surface stream ([80], and references therein). Although surface water was often not present around this piezometer, communities were less variable in composition over time there, than closer to the surface (i.e., in the shallow piezometer). This was expected since in the upper layer of the hyporheic the aquifer level was very variable, and water level lowered, ultimately leading to the desiccation of the hyporheic sediment. Conditions were therefore very selective for all types of fauna, requiring adaptations and ability to move within the sediment following the recession of the water table, whereas at the depth reached by the deep piezometer water was rarely absent. However, there were in general less individuals and taxa in the deeper samples. This trend was stronger for stygoxenes and stygophiles, as the environmental conditions became more selective for these surface-related taxa with increasing depth (especially the decrease in oxygen concentration and percentage saturation). Unexpectedly, a similar trend was recorded for stygobites.

Differences between depths were due to species loss with increasing depth for all faunistic groups, and more strongly so for stygoxenes, with the assemblages at the deep site being predominantly a subset of those recorded at the shallow intermittent site. The use of the hyporheic zone was taxon-specific: only Baetis sp. 1–2 mm among the surface-related taxa was characteristic of this site. Baetis sp. was one of the taxa that increased in the hyporheic zone after 15 to 24 h of surface water drying in the experiment of [75], confirming the ability of this mayfly to respond to drought by seeking refuge into the hyporheic zone. We also recorded some taxa exclusively at int-3 when water was present also at int-1, i.e., taxa that preferentially migrated at lower depths; surprisingly, these were all stygoxenes, represented by insect larvae (Chironomidae 2–3 mm) and nymphs (Leptophlebidae 3–4 mm), and one harpacticoid copepod (Bryocamptus (Rheocamptus) zschokkei).

During the 2017–2018 supraseasonal drought, when the aquifer level was lower than 1 m from the surface for about 7 months, stygoxenes and stygophiles moved to a greater depth and increased in abundances with the progressing of the drought, whereas the proportion and abundance of stygobites remained similar to that of non-drought periods. Similar results were observed by [72], who recorded an increase in the number of benthic taxa and the proportion of benthos within the hyporheic zone during the latter stages of a supraseasonal drought, although at a depth of 20 cm from the surface.

Taken together, these results support the hypothesis 1: although for stygophiles and especially for stygoxenes the conditions of the deeper aquifer are rather selective, the progressive lowering of the aquifer during the supraseasonal drought drives them to disperse to a depth of 3 m from the surface (probably with both passive and active mechanisms). As reviewed by [30], most of the assessments of the hyporheic zone as a refuge for benthic taxa have been so far based on samples collected from depths less than 1 m, with only few instances of deeper samples, with a maximum recorded of 2.6 m. Studies directed to the permanent and occasional hyporheos (i.e., stygophile taxa) and the effects of drought and low flows on these faunistic components, are still scarce (see review in [28]), but recently [78] recorded taxon-specific migration of temporary hyporheos into the hyporheic zone (with a maximum sampling depth of 90 cm) within 24 h from sediment drying. Our study therefore represents one of the few cases that demonstrated the persistence of benthic (stygoxene) and stygophile taxa in the hyporheic at a depth from less than 1 to 3 m for several months during the supraseasonal drought.

Stygobites living in the deep hyporheic habitat were less impacted by drought, as a large component of this assemblage occupies this habitat permanently. However, there were in general less individuals and taxa of this group in the deeper samples, possibly due to the impact of droughts. Stygobitic taxa exploited the more superficial hyporheic when it was saturated by water. In fact, the shallow aquifers are areas of marked fluxes of nutrients, and organic carbon concentrations decline with depth and distance from the surface [81]. Nutrient and carbon fluxes in the hyporheic zone are strictly depending from the hydrological exchange between surface and subsurface water; when the press disturbance caused by drying occurs, subsurface macro- and meiofauna might become deprived of food.

4.2. Hp 2: Hydrological Connectivity and Aquifer Flow Dynamics Drive Faunal Patterns

The second hypothesis we tested was that the connection with surface water and the direction of the vertical aquifer flow controlled the faunistic composition at different depths. We therefore analyzed how the community metrics were related to the hydrological phases, with a particular focus to the first rewetting, after the longest supraseasonal drought. Hp 2 appears to be supported, because the connection with surface water and the direction of the aquifer movement determined variations in taxonomic and functional diversity over space and depth; the different assemblages recovered fast in abundance after the very long supraseasonal drought.

The crucial role of the direction and strength of hydrological exchange in shaping hyporheic communities has been widely investigated and reviewed in detail by [30]; contrasting communities are usually reported as characterizing upwelling and downwelling zones, with stygobitic taxa (Copepoda and other microcrustaceans, Amphipoda, Isopoda) taxa dominant in upwelling zones, stygoxenes and stygophiles in downwelling zones ([28], and references therein). Our study reach was characterized by a quick rainfall infiltration rate and, as a consequence, the groundwater table alternatively switched between stationary and rewetting phases during wet periods with sustained precipitation, and recession phases when rainfall events were less frequent and intense. If animals move downwards following the recession of the aquifer, and conversely move upwards towards the surface during rewetting, the aquifer oscillation rate would be relevant. In fact, the animals moving within the saturated sediment would be stranded during the recession phases if the rate is faster than their speed. Unfortunately, field measurements of rates and velocity of invertebrates’ movements into the hyporheic zone from the benthic layer and vice versa are scarce. The measured recession rates in our study (maximum values of 0.016 cm min⁻¹ during rewetting, and 0.014 during recession) are well in the range of the only existing (to our knowledge) measurement of behavioral vertical movements of meiofaunal invertebrates [82] in flume experiments, where copepods and chironomids moved downwards at a rate of 0.09–1.6 and 0.47–1.3 cm min⁻¹, respectively. In our study the stygophile and stygobite assemblages were mainly composed of Chironomidae and Copepoda (Bryocamptus (Arcticocamptus) cuspidatus, Acanthocyclops cf magistridussarti, Diacyclops antrincola, Speocyclops cf franciscoloi). Although we do not have data on the crawling abilities of the remaining abundant taxa (Hydracarina and Oligochaeta), we can assume that also the organisms in these two assemblages were able to move both actively and passively following the movements of the aquifer, rather than remaining stranded in the drying sediment and die. Organisms of the meiofaunal size are reported to be highly susceptible to passive drift, and their distribution is influenced by streamflow at scales ranging from 10 s to 100 s of meters [83]. The analysis of the effect of the first rewetting in the shallow intermittent site after the first and longest drought (lasted from July 2017 to January 2018), showed a fast recovery in abundance of stygobites, stygophiles, stygoxenes. The aquifer oscillation rate at the onset of this main rewetting was 24.2 cm day⁻¹, for 1.9 m level rise; as a consequence, fauna should have taken a maximum of 8 days to reach the 3 m depth. Unfortunately, data on recolonization rates of the hyporheic zone in natural conditions are scarce, reporting maximum densities reached in a range of 1–3 days, but for a maximum depth of 20 cm [43,84,85,86], i.e., much shallower than our sites.

The communities did not converge in composition under hydrological stress (during periods of low flow and drought) as a result of taxon-specific responses, but converged during no-stress periods (high flows). Nonetheless, our results suggest that faunistic assemblages indeed responded to the different hydrological phases. Water was near the surface (less than 1 m depth) for less than 20% of the entire sampling period, and the water permanence at the shallow intermittent piezometer was concentrated over a short period. Thus, the aquifer level, more than the oscillation rate, was important in determining invertebrates survival in the hyporheic zone: the lower the level for 2 to 5 days, the higher the abundance of stygoxenes and stygophiles, indicating that these taxa move quickly downwards when the aquifer level lowers (during droughts and recession). Conversely, the abundance of stygobitic taxa, which permanently dwell in the aquifer, was not related to the aquifer level, suggesting that these assemblages find a suitable habitat at lower depths. The composition of the stygobitic assemblages differed between some phases of hydrological exchange calculated over the short time span of 2 days: if the aquifer was lowering and there was surface water feeding the aquifer or not, or if there was no surface water but the aquifer was lowering or raising (in the latter case, either lateral flow or subterranean flow was locally recharging the aquifer). In the first case, if the recession of the water table increases the downwards movements of stygobites by active or passive movements, the presence or absence of surface water recharging the aquifer determines the input of dissolved oxygen, organic matter and nutrients. Where subsurface water is enriched in nutrients and organic matter, due to intense exchange with the surface, stygobites may not be able to compete with faster-growing stygophile or even stygoxene fauna [76]. In the second case, the raising level of the aquifer was driven by lateral or deep upwelling flow, which would promote the passive dispersal of stygobites, as shown by the highest abundances recorded when comparing to the phase of aquifer lowering. Active movements of stygobitic fauna within the hyporheic zone probably occurred as well. [44] showed that permanent meiofauna moved actively in the vertical (downwards within the interstitial habitat), downstream and upstream in the shallow hyporheic of a mountain gravel stream.

The hyporheic habitat is patchy at the fine scale [40,87] as a result of the complex responses of organisms to interstitial water velocity, sediment pore size, organic matter content, dissolved oxygen concentration and other environmental parameters as well as biological interactions [87,88,89]. As a result that only a few measurements of the physico-chemical variables were taken in the intermittent sites, we could not properly evaluate the effects of changes in temperature and oxygen content. However, the values of these variables did not appear to differ much among sites. Given the broad oscillations of the aquifer level and the duration of the unpredictable supraseasonal drought, we assume that hydrology was the overwhelming driver of the observed faunistic patterns, as reported in numerous other studies [30].

4.3. Hp 3: the Hyporheic Use Increases the Resilience of the System, with Faunistic-Specific Resilience Responses

Following the conceptual models of [24,90], the analysis of nestedness and turnover provides indications of the mechanisms structuring invertebrate communities, i.e., invertebrate communities structured by resilience are nested, whereas communities structured by resistance have high rates of turnover. When drying is unpredictable, nestedness is replaced by turnover along a gradient of increasing disturbance [24,90]. We extended this interpretative model of beta diversity partitioning to the hyporheic habitat under unpredictable intermittence, focusing on the responses of different faunistic groups, with an increasing stygobization level. We capitalized on the work of [51], who investigated the mechanisms and rate of recovery for macroinvertebrates in the same sites investigated in the present work, when surface flow resumed after the supraseasonal drought of July 2017–January 2018. Reference [51] found that benthic macroinvertebrate communities recovery was driven by resilience rather than resistance. In fact, drought markedly reduced the diversity and density of macroinvertebrates, and passive drift from upstream reaches was the most probable source of post-drought recolonizers. Recovery by recolonization (recruitment) depends on species’ abilities to move through space and time from adjacent habitat sources either by active dispersal or by resting stages [91]. In the case of the intermittent reaches in the Po River, investigated in the present work, in-stream refugia had a negligible contribution, because the intensity and length of the drought caused the disappearance of surface pools [51]. Drought survival in the hyporheic zone is primarily linked to the permeability of the substratum, the severity of disturbance (magnitude, duration and timing) and the resistance/resilience of individual taxa [30]. In our study communities became increasingly taxonomically dissimilar with increasing disturbance by drought. For stygoxenes and stygophiles, the filtering effect (i.e., nestedness) was the main driver of beta diversity. Resilience was the principal persistence mechanism for stygophiles and stygoxenes along the intermittence/harshness gradient, as taxa-poor communities at the most temporary reaches were nested-subsets of richer communities found at the least temporary and perennial reaches. Stygobites, on the other hand, are organisms typically adapted to inhabit the subsurface environment, do not occur in surface water, follow the aquifer while it is receding (during droughts) and remain in the hyporheic when this is reconnected with surface water. For this assemblage, there was no filtering effect due to intermittence at the shallow site, whereas, unexpectedly, there was a loss of taxa from the shallow to the deep hyporheic, suggesting a possible impact of drought on this assemblage.

In a following step, we analyzed which traits characterize the three faunistic assemblages along the gradient. The results of the species-traits analysis further strengthened the outcomes of beta diversity partitioning, although we used a different taxonomic and faunistic classification for this set of analysis, in order to be able to apply the resilience/resistance trait approach proposed by [24]. Resilience traits were always dominant along the intermittence gradient and during all hydrological phases, for both benthic and hyporheic assemblages. For both benthos and hyporheos, these types of traits were proportionally more abundant than resistance traits. Resilience traits were more abundant during low flows and at the shallow intermittent site, i.e., when benthic invertebrates would enter the hyporheic zone as conditions worsen in the surface stream, and hyporheic taxa might use the shallow hyporheic zone during adverse conditions in the groundwater environment [75], and at the depth most affected by drought and aquifer recessions. The hyporheos was more resilient and less resistant than benthos, as expected as these assemblages are permanently living in the interstitial spaces, and can passively follow the recession of the water table [81], thus persisting in water-saturated layers.

The third hypothesis i.e., that the hyporheic act as a refuge increasing the resilience of the system to intermittence, but resilience decreases with increasing degree of specialization to life in the hyporheic was therefore supported. In fact, intermittence acted as a filter for stygoxenes and stygophiles, and not for stygobites, and all communities collected in the hyporheic were mainly resilient to intermittence. The lack or paucity of resistance traits, however, implies that, where water in the hyporheic interstitial space is not available, the hyporheic zone fails as a refuge [92].

5. Conclusions

Alpine streams are currently facing a “Mediterraneization” process, with the increase of flow intermittence, because of the combined effects of climate change and anthropogenic pressures, e.g., intensification of water abstraction and land-use alteration. Almost no data are available at present about the importance of the hyporheic zone in increasing the resilience of previously-perennial mountain streams facing intermittence. As reviewed by [28], the hyporheic zone is an important component in the suite of refuges that facilitate community resilience to disturbance events. However, invertebrate colonization following the resumption of flow remains a poorly-known process [24,30], even more so when the source of recolonizers is the hyporheic zone. In fact, the role of vertical active and passive migrations into and from the hyporheic zone is still poorly known [28], although vertical movement of benthic and hyporheic taxa into the sediment do occur even in absence of a disturbance [44].

Although the intermittent reach in our study belongs to the “downstream drying” pattern discussed by [19], and ref. [51] suggested passive drift from upstream reaches was the most probable source of recolonizers for post-drought recovery of surface benthic taxa in the study area, in the present work we were able to assess the role of the hyporheic zone as a possible source of recolonizes after flow resumption, or after phases of stress caused by the strong withdrawal of surface water. This role extends also to the off-channel hyporheic zone along the river banks (adjacent to our deep piezometer). Ref. [93] challenged the notion that drift is the primary source of recolonization in intermittent alluvial river, and in a manipulative experiment showed that colonization from the hyporheic zone following rewetting was the primary process promoting benthic communities persistence. We also confirmed that refuge use is taxon-specific, depending on a range of morphological, behavioral and physiological traits [28].

We are aware of the spatial limitations of our study, nonetheless, our results provide a description of the surface-hyporheic faunal dynamics over one very long, and two shorter droughts in river reaches which historically were characterized by a perennial flow regime. The hyporheic refuge use is patch-specific in heterogeneous habitats such as low-order gravel streams, and gathering information on short-scale dynamics is relevant because localized refugial hot-spots may support enough individuals for subsequent recruitment and recolonization of the surface [94]. A better knowledge of the role of the hyporheic zone in promoting resilience and recovery is important in guiding conservation and management decisions on a local and regional scale. This becomes more relevant, as the prolonged drying caused by climate change, exacerbated by water withdrawals for multiple uses, often exceed the thresholds of native species’ adaptations and can change community structure and ecological processes [95].

Studies directed to the permanent and occasional hyporheos (i.e., stygophile taxa) or groundwater fauna (stygobionts), and the effects of drought and low flows on these faunistic components, are still scarce (see review in [28]). [75] warned against making inferences about hyporheic communities based on benthic sampling, and indicated a pressing need for research analyzing benthic and hyporheic communities simultaneously. We provided a contribution by extending the assessment of the effects of droughts to the stygobitic assemblages, which resulted to be partly affected by the supraseasonal droughts. Stygobites provide important ecosystem services such water purification, bioremediation and water infiltration [26]. However, in many areas of the world groundwater communities remain poorly studied and they have never (to our knowledge) been addressed as at risk in newly-intermittent streams. As a result that fully successful river rehabilitation must include restoration of vertical linkages between the river and its shallow groundwater aquifers [95], disentangling the effects of intermittence on the different faunistic component in the hyporheic zone can help guiding effective, holistic river protection and restoration measures, especially in river systems with temporary reaches.