The Trebbia River is a 120 km-long tributary of the Po River with a mean annual discharge of nearly 20 m³ s^-1. Water is diverted for irrigation at ~20 km upstream from the confluence with the Po River and the environmental flow downstream of the last withdrawal is set at 1.5 m³ s^-1. In July and August, the median discharge upstream of the water withdrawal is 3.9 m³ s^-1, with an interquartile range of 2.3 and 5.7 m³ s^-1, respectively (data from the 2003–2015 period). The environmental flow does not generally suffice to feed the river up to the confluence with the Po River and the last ~5 km upstream of the confluence generally dry out in summer.

Six sampling sites from three different Trebbia River sections were selected (Fig. 1): i) upstream of the water abstraction - perennial (St1, St2); ii) downstream of the water abstraction - perennial (St3, St4); iii) downstream of the water abstraction – intermittent (St5, St6). Altitude ranges from 162 m a.s.l. (St1) to 49 m a.s.l. (St6). Intermittent sites experienced a dry period of nearly 2 months in summer 2017, during which the streambed completely dried. Macroinvertebrates were collected in December 2017, 2 months after flow resumption, and at the end of June/beginning of July (hereafter July) 2018, before the dry period. At each site, macroinvertebrates were collected at 10 locations in the main channel with a Surber net (0.05 m² frame area). Surber replicates were pooled together and organisms were sorted directly in the field and preserved in absolute ethanol. Ethanol was replaced 1 day after collection and samples were preserved at 4 °C. Organisms were morphologically identified at the lowest possible taxonomic level, usually genus or family, using a stereomicroscope according to Campaioli et al. (1994, 1999), Tachet (2010) and Miller and Bergsten (2016).

Figure 1.

Map of the study area.

Metabarcoding was performed according to Elbrecht and Leese (2015). Briefly, each sample was dried overnight at room temperature and then ground with an IKA Ultra Turrax Tube Drive Control (Staufen, Germany) for 30 min at 4000 rpm. For each sample, only one part of the body of larger organisms (e.g. legs) was retained to increase the probability of detecting smaller organisms and rare species. DNA was extracted according to a modified salt extraction protocol (Weiss and Leese 2016). RNA was digested with RNase A to improve downstream analysis steps. DNA was purified using the NucleoSpin Gel and the PCR clean-up kit (Macherey and Nagel, Düren, Germany). DNA sample concentration was equilibrated to 25 ng μl^-1 and the first PCR step (20 cycles) was performed using the BF2+BR2 primers according to Elbrecht and Leese (2017) with the Qiagen Multiplex PCR Plus kit. Two replicates per sample were analysed to increase accuracy. Samples were multiplexed following the tagging strategy described in Elbrecht and Steinke (2019) in the second PCR step (15 cycles). The final library was prepared by pooling an equal amount of amplicons for each sample. The library was purified using 0.76× SPRI select (Beckman Coulter, Krefeld, Germany) and sequenced with an Illumina HiSeq Rapid Run 2 × 250 bp as part of a library containing nearly 280 samples at the Eurofins Genomics (Ebersberg, Germany). Illumina reads have been deposited in the SRA database under the accession number PRJNA655393.

One of the most abundant species recovered by morphological identification, Oligoneuriella rhenana, was not detected by metabarcoding neither at species nor genus when a 97–95% cut-off was used for species/genus assignments. Single O. rhenana individuals were thus barcoded to check if any problems arose in the metabarcoding results (Suppl. material 2: Table S1, PCR settings upon request). All 24 O. rhenana barcode sequences were deposited in BOLD (sequence pages: TRBB001-20 - TRBB024-20).

