In Alpine rivers, rapid melting of glaciers will lead to changes in aquatic ecosystems, yet their current state remains poorly known. Benthic macroinvertebrates have been a primary focus of monitoring such ecosystems because they are highly sensitive to environmental disturbances. Traditionally, monitoring involves physical sampling such as using the kick-net method. Yet, such monitoring can be challenging or unfeasible due to harsh climatic conditions. Environmental DNA (eDNA) has been proposed as an alternative approach, as samples can be obtained without direct access to the river. However, the year-round feasibility of this technique, especially during winter months, and differences relative to kick-net sampling have not yet been fully assessed. Here, we compare kick-net and eDNA sampling as year-round tools to monitor alpine benthic communities. We sampled these communities every two months with both methods and examined the results in relation to the physicochemical parameters of the river (sediment load, temperature and conductivity). We found a significant correlation between the proportion of kick-net individuals and their detectability using eDNA, but this correlation was not as clear when using the relative read abundance of sequences instead of the binary detection of taxa. As for kick-net data, we did not find a direct link between sediment load and a decrease on the number of individuals recovered, but we found a strong negative correlation with temperature, indicating an increase of individuals with colder water temperatures. Overall, our results suggest that eDNA may serve as a year-round alternative to kick-net sampling, while also highlighting the need for further research to elucidate technique-specific biases among macroinvertebrate groups.

Study site

The Val Roseg is located in the south-eastern part of the Swiss Alps, in the canton of Graubünden, on the border with Italy (Figure 1). The Roseg river originates on the Bernina massif at 2160 m a.s.l, from the natural proglacial lake called Lej da Vadret (or Lake Roseg), which collects the meltwater discharge of the Roseg glacier and other minor glaciers (Figure 1). The Lej da Vadret was formed naturally in the mid-20th century due to the damming action of the central moraine, impacting the distribution patterns of organic matter and the benthic community in the river, due to its retention capacity of sediment loadings and the consequent decrease in downstream turbidity (Brittain & Milner, 2001; Burgherr et al., 2002; Tockner et al., 2002). The Roseg River is a second-order tributary of the Inn River, which is in turn tributary of the Danube River, and covers a distance of 11.3 km with a difference in altitude of 405 m, at the end of which joins the Bernina River (Ova da Bernina), before finally flowing into the Inn River (Burgherr & Ward, 2000; Tockner et al., 2002). The Roseg River drains a catchment area of 66.5 km², with an altitude ranging from 1755 m a.s.l (confluence with the Bernina River) to 4049 m a.s.l (Piz Bernina). Precipitation data for 1951 to 1980 indicate a mean value of 1600 mm/year, half of which fell as snow, while the annual mean flow is 2.76 m³/s (over the period 1955-1997) (Burgherr & Ward, 2001; Tockner et al., 2002).

The Roseg river is mainly fed by two glaciers, and therefore presents a typical glacial hydrological regime, with a peak flow in summer and a strong daily variability of discharge. In winter, river flow depends mainly on groundwater-fed contributions and on a constant minimum flow provided by the lake. The floodplain of the Val Roseg has been described in detail (Malard et al., 1999; Tockner et al., 1997). It lies between 1995 and 2060 m a.s.l, is between 100 and 500 metres wide and covers an area of 0.67 km² (Malard et al., 1999; Zah et al., 2001). Along the floodplain, the river takes on a very complex braided structure, with a wide variety of channels with very different physical and chemical characteristics, forming a heterogeneous network of aquatic habitats (Figure 1) (Tockner et al., 1997). In summer, despite the intense daily discharge variation due to glacier melt, the Roseg floodplain is exceptionally stable, and some channels have remained unchanged for more than 50 years. In winter, the volume of water flowing through the valley decreases considerably, leading to the drying out of many channels in the floodplain (Malard et al., 1999).

Sampling locations and experimental design

