Environmental DNA simultaneously informs hydrological and biodiversity characterization of an Alpine catchment

Alpine streams are particularly valuable for downstream water resources and of high ecological relevance; however, a detailed understanding of water storage and release in such heterogeneous environments is often still lacking. Observations of naturally occurring tracers, such as stable isotopes of water or electrical conductivity, are frequently used to track and explain hydrologic patterns and processes. Importantly, some of these hydrologic processes also create microhabitat variations in Alpine aquatic systems, each inhabited by characteristic organismal communities. The inclusion of such ecological diversity in a hydrologic assessment of an Alpine system may improve our understanding of hydrologic flows while also delivering biological information. Recently, the application of environmental DNA (eDNA) to assess biological diversity in water and connected habitats has gained popularity in the field of aquatic ecology. A few of these studies have started to link aquatic diversity with hydrologic processes but hitherto never in an Alpine system. Here, we collected water from an Alpine catchment in Switzerland and compared the genetic information of eukaryotic organisms conveyed by eDNA with the hydrologic information conveyed by naturally occurring hydrologic tracers. Between March and September 2017, we sampled water at multiple time points at 10 sites distributed over the 13.4 km² Vallon de Nant catchment (Switzerland). The sites corresponded to three different water types and habitats, namely low-flow or ephemeral tributaries, groundwater-fed springs, and the main channel receiving water from both previous mentioned water types.

Accompanying observations of typical physicochemical hydrologic characteristics with eDNA revealed that in the main channel and in the tributaries, the biological richness increases according to the change in streamflow, $\mathrm{d}q/\mathrm{d}t$, whereas, in contrast, the richness in springs increased in correlation with electrical conductivity. At the catchment scale, our results suggest that transport of additional, and probably terrestrial, DNA into water storage or flow compartments occurs with increasing streamflow. Such processes include overbank flow, stream network expansion, and hyporheic exchange. In general, our results highlight the importance of considering the at-site sampling habitat in combination with upstream connected habitats to understand how streams integrate eDNA over a catchment and to interpret spatially distributed eDNA samples, both for hydrologic and biodiversity assessments. At the intersection of two disciplines, our study provides complementary knowledge gains and identifies the next steps to be addressed for using eDNA to achieve complementary insights into Alpine water sources. Finally, we provide recommendations for future observation of eDNA in Alpine stream ecosystems.

1 Introduction

Alpine environments are often considered to be water towers for lowland areas (Viviroli et al., 2007) and hotspots for biodiversity (Körner, 2002) and will likely be disproportionately affected by climate change compared to other areas (Jacobsen et al., 2012; Grabherr, 2009). The complex topography of mountain areas results in a tremendous variation in the physical environment in terms of elevation, slope, and exposure, which subsequently creates large gradients of received radiation and thus extreme variations in temperature at small scales. Subsequently, the topography and temperature gradients translate into a landscape with an enhanced heterogeneity of potential water flow paths. These paths are additionally governed by highly seasonal incoming water (rain or snow), which is released as a seasonal alternation of runoff from rain, stored snow, and snowmelt. Consequently, channel networks in Alpine areas are known to regularly expand and contract, leaving a portion of streams ephemeral (Godsey and Kirchner, 2014; van Meerveld et al., 2019).

Despite the importance of mountain water resources, hydrologic processes are still often poorly characterized in mountain catchments. From a hydrologic as well as from an ecological viewpoint, the central phenomenon determining habitat availability is the origin of stream water, for example, whether streams are fed by seasonal snowmelt or by recent precipitation (Nolin, 2012; Cochand et al., 2019; McGuire and McDonnell, 2015). In headwaters of mountain environments, the relevant hydrologic mechanisms are difficult to observe because they occur out of direct reach, for example in the atmosphere or the subsurface (McDonnell et al., 2007). There is an ongoing search for hydrologic tracers to decompose river water into its most recent storage compartment (e.g., glacier, snow, soil, groundwater) and to identify when and where it entered the stream network (Abbott et al., 2016; Blume and Van Meerveld, 2015; Williams et al., 2009; Mosquera et al., 2018).

Spatial and temporal variations in the physicochemical properties of water in Alpine catchments, such as in temperature, discharge, stable isotopes, and electrical conductivity (EC), offer clues to identify water flow from glacier, snow, soil, and groundwater and its associated timescales (Zuecco et al., 2019; Tetzlaff et al., 2015; Penna et al., 2014; Klaus and McDonnell, 2013; Beria et al., 2018). For example, EC can relatively consistently discriminate snowmelt, rain, and potentially glacier melt, which all contain few solutes due to their little contact with rock and soil surfaces, from groundwater, which typically has much higher levels of solutes due to extended contact with surfaces (Williams et al., 2006; Cochand et al., 2019; Kobierska et al., 2015). On the other hand, water temperature cannot be used as a conservative tracer because it quickly changes along a water course, but it can provide insights on separating surface water (having a strong diel variation) and influxes of water from the subsurface (having more constant temperatures; Constantz, 1998; Comola et al., 2015; Hoehn and Cirpka, 2006; Westhoff et al., 2007; Selker et al., 2006). Finally, stable isotopes of water, δ¹⁸O and δ²H, have classically been used as tracers of source water, flow paths, and precipitation contributions. Precipitation becomes depleted, or has lower values, when it falls at higher elevations and in lower temperatures, especially when it falls as snow. This natural variation in isotopic ratio in precipitation according to season and elevation offers clues as to when and where it entered the system. However, in high Alpine contexts, it often falls short in offering additional insights into dominant flow paths because the isotopic ratio of all potential water sources is dominated by the precipitation phase of incoming water (rain versus snow; Lyon et al., 2018; Beria et al., 2018; Penna et al., 2018). Landwehr and Coplen (2006) improved upon the direct observation of isotopes with the line-conditioned excess (lc-ex), a metric that is equal to the residual between two isotopes of water, δ¹⁸O and δ²H. Higher values of lc-ex thus indicate when water has experienced losses to evaporation or other fractionation, allowing it to be used as an additional tracer alongside the ratios of oxygen and hydrogen isotopes individually.

Similar to their high variations in water proprieties, Alpine catchments are typically dominated by high fluctuations in discharge and corresponding changes in the extent of the stream network (Godsey and Kirchner, 2014). The change in discharge at a catchment's outlet over time, $\mathrm{d}q/\mathrm{d}t$, can thus serve as a proxy for stream network recession and expansion (Biswal and Marani, 2010; Mutzner et al., 2013). When $\mathrm{d}q/\mathrm{d}t$ is negative, it indicates how the upstream stream network contracts and the rate at which stored water, whether in soil or as groundwater, continues to flow into the stream. When it is positive, it indicates network expansion after a precipitation input event, whether rainfall or snowmelt, which has led to subsurface saturation, overland flow, and fast subsurface flow.