The metabarcoding data were analysed with JAMP, version 0.69 (https://github.com/VascoElbrecht/JAMP), of the R software (R Core Team 2017). JAMP is an R-based metabarcoding pipeline that mostly relies on USEARCH 11.0.667 (Edgar 2010) and VSEARCH 2.10.4 (Rognes et al. 2016). First, the sequencing files were demultiplexed according to sample tags and quality-checked using FastQC. The paired-end reads were firstly merged and reverse complements were built whenever needed. Primers were trimmed with Cutadapt (Martin 2011) and sequences were filtered using a maximum expected error (max ee) filtering of 0.5 (Edgar and Flyvbjerg 2015), after which only those ± 10 bp of the expected length (421 bp) were retained. Given differences in sequencing depth, the number of reads of each sample was subset to the lowest number of reads found in our samples (547,000, site 1, 2^nd replicate). OTU clustering was done by USEARCH with a clustering threshold of 97% similarity and reads, including singletons, matching the OTUs to generate an OTU table. The OTUs with a read abundance below 0.01% per sample were discarded, while only the OTUs present in both replicates were retained for the downstream data analysis. Finally, a taxonomic assignment was performed using BOLD (Ratnasingham and Hebert 2007), accessed on 07-05-2019. Besides OTUs, haplotypes were inferred by following the denoising approach as implemented in Elbrecht et al. (2018) using the default settings.

Pairwise correlations (Spearman’s rank) were calculated to test if alpha diversity patterns were consistent amongst the identification levels (family, mixed, OTU and haplotype) and methods (morphology, metabarcoding). Organisms were identified to a mixed level for both identification methods and resulted in taxa being present at multiple levels of the taxonomic hierarchy. In our work, mixed identification levels were present because of immature/damaged specimens for morphological identification and because related sequences were assigned to different taxa with metabarcoding. This inconsistency can affect the overall result of the study because of ambiguous parent-child pairs being present (e.g. Baetis and Baetidae). We decided not to resolve ambiguous parent-child pairs because different methods can lead to different results (Cuffney et al. 2007). Taxonomic information was managed with the package ‘biomonitoR’ (https://github.com/alexology/biomonitoR).

Non-metric multidimensional scaling (nMDS) was performed with the metaMDS function of the R package ‘vegan’ (Oksanen et al. 2019) using different combinations of identification method and taxonomic level: i) mixed levels and family from the morphological identification using both abundance and presence/absence data; ii) mixed level and family from metabarcoding presence–absence data; iii) OTUs presence–absence data; iv) data on intraspecific diversity, i.e. haplotype presence–absence data. Jaccard and Bray-Curtis dissimilarities were used for presence–absence and abundance data, respectively. nMDS was followed by a PROTEST analysis which allows to test the non-randomness of the correlation between two ordinations that is bound between 0 (no fit) and 1 (perfect fit) (Peres-Neto and Jackson 2001). It allows to test if different taxonomic levels or identification methods provide similar community structure information. The turnover and nestedness components were then calculated to test if the estimated beta diversity components were affected by TS. Beta diversity partitioning was performed according to the approach proposed by Baselga (2010), where nestedness is calculated as the difference between total dissimilarity, estimated by the Sørensen index, and turnover, estimated by the Jaccard index. The contribution of turnover and nestedness was then calculated as percentages of total dissimilarity. The correlation of turnover percentage amongst taxonomic levels and identification methods was then assessed by a Mantel test based on Spearman’s rank correlations. The Mantel tests obtained from nestedness and turnover contributions were equal because both were calculated as the percentage of the total dissimilarity. The probability values for alpha diversity, PROTEST and the Mantel test were corrected for multiple comparisons with the Holm correction. The effect of phylogenetic structure on the partitioning of beta-diversity was evaluated by partitioning the PhyloSor index (Bryant et al. 2008) in turnover and nestedness components according to Leprieur et al. (2012). At first, OTU and haplotype sequences were aligned using MUSCLE (Edgar 2004). Maximum likelihood trees were constructed from OTU and haplotype sequences using RAxML version 8 (Stamatakis 2014) called inside R with the ‘ips’ package (Heibl 2008). Briefly, we used the general time reversible model (GTRCAT) and used the automated, rapid bootstrapping analysis (autoMRE) to search for the bestscoring ML tree in one program run. Beta diversity partitioning analyses were performed with the package betapart’ (Baselga et al. 2018). The package ggplot2 (Wickham 2009) was used to plot the results.