Sampling sites were selected to cover the floodplain, the tributaries and the two main sources of the Roseg river, i.e. Tschierva glacier and Roseg lake (Figure 2). Samples were collected every two months between October 2021 and August 2022, both included. The experimental design consisted of sampling the junctions between the tributaries and the main river course, in order to study the relative contribution of the tributaries on the main river, as part of a larger monitoring project (Barchi & Lane, 2022). We selected the three main tributaries, i.e. T1, T2 and T3, but other tributaries feed the floodplain, which we did not sample for logistic reasons (Figure 2). We added three sampling points away from the junctions, i.e. glacier site, in the middle of the floodplain, i.e. main site, and at the end of the floodplain, i.e. main end (Figure 2). Samples were labelled as part of the main course or as tributaries. The two sampling methodologies were used at each point, i.e., kick-net sampling and eDNA filtering.

In a few cases, we were unable to collect samples (see Figure 2 for details). For eDNA samples, the three samples from T3 junction in October were not sampled for eDNA. Glacier and Main samples in December were not sampled due to the impossibility of accessing them. Sample T1 in February was not collected due to the scarcity of the flow, which impeded filtering any water. For kick-net samples, T1, Glacier and Main in December and T1 and T1 up in February were not sampled for kick-net individuals, either due to scarcity of flow or impossibility of access.

Kick-net sampling

For each sampling site, macroinvertebrates were collected using the kick-sampling method, following the standard procedure described by (Stark, 2001), with duplicates. A 25 cm x 25 cm kick-net with mesh gaps of 500 um was held downstream of the point to be sampled. The substrate material upstream of the net (0.35 m², i.e. 30-40 cm x 100 cm surface area) was disturbed with the aid of a foot for the duration of approximately 30 seconds, allowing the raised material to flow independently into the net, carried by the current (Hauer, 2007). After roughly removing the larger grains and any plant matter unintentionally collected, the contents of the net were transferred to bottles filled with 70% EtOH (95% in winter) and refrigerated. Macroinvertebrates were later sorted using a stereoscopic microscope, counted and identified to the family level by means of the identification key drawn up by Tachet et al. (2000), and with the help of the online portal Perla (DREAL, 2021).

In all cases, the downstream site was sampled first and the upstream site second, to avoid influencing the results or leading to eDNA cross contamination, i.e. T down first. Important differences in habitat and microhabitat conditions can occur at the level of the individual sampling site, yet the two sets of data were collected by selecting habitats that were as similar as possible, specifically favouring riffle habitats. This was done to maximise the comparability of results between different sites, minimising data variability (Stark, 2001). Given strong seasonal fluctuation in water discharge, and thus the variation in wetted width and water velocity, priority was given to keeping microhabitat conditions as constant as possible throughout the year and as consistent as possible between different sites. Kick-net communities were then transformed in two ways: to proportion of individuals and to scaled number of individuals of each family.

eDNA sampling

eDNA water samples were taken in triplicates. For each replicate, 30 L of water was pumped through an eDNA capsule (0.45 µm, Waterra, Canada) to retain DNA from aquatic organisms in the filter within. A total of 90 L of water was pumped for each sampling location. To preserve the DNA until DNA extraction, we filled the capsules with 40 mL of Nucleic Acid Preservation buffer (NAP buffer described in Camacho‐Sanchez et al., 2013). Samples were preserved at room temperature until extraction.

Environmental variables

We measured temperature (°C), electrical conductivity (μS/cm) and suspended sediment load (mg/L) at each sampling point. These variables are often used to assess the glacial influence on the development pattern of the macroinvertebrate community, and are often associated with the environmental harshness gradient that describes glacial river habitats  (Castella et al., 2001; Dangles, O., 2012; Ilg & Castella, 2006; Lods-Crozet et al., 2001; Niedrist & Füreder, 2018).

