Fresh insights into Mediterranean biodiversity: environmental DNA reveals spatio-temporal patterns of stream invertebrate communities on Sicily

The Mediterranean region with its islands is among the top biodiversity hotspots. It houses numerous freshwater taxa with a high rate of endemism, but is heavily impacted by anthropogenic pressures and global climate change. To conserve biodiversity, reliable data on species and genetic diversity are needed especially for the scarcely known insular freshwater ecosystems. Environmental DNA (eDNA) metabarcoding provides a straight-forward opportunity to assess aquatic biodiversity. Therefore, we conducted the first eDNA metabarcoding study in one stream catchment on Sicily. Specifically, we aimed to (i) investigate spatial diversity patterns of macroinvertebrate communities, (ii) assess seasonal changes (autumn and winter), and (iii) check if dispersal barriers can be identified. Water samples were taken at 27 different sites in two seasons and eDNA metabarcoding was performed using a fragment of the mitochondrial cytochrome c oxidase subunit I gene as a marker. In total, we detected 98 macroinvertebrate species, including 28 taxa potentially new to Sicily. Exact sequence variant and species composition data showed that diversity differed between seasons with less taxa detected in winter. We also detected a dispersal barrier, which had a stronger effect in autumn. Our findings show that eDNA metabarcoding provides valuable information on Sicilian freshwater biodiversity. We therefore encourage its application for understudied regions to better understand the state and dynamics of freshwater biodiversity.