With the morphological approach 4050 organisms belonging to 69 taxa were identified, while 132 taxa, 229 OTUs and 513 haplotypes were detected by metabarcoding and bioinformatic protocols. According to morphological identification, the most abundant order was Ephemeroptera (1499 organisms), followed by Diptera (1157), Trichoptera (738), Trombidiformes (261), Plecoptera (165), Coleoptera (132), Odonata (56), Heteroptera (12) and other groups (30). Diptera and Ephemeroptera were the taxonomic groups that benefited most from the increase in the identification level in terms of richness (Fig. 2). Moreover, the abundance of the most abundant families of Diptera, Ephemeroptera, Trichoptera and Coleoptera was related to haplotype richness (Suppl. material 2: Fig. S1).

Figure 2.

Richness from family to haplotypes for the main macroinvertebrate groups found in this study.

Fewer taxa than the number of OTUs were found because several OTUs were assigned to the same taxon, while others were not assigned at all. Table 1 reports the number of taxa found at the different taxonomic levels by both the morphological and metabarcoding approaches. Metabarcoding missed some taxa identified by the morphological approach at both the family and genus level. Families Blephariceridae, Ceratopogonidae, Dugesiidae, Philopotamidae, Physidae and Scirtidae were detected only by the morphological approach, while Hydracarina families Lebertiidae, Mideopsidae and Tetrastemmatidae were detected only by metabarcoding. The mean family richness obtained by the morphological approach (18.2 ± 7.0, mean ± standard deviation) was not significantly different (paired t-test, t-value = -1.63, df = 11, p-value = 0.13) from that inferred via metabarcoding (17.2 ± 5.8) (Fig. 3). The overlap between the morphological and metabarcoding results, calculated as the ratio between the number of shared families and the number of total families per site, was 84.2% on average, with a minimum of 68.4% and a maximum of 100%.

Figure 3.

Family (a–b) and mixed taxonomic level (lowest possible level; c–d) richness for the morphological identification and metabarcoding in December 2017 (shortly after the dry season) and July 2018 (shortly before the next dry season). OTU and haplotype richness are also shown (e–f).

The genera Esolus (Coleoptera), Wormaldia (Trichoptera), Dugesia (Turbellaria) and Dinocras (Plecoptera) were not detected by metabarcoding. The sequences related to these genera were probably assigned to coarser taxonomic levels as our OTU table comprised unassigned Insecta, Coleoptera, Trichoptera, Turbellaria, Dryopidae and Perlidae (e.g. similarity to reference sequences > 10%). Metabarcoding originally identified the trichopteran species Psychomyia pusilla as Rhyacophila pusilla, which is only known from China (Malicky and Changhai 2002). This assignment proved to be a misclassification problem in the BOLD database and was corrected (records HMKKT101-10 and HMKKT369-10 on the BOLD website). The only mollusc taxon collected in the Trebbia River (Physa) was not detected by metabarcoding, probably because of a primer bias (Elbrecht and Leese 2015) and Turbellaria were assigned only at the class level, probably because suitable reference data were lacking.

Barcoding of single Oligoneuriella specimens confirmed both the morphological and metabarcoding identifications. For Italy, the family Oligoneuriidae is represented by the species Oligoneuriella rhenana. The best match of our specimens to the O. rhenana sequences found in BOLD was 91.79% similarity.

A high rank correlation between the morphological- and metabarcoding-based alpha diversity estimates was found for both the mixed (0.88, p < 0.01) and family (0.93, p < 0.001) level. The high correlation for the mixed taxonomic level was present despite the different taxonomic levels used with the two identification methods. In metabarcoding, a high correlation was found for the mixed taxonomic level with OTUs (0.80, p < 0.05) and haplotypes (0.80, p < 0.05). Non-significant correlations were observed for the family level identification with OTUs (0.51, p > 0.05) and haplotypes (0.50, p > 0.05). All the results are reported in Table 2.