Temperature and electrical conductivity were measured with a multi-parameter portable metre (WTW MultiLine Multi 3630 IDS), calibrating the respective probes daily before taking measurements. For the sediment load, a 500 mL sample of water was collected at each sampling site. In the laboratory, the water was filtered through PES 0.22m membrane filters with the help of a pump filtration system. The filters were then dried completely in the oven at 105°C (for approximately ~2h30) and the dry sediment mass derived with an analytical balance.

eDNA extraction

Extractions were performed in a pre-PCR laboratory dedicated to low DNA-content analyses. In brief, capsules were shaken thoroughly for 2 min to maximise DNA yield from the filter. The content of the capsule was then transferred to 50 ml Falcon tubes. All the material for this procedure was sterilised between samples using bleach and Mili-Q water, in this order. Tubes were then centrifuged at 7500 rpm for 1h at 4°C. Supernatants were discarded and 1440 μl of ATL buffer from the DNeasy Blood & Tissue Extraction Kit (Qiagen) and 40 μl of proteinase K were added. After vortexing and resuspending the pellet, the mixture was transferred to 2 ml Eppendorf tubes and incubated for 2h at 56°C. The DNA extraction was pursued following the NucleoSpin Soil Kit protocol (Macherey-Nagel). Elution was done with 2 x 50 μl of SE buffer. Negative controls were included at all steps. The extractions were tested for inhibitors with qPCR to identify the best dilution fold and PCR cycles. qPCR reagents and conditions were the same as in PCR amplification (see below), with the addition of SybrGreen (Thermo Fisher Scientific). Based on these results, and follow-up PCR tests (data not shown), all samples were subsequently diluted 10-fold before 55 cycles of PCR amplification.

DNA metabarcoding

Each biological replicate was amplified in triplicates (3 biological replicates, i.e., 9 PCR replicates per sampling location and time). Each PCR reaction was done using barcoded (tagged) primer pairs, i.e. with unique 8bp sequences attached to the 5’ end of both forward and reverse primers, to label and trace back sequences to each unique PCR replicate. DNA extracts were amplified using a generalist arthropod primer pair (Arth02, (Taberlet et al., 2018)) targeting the 16S mitochondrial rDNA (76-168 bp) of arthropod species. Thermocycling conditions were as follows: denaturation at 95°C for 10 min, followed by 55 cycles of 30 s at 95°C, 30 s at 49°C and 1 min at 72°C, with a final elongation step of 7 min at 72°C. Each 96-well PCR plate contained 20 biological replicates, each amplified in triplicate (n=60), while the remaining wells were used for blanks, negative extraction controls, negative PCR controls, and positive controls (DNA assemblies of a mix of eight arthropod species, see Suppl. Table 1). Blanks correspond to empty wells and allow to estimate the proportion of tag switches (i.e., false combination of tags, generating chimeric sequences) occurring during the sequencing process (Schnell et al., 2015). Amplification success was confirmed on a 1.5% agarose gel. Then, PCR products from each PCR plate were subsequently pooled and purified separately using a MinElute PCR Purification Kit (Qiagen). Purified PCR products were quantified using a Qubit Fluorometer (Life Technology Corporation). Libraries were prepared following the Tagsteady protocol as in Carøe & Bohmann (2020), which reduces tag jumps. Nine libraries were constructed, each containing the amplicons from 3 PCR plates pooled together. Each library contained a unique index combination to ensure the combination of barcoded primers and library indexes was unique for each PCR replicate. Final libraries were quantified (KAPA Library Quantification Kit, Roche), normalised and pooled before 150 paired-end sequencing on an Illumina MiniSeq sequencing system with a High Output Kit (Illumina, San Diego, CA, USA).

Bioinformatic data analyses

