Assessing different components of diversity across a river network using eDNA

Assessing individual components of biodiversity, such as local or regional taxon richness, and differences in community composition is a long‐standing challenge in ecology. It is especially relevant in spatially structured and diverse ecosystems. Environmental DNA (eDNA) has been suggested as a novel technique to detect taxa and therefore may allow to accurately measure biodiversity. However, we do not yet fully understand the comparability of eDNA‐based assessments to classical mor‐ phological approaches. We assessed may‐, stone‐, and caddisfly genera with two contemporary methods, namely eDNA sampling followed by molecular identifica‐ tion and kicknet sampling followed by morphological identification. We sampled 61 sites distributed over a large river network, allowing a comparison of various diversity measures from the catchment to site levels and providing insights into how these measures relate to network properties. We extended our data with historical mor‐ phological records of total diversity at the catchment level. At the catchment scale, identification based on eDNA and kicknet samples detected similar proportions of the overall and cumulative historically documented richness (gamma diversity), 42% and 46%, respectively. We detected a good overlap (62%) between genera identified from eDNA and kicknet samples at the regional scale. At the local scale, we found highly congruent values of local taxon richness (alpha diversity) between eDNA and kicknet samples. Richness of eDNA was positively related to discharge, a descriptor of network position, while kicknet was not. Beta diversity, a measure of dissimilarity between sites, was comparable for the two contemporary methods and is driven by species replacement and not by nestedness. Although eDNA approaches are still in their infancy and optimization regarding sampling design and laboratory work is still needed, our results indicate that it can capture different components of diversity, proving its potential utility as a new tool for large sampling campaigns across hitherto understudied complete river catchments.