The community structure assessed with nMDS was similar at all the taxonomic levels investigated by the morphological and metabarcoding approaches (Fig. 4, Suppl. material 2: Fig. S2). The visual similarity between the nMDS ordinations was also supported by the PROTEST analyses, which revealed significant comparisons for all the analysed pairs (Table 3). However, while the correlation was greater than 0.88 for different Linnaean taxonomic levels (range 0.88–0.98), the correlation between the OTUs/haplotypes and the Linnaean taxonomic levels dropped to range from 0.76 to 0.86. These findings indicate changes when using the highest, even intraspecific (haplotype), resolution. The July and December samples formed two separate clusters in all nMDS analyses. Amongst the December samples, intermittent sites St5 and St6 deviated from the bulk of perennial sites. This pattern was consistent amongst all the tested taxonomic levels and also between these levels and the haplotype level, but was less distinct for the ordinations obtained with morphological data, where St5 was less separated from perennial sites (Figure 4a, b). Distance amongst the July samples on the ordination space was shorter than for the December ones and both perennial and intermittent sites tended to cluster together. Once again, the nMDS performed on OTU/haplotype data provided a slightly different interpretation, because intermittent St6 was very much isolated from the other sites (Figure 4c, d). The metaMDS function used to perform nMDS rotates the configuration to maximise the variance of points on the first dimension, so it was possible to directly compare the patterns highlighted by different ordinations. When we took the investigated sites as a continuum along the river, the nMDS performed on the Linnaean taxonomic levels showed a gradient orientated in a bottom-left to top-right direction. This was not the case for the nMDS performed on OTUs, where the gradient was orientated in a bottom-up direction. This scenario highlights a temporal trend on the first axis and a spatial trend on the second one. The nMDS performed on the abundance data behaved similarly to that performed with the presence–absence data. However, the nMDS performed at the family level with the morphological approach identified a subgroup formed by the perennial sites located upstream of the water abstraction (St1 and St2) in July, while the others did not.

Figure 4.

nMDS ordination plots at different identification levels: family level for morphology (a) and metabarcoding (b) data, OTUs (c) and haplotypes (d). Stress was lower than 10% for all analyses. Numbers refer to sites: 1–2 upstream perennial, 3–4 downstream perennial and 5–6 downstream intermittent.

Lowering the taxonomic level had a huge impact on beta diversity, which increased on average when taxonomic resolution rose (Fig. 5, Suppl. material 2: Fig. S3). The percentages of the nestedness and turnover components of beta diversity reversed by going from a higher (family) to a lower (haplotype) identification level, at least for the December samples. When looking at the beta diversity partitioning calculated from the OTUs/haplotype data between sites and seasons, turnover came over as the most important component (Fig. 5, Suppl. material 2: Fig. S3). On the contrary, the contribution of nestedness varied across seasons. In December, nestedness was a consistent component for St6 (Fig. 5c) and for St4 and St5 to a lesser extent while in July it was generally negligible (Fig. 5d, Suppl. material 2: Fig. S3). The correlation of nestedness at different taxonomic levels decreased with increasing the distance amongst them (Table 4). Regarding the metabarcoding-generated data, good rank correlations were found between OTUs and the mixed taxonomic level (0.84, p-value < 0.01) and between the mixed taxonomic and family levels (0.76, p-value < 0.01), while the other correlations were always below 0.6 despite being significant. The inclusion of phylogenetic structure into the beta diversity partitioning analysis increased the amount of nestedness at the expenses of turnover (Suppl. material 2: Figs S4, S5).

Figure 5.

Values of the turnover and nestedness components for St1 and St6 in both December and July and for different identification levels (family, taxa = mixed taxonomic level) and methods (mor = morphological identification, meta = metabarcoding). The height of bars represents the average of the total dissimilarity between the target site and other sites, while the error bars the standard deviation.

The agreement between morphological and metabarcoding identifications is an important step towards using metabarcoding in community ecology. In our work, metabarcoding detected most of the taxa identified with morphological identification. Similar results were found by Elbrecht et al. (2017), Emilson et al. (2017) and Kuntke et al. (2020) for freshwater macroinvertebrates. The mismatches recorded between both methods may be due to the lack of reference sequences and misclassifications in BOLD and the primer bias. The lack of reference sequences in BOLD probably had the strongest impact on the results of this study because 29 OTUs were assigned to higher taxonomic levels (e.g. order or class). Taxa with no reference database entry cannot be assigned to species with molecular methods. The incompleteness and lack of accuracy in the reference databases have been advocated as the most important problem leading to mismatches between morphological and molecular method results (Carew et al. 2017, Pawlowski et al. 2018, McGee et al. 2019). Similar problems were found by Andújar et al. (2018) and Elbrecht et al. (2017) for macroinvertebrates. Incomplete reference databases have also been reported by several authors for different biotic groups, such as littoral hard bottom communities (Wangensteen et al. 2018), macrophytes (Kuzmina et al. 2018) and ichthyoplankton (Hubert et al. 2015). In our work, we likely deflated the rate of families detected by metabarcoding because we did not trust hits below 90% similarity, a very conservative approach in terms of assignment. The problem of unassigned taxa is expected to improve in the future due to ongoing barcoding activities and because molecular methods are being increasingly used in ecological studies (Curry et al. 2018, Weigand et al. 2019).