In the highly variable Alpine environment, water flow and its properties can mean proliferation or eradication of habitats; thus it determines persistence of biodiversity. Purely disciplinary views to address current global changes and the consequences for habitats and diversity are not sufficient to address this challenge. Future indicators that simultaneously offer information about the hydrological and ecological dimensions of the environment may be critical to get a coherent and complete picture of aquatic systems. Such novel methods should complement existing tools by either improving their accuracy or by revealing additional, complementary insights into a system, such as complementing hydrologic information with biological information. Alpine aquatic systems are especially suited for such an interdisciplinary approach, because the high variability and complexity of physicochemical properties and dynamic flow paths is known to be matched by equally complex and diverse biological habitats that are inhabited by highly specialized organismal communities (Brown et al., 2003; Milner and Petts, 1994; Ward, 1994). For example, midges of the genus Diamesa are characteristic or sometimes even the only inhabitants of rivers dominated by glacial input, resisting the cold temperatures and high turbidity (Milner and Petts, 1994). Contrarily, plecopterans of the genus Protonemoura or the species Leuctra nigra are typical inhabitants of Alpine springs (Staudacher and Füreder, 2007; Hahn, 2000). The same is true for microorganisms, where higher abundances of α-Proteobacteria were found in glacial streams than in other aquatic systems, likely reflecting the trophic status (i.e., the primary productivity) of the habitat (Battin et al., 2004). Finally, algae of the genus Chamaesiphon and the species Hydrurus foetidus are indicative of glacier-dominated sites (Hieber et al., 2001). Consequently, the presence and drift of biological organisms are expected to be not only of high ecological relevance, but also to have the potential to trace connectivity of the stream network, a hypothesis which has been confirmed by some studies on bacteria (Pfister et al., 2009; Lan et al., 2019) and diatoms (Wang et al., 2017).

In recent years, the use of environmental DNA (eDNA) has been shown to be highly powerful for monitoring organisms in aquatic ecosystems (Valentini et al., 2016; Deiner et al., 2015; Bohmann et al., 2014; Altermatt et al., 2020). All organisms, ranging from microbes to vertebrates, leave traces of their genetic material in the environment. This DNA can be collected in environmental samples and is referred to as “environmental DNA” (Ficetola et al., 2008; Taberlet et al., 2012). With a single sampling technique, eDNA samples can detect communities of a broad taxonomic scale by using a high-throughput sequencing approach of a barcoding region, so-called metabarcoding. In an ideal case, the barcoding region allows for identification at the level of species. In stream systems, several studies have shown that eDNA can travel for distances of tens to hundreds of kilometers (e.g., Deiner and Altermatt, 2014; Pont et al., 2018). Thus, stream networks congregate eDNA from upstream areas (Deiner et al., 2016), highlighting its potential to derive flow path information at the catchment level as well as to reconstruct diversity patterns using hydrologic models (Carraro et al., 2020a). Recent work from lowland streams and their wastewater inflows also showed that mixing of different water sources can be traced to a high level of reliability with eDNA (Mansfeldt et al., 2020). While artificially introduced DNA attached to particles has already been used as a hydrologic tracer (Dahlke et al., 2015; McNew et al., 2018; Foppen et al., 2013), naturally occurring DNA has just started to enter into the repertoire of hydrologic tracers (Good et al., 2018; Carraro et al., 2018, 2020a, b). Particularly, it may complement existing physical tracers by not only indicating possible water sources, but also indicating their organismal communities, which are regularly used as a water quality component.

In this study, we tested whether the inclusion of eDNA in a hydrologic assessment of an Alpine system improves our understanding of hydrologic flows while simultaneously delivering biological information. We also explored how this would benefit biologists seeking to quantify Alpine diversity by providing clear recommendations regarding where and when to sample eDNA in river networks for assessments of diversity. To do so, we repeatedly sampled water in an Alpine catchment from spring to summer. We selected 10 sites corresponding to three different water types: low-flow streams that are likely ephemeral, fed by snowmelt, glacier melt, or rain (tributaries); groundwater emerging from the subsurface (springs); and the higher flow stream fed by the other two (main channel). Throughout this study, we will refer to these three water types to encompass the characteristic flows and sources of this catchment and their corresponding aquatic habitats. Firstly, we hypothesized that different Alpine water types carry significantly different eDNA signals that can be used to discriminate between them. Secondly, we hypothesized that hydrologic variability, i.e., the change in streamflow and subsurface water flow as observed through various physicochemical indicators, drives the temporal and spatial eDNA signal. Thirdly, we expected the sum of eDNA from upstream sampling points to shape the eDNA signal at the most downstream catchment outlet. Our study provides an initial assessment of the spatial and temporal heterogeneity of eDNA in an Alpine stream system. We give concrete recommendations for sampling eDNA in this context and identify opportunities but also challenges to its possible use for gaining additional hydrologic insights.

2 Methods

We monitored eDNA in an Alpine catchment in Switzerland in parallel with hydrometeorological observations. On 11 field days in 2017 (for exact dates see Table S1 in the Supplement), we sampled eDNA and simultaneously observed the EC, water temperature, and stable isotope ratios of the water (δ¹⁸O, δ²H). In addition, we continuously monitored discharge at the catchment outlet and meteorological parameters at four stations distributed across the catchment (Fig. 1). We first describe the field site and instrumentation for hydrologic observations. Then, we explain the eDNA sampling and laboratory procedures before describing the data analysis and statistics used.

2.1 Site description

The Avançon is the main stream in the Vallon de Nant catchment (Canton of Vaud, Switzerland; Fig. 1). The 13.4 km² catchment ranges from 1200 m a.s.l. (as defined by a gauging station at 46.25301^∘ N, 7.10954^∘ E) to 3051 m a.s.l. (Le Muveran). It has a dominant north exposition, with almost no direct sun during winter in its upper parts, allowing a small glacier (Glacier de Martinets, 0.36 km²; GLAMOS, 1881–2018) to persist at relatively low elevation. Permafrost is likely present above 2400 m a.s.l. (Giaccone et al., 2019). The streamflow regime is of the nivo-pluvial type (Aschwanden and Weingartner, 1985), with a monthly streamflow peak in June. The stream network in the catchment consists of streams that have been classified as first-, second-, and fourth-order streams according to the Strahler terminology. Correspondingly they have been characterized as having low (<0.05 m³ s⁻¹ and occasionally ephemeral) and medium annual mean discharge (0.05–1 m³ s⁻¹), which we will refer to as tributaries and main channel, respectively; see Schaffner et al. (2013), Pfaundler (2005), and Table S2 in the Supplement. The stream network is fed by melt from the seasonal snowpack including avalanche deposits, which can persist into the summer months, and by drainage from the soil (abundant unconsolidated alluvium) and from deeper groundwater via springs. Permafrost and glacier melt enters the catchment in the upper tributaries.