| INTRODUC TI ON
Quantifying biodiversity accurately is a long-standing challenge of primary importance in ecology (Dornelas et al., 2013;Gotelli & Colwell, 2011;Whittaker, 1972). On the one hand, there is a fundamental interest to understand the distribution of diversity in time and space and the mechanistic drivers from regional to local scales (e.g., Gaston, 2000;Gotelli & Colwell, 2011;Koleff, Gaston, & Lennon, 2003). On the other hand, the ecological state and functioning of ecosystems is often inherently linked to biodiversity (Pennekamp et al., 2018) and the current loss of biodiversity has potentially large negative consequences on the functions and services of ecosystems (Chapin et al., 1998;Isbell et al., 2017). This is especially relevant for freshwater habitats because they provide crucial ecosystem services such as drinking water, food security, or recreational value to humanity (Cardinale et al., 2012;Dudgeon et al., 2006;Postel & Carpenter, 1997).
In river systems, may-, stone-, and caddisflies (Ephemeroptera, Plecoptera, and Trichoptera; thereafter abbreviated as EPT) are often used as indicators of water quality due to their sensitivity to environmental change, their different preferences for ecological niches, and their relatively well-known taxonomy (Schmidt-Kloiber & Hering, 2015). The presence or absence of certain EPT taxa or their overall richness is highly informative and can be tightly linked to habitat quality (Resh & Rosenberg, 1993). Importantly, they not only describe the current state of a water body, like a single water chemistry samples, but also integrate its changes over time and therefore also inform on the long-term status of the water body. Thus, EPT are at the heart of many freshwater quality assessments around the world and are included in many regulatory frameworks, such as the Water Framework Directive (Directive 2000/60/EC, but see also Borja, Miles, Occhipinti-Ambrogi, & Berg, 2009)
Classically, EPT are collected with a standardized kicknet method (Barbour et al., 1999). Taxa are then identified under a dissecting microscope, which is time-consuming and therefore costly. Taxon richness is the most fundamental approach to estimate biodiversity and is still widely used. Even though the number of taxa is a convincingly intuitive proxy of biodiversity and the basis of many fundamental concepts in ecology, it is a difficult variable to measure accurately (Gotelli & Colwell, 2011;Purvis & Hector, 2000). Diversity can be further divided into different metrics, such as local richness, regional richness, and between-site dissimilarity (also known as alpha, gamma, and beta diversity, respectively). Importantly, the latter metric also includes information on taxonomic identity. The emerging technique of environmental DNA (eDNA) metabarcoding is expected to become a complementary or even replacement method of classical morphological identification (Deiner et al., 2017;Hajibabaei, Shokralla, Zhou, Singer, & Baird, 2011;Lawson Handley, 2015). About a decade ago, the first study was published on the detection of species through DNA in an environmental sample (Ficetola, Miaud, Pompanon, & Taberlet, 2008). With the implementation of high-throughput sequencing, which allowed not only detection of single species but also detection of whole communities, eDNA metabarcoding has been proposed to revolutionize biodiversity assessments (Bohmann et al., 2014;Lawson Handley, 2015;Shokralla, Spall, Gibson, & Hajibabaei, 2012). As with every novel technique, eDNA metabarcoding creates new opportunities but also challenges, especially in terms of recognizing what information it provides and how it compares to previously implemented and established methodologies (Blackman et al., 2019).
Comparisons between eDNA and traditional morphological methods have hitherto mostly focused on either local or regional richness comparisons. In river systems, previous studies have generally detected higher taxon richness with eDNA than kicknet approaches Olds et al., 2016;Valentini et al., 2016). Those results were likely due to eDNA at one location integrating taxon information from upstream reaches via downstream transportation of eDNA (Deiner & Altermatt, 2014;Deiner, Fronhofer, Mächler, Walser, & Altermatt, 2016;Pont et al., 2018). This suggests that traditional techniques represent a more accurate local estimate, while eDNA integrates information across space. In that context, any comparison of the two techniques will be influenced by the spatial scale of the study. However, we still do not know whether findings are directly comparable, complementary, or different across different spatial scales and across different components of biodiversity, due to the specific spatial properties of kicknet being a local sampling technique versus eDNA reflecting more an integrated sample. A detailed understanding is needed to make decisions on how to sample biodiversity in complex landscapes and to compare both local and regional measures. This is particularly relevant in river landscapes, where the typical underlying dendritic network structure is known to affect biodiversity (Altermatt, 2013;Altermatt et al., 2013;Carrara, Altermatt, Rodriguez-Iturbe, & Rinaldo, 2012;Harvey et al., 2018;Tonkin et al., 2018).
In our study, we compared different measures of biodiversity of EPT sampled traditionally (i.e., by kicknet) or by eDNA, using a spatially structured approach that representatively covered a river network in a 740-km 2 catchment. We analyzed how these two different approaches capture the facets of biodiversity at the level of alpha (local site), beta (between sites), and gamma diversity (catchment level). Gamma diversity information was supplemented with all historically available data. We hypothesize that (a) revealed gamma diversity would be similar in the two K E Y W O R D S dendritic networks, Ephemeroptera, freshwater biodiversity, metabarcoding, Plecoptera, Trichoptera contemporary methods; (b) alpha diversity estimated by eDNA is higher compared to kicknet sampling and morphological analysis; and (c) beta diversity between sites is driven by nestedness for eDNA samples but should be less pronounced for kicknet samples. The latter two hypotheses are expected to be driven by the transport of eDNA in the river system. Given that taxon richness, community composition, and taxon identity are among the most commonly studied biodiversity variables, we discuss the design of biodiversity monitoring with eDNA in dendritic river networks.