Lack of reference sequences from BOLD is more likely for taxa with complex diagnostic characters (e.g. Hydracarina or Turbellaria). However, our results showed that some well-studied insect taxa were not assigned to their respective higher taxonomic level. In this context, wider divergence, i.e. larger than 10%, is reported within mayfly genera (Lucentini et al. 2011, Morinière et al. 2017). Therefore, the presence of morphologically cryptic taxa with great phylogenetic divergence and the use of fixed thresholds to define taxonomic levels (Brown et al. 2015), could explain our findings. In our work, some morphologically identified taxa included potential cryptic species (Jackson et al. 2014), as seen for instance for the mayfly O. rhenana. This species, the only Oligoneuriidae species reported for Italy, was amongst the most abundant ephemeropteran in our samples in July. However, this was detected only by metabarcoding at the family (OTU_112) and order (OTU_32) levels. The barcoding of O. rhenana individuals confirmed the metabarcoding results, indicating the presence of a potential cryptic species. Furthermore, the employed fixed threshold could have impeded the assignment of Oligoneuriella OTUs at the species or genus level as the mean intra- and interspecies distances of the O. rhenana sequences retrieved from BOLD were 12% and 15%, respectively. The herein applied use of a 95% threshold to assign genera and of 90% to assign families might not suffice, at least for Oligoneuriella/Oligoneuriidae. Despite more research being necessary to define the taxonomic status of this taxon (e.g. by analysing nuclear markers, morphological identification based on adult specimens, comparisons to morphology and genetic markers of type material), the implications of these findings are worth discussing. Metabarcoding can in fact help to identify inconsistencies in the current morphological taxonomy and to improve our freshwater fauna knowledge (Cardoni et al. 2015). Detecting cryptic species is also important, because they may have different environmental requirements (Eisenring et al. 2016, Leys et al. 2016), contrasting responses to multiple stressors (Sturmbauer et al. 1999, Macher et al. 2016, Beermann et al. 2018) and might provide valuable information on biogeographical processes.

We demonstrated that the increased taxonomic resolution provided by metabarcoding, even down to intraspecific diversity (haplotypes), is useful for interpreting macroinvertebrate community patterns as highlighted by alpha diversity, nMDS and, in particular, turnover/nestedness results. We found good concordance amongst the alpha diversity estimates for both the morphology and metabarcoding identification methods at the family and mixed levels (down to species, whenever possible). Concordance declined when comparing OTUs, or even haplotypes, with higher taxonomic levels. In this context, we can expect good concordance between different identification levels when the lower to higher taxonomic level ratio is low (Gaston 2000a). This happens, for example, when the number of species within each family is small and thus the ratio approaches 1. The use of molecular data is, in turn, expected to exacerbate this situation given the increased ability to assess biodiversity as highlighted in our work, where the number of entities to be studied (richness) increased by a factor of 10 when moving from the family level to OTUs/haplotypes.