Figure 1 Map of Vallon de Nant showing sampling sites and stream network. Main channel sites are marked by red squares, tributaries by green triangles, and springs by blue circles; in all three cases, shading lightens with elevation. The two letter codes are used to identify the sites (see also Table S2). The main channel is the Avançon de Nant, which is shown in dark cyan. Predominately ephemeral tributaries of various origins are shown in turquoise. The remaining Glacier de Martinet is shown in light blue. Yellow stars with a red dot in their center indicate where meteorological stations also corresponded to points of rain collection for isotope analysis. Discharge observation occurred at the outlet, which is indicated by the MR/ER sampling site icon. The location of the catchment on an outline of Switzerland with its main hydrologic network is shown to the lower right. On the right are pictures of the main channel (MR/ER, top), a tributary (ST, middle), and a spring (BS, bottom). Colored circles show sampling points, and arrows show direction of water flow. Elevation and river network data courtesy of the Swiss Federal Office of Topography (geodata © swisstopo (DV084371)).

2.2 Hydrologic parameters

2.2.1 Instrumentation and field observation

Water level at the outlet (Fig. 1) was measured using a VEGAPULS WL-61 optical height gauge (VEGA, Schiltach, Germany) over a geophone weir, and discharge was determined using an established rating curve (Ceperley et al., 2018). A network of four distributed meteorological stations (at 1253, 1500, 1780, and 2100 m a.s.l.; Fig. 1) included precipitation and air temperature, mainly measured with a Lufft-WS 300/400 (Fellbach, Germany), and solar radiation, measured with an Apogee incoming radiometer (SP230, Logan, Utah, USA), mounted in an energetically autonomous wireless network with a GPRS connection (Sensorscope Sarl, Lausanne, Switzerland; Michelon, 2017; Ingelrest et al., 2010). At a single point (1253 m a.s.l.), a MADD tipping bucket rain gauge provided back-up measurements (Fallot, 2013, Yverdon, Switzerland,). Data from the Swiss automatic meteorological station network (MeteoSwiss, 2019) are used to fill gaps of the local observations (12 missing days in 2017). Water temperature was recorded every 15 min or more frequently at all but one sampling site (TT; see Fig. 1) with Hobo temperature/light pendants and one temperature/conductivity logger for parts of the study period (Supplement, Table S2; Onset Computer Corporation, Bourne, MA, USA). One WTW-Tetracon 325 (Xylem Analytics, Germany) was installed, immersed at the outlet for the entire year, and logged temperature every minute.

At each point, temperature and EC were measured simultaneously with water sampling using a WTW probe (multi-3510 with a IDS-tetracon-925, Xylem Analytics, Germany). When temperature was not measured during sampling, it was substituted by an observation from the continuous water temperature loggers. When EC was not measured in the field, it was post-analyzed using a glass 6 mm probe in the laboratory (Jenway 4510, Staffordshire, UK).

A total of 135 composite rain samples were collected using funnels flowing into insulated bags at three locations corresponding to the rain gauges (1253, 1500, and 2100 m a.s.l.; see Fig. 1) and emptied approximately weekly or biweekly during the summer seasons between June 2016 and November 2018 (3). In total between February 2016 and April 2018, 199 snow samples were collected in a distributed manner across the catchment, with between 3 and 26 samples per month in the winter, and included the entire snowpack depth when possible (see additional details in Table S3).

The snow-covered area was calculated as a ratio from 0 to 1 over the non-forested part of the catchment using 21 , images and 24 Sentinel-2 images between August 2016 and December 2017 (Michelon et al., 2018). Time steps between the images ranged from 0 to 43 d, averaging 11 d. A linear interpolation of the combined time series was used to provide a fractional area covered with snow per sampling date. Satellite images were taken on average 7 d before or after sampling, though some were on sampling days and others were up to 20 d prior or after.

2.2.2 Flow calculations

Both baseflow and $\mathrm{d}q/\mathrm{d}t$ were calculated from the discharge data at the outlet at the highest recorded streamflow resolution (10 min). Baseflow was defined as the lowest value in a 10 d moving window. On the other hand, the change in discharge over time, $\mathrm{d}q/\mathrm{d}t$, was determined as the slope of the linear regression fitted on a 48 h window preceding the moment that each sample was collected. This 2 d window was chosen to reduce sensitivity to the noise of streamflow records, but it also corresponds to the persistence of eDNA in the environment while also preventing the domination of the diurnal streamflow cycle. For this study, $\mathrm{d}q/\mathrm{d}t$ is a proxy for how the stream network contracts and expands in combination with overland flow. Rather than assuming a particular scaling relationship, we took the $\mathrm{d}q/\mathrm{d}t$ at the outlet to be a good proxy for the upstream network contraction, which is justifiable as significant, and rapid increases in discharge correspond to precipitation events covering large parts of the catchment (i.e., are not concentrated on some small sub-catchment) in the Vallon de Nant (Michelon et al., 2020). Detailed results of a sensitivity analysis regarding the optimal time window length are reported in the Supplement (Fig. S1).

2.2.3 Stable isotopes of water

Water was sampled at the 10 stream sites directly following the sampling of eDNA using 12 mL amber screw vials with a solid polypropylene cap and a silicone rubber/PTFE septa (BGB Analytik, Boeckten, Switzerland) at the sampling points (Fig. 1 and Supplement Table S3). Water was analyzed for stable isotope compositions of oxygen and hydrogen using a Picarro Wavelength-Scanned Cavity Ring-Down Spectrometer (WS-CRDS) 2140-i (Santa Clara, California, USA). Samples were injected a minimum of six times, and the last three measurements were averaged to calculate a raw value. Values were corrected according to a standard curve determined with internal standards, which are regularly calibrated against international standards VSMOW (Vienna Standard Mean Ocean Water) and SLAP (Standard Light Antarctic Precipitation) of the IAEA (International Atomic Energy Agency). Delta units of isotope compositions are reported (Coplen, 1994).

All available samples of precipitation were used to determine a local meteoric water line (LMWL; Fig. S2), which was used to determine the line-conditioned excess (lc-ex, L), or offset from the LMWL:

$$\begin{array}{}\text{(1)}& L={\mathit{\delta }}^{\mathrm{2}}\mathrm{H}-\mathrm{7.82}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{10.47}.\end{array}$$

A positive value of lc-ex falls above the LMWL and a negative one below (Landwehr and Coplen, 2006). The value of lc-ex reflects a combination of the source water and the physical processes that have occurred since precipitation, with more negative values indicating more evaporation and condensation (Sprenger et al., 2016).

2.3 Environmental DNA

2.3.1 Sampling and laboratory procedure

eDNA was sampled by filtering 250 mL of stream water on each of four filter replicates per site using GF/F filters (25 mm diameter, 0.7 µm, Whatman International Ltd., Maidstone, UK) directly in the field at the 10 sites immediately before stable isotopes were sampled. Following filtration, we dried the filter by squeezing air through it, rolled the filter with tweezers, and stored the filter in a 1.5 mL tube following the description of Mächler et al. (2018). If it was not possible to collect 250 mL on a single filter, we used more filters until we reached the goal of 1 L filtered water per site. The tubes were transported on ice and were then stored at −20 ^∘C. We implemented a negative filter control consisting of 1 L Milli-Q water that was previously treated with UVC light, before beginning sample collection in the field, to verify a DNA-free status of materials used. Sites were sampled in all cases except during extreme avalanche conditions, snow cover of the sampling site, or when the sampling site was dry.