The bioinformatic processing of the raw sequence output and first filtering was done using the OBITools software (Boyer et al., 2016). Forward and reverse reads were assembled with a minimum quality score of 40. The aligned sequences were assigned to samples based on unique tags combinations. Assigned sequences were then de-replicated, retaining only unique sequences. All sequences with less than 100 reads per library were discarded as well as those not fitting the range of metabarcode lengths. This was followed by two different clustering methods. First, pairwise dissimilarities between reads were computed and lesser abundant sequences with single nucleotide dissimilarity were clustered into the most abundant ones. Second, we used the Sumaclust algorithm (Mercier C, 2013) to further refine the resulting clusters based on a sequence similarity of 97%. It uses the same clustering algorithm as UCLUST (Prasad, D.V., 2015) and it is mainly used to identify erroneous sequences produced during amplification and sequencing, derived from its main (centroid) sequence. Sequences were taxonomically assigned to taxa with a database for Arth02 generated using the EMBL database (European Molecular Biology Laboratory, Release 143).

Further data cleaning and filtering was done in R (version 4.0.2) using the metabaR package (Zinger et al., 2021). Sequences that were more abundant in extraction and PCR controls than in samples were considered as contamination and removed. Operational taxonomic units (OTUs) with similarity to the reference sequence lower than 97% were removed from the dataset. Removal of tag-leaked sequences was done independently for each library, using the function tagjumpslayer from the same package, to reduce the number of false positives. It is based on the principle that abundant OTUs are more likely to be detected in low abundance in samples where they are not supposed to be compared to rare OTUs. Different thresholds are tested to evaluate how the filtration procedure affects the dataset, and each OTUs abundance is corrected accordingly when their relative abundance over the entire dataset is below this threshold.

This approach allowed us to discard single OTUs instead of whole PCR replicates. However, PCR replicates with too small reads count were also discarded. Taxonomic assignment was refined by manually inspecting retained OTU with more than 100 reads count using the BLAST algorithm of GenBank. For cases with multiple candidate species, the geographic range was considered to select the most probable species or the OTU was assigned at a higher taxonomic level. If more than one species was equally matching an OTU, we assigned it to a higher taxonomic level. All unique sequences and their taxonomic assignation are available (Suppl. Data). Finally, singletons were removed. Remaining PCR replicates were merged by sample location and sequence, keeping the mean relative read abundance (RRA) and the Presence/Absence (P/A) for each OTU. We also grouped the OTUs by species, family and order.

Statistical analyses

All statistical analyses were carried out in R (Version 4.0.2). First, to get an overview of how environmental variables varied through time and space, we calculated a principal component analysis (PCA) and visualised their differences for all sampled locations and sampling months (hereafter referred to as "season"). Second, in order to quantify the similarity between the two methodologies, we modelled the likelihood of eDNA detection (Presence/Absence, i.e. P/A) based on the proportion of kick-net sampled individuals (Binomial distribution, Model 1) and the sampling season. This was also done using the RRA dataset (Beta distribution, Model 2).

P/A      ~          Prop. Kick-net individuals      *         Season                         Model 1

RRA     ~          Prop. Kick-net individuals     *          Season                        Model 2

We hypothesised that more abundant taxa detected through kick-net sampling would be more likely to be detected through eDNA. Model 1 and 2 included only eDNA data from taxa which had been found using kick-net sampling. Comparisons would have been misleading if taxa only recorded through eDNA were also included, due to the greater diversity of taxa detected with eDNA and because it also detects species which are not present in the river.

Third, we focused only on the key water quality indicator arthropod groups (i.e. Ephemeroptera, Tricoptera, Plecoptera and Diptera) to further study the correlation between these two methodologies, and assess how eDNA quantifies them compared to Kick-net sampling, and what is the potential variability across seasons. We modelled using a linear mixed effect model using the eDNA data at order level against the proportion of kick-net individuals (Model 3).

RRA     ~ Prop. Kick-net individuals   * Arthropod order * Season + (1 | Sample)       Model 3

Fourth, we investigated to what extent the measured river water properties influenced the abundance of kick-net communities (Model 4). As predictor we used the scaled number of individuals, which followed a Normal distribution. We tested each environmental variable in a separate model, i.e. Conductivity, Sediment load and Temperature. All modelling was done using the glmmTMB package.

Scaled Sum of Individuals       ~          Environmental variable                      Model 4