Our results showed that the correlation amongst the ordinations performed with the data at different Linnean taxonomic levels for both morphology and metabarcoding was high (> 0.88). Good concordance between community patterns at different taxonomic levels has been found for arthropods (Timms et al. 2013), phytoplankton (Carneiro et al. 2010) and spiders (Cardoso et al. 2004). For macroinvertebrates, Slimani et al. (2019) reported a marked concordance of the community structure at the family and genus levels with those obtained at species level for different insect orders. Ordinations were also consistent when comparing abundance and presence-absence data, despite minor differences arising in ecological interpretations. In our work, the good correlation between presence–absence and abundance ordination was likely due to the transformation of the abundances made prior to running the nMDS. Transformation of abundances is the default option in the vegan package when wide abundance variability is present and is generally performed to reduce the effects of abundant taxa on the ordination result. In line with this, Heino (2008) found good concordance between abundance and presence-absence matrices only after data transformation. Despite the good concordance with the Linnaean taxonomic levels, the ordinations performed with OTUs/haplotypes were effective in tracking temporal and spatial trends in community structure and in detecting the intermittency effect in December. The clearer clustering of temporal trends with metabarcoding was probably due to the increased taxonomic resolution of OTUs, which led to more taxa being exclusive to only the first or second season.

Morphological identification and metabarcoding gave contrasting responses when looking at beta diversity partitioning. This was especially true in December, when turnover turned out to be the most important fraction for OTU/haplotype level, while higher taxonomic levels highlighted the opposite pattern. In our work, not all organisms were identified at species level by either morphological identification or metabarcoding. This is a potential limitation of our study, as beta diversity dynamics at small spatial scale should be better highlighted by using finer taxonomic levels. The use of OTUs or inferred haplotypes allowed us to circumvent this problem by disengaging ecological patterns from the Linnean taxonomy. In this context, recent findings have shown that inferences in beta diversity dynamics based on genetic methods, performed better than those based on morphological identifications (Gill et al. 2012, Vodă et al. 2015, Bringloe et al. 2016). Including the phylogenetic information into the beta diversity partitioning analysis led to a decrease of the turnover percentage, meaning that OTUs and haplotypes found at intermittent sites were phylogenetically related to those of perennial sites. However, the inclusion of phylogenetic information lacks a direct link with desiccation, because strictly related species or genera can have different responses to drying (Suemoto et al. 2004, Gusev et al. 2014).

The increasing availability of data at higher taxonomic resolutions poses questions about the taxonomic level needed to study key ecological processes. At first, our results showed how an increased resolution was beneficial for the most abundant orders found with morphological identification. This was particularly evident for haplotypes, to which Diptera and Ephemeroptera contributed 82%. Haplotype richness could thus reflect both the abundance of organisms of a target group, the presence of bottleneck or founder effects, as well as the species-specific variability of the COI fragment tested in this study. Moreover, the use of OTUs or haplotypes instead of the family level data resulted in different beta diversity partitioning results. However, a finer taxonomic resolution might not always be beneficial for inferring ecological patterns from monitoring data. For example, Martin et al. (2016) demonstrated that the explanatory power of metacommunity models increased with coarse taxonomic resolution because of the assumed ecological interchangeability of species within genera or families. Moreover, the increase in turnover at the expense of nestedness when finer taxonomic levels were used could be related to the better chance of finding rarer taxa when increasing alpha diversity, rather than to structural patterns governed by flow intermittency. Additionally, the effect of a low sampling effort increases with increasing taxonomic resolution, which advises caution when interpreting nestedness/turnover data based on haplotypes inferred from very few individuals. Thus, the ideal taxonomic level for studying macroinvertebrate dynamics in river networks obviously depends on the research question. However, we argue that emerging genetic techniques have the potential to allow a more complete description of the biological connectivity amongst sites, which will refine our knowledge of species-specific autecological preferences. The advantage of genetic data is that a lower resolution can be easily obtained by lumping taxa into higher ranks (family, order etc.).

Recolonisation patterns of macroinvertebrates in IRES are driven by different factors, amongst which distance of a given intermittent site from the nearest perennial site and length and predictability of the dry period are the most important (Datry et al. 2017). After a dry phase, reduced alpha diversity at different taxonomic levels was evidenced for intermittent sites (Fig. 2). This trend was not observed in July, which denotes the resilience of the macroinvertebrate community to lack of water. Similarly to our findings, the recovery of the macroinvertebrate community to pre-drought levels was found by Skoulikidis et al. (2011) in an artificial intermittent stream in Greece. Datry (2012) reported similar results in the Albarine River, although the degree of recovery depended on the distance from the nearest perennial source.