We extracted the DNA from the filters after all samples had been collected. We used the DNeasy^® Blood and Tissue kit (Qiagen, Hilden, Germany) following the protocol for animal tissue besides a few changes (see Supplement, Sect. S3.1). In order to reduce biases due to lab procedure, we extracted the samples in a random order. During extraction, we combined two of the four filter replicates, resulting in two extractions per sampling event at the individual sites. We included a negative extraction control, containing a filter previously treated with UVC light, for each batch of extractions, resulting in eight extraction controls. The extracted DNA was stored at −20 ^∘C until further processing.

Generally, eDNA exists in low concentrations in the environment. Thus, to sequence the DNA region of interest, it first needs to be amplified. In our study, we targeted a subsection of a previously identified barcoding gene region for eukaryote species (Hebert et al., 2003), a 313 base pair fragment of the cytochrome oxidase I gene (COI; Geller et al., 2013; Leray et al., 2013, Table S4). We used the Illumina dual-barcoded two-step protocol to amplify and sequence our samples consisting of polymerase chain reactions (PCR), replicates to minimize stochasticity effects, implementation of additional controls, indexing, and pooling (see Supplement, Sect. S3.2, for individual steps).

eDNA laboratory methods are optimized for the detection of low amounts of DNA and thus very sensitive to contamination. In order to reduce false positives, we followed the previously described laboratory precaution for work with eDNA (Deiner and Altermatt, 2014; Deiner et al., 2015; Mächler et al., 2015). Reused field material, such as filter housings and syringes, was soaked 40 min in 2.5 % sodium hypochlorite (i.e., bleach), rinsed with deionized water, and treated with UVC light prior to reuse in the field. Furthermore, four types of controls were sampled during the field (filtration) and laboratory procedure (extraction, PCR negative and PCR positive; see Supplement, Sect. S3.3), which were run alongside the eDNA samples during the laboratory workflow and are later used for data cleaning.

2.3.2 Molecular data processing

The main goal of the bioinformatic analysis is to remove errors due to sequencing techniques and regroup sequences into operational taxonomic units (OTUs, operational taxonomic units). High-throughput techniques produce errors during sequencing, and in order to increase appropriate sequence identification, the resulting reads are cleaned and undergo an initial check for quality (see Supplement, Sect. S3.4). We determined variations of the amplified DNA sequences and reduced the influence of sequencing errors with UNOISE3 (usearch v10.0.240; Edgar, 2016). The detected sequence variants were clustered additionally at 99 % sequence identity to reduce diversity and to account for possible amplification errors. This clustering of the sequences results in so-called ZOTUs (i.e., zero-radius operational taxonomic units). In this sense, a ZOTU is a cluster of DNA sequences that are very similar and can be seen as a rough proxy for a species. As a final step, the ZOTUs were assigned to a taxonomic name if possible by comparing it to taxonomic databases (blast 2.3.0 and usearch v10.0.240, tax filter=0.9). However, databases are still highly incomplete (Weigand et al., 2019), and thus matches cannot be expected at a high rate.

The three types of negative controls produced during the field and laboratory procedure were used to remove noise caused by possible contamination. In a first step, we addressed possible contamination following Evans et al. (2017): for each ZOTU found in one of our negative controls, we calculated a relative frequency by dividing the sum of reads for an individual ZOTU by the sum of all ZOTUs in the negative controls, which was then used as a threshold to clean the field samples (see Supplement Sect. S5). Any ZOTU with a frequency below the threshold was removed from further analysis. We then merged the two eDNA replicates for each sampling point and removed any samples that were below 20 000 reads (2 out of 107 samples), likely originating from errors in the field or laboratory work.

The number of species, also called richness, is a simple measure of biodiversity in ecology. However, richness is strongly affected by the sampling effort (i.e., the higher the sampling effort, the more species that will be found). Therefore, in order to compare diversity measures between sites, data were rarefied, a method to standardize effort (Simberloff, 1978). The number of reads per sample is used as a proxy for sampling effort. Consequently, we rarefied the data to the sequencing depth, or number of reads, of the lowest field sample (26 759 reads). The Illumina MiSeq run resulted in 12.5 million reads and 9858 ZOTUs after the bioinformatics pipeline, of which 2.9 million reads and 9635 ZOTUs remained after removal of contamination and rarefaction (see Table S5 in the Supplement for detailed information).

2.4 Data analysis and statistics

We used multivariate analyses to answer three specific questions: (i) whether the specific composition of eDNA (Sect. 2.4.1) could differentiate between the three water types; (ii) whether hydrologic variability (according to the aforementioned indicators) correlated with the eDNA richness (Sect. 2.4.2); and (iii) whether the main channel's eDNA variation could be used to separate the flow contribution of springs and tributaries (Sect. 2.4.3). All statistical analyses following the initial test were performed using R (R Core Team, 2018, version 3.5.3) and the R package phyloseq (McMurdie and Holmes, 2013, version 1.24.2).

2.4.1 Differentiation of water type by eDNA