Introduction
Freshwater ecosystems are hotspots of biodiversity (Strayer & Dudgeon, 2010). Even though they only comprise roughly 1% of global land surface area (Dudgeon et al., 2006;Strayer & Dudgeon, 2010), freshwater ecosystems host nearly 9.5% of Earth's described animal species (Balian et al., 2007). However, degradation of freshwater ecosystems is a global phenomenon, affecting in particular riverine habitats (Vörösmarty et al., 2010;Grill et al., 2019). Anthropogenic pressures are superimposed by the effects of global climate change, in particular heat waves, floods and droughts (Arnell, 1999;Vörösmarty et al., 2010). As a consequence, freshwater faunal biodiversity is declining, making rivers and lakes the most endangered ecosystems in the world (Dudgeon et al., 2006;Reid et al., 2019;Almond et al., 2020). This holds true in particular for the Mediterranean region (De Figueroa et al., 2013), which is considered as one of the top biodiversity hotspots in the world (Médail & Quézel, 1999;Myers et al., 2000). Here, diversity is especially high in freshwater ecosystems, housing more than 6% of the global freshwater biodiversity, with particularly high numbers of endemic and rare taxa found on Mediterranean islands (De Figueroa et al., 2013). Yet, islands are also among the most threatened ecosystems within the region, with a predominance of non-perennial rivers and streams particularly vulnerable to anthropogenic impacts (Hopkins, 2002;Skoulikidis et al., 2017).
Freshwater invertebrates, and in particular aquatic macroinvertebrates, are biological key components known to shape local communities, and are prime indicators of water quality. This also seems to be true in the streams on Mediterranean islands, where local macroinvertebrates are recognised as bioindicators of ecological disturbance (Feio et al., 2014;García et al., 2014;Erba et al., 2015). However, evidence about drivers influencing the observed insular aquatic macroinvertebrate diversity over time is scarce and restricted to a limited area (Garcia et al., 2017;Lobera et al., 2019). Studying their local spatio-temporal patterns of distribution is crucial for understanding the state and dynamics of freshwater ecosystems (Skoulikidis et al., 2017). Given that species diversity is high but taxonomic expertise, in particular determination keys, is often not available, biological surveys and bioindication often rely on data from higher taxonomic levels only, typically family level. While this level increases reliability of the data, it also misses out the true diversity. DNA metabarcoding provides new opportunities to study diversity of freshwater communities (Carew et al., 2013;Elbrecht et al., 2017). In particular, DNA metabarcoding of environmental DNA (eDNA) provides a great tool for studying aquatic biodiversity in a non-invasive and time-efficient way (Rees et al., 2014;Goldberg et al., 2015;Deiner et al., 2016). While most aquatic eDNA surveys have focused on vertebrate species (Hänfling et al., 2016;Harper et al., 2018;Closek et al., 2019), recent studies have shown that a lot of information can also be obtained for invertebrate species (Macher et al., 2018;Mächler et al., 2019;Leese et al., 2021), despite the problem that the template molecules from the target taxa could be 'watered down' (Hajibabaei et al., 2019). Getting a broad picture of local macroinvertebrate communities from multiple fine-scale localities and from different time points within a single river system without investing a large amount of time in taxonomic identification could be of paramount importance for a better understanding of the observed diversity patterns as well as for its management (Stubbington et al., 2018).
Knowledge about the diversity and distribution of freshwater macroinvertebrates in the Mediterranean region is constantly growing, both due to multinational (e.g. Fauna Europaea, GBIF) and local initiatives (Ruffo & Stoch, 2006), but is still missing for many regions. While the accessibility of reference DNA barcodes for macroinvertebrates is growing, there is still a significant gap present (Curry et al., 2018;Weigand et al., 2019). This holds true in particular for Mediterranean countries for which operational routine monitoring taxa lists often use higher taxonomic levels than species. Thus, given the incompleteness of the reference databases, taxonomic assignment of obtained DNA information can be troublesome. In taxonomic groups with even scarcer reference DNA sequences available, like diatoms and other unicellular eukaryotes or meiofauna, where the majority of obtained information cannot be assigned to species level, there is a notion of using molecular operational taxonomic units (OTUs) (Westcott & Schloss, 2015;Feio et al., 2020) or exact sequence variants (ESVs) (Tapolczai et al., 2019(Tapolczai et al., , 2021 for analyses of diversity patterns. By often providing more ecologically relevant information (Zizka et al., 2020;Tapolczai et al., 2021) and higher reusability compared to OTUs (Schmidt et al., 2015), ESV-based approaches could provide a promising solution for studying remote and understudied regions, where much of the observed molecular diversity cannot be assigned to species level, yet.
In this study, we conducted the first eDNA survey for a stream ecosystem on Sicily. Sicily is the largest Mediterranean island and is located in the central part of the Mediterranean basin. The freshwater macroinvertebrate fauna is relatively well documented with approximately 1300 species being reported from the insular freshwater ecosystems (Ruffo & Stoch, 2006). The freshwater ecosystems on the island exhibit a comparatively high level of local endemism (Stoch, 2000). In some groups, like Malacostraca, even more than 50% of the taxa are endemic (Hupało et al., 2021). However, still very little is known about Sicilian freshwater macroinvertebrate communities and their local spatio-temporal distribution and dynamics. By focussing on a single river system in the central part of Sicily, we aimed to provide insights into (i) the general biodiversity patterns of macroinvertebrate communities, (ii) the turnover of local macroinvertebrate communities between autumn and winter, and (iii) whether dispersal in those communities is hindered by a barrier in the riverscape. Additionally, our aim was to compare patterns inferred from a species-based to an ESV-based approach. With that, we aimed to evaluate if eDNA metabarcoding data can provide valuable information for assessment and monitoring of Mediterranean freshwater biodiversity, even when the DNA reference databases are incomplete.

Field sampling
Water samples were collected during two sampling campaigns conducted in September 2019 (autumn) and January 2020 (winter) in the Fosso del Tempio river system on Sicily. The studied river system belongs to the basin of Fiume dei Monaci, which covers approximately 590 km 2 and belongs to the Fiume Simeto catchment. The studied river originates from the slopes of Monte Moliano e Montagna on the border of the territory of the Municipalities of Aidone and Piazza Armerina. The main stretch of Fosso del Tempio is approximately 30 km long and after the first main confluence, it changes its name to Fosso Pietrarossa. To our knowledge there is no information regarding the average discharge of Fosso del Tempio.
A total of 27 sites were visited, with 18 being sampled both seasons, partially due to intermittency of parts of the studied river system, resulting in two sites (A6, E9) being completely dry during autumn sampling ( Fig. 1, Table S1). To more broadly assess regional diversity, nine additional sites from neighbouring catchments were sampled at least in one season (Fig. 1). The specific sampling sites were initially chosen for another study focusing on cryptic species coexistence of amphipod taxa formerly discovered at one of the sites (E3). For sampling, a sterile 1 l bottle was used. Depending on the size of the stream, the water was taken either from the water surface at a single spot (1 l), or from each side of the stream (2 9 0.5 l). Water samples were taken against the flow direction wearing single-use nitrile gloves. Subsequently, samples were filtered on-site using a vacuum pump (Druck-Vakuumpumpe Typ:00A7.400 12 V DC, Industrievertretung Neubauer) and a nitrocellulose filter (0.45 lm pore size, VWR, 513-1454). At each filtration location, a field blank (left open on the equipment for 30 s) was included and further processed with the other samples. If weather conditions did not allow for on-site filtration, samples were stored in a cooling box (4°C) and filtered later during the same day (maximum 5 h after sampling). If the filter clogged during filtration, a second or third filter was used (Table S1). After filtration, the filters were folded with sterile tweezers and transferred into 1.5 ml Eppendorf tubes filled with 96% technical ethanol. Filters were immediately put into transportable ice boxes (\ 0°C) and within a few hours after sampling stored at -20°C until DNA extraction.

DNA extraction and purification
All laboratory work was conducted under sterile conditions in a specialised laboratory ('eDNA lab') that is only used for work on eDNA, regularly sterilised using UV light and where the air is filtered with a HEPA 13 filter. Full body protective clothing was used (disposable overalls, facemasks, gloves, shoe covers) and the laboratory was irradiated with UV light prior to all working steps. Handling of the samples (DNA extraction; PCR) was conducted under independent UV hoods. First, the filters were taken from the ethanol and laid out in sterile petri dishes to dry for 15 to 20 h. The petri dishes were covered with aluminium foil to protect them from UV irradiation. After drying, the filters were ripped into 6 to 7 pieces and transferred into 2 ml Eppendorf tubes containing 610 ll of lysis buffer (600 ll TNES and 10 ll Proteinase K (10 mg/ml)). DNA extraction was conducted using a modified salt extraction protocol as described in Weiss and Leese (2016). Around 10 to 18 samples were extracted together. The field blanks were randomly distributed among sample batches and treated like normal samples in extractions.
Prior to DNA purification, 1 ll RNase (10 mg/ml, Thermoscientific, EN0531) was added to each sample and incubated at 37°C for 30 min to remove RNA. After incubation, DNA was purified using the MinE-luteÒ Reaction Cleanup Kit (QIAGEN). All eDNA samples, where more than one filter had been used (Table S1), were pooled during purification by loading them onto the same column. After DNA purification, 68 DNA samples were retained for further processing (including negative filters). Purified DNA was resuspended in 20 ll nuclease-free H 2 O and stored at -20°C. The DNA concentration was measured using the Qubit 2.0 (dsDNA High Sensitivity Array, Thermo Fisher Scientific, Beverly, USA).

PCR amplification
A 205 bp long fragment of the cytochrome c oxidase subunit 1 (COI) mitochondrial gene was amplified using the fwh2 primer set (fwhF2 ? fwhR2n)  and the Multiplex PCR Plus Kit (100) (QIAGEN). A two-step PCR protocol was employed, the first step with the untagged fwh2 primers and the second step with indexing primers (Leese et al., 2021).
To account for PCR stochasticity and to increase consistency of the results, two PCR replicates were processed for each sample (Zizka et al., 2019). In the first step, 1 ll of DNA was amplified in a 25 ll reaction (12.5 ll of the Multiplex Mastermix, 2.5 ll Color Dye, 0.25 ll fwhF2 (10 lM), 0.25 ll fwhR2n (10 lM) and 8.5 ll molecular grade H 2 O). The cycling conditions for the first step were split into two parts, starting with 5 min at 95°C for initial denaturation, followed by 10 cycles of a temperature step-down PCR (-1°C per cycle) lasting 90 s going from 68 to 58°C, followed by final elongation for 30 s at 72°C. The second part consisted of 25 cycles of 95°C for 30 s, 58°C for 90 s, and 72°for 30 s. Final elongation was conducted for 10 min at 68°C.
In the second PCR step, indexing primers were added. Each sample and each replicate were given a unique primer combination via a 8-bp Illumina index added to the PCR amplification primer to bioinformatically assign sequences to their original sample after sequencing (Table S2). The second PCR step was conducted with 1 ll of the product from the first PCR step, 12.5 ll of the Multiplex Mastermix, 2.5 ll Coral Color Dye, 0.25 ll forward indexing primer (10 lM), with the potential in-stream barrier in autumn and in winter 0.25 ll reverse indexing primer (10 lM) and 8.5 ll PCR H 2 O. Cycling conditions were as follows: 95°C for 5 min for initial denaturation, then 20 cycles of 95°C at 30 s and 72°C for 120 s with a final elongation at 68°C for 10 min.
The PCR success was checked using a 1% agarose gel and the DNA concentration was measured with both the Qubit (2.0) High Sensitivity kit and the Fragment Analyzer TM Automated CE System (Advanced Analytical Technologies GmbH) using the NGS Standard Sensitivity Kit, before equimolar pooling (SequalPrep Normalization Plate, Applied Biosystems, Foster City, CA, USA). Residual primers were removed with a left sided SPRIselect size selection (Beckman Coulter, ratio: 0.76x). The final library was sent to Macrogen Korea and paired-end sequenced on two HiSeq X Illumina lanes (read length 2 9 150 bp). To increase sequence diversity, 5% PhiX was added by the sequencing company.

Bioinformatic analyses
Raw reads for the two libraries were received as demultiplexed fastq files. The sequencing data was processed using a pre-release version of the graphicaluser interface pipeline MetaProcessor (available at https://github.com/TillMacher/MetaProcessor). First, the quality of the raw reads was checked using FastQC (Andrews, 2010). Subsequently, samples were renamed using a custom python script (Metaprocessor: rename_samples.py). Paired-end reads were merged using VSEARCH version 2.11.1 (Rognes et al., 2016), allowing for 25% differences between merged pairs and a minimum overlap of 20 bp. Afterwards, primers were trimmed using cutadapt version 2.8 (Martin, 2011), using the linked adapter option without anchoring. Reads were then filtered by length (195-215 bp threshold for fwh2 target fragment) and by maximum expected error (maxee = 1), using VSEARCH. The filtered reads were dereplicated with singletons and chimeras removed with VSEARCH. All reads were then pooled using a custom python script and globally dereplicated (MetaProcessor: v_derep_singletons_uchime.py). Sequences were denoised into ESVs, using VSEARCH's '-clus-ter_unoise' function, with a minimum size of 8 and an alpha value of 2. The ESVs were remapped (use-arch_global function, 100% similarity) to the individual sample files to create the read table. Taxonomic assignment of ESVs was conducted using BOLDigger version 1.1.10 (Buchner & Leese, 2020) and the Barcode of Life data system (BOLD) database (Ratnasingham & Hebert, 2007). The option ''JAMP filter'' was applied to extract the final taxonomy table.
Downstream processing of the dataset was conducted using TaxonTableTools (TTT) version 1.3.0 (Macher et al., 2021a). First, the full taxonomy table was processed for replicate consistency, where only ESVs present in both replicates were retained by using 'merge replicates' option. Afterwards, the resulting TaXon table was filtered using a read-based filter with 0.0005% threshold filtering. Using such a low threshold value was possible due to the high sequencing depth obtained combined with the dual indexing strategy (same index in i7 and i5 read) and lack of reads in negative controls, ensuring high reliability of data with preserving the records of rare taxa observed. The final filtering step included removal of taxonomically unmatched ESVs. General patterns of diversity were studied based on read proportion pie charts generated on phylum and order level. The resulting charts were superimposed on circular neighbor-joining phylogenetic trees based on p-distance generated with MEGA7 software (Kumar et al., 2016). To investigate seasonal patterns, the TaXon tables including autumn and winter samples were compared by generating Venn diagrams. In order to evaluate the effect of the potential in-stream barrier, the TaXon tables were first converted to presenceabsence (p/a) tables and subsequently, beta-diversity heatmaps and 3D Principal Coordinate Analysis (PCoA) charts were generated and the respective ANOSIM (analysis of similarities) R values were calculated. In principle, the ANOSIM is an analysis of variance using the average measures of dissimilarity (in this case Jaccard distances) in species composition, comparing them between and within given samples. Resulting test statistic value R usually ranges between 0 and 1, where values close to 1 generally support the sample dissimilarity contrary to values close to 0 which indicate no difference (Clarke, 1993;Chapman & Underwood, 1999). Finally, to visualise species distribution, Site Occupancy plots and Parallel Categories (ParCat) diagrams were produced.

Datasets for analyses
The final dataset comprising all ESVs retained after quality filtering ('All ESVs' dataset) was subsequently filtered for ESVs assigned to benthic macroinvertebrates (BMI) which are recognised as freshwater taxa, including ones which are considered amphibiotic with a larval stage being confined to freshwater habitats ('BMI ESVs' dataset). Taxonomic filtering was performed according to taxa information stored in the Global Biodiversity Information Facility (GBIF; https://www.gbif.org) and World Register of Marine Species (WoRMS; Horton et al. 2021) databases. Freshwater macroinvertebrate taxa identified as potentially new for Sicily were cross-validated with the publicly available checklist of the Italian fauna (Latella et al., 2007) as well as species' distribution maps in GBIF database (which, however, may also contain some unvalidated entries). Out of all ESVs assigned to benthic macroinvertebrates, ESVs belonging to the same species were grouped together and represented by a single entry on species level for the final dataset ('BMI species' dataset). Furthermore, depending on the research question, the number of sampling sites for certain analyses varied. Thus, for analyses of the seasonal differences, only the sites, which were sampled both seasons were used (= 18 sites; Fig. 1). For analyses of possible barrier effects, only the sites from the Fosso del Tempio systems were used (highlighted in green and violet in Fig. 1), regardless if the site was visited both seasons or not, since the data for the barrier effect were analysed separately for each season. The downstream dataset consisted of nine sites of the Fosso del Tempio system, which were sampled in both seasons. The upstream dataset varied between seasons with seven sites being included in autumn and six in winter. There were five core upstream sites included in both seasons with three sites being sampled only in one season (A11, A18 in autumn and E9 in winter; Fig. 1). The downstream site A6, sampled only in winter, was excluded from the dataset given very low observed ESV/species numbers, pointing to potential strong effects of nonperennial flow, which could introduce a strong bias in the barrier analysis. All filtering steps described above were conducted using taxon-based and sample-based filtering options in TTT software.

Results
We obtained 445,683,202 raw read pairs, which also included reads of a second library with a different primer that was run in parallel on the same sequencing run. After paired-end merging, 232,498,448 reads were retrieved from fwhF2 ? fwhR2n primer pair.
After the final quality filtering step 120,167,492 raw reads remained (31,900 reads assigned to negative controls), which were denoised into 17,859 ESVs out of which 17,409 could be taxonomically assigned. The subsequent merging of the PCR replicates retained 110,605,585 reads (no reads present in negative controls). The final read table (corresponding to 'All ESVs' dataset; Table S3) obtained after 0.0005% threshold filtering and removal of taxonomically unmatched ESVs, consisted of 110,148,158 reads, which were assigned to 7331 ESVs. Only around 6.5% (= 474) of all ESVs could be assigned to species level, resulting in 340 unique species. Taxonomic filtering for freshwater macroinvertebrates resulted in a read table consisting of 14,868,364 reads, assigned to 466 ESVs ('BMI ESVs' dataset; Table S4). Nearly 36% of macroinvertebrate ESVs (= 167) could be assigned to species level with 98 unique freshwater macroinvertebrate species retrieved.

General patterns of diversity
The ESVs belonging to the 'All ESVs' dataset could be assigned to 29 phyla, with the highest proportion of reads as well as the highest number of ESVs being assigned to Arthropoda (Fig. 2). Apart from Arthropoda, high numbers of ESVs were retrieved for Heterokontophyta, Ochrophyta, and Bacillariophyta. Altogether, above 40% of reads were assigned to nonmetazoan taxa (Fig. 2). The majority of all species were assigned to Arthropoda, with Heterokontophyta and Annelida being the only other phyla with more than 10 species retrieved. Less than 0.5% of all reads consisted of 134 ESVs belonging to 15 phyla, both metazoan and non-metazoan.
The ESVs of the 'BMI ESVs' dataset were assigned to 14 orders, with nearly 65% of reads belonging to dipteran ESVs, which were also characterised by highest ESV and species diversity (Fig. 2). Next most read-abundant macroinvertebrate orders were amphipods and haplotaxidan annelids. More than 4% of reads were assigned to 102 poriferan ESVs, none of which could be assigned to species level with four being assigned to the orders Bubarida and Desmacellida, respectively (Table S4). Nearly 10% of all macroinvertebrate ESVs (22 species of five orders) comprised less than 1% of the reads.

Temporal turnover
The number of ESVs or species detected was higher in autumn than in winter regardless of the dataset, with a minor fraction of shared species/ESVs (Fig. 3). In the 'All ESVs' dataset, 4698 ESVs were retrieved in autumn and 2554 in winter, with 1288 ESVs being Fig. 2 Read proportion charts with taxonomic composition, the number of ESVs assigned to a certain taxa group (phylum level for all ESVs; order level for freshwater macroinvertebrates*) and the number of unique species assigned to ESVs. Upper graph shows the read proportions for the entire dataset and lower graph read proportions in the freshwater macroinvertebrate dataset. *Since most ESVs assigned to Porifera could not be assigned to order level, they were kept to phylum level respectively shared between seasons. In the 'BMI ESVs' dataset, 337 ESVs and 170 ESVs were retrieved in autumn and winter, respectively, with 118 shared between seasons. In the 'BMI species' dataset, 72 species were detected in autumn and 37 in winter, with 29 species found in both seasons. Less than half of the ESVs and macroinvertebrate species were found in winter compared to autumn, regardless of the dataset composition. The proportion of shared entities varied between approaches. In general, the percentage of shared diversity was inversely proportional to the broadness of the dataset, regardless of season analysed. In autumn, the proportion of diversity shared varied between 27.4% for 'All ESVs' to 40.3% for 'BMI species', whereas in winter it ranged from 50.4% for 'All ESVs' to 78.4% for 'BMI species'.

Barrier effect
We found strong community structuring across the barrier in autumn for all three datasets analysed (Fig. 4). ANOSIM R values (measure of average dissimilarity) were all highly significant (P \ 0.001) and ranged from 0.427 in 'All ESVs' dataset to 0.593 in 'BMI ESVs' dataset. In contrast, no significant differences were found for any of the datasets in winter, where ANOSIM R values varied from 0.088 in 'BMI species' dataset up to 0.157 in 'BMI ESVs' dataset (P \ 0.05). This pattern was also reflected by the average Jaccard distances obtained per site (Fig. 4A, Table S5) as they were generally higher between sites across the barrier than on each side of the barrier regardless of season, except for upstream sites during winter. Considering individual sites, the Jaccard distances were on average higher across the barrier than within a stream section, with exception of two sites in autumn (A6F in 'BMI ESVs' and A1 in 'All ESVs' datasets) and 8 sites in winter (E3, E4, E5, A6F, E8, E9, A7 and A15 across various datasets; Table S5). However, the differences in the Jaccard distances between stream sections compared to sites within a stream section were more pronounced in autumn compared to winter. These results are further supported by the PCoA analyses performed on all sites within a stream section. In autumn clear separation between samples from downstream and upstream can be observed, when in winter there was no clear distinction between samples from the opposite sides of the in-stream barrier (Fig. 4B). Moreover, the direct comparison of sites closest to the in-stream barrier (E5 from the downstream and E8 from the upstream) indicate the proportion of diversity shared between those sites varies between seasons. In autumn, Jaccard distances ranged from 0.89 to 0.91, whereas in winter the values were between 0.67 and 0.76. The value of 1 obtained for those sites when using 'BMI species' dataset in winter should be treated with caution, since it is a result of no ESVs assigned to species level on site E5. A similar pattern can be observed in the proportion of shared diversity which ranges from 9% in case of 'BMI ESVs' to 10.5% for 'BMI species' whereas in winter it ranges from 23.5% in case of 'BMI ESVs' to 32.6% in 'All ESVs', with the exception of BMI species where none of seven ESV could be assigned to species level in site E5 (Fig. S1).

Macroinvertebrate species composition and distribution
Apart from looking at the diversity statistics, we also analysed the observed macroinvertebrate species composition as well as the seasonal turnover, including the potential effect of the barrier on species diversity. The distribution of the 98 freshwater macroinvertebrate species reported in the 'BMI species' dataset differed between seasons (Fig. 5A)   Fig. 3 Venn diagrams showing number of species or ESVs detected in the autumn and winter samples as well as the fraction of shared species/ESVs between seasons for the three datasets analysed and stream sections (Fig. 5B, C). Overall 32 of the 98 freshwater species were shared between Fosso del Tempio and neighbouring streams, 52 species were exclusive to Fosso del Tempio and 14 to the neighbouring streams (Fig. S2). We identified 28 species potentially new for Sicily: 23 dipteran, 3 annelid, 1 versions available as html files in Supplementary Material. The upper graphs (orange box) present the results for autumn, the lower graphs (blue box) present the results for winter. The color code corresponds to the one presented in Fig. 1: violet = downstream sites, green = upstream sites caddisfly and 1 water beetle species, respectively (Table S6). Ten of them were found exclusively in the neighbouring catchments.

Temporal turnover
Regarding seasonality (Fig. 5A), there were 72 species found in autumn and 37 species found in winter with 29 retrieved in both seasons. Nearly half of Fig. 5 Macroinvertebrate species' site occupancy plots (left) and ParCat distribution plots (right) for the 'BMI species' dataset. A Species found in the two sampling seasons, B species sorted after stream section divided by the barrierautumn, C species sorted after stream section divided by the barrier-winter. The species' names marked in red indicate taxa potentially new for Sicily all species found were dipterans. Over 40% of all species found in each season were only retrieved from a single site with mayfly Cleon dipterum being the most widespread species in autumn and Echinogammarus sicilianus being retrieved from most sites in winter. Notably both amphipod species and the only snail species were found in both seasons along with nearly all caddisfly species (Table S7). On the contrary, most dragonflies and water beetles were observed in autumn, whereas two mayflies and a freshwater shrimp were retrieved only in winter.

Barrier effect
In autumn, 55 of the total 77 species were found in the downstream and 45 in the upstream section relative to the barrier, with 23 being shared between the sections (Fig. 5B). In the downstream sites eDNA detections of two species were retrieved from all sites: Cleon dipterum and Echinogammarus adipatus, whereas in the upstream Echinogammarus sicilianus and the chironomid Parametriocnemus stylatus were recorded from each site. Over 50% of the species detections were from single sites only. In winter, 28 of the total 33 species detected via eDNA in the Fosso del Tempio river system were found downstream of the barrier and 13 upstream with 8 shared between sections (Fig. 5C). Similarly to autumn, Echinogammarus sicilianus was retrieved from all sites in the upstream, whereas no species was recorded from all downstream sites with Simulium intermedium recorded from all sites but one. Again, in both stream sections, the majority of species was recorded from single sites only. Interestingly, regardless of the season, most caddisfly species, along with the only hemipteran, were found in the upstream sections only, whereas most dragonflies and water beetles were retrieved downstream of the barrier (Table S7). Notably both amphipod species were found on both sides of the barrier in both seasons.

Discussion
Observed diversity-freshwater macroinvertebrates and beyond By using eDNA metabarcoding, we detected the presence of 98 freshwater macroinvertebrate species in the Fosso del Tempio river system and its neighbouring catchment. The detected diversity counts for approximately 7.5% of all freshwater macroinvertebrate species known from the island (Ruffo & Stoch, 2006;Latella et al., 2007). To the best of our knowledge, no freshwater macroinvertebrate biomonitoring data for Fosso del Tempio exist and thus, it is hard to estimate what fraction of the local diversity we were able to retrieve. However, given that a relatively small portion of Sicily's hydrological network was covered with this study, the number of macroinvertebrate species retrieved is quite high. Moreover, taking into consideration that nearly 65% of obtained diversity could not be assigned to species level, the detected diversity is even higher. The majority of all freshwater macroinvertebrate species detected are arthropods, which were also the most dominant group in terms of read abundance. Among arthropods, the most dominant order in terms of both read abundance and diversity was Diptera. This result was to be expected since nearly half of all freshwater macroinvertebrates reported from Sicily belong to this order (Ruffo and Stoch 2006;Latella et al., 2007). It is also in agreement with other eDNA-based bioassessment of freshwater macroinvertebrates, where dipteran representatives are among the most dominant groups (Fernández et al., 2019;Leese et al., 2021). Within Diptera, we have detected 21 species belonging to the ecologically important family Chironomidae, which along with 12 species of 'so-called' EPT (Ephemeroptera, Plecoptera, Trichoptera) taxa, could serve as a basis for estimating the ecological state of the freshwater ecosystem in Fosso del Tempio river system (Wallace et al., 1996). Notably, we have also detected the presence of a single freshwater snail, the highly invasive snail Physella (Physa) acuta, known from Sicily already since the mid-nineteenth century (Vinarski, 2017). Since it was also the only snail species found, we have considered the possibility of primer bias towards inadequate amplification of molluscan DNA. However, after evaluating the primer binding for ten out of approximately 23 freshwater snail species known from Sicily (Ruffo & Stoch, 2006), for which DNA information was available, we excluded that possibility since all of them could be hypothetically amplified with used primers (unpublished data). Absence of other snail taxa could be also resulting from low amounts of DNA shed in the environment either due to hard shell structure also possibly combined with very low numbers of individuals, which largely influence the detectability of an eDNA signal (Harrison et al., 2019;Stewart 2019). However, given we were still able to detect considerable diversity of taxa that bear hard carapax like crustaceans or water beetles, it could be rather due to the low number of individuals present. On the other hand, the observed low gastropod diversity could possibly indicate that invasive P. acuta have outcompeted other native snail species as already documented in other parts of the world (Dobson 2004;Zukowski and Walker 2009). This hypothesis is difficult to verify though, since no data on historical freshwater snail diversity in Fosso del Tempio is available. However, since P. acuta was only detected in relatively few samples, it might indicate either low abundance of snail species or low amount of DNA shed to the water.
In some cases, we also noticed a high level of intraspecific diversity, particularly in three species where more than 10 ESVs were ascribed to a single taxon namely a chironomid Cricotopus bicinctus, a simulid Simulium rubzovianum and an amphipod Echinogammarus sicilianus. Those findings confirm prior observation in all three taxa indicating already high degree of intraspecific diversity and in case of E. sicilianus the presence of potential cryptic diversity (Sinclair and Gresens 2008;Sari et al., 2012;Hupało et al., 2021). Interestingly, we have detected 28 species that to our knowledge, were not reported from Sicily before, indicating they might represent taxa potentially new for Sicily. Most of them belong to dipteran families Chironomidae, Simuliidae, Culicidae, Limoniidae, Dixidae, Empididae, Tipulidae and Muscidae with all species known to have an aquatic larval stage. Most of those species are known to occur in the Mediterranean basin indicating they might represent formerly overlooked valid records. However, five of them: Culex inconspicuosus, Gonomyia tenella, Limnophora olympiae, Simulium ruficorne and Stempellinella ciliaris were not reported from vicinity of Mediterranean basin with known occurrence data coming from Central and Northern Europe as well as Central and Southern Africa, indicating that those might represent false positive signals. Those false positives might result from mistakes present in public databases including incorrect synonyms and wrong sequence entries among others (Pentinsaari et al., 2020), which could contribute to observed records despite BOLDigger strategy to use the most common species across the most frequent of the most similar identification hits to enhance reliability of the record (Buchner and Leese 2020). Regarding other taxa, which might be new records for Sicily, there are three aquatic oligochaetes (Dero furcata, Nais christinae, Nais elinguis), all reported from the Mediterranean basin, one elmid beetle Limnius volckmari also occurring in the Mediterranean with a population on Corsica and one caddisfly Hydropsyche instabilis reported also from southern Italy. Although all of those occurrences might represent valid new records, one should treat them with caution, since there were no actual specimens observed or collected from the sites, indicating that these require further investigation.
Our target organism group-freshwater macroinvertebrates-comprised 466 ESVs, which reflected about 6.5% of all ESV diversity obtained in this study. Although the vast majority of the remaining 6865 ESVs could not be assigned to species level, they still provide valuable information about the ecosystem surrounding the Fosso del Tempio system. We have retrieved traces of multiple entities of diatoms, other unicellular algae as well as heterotrophic and mixotrophic protists and fungi. Although only some could be assigned to species level, most of them could be resolved to higher taxonomic levels only, like phyla or classes. Since COI is suggested as barcoding region primarily for metazoans, in many cases the taxonomic assignment could probably be resolved better if other, group-specific molecular markers would be used e.g. rbcl, 18S rDNA or ITS (Evans et al., 2007;Pawlowski et al., 2012;Schoch et al., 2012). We have also retrieved DNA traces of several terrestrial arthropod species including 27 dipterans, 26 butterflies, 21 true bugs, 20 beetles, 12 hymenopterans, seven spiders, five orthopterans, five collembolans and a single mantid, phasmid, lacewing, bark lice and isopod. There was also a significant amount of information deriving from freshwater organisms, which we did not consider as our target group. We have also observed DNA signals coming from planktonic biota including five rotifers and two Hydra species. Moreover, the DNA trace of common carp (Cyprinus carpio) was retrieved from a single site. Interestingly, C. carpio established stable wild populations in inland waters of Sicily and is considered parautochthonous since being introduced supposedly by Romans (Marrone & Naselli-Flores, 2015). It is likely representing a valid record since results obtained using fish-specific 12S marker (Taberlet et al., 2018) from the same sites confirmed its presence (Hupało et al., unpublished results). Altogether, the information included in the non-target part of the dataset might provide valuable information about the surrounding ecosystem. More and more studies indicate the significance and validity of information contained in so-called bycatch derived from aquatic eDNA data (Macher et al., 2021b;Mariani et al., 2021). Thus, including and sharing the entire eDNA sequence dataset could be of paramount importance for further studies from the region and effort should be made towards their storage and public availability (Berry et al., 2020;Makiola et al., 2020). On the other hand, we have also retrieved a number of ESVs that represent obvious false positive records, e.g. including ones that were taxonomically assigned to Echinodermata, a phylum that is known to occur only in marine waters. However, given the low level of genetic similarity against the reference data i.e. even going below 80%, one can likely treat this as misidentification rather than a contamination, especially considering that no DNA signals were retrieved from negative controls. Regardless, one still should treat the records mentioned above with caution, because even though many of them could represent valid records, some of them might be a result of misidentification due to faultiness and incompleteness of reference databases.

Insight into seasonality and barrier effects
Our results obtained with eDNA-based information provide an insight into the dynamics of the Fosso del Tempio river system. By conducting sampling campaigns in both autumn and winter seasons, we were able to investigate seasonal patterns of freshwater macroinvertebrate diversity. We observed more than twice as many macroinvertebrate species in autumn compared to winter. This result is to some degree unexpected. Even though differences in species composition in the Mediterranean region could be partially attributed to a high degree of water level fluctuation observed between arid summer and wet winter, the species richness is expected to remain similar between seasons (Gasith & Resh, 1999;Bonada & Resh, 2013). Given that samples were collected in September, at the end of dry season in Sicily and in January, which marks approximately the middle of wet season, the habitat and resource availability for aquatic biota are very different. This should reflect the changes in species' assemblages and dominances, however the majority of species should be evolutionarily adapted to those varying conditions, resulting in certain similarities in species diversity between seasons (Gasith & Resh, 1999). One of the reasons behind observed differences between the seasons could be visibly higher water levels during winter sampling. This could also have an indirect effect on lower diversity observed as the DNA concentration gets more diluted and thus, less taxa are being retrieved. That being said, with a relatively high number of species shared between seasons, one may assume that at least some of the taxa reported only from autumn could in fact also occur on sites in winter and could have been retrieved if a higher water volume would have been sampled. On the other hand, there were eight species that were only retrieved from winter samples. Although in some cases like Eukieferella claripennis, Limnius volckmari and Tipula lateralis, it seems to be the case of rarity of those taxa in autumn since they were only retrieved locally, other cases might show some true seasonal patterns. In case of mayflies Baetis pavidus and Caenis pusilla a possible niche overlap with Caenis luctuosa mediated with supposed seasonal fluctuations might be the reason behind observed species presence in close-by sites in autumn and winter. The same could be true for the water beetle Agabus brunneus 'replacing' Agabus bipustulatus on site A12 in winter compared to autumn. Interestingly, we have also retrieved an eDNA signal from the freshwater shrimp Atyaephyra desmarestii, native to the Mediterranean region. The retrieved information from winter eDNA sampling corresponds with the observations made in the field where numerous specimens were observed only during winter sampling with no specimen found in autumn. Although the species is known to have strong seasonal fluctuations (Dhaouadi-Hassen & Boumaiza, 2009;Fidalgo et al., 2015;Schoolmann et al., 2015), the species' abundance should be higher in autumn compared to winter. Since the opposite is the case here, this puzzling observation likely requires further investigation.
Our data also indicate that the concrete weir present in the stream course might have an effect on dispersal of resident aquatic taxa. There are numerous studies confirming that weirs may act as barriers for macroinvertebrate dispersal significantly reducing both connectivity and observed species diversity on the opposite sites of the barrier (Poff & Zimmerman, 2010;Brooks et al., 2018). The effect seems to be even more profound when there is a high degree of water level fluctuations, which is often the case in Mediterranean streams. Here, we observed a significant change of water level and flow velocity between dry and wet seasons. During sampling conducted in September, we have seen virtually no water going over the weir, whereas in winter with the higher water level, the change in amount of water was clearly visible, although flow was still disrupted by the barrier. This effect seems to be reflected by the results of similarity-based analyses with a clear distinction between the observed diversity in the downstream and upstream samples during dry season compared to wet season, where no visible difference was observed. This is also reflected in the proportion of species that were found only upstream catchment between seasons with a much lower percentage of species uniquely found in the upstream in winter compared to autumn. A similar pattern was found when comparing diversity observed on two border sites separated by the barrier, where a comparatively higher proportion of diversity was shared between those sites in winter compared to autumn. Although, any species-related comparisons between the seasons have to be taken with caution due to various reasons described above, the observed effect of an in-stream weir could have significant implications for observed diversity in the Fosso del Tempio system. Based on our results, the weir seems to hinder the connectivity between upstream and downstream parts of the system with a portion of taxa being found on both sides of the barrier. Regardless of season, the majority of taxa shared between stream sections are the hemilimnic ones with flying adult stages. Increased dispersal ability in flying insects with aquatic larvae has been proven to be important in terms of observed patterns of diversity in a stream ecosystem with a high degree of similarity between river sections regardless of any possible barriers (Hughes et al., 2009). Similarly, the opposite is true for the species with lower dispersal abilities with more fractionated and restricted patterns of genetic diversity. In such cases, the in-stream barriers might lead to genetic differentiation leading to strong population structuring and in some cases, to allopatric speciation. We have observed a similar case in our data, where certain ESVs of Echinogammarus sicilianus, a freshwater amphipod with limited dispersal abilities, are present only in the downstream or upstream section of the Fosso del Tempio, regardless of the season. This finding indicates that the intraspecific diversity observed might, at least partially, be resulting from the effect of the in-stream weir.

Different approaches, similar results-ESV vs. species comparability
The results regarding seasonality patterns and the effect of the in-stream barrier obtained with macroinvertebrate species-based approach were similar to the ones obtained with all macroinvertebrate ESVs including more than 60% of all macroinvertebrate ESVs which could not be assigned to species level. Although notable progress in DNA barcoding of aquatic biota has been made in recent years, there is still a significant amount of data missing in reference databases (Weigand et al., 2019). The difference in freshwater macroinvertebrate barcode coverage varies between different regions due to unbalanced efforts taken in filling the missing gaps. Our results indicate that Sicilian freshwater macroinvertebrates seem to be relatively underrepresented with only 35% of our dataset being assigned to species level, compared to the European average of approximately 64.5% freshwater macroinvertebrate species having at least a single barcode publicly available (Weigand et al., 2019). In those cases, ESV-based approaches could provide an alternative solution for studying freshwater ecosystem dynamics. This seems to be also of particular importance for research based on organisms where taxonomy is complex and not fully resolved like diatoms where the use of ESVs has been successfully implemented (Tapolczai et al., 2019). Moreover, ESVs provide a finer level of detail to observed diversity by providing an insight also into intraspecific diversity, which is discarded when species or clusters serving as species proxy (e.g. OTUs) are considered (Zizka et al., 2020;Tapolczai et al., 2021). This could be of paramount importance when analysing phylogeography or population genetics of observed taxa (Turon et al., 2020;Antich et al., 2021).
Since the primers used in this study amplify a broad range of taxa groups and given that macroinvertebrate diversity represented only a fraction of all retrieved ESVs, we have also decided to look at the results obtained with an approach where all generated diversity will be taken into consideration. Surprisingly, despite taxonomic broadness of the dataset, including also a high portion of non-freshwater biota, the results obtained were similar to those obtained only with macroinvertebrate data. This finding seems particularly interesting pointing out the high sensitivity and reproducibility of the generated results. Even though the results deriving from this general approach should be looked at with increased awareness, they still can inform about community connectivity or seasonality patterns in a stream ecosystem.
Implementation of (e)DNA-based results into routine bioassessments Regular biomonitoring of freshwater ecosystems towards improving water quality is now a standard procedure taking place in numerous countries worldwide. In Europe, the achievement of 'good ecological status' of the surface waters is regulated by the European Union's Water Framework Directive (WFD, Directive 2000/60/EC), which legally imposes the continuous bioassessment of several key organism groups being 'Biological Quality Elements' (BQEs), including, among others, benthic macroinvertebrates. Even though the diversity estimates within WFD are based on classical morpho-taxonomic approaches, eDNA metabarcoding has been shown to have great potential for being a powerful complementary tool for routine biomonitoring (Brantschen et al. 2021, Pawlowski et al., 2018. It often shows comparable or even higher levels of precision compared to traditional morphological assessments Leese et al., 2021). In principle, the eDNA-based results obtained in this study could provide valuable information about the studied river system. However, there are certain limitations to its direct applicability. Firstly, to our knowledge, there are no past or present biomonitoring data available for regional macroinvertebrate diversity in the Fosso del Tempio system and it is therefore difficult to place the obtained results in a larger context or make any comparisons with prior assessments. Another issue is the incompleteness of the DNA reference databases, hindering species level assignment. Finally, for implementing eDNA-based methods in routine biomonitoring, further validation and standardisation including modification of the assessment indices is needed (Hering et al., 2018). Given the constant improvement of reference libraries combined with the growing interest about Mediterranean freshwater biodiversity, there is a potential for using the generated eDNA results in future bioassessments.

Conclusions
Our results provide a first insight into freshwater diversity and community dynamics of the Fosso del Tempio river system on Sicily, used as an exemplary Mediterranean insular freshwater ecosystem. We showcase the potential that comes with using eDNA metabarcoding for studying freshwater ecosystems, even in understudied regions with still significant portions of genetic information missing like the Mediterranean islands. Based on our findings, an ESV-based approach provides a promising, highly reproducible approach for studying freshwater community patterns and dynamics even, or in particular, in highly understudied regions. By providing reproducible taxonomic units, data could be easily reused in the future with new regional biodiversity data added and with increasing completeness of local reference databases. Hence, we highlight the need for a unified digital storage and access solution, which could ensure the public availability of the data at both regional and international level. Specifically for regions with non-perennial flows, hydrological differences between seasons need to be considered and incorporated in bioassessment planning, sampling and interpretation. We stress the importance of proper sampling design taking into consideration seasonality and equal representation of sites across the entire river system to maximise the representation of the regional community and its dynamics.

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