nMDS results supported conclusions drawn from alpha diversity patterns, which revealed greater dispersion for December (post-drying) samples than for July (pre-drying) ones. The greater dispersion of December samples was driven by intermittent sites, which could host a nested subset of communities at perennial sites or a completely different pool of taxa. Community nestedness is expected to be dominant along flow intermittency gradients in IRES (Datry et al. 2014). In this context, our results revealed the opposite pattern for both sampling seasons. Given the distance from the nearest perennial source and the absence of supra-seasonal dry periods, the observed pattern was expected in July as 9 months had elapsed since flow resumption, which thus allowed complete recolonisation or emergence from riverbed by macroinvertebrates at intermittent sites. Recovery of the macroinvertebrate community to the pre-drying conditions is relatively fast in IRES (Acuña et al. 2005, Datry 2012). On the contrary, the predominance of turnover over nestedness upon rewetting was not expected and was clearly detectable only with the metabarcoding data. Recolonisation after rewetting was herein driven by a mixture of physiological adaptations to drying, dispersal abilities and behavioural adaptations. Organisms with desiccation resistance stages or with excellent dispersal abilities are favoured when the distance between perennial and intermittent sites is long, whereas dispersal of organisms is not limited when they are in close proximity (Bogan et al. 2013). Dispersal limitation was likely a minor factor in our study because of the short intermittent stretch and the maximum distance of 2.5 km to the nearest perennial source. After rewetting, however, 26% of OTUs were exclusive to intermittent sites, which suggest recolonisation from the Po River or from organisms with physiological or behavioural adaptations to intermittency. Intermittent sites hosted 16 exclusive OTUs, with eight OTUs of Chironomidae, three OTUs of Simuliidae, one OTU of Trombidiformes, one OTU of Calopteryx splendens and one OTU of Cloeon dipterum. Chironomidae can survive dry periods in desiccation-resistant stages (Frouz et al. 2003) or by using the hyporheic zone as a refugium (Stubbington and Datry 2013). Amongst the Chironomidae species exclusive for intermittent sites after re-wetting, Cricotopus bicinctus has often been found in both lentic and lotic temporary water bodies (Boix et al. 2001). Simuliidae larvae and their eggs are generally sensitive to desiccation (Adler and McCreadie 2002) and the presence of OTUs exclusive of intermittent sites is probably due to either the dispersal of adults and larvae or their resistance in the hyporheic zone (Lencioni and Spitale 2015). However, see Bogan et al. (2013) for exceptions. Trombidiformes present resistant stages, as rehydration experiments have clearly shown (Stubbington and Datry 2013), but also a high drift propensity (Imbert and Perry 2000, Fenoglio et al. 2004). As C. splendens is not resistant to desiccation, its presence at intermittent sites was probably due to drift from upstream sites or from eggs hatched from late oviposition (Gallesi and Sacchi 2019). C. dipterum was present only at site 6 in December. Thus, its presence could be attributed to the spatial proximity of this site to the Po River confluence or to other perennial water sources in the nearby surroundings (e.g. drainage channels). C. dipterum is ovoviviparous and lays eggs containing mature embryos that hatch immediately after oviposition, which thus allows impacted sites to be rapidly recolonised (Degrange 1959).

Despite studying only few sites and limited repeated sampling over time, our work provided interesting insights into the use of metabarcoding to study ecological processes in IRES. At first, community metabarcoding yielded reliable information about taxonomic composition at all sites. In the Trebbia River, community dynamics inferred from metabarcoding data can be used alternatively to those derived from morphological identification when family/mixed identification levels and presence–absence information are needed. Second, we demonstrated that an increase in the identification level only moderately affects alpha diversity and community structure patterns, while it strongly affects beta diversity partitioning. The latter showed a switch from the nestedness to the turnover component moving from family to OTUs/haplotypes. Finally, an increased taxonomic resolution to species level using metabarcoding data was beneficial to interpret the observed patterns. Our work highlights how inferences drawn from highly resolved identification levels, even down to haplotypes, can provide a more comprehensive picture about the responses of biological communities to flow intermittency. This is a key point for IRES, where detailed information about recolonisation processes are needed to detect ongoing dynamics and for planning optimal conservation strategies.