| MATERIAL S AND ME THODS
We studied a river network in a 740-km 2 catchment containing 61 sampling sites in the upper river Thur in northeastern Switzerland ( Figure 1, Table S1). The catchment comprises three main river stems: Thur, Glatt, and Necker, the latter two draining into the Thur. In a sampling campaign conducted from 11 June to 22 June 2016, we collected eDNA and benthic invertebrate kicknet samples, with the two sampling methods performed at each site within a two-day window. To characterize the position of sites in the network, we extracted stream order, catchment area, and the mean annual discharge data for each site from existing databases (BAFU, 2013(BAFU, , 2014Pfaundler & Schoenenberger, 2013).

| Historical data
We obtained long-term biodiversity data on gamma diversity in our catchment from the Centre Suisse de la Cartographie de la Faune (CSCF). These data include all EPT species ever recorded in the F I G U R E 1 Field sites in the river Thur catchment in northwestern Switzerland (insert at bottom right). Colors are coding for stream order (first-to seventh-order streams: blue, cyan, green, yellow, orange, orange-red, and red, respectively). Black lines indicate the three major subcatchments: Thur, Necker, and Glatt. Data source: swisstopo: VECTOR200 (2017), DEM25 (2003), SWISSTLM3D (2018); BAFU: EZG (2012); Bundesamt für Landestopographie (Art.30 Geo IV): 5,704,000,000, reproduced by permission of swisstopo/ JA100119 whole catchment over the time period 1981-2016. The data are of various sampling origins, but of high quality, thus giving a highly robust and reliable cumulative estimate of EPT genus richness (and the respective EPT genus identity) at the whole catchment scale. The cumulative data consist of 3,467 individual records based on observations of the species in an area of the catchment, which we then converted into genus richness.

| Contemporary kicknet data
We collected benthic macroinvertebrates based on three-minute kicknet sampling applied to three microhabitats present at a given site (Barbour et al., 1999). Leaves and debris were removed from the sample, and the remaining material was pooled and stored in 96% molecular grade ethanol. In the laboratory, all EPT individuals were identified with a microscope to species level. A few taxa, only present as early instar larva or containing cryptic species, were grouped in predefined complexes, subsequently treated at the genus level.
We could not assess EPT taxa at one site due to the loss of a sample.

| eDNA filtration in the field
At each site, we sampled three times 250 ml of river water, each on a separate GF/F filter (pore size 0.7 μm, Whatman International Ltd., Maidstone, UK; for more details, see supplementary materials and Mächler, Osathanunkul, & Altermatt, 2018). We collected eDNA samples about 5-10 m upstream of the kicknet sampling to minimize cross-contamination. We directly sampled water with a disposable syringe out of the water body; for more detailed sampling procedure, see also Mächler et al. (2018). Samples were stored in a Styrofoam box equipped with cooling elements until we came back from the field (no longer than 9 hr). Thereafter, the samples were stored at −20°C until further processing. On each field day, we performed a replicated filter control (FC) that consisted of 250 ml of nanopore water previously treated with UVC light and sealed in the clean laboratory facilities. We filtered the FC in the field before any sampling site was visited in order to check whether the reused material (i.e., filter housings and syringes) was clean. In total, we generated 11 such filter controls, each consisting of three replicates and thus resulting in 33 filter controls. Zero-radius operational taxonomic units (thereafter called ZOTUs) found in at least two replicates of a filter control from the same date were removed from all samples for the further analysis. Further information on filtration, eDNA facilities, and material preparation can be found in the supplementary file.

| Extraction and library preparation
Detailed information about the extraction and library preparation can be found in the supplementary file. In short, we used the DNeasy Blood and Tissue Kit (Qiagen GmbH) to extract the DNA from eDNA samples in a randomized order, also including extraction controls (EC) and filter controls (FC). Finally, DNA was eluted in 75 µl AE buffer at the end of the extraction and cleaned with One Step PCR Inhibitor Removal Kit (Zymo Research). We used the Illumina MiSeq dual-barcoded two-step PCR amplicon sequencing protocol.
We targeted the short COI barcoding region specified by Leray et al. (2013) and Geller, Meyer, Parker, & Hawk (2013). The first PCR was performed with primers that contained an Illumina adaptor-specific tail, a heterogeneity spacer, and the amplicon target site (see Table   S2). On each of the 96-well PCR plates, we included a negative (

| Bioinformatic data processing
After the two successful Illumina MiSeq runs, the data were demultiplexed and the quality of the reads was checked with FastQC (Andrews, 2015). Raw reads were end-trimmed (usearch, version 10.0.240, R1:30nt, R2:50nt) and merged with an overlap of min 15 bp max 300 bp (Flash, version 1.2.11). Next, the primer sites were removed (full length, no mismatch allowed (cutadapt, version 1.12)), and thereafter, the data were quality-filtered (prinseq-lite, version 0.20.4) using the following parameters: size range (100-500), GC range (30-70), mean quality (20), and low complexity filter dust (30). In a next step, UNOISE3 (usearch, version 10.0.240) was used to determine amplicon sequence variants (ZOTUs). UNOISE3 has a build-in error correction to reduce the influence of sequencing errors (Edgar, 2016). An additional clustering at 99% sequence identity was performed to reduce sequence diversity and to account for possible amplification errors in the first PCR. The resulting ZOTUs were checked for stop codons using the invertebrate mitochondrial code, to ensure an intact open reading frame. This resulted in 27 M reads corresponding to 11,313 ZOTUs (Table S3). In a first step, all COI-related sequences were downloaded from NCBI. In a next step, the ZOTUs were blasted against the NCBI COI collection and the top five best blast hits were extracted. We used the R packages "taxize" (Chamberlain & Szöcs, 2013, version 0.9.7) and "rentrez" (Winter, 2017, version 1.2.2) to obtain taxonomic labels based on the accession numbers. The fasta headers of the selected COI sequences were modified with the taxonomic labels in order to index the database. As a final step, the ZOTUs were assigned to taxa using Sintax (usearch, version 10.0.240) and the NCBI COI-based reference. The annotations of taxonomic problematic or missing predictions were verified later using Sintax and the MIDORI references (Machida, Leray, Ho, & Knowlton, 2017).

| Data analysis
Our analysis was based on the following strategy: First, we confirmed that the two Illumina runs could be combined, and second, we cleaned the sequencing data and selected only ZOTUs that were assigned to an EPT order (see detailed information in the supplementary file). For individual sites, ZOTUs were only counted if they were present in at least two of the three independent replicates ( Figure S1), which is a highly stringent assumption and conservative with respect to detection of taxa. Thereafter, we were able to analyze various diversity measures to identify congruence and differences between eDNA and kicknet sampling approaches: (a) gamma, (b) alpha, (c) beta diversity, and (d) between-method community dissimilarity. All statistical analyses were performed using R (R Development Core Team 2008, version 3.5.2).

| Gamma diversity
We used the R package "venneuler" (Wilkinson & Urbanek, 2011, version 1.1-0) to draw Venn diagrams for gamma richness of historic, eDNA, and kicknet data. We used the R package "vegan" (Oksanen et al., 2011, version 2.5-4) to calculate two different taxon accumulation curves of the two sampling methods over the whole catchment: the first one with "random" method and the second in "collector" method where we ordered the samples from down-to upstream.

| Alpha diversity
To compare diversity estimates delivered by the two methods, we used the R packages "phyloseq" (McMurdie & Holmes, 2013, version 1.24.2) to calculate richness. With Pearson's correlation test, we identified whether richness measures at sampling sites correlate between the two contemporary methods. We then tested whether genus richness of the two methods increased with discharge, using a linear model.

| Beta diversity
We used Sørensen dissimilarity as a measure of beta diversity, based on the presence/absence data. With the R package "betapart" (Baselga, Orme, Villeger, De Bortoli, & Leprieur, 2017, version 1.5.1), we calculated Sørensen dissimilarity, which can be further partitioned into nestedness and turnover components, allowing us to distinguish dissimilarity arising through embedment (i.e., loss) or replacement of species (Baselga, 2010;Baselga et al., 2017). To observe how beta diversity measures relate to stream distance, we extracted distances between sites and nodes (junctions where two rivers join) from GIS data (swisstopo) and added these as edge weights. We then constructed the adjacency matrix representing our fluvial network consisting of edges and vertices with the R package "igraph" (Csardi & Nepusz, 2006, version 1.2.4) to calculate stream distance between flow-connected sampling points. We used linear models to test whether distance and the difference in stream order of the compared sites explain patterns in the beta diversity measures. We compared models including a null model (by fitting a constant to test the null hypothesis that the slope is zero), models containing only one of the explanatory variables, or both variables (with and without interaction), and performed model averaging in order to calculate relative variable importance (R package "MuMIn", Barton, 2009, version 1.43.6). We additionally calculated Sørensen dissimilarity and checked for differences among stream orders.
First, we tested for heterogeneous variance with a Bartlett test, and if it was significant, we performed a Kruskal-Wallis test to identify whether there was at least one difference in means. If this was true, then we followed a multiple mean comparison post hoc test on rank sums (R package "pgirmess", Giraudoux, 2018, version 1.6.9).

| Between-method community dissimilarity
We also calculated Sørensen dissimilarity and its two components, nestedness and turnover, between eDNA and kicknet samples for each individual site to detect discrepancies in community composition between the two contemporary methods.

| Gamma diversity
At the regional scale (i.e., the whole catchment level), 96 different EPT genera were historically documented. We found 47 EPT genera with our kicknet samples and 42 genera with our eDNA samples, reflecting 46% and 42% of the historically established taxon richness, respectively ( Figure 2 and Table S4). A high proportion of these sampled taxa were already present in the historic records (94% and 95%, respectively). Thirty-six of these EPT genera were detected with both the kicknet method and the eDNA method, reflecting a 62% overlap of the two methods. Both methods also detected a similar additional proportion of taxa found by one method only, but present in the historic dataset ( Figure 2). Finally, we found four genera with eDNA, kicknet sampling, or both approaches that were not previously listed in the historic data. All of these four genera occur at the border of our studied catchment, and due to more recent range expansions rare, single appearances within the catchment are possible.
Separate taxon accumulation curves (accumulating gamma richness with number of sites included) for the kicknet and the eDNA methods were qualitatively similar ( Figure S2A), based on visual comparison of the curves and their 95% confidence band. This pattern remained similar even if sites were accumulated in the order from the most downstream to the most upstream sites ( Figure S2B).

| Beta diversity
Overall, we found comparable Sørensen dissimilarity between flow-connected sites for eDNA and kicknet. For both methods, the turnover component contributed more to the dissimilarity than nestedness-related components ( Figure 5; see Table S5 for information on detected genera per site), which was more pronounced for eDNA than kicknet. The analysis for the relative importance of variables showed that differences in stream order were generally more important for all beta diversity measures; however, for turnover of both eDNA and kicknet, pairwise distance between sites showed only partially lower importance compared to the differences in stream order (Table S6). Both methods showed significant differences in mean Sørensen dissimilarity among stream orders, and we found significant group differences for the post hoc mean comparisons in stream orders ( Figure 6 and Table S7).

| Between-method community dissimilarity
We found an intermediate discrepancy in the community compositions described by the two contemporary methods (M = 0.48, SD = 0.14, Figure S3), and partitioning this difference into nestedness and turnover components indicates that turnover contributes more to the differences in detected EPT genera.

| D ISCUSS I ON
By using a spatially structured approach covering a major river network, we found similar patterns of diversity measures of EPT sampled with kicknet and eDNA. Using a unique historically assembled overview of cumulative, "true" gamma diversity within the study region, Stream order we were able to put our contemporary samples in a historic context.
We found a quantitatively similar overlap between each contemporary approach and historic gamma diversity, in accordance with other studies comparing eDNA with long-term data in freshwater systems (e.g., Hänfling et al., 2016;Valentini et al., 2016). Given that the his-  across hitherto understudied complete river catchments may be a promising avenue, since minimally trained people without any taxonomic expertise can collect great parts of regional or even landscape richness in a rapid time frame with this method.
We also found a reasonable congruency in local alpha diversity of the two contemporary methods and identified a new dependency of eDNA estimates on discharge level. Overall, local richness estimates of eDNA and kicknet sampling are highly comparable. This comparability may be strengthened by our highly stringent inclusion criteria for eDNA estimates, which are more conservative than those used in previous studies that detected higher richness with eDNA compared to traditional methods (e.g., Deiner et al., 2016;Valentini et al., 2016). In addition, we used a COI barcoding primer targeting a broad taxonomic range, which may also have a lower detection rate for specific taxonomic groups. In near future, the forthcoming design and use of EPT-specific primers or the completion of EPT sequence references should reduce the current drawbacks of eDNA methods regarding taxonomic identification. For eDNA, we found a positive relationship between richness and discharge, which was not the case for kicknet data and could be attributed to eDNA integrating biodiversity detection across space due to downstream transport (Deiner & Altermatt, 2014;Li et al., 2018;Pont et al., 2018). Alternatively, richness indeed increased with discharge (i.e., downstream), but the slightly suboptimal sampling period (summer) and the simplified kicknet protocol we used was an inappropriate technique for sampling larger streams, and thus obscured the relationship. Both aspects reduce detection of genera with the traditional methods, as some species might be missed due to reduced sampling effort or immature larval stages hampering identification.
Surprisingly, the 7th stream order sites diverge from this pattern, potentially because we could only access part of these wide river cross sections, and sampling was restricted to the river edge where less mixing occurs, where more extensive sampling may be needed with eDNA (Bylemans et al., 2018). Also, recent studies indicate that eDNA is not evenly distributed in the water column (Macher & Leese, 2017), and it is still unclear how spatial variance in structures, such as riffles and ponds, is affecting the mixing of the water column and thus the equal detection of eDNA along the water column. We speculate that only in smaller streams (1st to 5th order), one or two samples from the edge or in the center adequately reflect the eDNA distribution across the river transect, while in larger rivers, multiple samples across the cross section may be recommended.
We found that pairwise beta diversity (Sørensen dissimilarity) at the local scale was similarly assessed by the two methods, which is congruent to findings by Li et al. (2018). Despite differences among sites in Sørensen dissimilarity, the turnover component (i.e., species replacement) contributes more to the dissimilarity than nestedness for both methods but is even more pronounced for eDNA than kicknet. Turnover indicates that species are replaced between the sites and implies that transportation of eDNA is not the main mechanisms driving differences; otherwise, nestedness would be expected to be stronger. Barnes & Turner (2016) presented many processes (e.g., degradation, re-suspension, or fragmentation) that influence eDNA in the environment and therefore challenge our mechanistic understanding. Thus, it remains unclear whether the detected differences stem from ecological or methodological variation.
Studies on diversity patterns have long focused on a linear view of streams. But it is now increasingly acknowledged that the underlying network structure plays a significant role in shaping species distributions (Carrara et al., 2012;Harvey & Altermatt, 2019;Harvey et al., 2018;Holyoak, Leibold, & Holt, 2005;Seymour, Fronhofer, & Altermatt, 2015;Tonkin et al., 2018). While experimental laboratory studies and field surveys have found general patterns of diversity distribution in river landscapes, we do not know whether eDNA will lead we detect differences between groups of stream orders for both methods. These differences are mainly between beta diversity of large stream orders and small stream orders, as expected by theory.
Overall, for local composition, we see some discrepancy of the two methods for community composition at specific sites.
This discrepancy is driven by turnover, indicating that detection of different species with the two methods differs and is not due to the integration of eDNA over distance. We see several mechanisms that could cause this difference between the two methods.
First, bias introduced in the laboratory process may have hampered the detection of species that were actually present (e.g., through primer bias (Elbrecht & Leese, 2017), PCR stochasticity , sequencing depth). The design of specific EPT primers may improve comparison and deliver even a better overlap than a universal eukaryotic primer as we used in this study.
Second, DNA shedding rates, densities, and activity can differ between genera (Bylemans et al., 2017;de Souza, Godwin, Renshaw, & Larson, 2016;Sassoubre, Yamahara, Gardner, Block, & Boehm, 2016) and affect the eDNA detection. Also, it remains unknown how habitat preferences of different genera affect detection due to limited mixing or flow of the preferred habitat. Third, it is possible that DNA got resuspended from sediments through the mixing of the water column or other disturbances Shogren et al., 2016Shogren et al., , 2017, resulting in a signal of locally extinct taxa. However, we have no indication from the historic data that genera detected only with eDNA have been absent in the catchment for a longer time.
The distribution of biodiversity and indicator species is of interest to many fields in ecology, from basic research and theory to applied projects of biodiversity conservation and ecosystem assessments. Our work identifies novel opportunities of eDNA as a reliable tool to detect biodiversity patterns, similar to traditional kicknet sampling in riverine networks. Our findings show high robustness of the method, allowing its use for rapid richness analyses and offering the potential to do quick assessments of biodiversity with untrained collectors for BioBlitz campaigns or citizen science projects. However, if the goal is to extend previous biomonitoring or biodiversity datasets with eDNA sampling, the spatial scale must be considered when designing sampling schemes to detect taxon richness.

ACK N OWLED G M ENTS
Data analyzed in this paper were generated in collaboration with the Genetic Diversity Centre (GDC), ETH Zurich. We would like to

CO N FLI C T O F I NTE R E S T
All authors declare that there is no conflict of interest regarding the publication of this article.

Sequencing data is publicly available on European Nucleotide
Archive under the study accession numbers (secondary accession number) PRJEB31920 (ERP114535) and PRJEB33506 (ERP116301).