In order to test the capacity of eDNA to provide insight into the flow path of the water, we assessed whether the eDNA composition varied according to sampled water type. In the study catchment, we identified three characteristic aquatic environments (tributaries, springs, and the main channel) representing unique habitats, each with corresponding eukaryotic communities. We first used a Kruskal–Wallis test to compare the mean ranks of six variables (δ¹⁸O, δ²H, lc-ex, electrical conductivity, water temperature, and rarefied ZOTU richness) to determine if there were statistically significant differences in richness between the seasonally spread sampling days of the different water types: the main channel (M), springs (S), and tributaries (T) (Kruskal and Wallis, 1952). Second, we tested whether we can distinguish eDNA samples from the three water types based on the identity of the detected eDNA sequences (i.e., ZOTUs' sequence). For this analysis, we used a non-metric multidimensional scaling (NMDS; Kruskal, 1964a) approach, an ordination method that represents pairwise comparisons, or rank orders, of sampling sites in a two-dimensional space. Sites clustered close together are more similar in terms of their ZOTUs than sites further apart, which enabled us to identify if communities detected with eDNA cluster according to water type or hydrologic parameters. We used weak ties (stress type 1), allowing ties to be broken where equal observed dissimilarities exist (Kruskal, 1964a). We calculated the stress of the NMDS according to Kruskal (1964b) as a measure of how well the configuration matches the observed data. The NMDS was plotted for the individual sampling days separately, to observe the development of differences over the sampling season.

Dissimilarity measures used for the NMDS are based on the sequence identities of the ZOTUs and were calculated to do a pairwise comparison between sites. We used an unweighted UniFrac (Lozupone and Knight, 2005) as a dissimilarity input to the NMDS, which is based on presence/absence and the phylogenetic tree (i.e., the evolutionary relationships among ZOTUs based on their sequence similarity) of the ZOTUs. We further fitted variables of interest (elevation, EC, water temperature, total daily solar radiation, baseflow, daily snow cover area, $\mathrm{d}q/\mathrm{d}t$, δ¹⁸O, and lc-ex) onto the NMDS ordination two-dimensional space with a commonly used method that maximized correlation between the ordination and the corresponding environmental variables (function envfit, R package vegan; Oksanen et al., 2007, version 2.5-4), to identify the directions towards which the vectors change most rapidly in the presented ordination space. Next, we focused on ZOTUs that have a taxonomic assignment and used the R package indicspecies to identify candidate taxa that are statistical significantly associated with a particular water type or combinations of those by using a permutation test (De Caceres et al., 2016).

2.4.2 Influence of hydrologic variability on eDNA diversity

Due to the spatiotemporal variability in water sources, Alpine catchments are a mosaic of different water types and habitats, indicated by a high range of physical and chemical processes (Ward et al., 1999; Tockner et al., 2000). We explored the influence of water source contributions on eDNA sampled over a season associated to two variables: the change in discharge at the catchment outlet, $\mathrm{d}q/\mathrm{d}t$ and subsurface exchange as revealed by changes in electrical conductivity (EC) in springs (Laudon, 1997). We used rarefied ZOTU richness, calculated for each individual sampling site and time point as a simple measure of eDNA diversity. We utilized a generalized linear mixed-effect model (GLMM; function glmer in the R package lme4; Bates et al., 2007) with ZOTU richness (R) as the response variable, change in discharge ($D=\mathrm{d}q/\mathrm{d}t$) and the water type (W) as fixed effects, and the sampling site (S) as a random effect. The univariate response variable Y=R was modeled with vectors X=[S] and $\mathit{Z}=[D,W,DW]$ of the explanatory variables associated with the fixed effect (Z) and the random effect (X). We assumed a Poisson distribution and the conditional mean $E\left(Y\right|\mathit{X},\mathit{Z})=\mathit{\mu }=g^{-\mathrm{1}}(\mathit{\eta })$ related to the linear predictor by a log-link function with inverse $h=g^{-\mathrm{1}}$ (Crawley, 2012; Coxe et al., 2009). The identity of the sampling sites was defined as a random effect in an intercept model, to take sampling site variability into account.

We compared the models with and without interaction of the two explanatory variables (D, W). The testing of the interaction allows us to identify if the intercept and slope of richness vary according to the water type, whereas a purely additive approach would only allow us to test whether the intercepts were different, while assuming the slopes were the same. We selected the best model based on a χ² test, classically used to compare model fits. Since the change in discharge, $\mathrm{d}q/\mathrm{d}t$, was negatively correlated with ZOTU richness in springs, we performed a second analysis where we replaced $\mathrm{d}q/\mathrm{d}t$ with EC (E), in the above model. EC is an indicator of the subsurface exchange that defines the water type of springs likely better than $\mathrm{d}q/\mathrm{d}t$. EC was scaled and centered for the model due to large eigenvalue ratios.

2.4.3 Separation of upstream water contribution

In order to examine the contribution of ZOTUs from upstream springs and tributaries to the richness in the main channel downstream, we calculated overlap in ZOTUs. We identified ZOTUs that were unique to samples from one water source, shared by only two, or all three water sources. If a ZOTU was only found in spring samples or in spring and main channel samples, but never in any tributary samples over the whole data set, the ZOTU origin was assigned to springs. If only found in tributaries or in tributaries and main channel samples but never in any spring sites, its origin was assigned to tributaries. The ZOTU's presence in the main channel did not influence this decision, since we assumed the main channel contains eDNA of both tributaries and springs, either because it may be shared due to the similar habitat or due to downstream transport of eDNA (Deiner and Altermatt, 2014; Pont et al., 2018).

Next, we partitioned the ZOTUs in main channel samples into those originating from springs and tributaries, providing a count of how much springs and tributaries, respectively, are contributing to the main channel. We then plotted the relative contributions separately against $\mathrm{d}q/\mathrm{d}t$ and EC, as these were identified as major drivers in our prior analysis of ZOTU richness (see Sects. 2.4.1 and 2.4.2). To confirm this relationship, we computed a generalized linear model (GLM) with the relative contribution (C) as a response variable, $\mathrm{d}q/\mathrm{d}t\left(D\right)$ and the origin of the ZOTU (O) as explanatory variables:

$$\begin{array}{}\text{(2)}& C={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}\cdot D+{\mathit{\beta }}_{\mathrm{2}}\cdot O+\mathit{ϵ},\end{array}$$

where ϵ comes from a quasi-binomial error distribution. This choice results from the fact that (i) binomial error distributions are appropriate for proportional variables; and (ii) our data set shows overdispersion, i.e., higher variability than what can be explained by the standard binomial error distribution. We compared the models with and without interaction of the two explanatory variables and selected the best model with a F test. As before, we also tested a model with EC (E) as n explanatory variable replacing D, as it might be more important to recover the input of springs:

$$\begin{array}{}\text{(3)}& C={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}\cdot E+{\mathit{\beta }}_{\mathrm{2}}\cdot O+\mathit{ϵ}.\end{array}$$

3 Results

3.1 Confirmation of water typology

Annual precipitation in 2017 was 1760 mm, which is close to the long-term average of 1920 mm over the period 1961–2017 (minimum 1470 mm, maximum 2600 mm; MeteoSwiss, 2011). Snowmelt peaked in early June, and snow cover was reduced to 9 % of the watershed area by 1 July (DOY 181) and reached 0 % on 18 August (DOY 229). Further details regarding the hydrometeorological observations can be found in the Supplement, Sect. S4. In general, the three water types (main channel, tributaries, and springs) are significantly different from each other according to their physicochemical characteristics. All comparisons of water temperature were significant (p<0.05) in mean rank difference according to the Kruskal–Wallis test (Table S6), making it the best discriminator of the variables besides eDNA (Table S7). Springs had the lowest and most stable water temperature; tributaries, including the uppermost point on the main channel, UR, showed high water temperature variability and were considerably warmer than springs; and the main channel had intermediary values. Overall, mean EC values were higher in springs than other observed water types and generally higher at lower elevations. The main channel showed mid EC values and a narrow temporal range of EC values, resulting from a mixing of sources (tributaries and springs) throughout the season. Two of the between-water-type comparisons of EC were significant (p<0.05) in mean rank difference according to the Kruskal–Wallis test: main channel to spring (M–S) and spring to tributary (S–T); see Table S6. Values of isotopes in the main channel were not consistently intermediary to those in springs and in the tributaries. We found significance (p<0.05) in the mean rank difference of δ¹⁸O between the main channel and springs (M–S) and the main channel and tributary (M–T) but not between springs and tributary (S–T). We found significance in mean difference of δ²H also in those two comparisons, M–S and M–T, and in one comparison of lc-ex, M–S.

Given the seasonality of the physical processes, the variables fluctuated over the observation period, though statistical confirmation of persistent seasonality was not possible with a single season of data. Enrichment and depletion of the stable isotope concentrations in streamflow are mainly determined by contributions of snow (including snowmelt) versus rainfall (Fig. 2). The alternation of sources is visible as a large range in isotopic ratios (Fig. S3) and is particularly apparent for two of the tributaries (ST and TT) and to a certain extent for the sampling point UR in the main channel.

Figure 2 Discharge (blue, upper), precipitation (red upper), total solar radiation (yellow, middle), snow cover area (cyan, middle), water temperature (purple, lower), and air temperature (green, lower) for the sampling season (2017). Meteorological variables were averaged over the catchment area, whereas water discharge and temperature were measured at the outlet (MR/ER). Shaded areas show the range of values of temperature within a day at the four stations (air, green) and at the outlet (water, purple) and the range of flow at the outlet for discharge and show the error of the snow-covered area estimated from satellite data. Sampling days are indicated as vertical black lines with the corresponding day of the year on the lower axis.

The analysis of the physicochemical tracers confirmed the hypothesis underlying the sampling design, namely that the different sites correspond to different water types with typical seasonal patterns of contributions, even if there was some overlap due to accumulations of different processes. More explicitly, these sites had different contributions from rain and snow inputs at different points in time, as evidenced by stable isotopes of water, on the one hand, and various ratios of surface to subsurface flow paths, on the other hand, as indicated by EC and temperature. Thus, each site had a unique mix of water sources and subsurface contributions, and some sites had characteristics of two types. For example, the uppermost point on the river, UR, seemed to correspond to the EC and temperature of the main channel sometimes but at other times to those of tributaries. Likewise the most upstream tributary, UT, seemed to have an EC approaching that of springs compared to the other tributaries, and the EC of one spring, BS, is consistently lower than that of the other springs (Fig. S3). Accordingly, we expect the seasonal and spatial biodiversity patterns observed at these sites to be at least partly explained by these different origins and water flow paths and, in turn, that eDNA observations might convey information on hydrologic processes.

Figure 3 Variation in δ¹⁸O (a), δ²H (b), lc-ex (c), EC (d), temperature (e), and rarefied ZOTU richness (f) for the individual sampling sites over the field season (given as day of the year). The sampling sites in the legend are clustered based on the water source and within the cluster ordered from low to high elevation. The figure legend corresponds exactly to that in Fig. 1, and the x axis corresponds exactly to that in Fig. 2.

3.2 Differentiation of water type by eDNA

There was a higher variability in the number of ZOTUs (i.e., the richness) found in main channel and tributary sites (i.e., the standard deviation was higher), while sites in springs were more similar (Fig. 3, Fig. S3). All comparisons of rarefied ZOTU richness, M–S, M–T, and S–T, were significant (p<0.05) in mean rank difference according to the Kruskal–Wallis test (Table S6).

Next, we focused on the analysis that included the sequence identity of the ZOTUs. The NMDS plot, which uses pairwise comparisons of multidimensional data to represent data in two-dimensional space, resulted in a stress value of 0.132, indicating that the reduced dimensions represent the data fairly well (Kruskal, 1964a). While spring sites clustered slightly apart from main channel sites, this was not the case for tributary sites that are located at the intersection between the main channel and springs (Fig. 4), indicating that the specific ZOTUs of the main channel are phylogenetic closer to the tributary-specific ZOTUs. The arrows in Fig. 4 indicate the direction of the fastest change in the environmental variables of EC, δ¹⁸O, water temperature, and lc-ex, respectively, which were all significant (Table S8), and their length is proportional to their correlation with the ordination of the genetic data (Oksanen et al., 2007). EC is the most important explanatory variable (i.e., it has the longest arrow in Fig. 4), followed by water temperature, lc-ex, and finally δ¹⁸O. The NMDS plots separated by sampling day (Fig. S5) show that in May and June (DOY 136 and 164), both tributaries and the main channel are more similar to each other compared to later in the season. The similarly low isotope values across all sites (δ¹⁸O=11.9 ‰, range 1.4 ‰ on DOY 164) allow us to attribute early season homogeneity in eDNA to the dominance of snowmelt, which flushes the system and reduces inter-site variability. In contrast, in the late season (DOY 220), there was the highest range of water isotope values (δ¹⁸O, 3.05 ‰) across all sites (Fig. 3), suggesting more mixing of water originating as different types of precipitation and subsequent different trajectories of storage before being released as streamflow. However, the persistent clustering in eDNA suggests that the flood event in early August (see Fig. S5, DOY 215) equalized the detected eDNA community of tributaries and the main channel for this day, while it was spread wider the sampling days before (DOY 207) and after (DOY 220).

A total of 3 % (286) of ZOTUs had a taxonomic assignment to genus level (Supplement Sect. S6), of which 14 ZOTUs were significantly associated with the main channel, 3 with springs, and 26 with tributaries according to the indicator analysis (Table S9). A total of 17 ZOTUs were associated with tributaries and main channel sites. Multiple ZOTUs assigned to the genus of Diamesa were associated with the main channel or both the main channel and tributaries, confirming traditional observations (Milner and Petts, 1994). It is a common challenge of eDNA metabarcoding for the COI gene that only a low proportion of ZOTUs have an associated taxonomic name (here 3 %, representing only 5 % of the reads), which largely results from a lack of reference data (Weigand et al., 2019). Thus, all further analysis were done at the level of ZOTUs, irrespective of whether they had a taxonomic name assigned to them or not, which is suggested to be the future avenue of metabarcoding (Apothéloz-Perret-Gentil et al., 2017; Cordier et al., 2018).

Figure 4 Community composition in the main channel, tributaries, and springs. Non-metric multidimensional scaling (NDMS) analysis of communities based on UniFrac dissimilarities (stress=0.132, indicating a fair representation of the data in the reduced dimensions). The axis represent the dimensions 1 (NMDS1) and 2 (NMDS2), respectively. Ellipses represent a 95 % confidence interval.

3.3 Influence of hydrologic variability on eDNA diversity

The link between environmental variability in the different water sources and detected richness was pronounced: both change of discharge and EC affected detected richness (Fig. 5a respectively b); and the best model included the interaction of the fixed effects for both the model with $\mathrm{d}q/\mathrm{d}t$ and the model with EC (Table S10). The first model, including $\mathrm{d}q/\mathrm{d}t$, revealed that the fixed effects and their interaction were significant (Table 1; see Table S11 for results of random effects). However, water source is a factorial variable, and the individual sources did not show the same response: the intercept (−0.637) and the slope (−2.055) for the spring samples are negative compared to the main channel, indicating an inverse effect of $\mathrm{d}q/\mathrm{d}t$ on the richness detected in springs. The slope for tributary sites is smaller compared to the main channel sites but still resulting in a positive relationship of richness and $\mathrm{d}q/\mathrm{d}t$. The second model, including EC instead of $\mathrm{d}q/\mathrm{d}t$ as a fixed effect, identified all fixed effects and their interaction to be significant, except the interaction of EC and tributary sites.

Table 1 Results from the generalized linear mixed model fit by maximum likelihood (Laplace approximation) for fixed effects; asterisks stand for the interaction. Note that the water source is a factor with three levels (main channel, spring, and tributary). The estimate of the slope (relationship) in the model is indicated.

Figure 5 Rarefied ZOTU richness plotted against $\mathrm{d}q/\mathrm{d}t$ (a) and EC (b). The sampling sites in the legend are clustered based on the water source and within the cluster ordered from low to high elevation. Solid lines indicate significant interactions of fixed effects, and the dashed line indicates significant relationship of the fixed effects and the response variable only (without observed significant interaction). Smoothing lines are based on the GLM function. Note that the scale of the x axis for EC (b) varies by water type.

3.4 Separation of upstream water contribution

Over the whole sampling campaign, 79.2 % of ZOTUs found in tributaries were recovered in the main channel, compared to 77.0 % found in springs and the main channel (Fig. 6). Over all ZOTUs, we found 1374 ZOTUs (14.3 %) that were detected only in spring sites but never tributaries and 3221 ZOTUs (33.5 %) that were exclusively assigned to tributary sites. The best generalized linear model for $\mathrm{d}q/\mathrm{d}t$ included the interaction with the origin of the ZOTUs (Table S12). Only the contribution from ZOTUs originating in tributaries to the main channel had a significant interaction with $\mathrm{d}q/\mathrm{d}t$, while no significant interactions were observed for ZOTUs originating from springs (Table 2, Fig. 6b). In contrast, the best model with EC, which was more correlated with richness in springs, did not include an interaction term (Table 2, Fig. 6c).

Figure 6 Venn diagram for the three sampled water sources; overlapping areas depict shared ZOTUs between water sources (a). Values indicate the numbers of ZOTUs in the given section, and the size of the circles is weighted by the number of ZOTUs. Relative contribution of ZOTUs from tributaries (green) or springs (blue) to samples in the main channel against $\mathrm{d}q/\mathrm{d}t$ (b) and EC (c). The solid line indicates significant interactions of fixed effects, the dashed lines indicate significant relationships of both fixed effects and response variable only (without significant interaction), and the dotted line indicates only the significant relationship of one fixed effect. Smoothing lines are based on the GLM function.

Table 2 GLM results for response variables based on a t test; asterisks stand for the interaction. Note that the origin of the ZOTU is a factor with two levels (spring and tributary).

4 Discussion

Our results highlight that eDNA can provide additional insights regarding the water type and origin when accompanied by observations of naturally occurring physicochemical tracers. Additionally, eDNA delivers information about time points of high connectivity between the sampling sites in a hydrologic highly diverse Alpine system. Thereby, it complements the classical observation of stable isotopes of water, water temperature, and EC. The measure of taxonomic richness assessed by eDNA showed highest values in tributaries, especially compared to springs, and showed less variation over time in tributaries than the main channel. Below, we further elaborate on what we learn about the hydrology of the study site from eDNA observations and on the relationship between eDNA and the two physical variables that were most strongly associated with it, namely the change in discharge $\mathrm{d}q/\mathrm{d}t$ and EC. Overall, in response to our own specific questions and hypotheses, we found that the eDNA signal not only provided unique biological information, but also partially discriminated water types, reflected hydrologic variability, and was determined by upstream contributions.

4.1 Differentiation of water type by eDNA

Our analysis showed that the eDNA composition of the three water types was indeed different but not to a level that made them entirely distinct. In fact, we always expect a portion of the eDNA signal that is non-informative on the water types, and this overlap can be explained by either shared species compositions due to ecological connectivity between sites (Pringle, 2003; Bracken and Croke, 2007) and/or by transport of eDNA between hydrologically connected sites, i.e., downstream in the main channel or laterally through groundwater exchange with the hyporheic zone (Bracken et al., 2013). Even through a hydrologic lens, these three water types are not entirely distinct; for example, groundwater contributes to the flow in the tributaries and main channel. The feasibility of using eDNA data to disentangle different hydrologic contributions in low-land streams, specifically inflows from wastewater treatment plants into natural streams, has recently been demonstrated (Mansfeldt et al., 2020). Thus, in a context of completely distinct water sources, source tracking with DNA from microorganisms is possible but may be more challenging in natural Alpine streams with overlapping signals (Wilhelm et al., 2013).

We identified EC and temperature as the two most important indicators for the community clustering over all samples. This corroborates previous, classical Alpine stream biodiversity work (e.g., Milner and Petts, 1994; Ward, 1994; Brown et al., 2003). Interestingly, phylogenetic community composition of the tributaries is between the main channel and the spring sites, represented in the intermediate position of the clustering in the NMDS. The transport of eDNA is not able to override these differences; otherwise we would expect the main channel to cluster between the tributary and the spring sites. By monitoring over a season, we identified time points of greater and lesser homogeneity among sites according to the eDNA community. This is consistent with the different hydrologic characteristics: while the main channel and tributaries resembled each other more (i.e., were more connected) on days with increased precipitation or snowmelt (e.g., DOY 136, 164, or 215), springs fed by groundwater were more stable over time and showed no reaction. eDNA provides promising new insights into the temporal evolution of the connectivity of the stream network system, a key to assessing dominant hydrologic flow paths.

Future studies might investigate how community composition could be exploited to enhance knowledge about ecosystem processes. However, our analysis was restricted to ZOTUs identity, as currently the assignment of species names is challenging, likely because the database coverage for reference sequences is still limited (Weigand et al., 2019) and might be even lower for high Alpine areas. While ZOTUs do not tell us anything about the identity of the organisms behind, they are sufficient to uncover shared versus unique signals from different water types. The advantage of the ZOTU level is that the analysis is not restricted to a subset of organisms that are commonly used in aquatic ecology, or those that are covered in the databases, but are using an overall “fingerprint” of organisms.

4.2 Influence of hydrologic variability on eDNA diversity

4.2.1 Change in discharge

Despite all variability, there was a positive relationship between the change in streamflow, $\mathrm{d}q/\mathrm{d}t$, and detected eDNA richness in the tributaries and main channel. High-discharge events resulted from either rainfall or snowmelt and cause expansion of the surface channel network, or in other words, increased overland and subsurface flow (Chorley, 1978), a well-established consequence (Biswal and Marani, 2010; Mutzner et al., 2013). Likely, this expansion put water in contact with additional surface areas, likely resulting in the wash-in of terrestrial-origin eDNA and causing an increase in the diversity signal. However, due to the lack of taxonomic information, the exact origin of these ZOTUs (i.e., terrestrial origin or not) cannot completely be identified. Conversely, when streamflow decreased, corresponding to an overall contraction of the network, reduced overbank flow, less exchange with the hyporheic zones, and diminution of channel bank erosion likely led to more characteristic DNA of the sampled water types in our record (Fig. 7). The tributaries can be considered to be the suppliers of the ecohydrologic system, gathering genetic diversity when they contain more water and hydrologic connectivity increases. The richness at spring sites decreased with increasing $\mathrm{d}q/\mathrm{d}t$, perhaps because contributions from precipitation diluted the eDNA signal. This is further supported by the observed positive relationship between $\mathrm{d}q/\mathrm{d}t$ and an increased number of ZOTUs originating from tributaries but not from springs.

Figure 7 Conceptual diagram of relation between $\mathrm{d}q/\mathrm{d}t$ and eDNA richness in the hydrologic system. When $\mathrm{d}q/\mathrm{d}t$ is negative (a), the limited hydrologic connectivity and the limited connectivity between the terrestrial and aquatic habitat led to a segregated, water source-specific eDNA record, whereas when it is increasing (b), exchange is enhanced between the stream and the hyporheic and riparian zones through hyporheic connectivity, overland flow, and groundwater, which leads to higher diversity of suspended eDNA at the outlet.

4.2.2 Electrical conductivity

Detected richness in springs was relatively stable through the seasons (as were stable isotopes and temperature); we even observed a slight negative relationship between richness and $\mathrm{d}q/\mathrm{d}t$. This is not surprising, as the springs have a relatively constant flow and stable surface area. However, EC increased over the season in springs, indicating elevated cumulative contact with subsurface substrate (Malard et al., 1999), which also positively correlated with richness. Our models suggest that, similarly to surface water, richness in springs could be increased by interaction with other habitats, in this case potentially subsurface areas (indicated by the higher EC). However, we found a relatively small contribution of eDNA from springs to the main channel, at least during our sampling season, and no relationship between EC and the contribution of ZOTUs to the main channel from springs. The quantity of the contribution from springs to the main channel, in terms of both eDNA and water, is most likely greater and thus more visible during winter baseflow conditions (Zhang et al., 2018), which were omitted from our study.

4.3 Separation of upstream water contribution

We hoped that eDNA would reveal where water flowed before emerging as surface water, such that it can be used to assess hydrologic connectivity. We tested this by comparing DNA from water types (tributaries and springs) with that of the main channel. There was a reasonably high level of recovery in the main channel, 77 %–79 %, in terms of ZOTUs, and a level, though minimal, of specificity to either springs or tributaries (14 %-33 %, respectively). Potentially, the spatial scale matters, as several studies have shown that eDNA can travel tens to hundreds of kilometers (e.g., Deiner and Altermatt, 2014; Pont et al., 2018). The transport distances in our study system can be assumed to have minor influence on the degradation of eDNA, since they are comparatively small, and accordingly, decay of eDNA can be neglected as a first approximation.

Figure 8 Discharge at outlet decomposed into water origin based on sampled eDNA on 11 d. Contribution of volume from springs is shown in blue, from tributaries in green, and from unknown sources in red. The annual hydrograph is shown for 2017; however the volume outside of the sampling period is not computed and is shown in grey. Vertical lines indicate days when eDNA was sampled.

To calculate the hydrologic source contributions over time, we decomposed the water from the main channel to the two origins based on detected eDNA traces. This has previously be done using stable isotopes of water for example by Ohlanders et al. (2013) to identify the glacial melt contribution to discharge. In order to use eDNA for this calculation, we made the following assumptions: (1) that the percent of the total ZOTUs found in water from the main channel, which were also found in either springs or tributaries, reflects the water flow from those channel types; (2) that the rarefied number of ZOTUs equally represented present ZOTUs as opposed to a bias introduced by the rarefaction step; (3) that each of those ZOTUs is composed of a similar number of species with a normally distributed amount of shed eDNA; and (4) that we can interpolate linearly to compute contributions of water types between sampled days (intervals ranging from 5–35 d). Because the volume of sampled water was consistent, we can estimate that there was a total contribution of 14 % and 50 % from springs and tributaries, respectively over the sampled season (Fig. 8). An annual total of almost 13.5 m³ flowed through the outlet in 2017, of which 9.9 m³ flowed during the eDNA sampling period, and of that, we estimate that 1.4 and 5 m³ were from the sampled or similar tributaries or springs, respectively. The remaining 3.6 m³ must be from other, unsampled sources. Our decomposition of the hydrograph is based on eDNA flow between sources and the main channel. Most hydrologic models base hydrograph decomposition on source area, runoff process, or input magnitudes without a validation method (Weiler et al., 2018). Our theoretical evaluation regarding an eDNA basis for hydrograph decomposition can now be compared to models with the same goal.

4.4 Recommendations and future steps

eDNA appears to offer high potential to not only provide insights into ecological diversity but also to inform interpretations of hydrological flow path connectivity. Future studies using eDNA as a supplemental indicator of hydrologic processes would benefit from more targeted sampling campaigns in terms of hydrologic conditions, including a systematic sampling of discharge recession periods or sampling before and after high-flow events generated by different hydrologic processes (e.g., high intensity precipitation events, snowmelt-induced floods, glacial melt). Although the construction of a complete species inventory was not feasible with our data, with the implementation of the indicator analysis, we did a first step towards the identification of key biota in naturally occurring eDNA in three characteristic Alpine water types. Future studies may increase the collection of water from unique origins (such as glaciers, permafrost, hyporheic zones, or swift channels) to identify characteristic biota that could be used as tracers and lead to a better understanding of their respective hydrologic contributions. Furthermore, future sampling could target single species known to be unique to certain water sources, water types, or processes and assess them for quantity. This would provide the data necessary for a mass balance of fluxes between points by tracing the concentration of eDNA. Importantly, however, it is currently highly debated whether this could be integrated into a metabarcoding approach (e.g., Deiner et al., 2017; Elbrecht and Leese, 2015) due to the influence of various steps in the laboratory procedures that prevent a true assessment of abundance.

In addition, our results expose the complexity of using a metabarcoding approach for riverine biodiversity assessments on a broad taxonomic scale; it is not possible to differentiate between confounding hydrological processes and cycles with opposite effects on diversity. We recommend that future assessments of species richness in Alpine stream networks base their sampling on two aspects, depending on their goals. If the goal is to detect the maximum diversity of a targeted catchment, then the detected diversity can be maximized with a single sampling campaign just after a rain or snowmelt event when the homogenization of the eDNA is highest, and we do not see differences in the species composition according to water type. In contrast, if the goal is to detect habitat-specific diversity or rare species, then sampling after an event resulting in high $\mathrm{d}q/\mathrm{d}t$ would not be recommended. In this case, the various habitats, or water sources, should be sampled when the environment is the most stable. Similarly, it is possible that EC could be used to guide sampling decisions to target specific flow conditions. Finally, future sampling campaigns should have a clear design to answer specific questions about connectivity, i.e., how the sampled eDNA is reflecting the hydrologic and geomorphic versus ecologic connectivity in space and time.

5 Conclusions

Our study is a first evaluation of the spatial and temporal heterogeneity of eDNA diversity in an Alpine stream system and how it relates to observed naturally occurring hydrologic tracers. The genetic data provided a new axis of information regarding the biological diversity, which in some locations and at certain moments in time corroborates the understanding of hydrologic processes gained from observation of electrical conductivity, stable isotopes of water, or water temperature while also bringing new insights regarding the connectivity of water flow. Most importantly, we identified that some genetic signals are associated with specific water types encompassing the characteristic water flow paths and sources of this catchment, especially springs. Furthermore, the genetic signals demonstrate seasonal patterns according to water type. Consequently, eDNA may have a high potential as an additional, complementary variable to be interpreted next to classical tracers in hydrologic studies. eDNA not only provides data that inform understanding of hydrologic processes but also provides insight into the ecological diversity and habitat characteristics. Combined hydro-ecological knowledge is essential for a coherent and complete picture of Alpine aquatic systems that can inform understanding and management of the physical and living environment.