{'fraction_non_alphanumeric': 0.035320156158740894, 'fraction_numerical': 0.040392241246810734, 'mean_word_length': 5.035621521335807, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decision making and best practices for taxonomy-free eDNA metabarcoding in biomonitoring using Hill numbers

Environmental DNA (eDNA) metabarcoding raises expectations for biomonitoring to cover organisms that have hitherto been neglected or excluded. To bypass current limitations in taxonomic assignments due to incomplete or erroneous reference data bases, taxonomic-free approaches are proposed for biomonitoring at the level of operational taxonomic unites (OTUs). However, this is challenging, because OTUs cannot be annotated and directly compared to classically derived data. The application of good stringency treatments to infer validity of OTUs and the clear understanding of the consequences to such treatments is thus especially relevant for biodiversity assessments. We investigated how common practices of stringency filtering affect diversity estimates based on Hill numbers derived from eDNA samples. We collected eDNA at 61 sites across a 740 km2 river catchment, reflecting a spatially realistic scenario in biomonitoring. After bioinformatic processing of the data, we studied how different stringency treatments affect conclusions with respect to biodiversity at the catchment and site levels. The applied stringency treatments were based on the consistent appearance of OTUs across filter replicates, a relative abundance cut-off and rarefaction. We detected large differences in diversity estimates when accounting for presence/absence only, such that the detected diversity at the catchment scale differed by an order of magnitude between the treatments. These differences disappeared between the stringency treatments with increasing weighting of the OTUs’ abundances. Our study demonstrated the usefulness of Hill numbers for comparisons between data sets with large differences in diversity, and suggests best practice for data stringency filtering for biomonitoring.


39
While initial metabarcoding approaches have mostly focused on mirroring or matching data on the 40 well-established indicator taxa (e.g., Douglas et al., 2012;Elbrecht, Vamos, Meissner, Aroviita, & Leese, 41 2017; Visco et al., 2015), (eDNA) metabarcoding allows us to go well-beyond these indicator groups and 42 has raised the expectation to improve the detection of taxonomic groups that have hitherto been ignored, 43 for example because morphological identification is difficult. Consequentially, eDNA metabarcoding may 44 be offering new opportunities to include currently underused groups of organisms for biomonitoring (e.g., 45 chironomids, Czechowski, Stevens, Madden, & Weinstein, 2020;or oligochaetes, Vivien, Apothéloz-Perret-46 Gentil, Pawlowski, Werner, & Ferrari, 2019) or even extend to groups that are abundant but currently 47 excluded in aquatic biomonitoring (e.g., rotifers, ciliates). In this context, an accurate taxonomic assign-48 ment of eDNA sequences is seen as a crucial but difficult step, because many of these indicator organisms 49 are not completely or inaccurately covered in reference databases (McGee,Robinson,& Hajibabaei,50 2019; Weigand et al., 2019). Also, taxonomic names change over time and entries in databases can differ 51 from currently accepted names. As an alternative, the use of taxonomic-free analysis of eDNA has been 52 proposed. This approach is based on the level of operational taxonomic units (OTU), and does not 53 depend/require any taxonomic association (Apothéloz-Perret-Gentil et al., 2017). An OTU is a cluster 54 of sequences grouped by sequence similarity (Blaxter et al., 2005) and is reflecting a pragmatic proxy for 55 species, but is not necessarily in line with classical taxonomy. Taxonomic-free approaches may facilitate 56 rapid implementation of eDNA in biomonitoring. They are, however, challenging because an OTU has 57 no independent verification of any taxonomic unit, unlike taxonomic names assigned to an OTU that 58 can be compared to classically derived data. Using OTUs thus requires good understanding of the data 59 analysis and the steps involved for validation. Of particular importance, especially for biomonitoring 60 and the assessment of biodiversity, is the decision and implementation of stringency treatments to infer 61 validity of OTUs. Multiple strategies to omit OTUs have been reported (see Deiner et al., 2017). Some 62 of these steps are also frequently implemented for eDNA metabarcoding analysis resulting in a taxa list, 63 but the final interpretation of the data is often made in comparison to previously detected taxa in the 64 study area (e.g., Deiner, Fronhofer, Mächler, Walser, & Altermatt, 2016;Hänfling et al., 2016;Mächler 65 et al., 2019). Unfortunately, implemented stringency treatments are often based on somewhat arbitrary 66 decisions but their effects on the downstream diversity analysis has been rarely investigated. Especially 67 in case of taxonomy-free eDNA approaches, we need better understanding of how our decisions during 68 the analysis are affecting the outcome.

69
After passing stringency selection, OTUs can be used to calculate and compare diversity and rich-70 ness measures. Virtually all biomonitorings and assessments of biodiversity are based on such data 71 giving numbers or divergence in diversity, and describing diversity at the regional, local or between-site 72 scale (gamma, alpha or beta diversity) to inform decision making processes (e.g., biodiversity monitoring Monitoring for the National forest system in USA, Manley & Van Horne, 2006). While this is com-75 mon practice, diversity is often calculated in non-compatible ways across its scales or studies, which is 76 severely limiting comparability between the diversity levels, as well as between studies, and thus chal-77 4 lenges interpretation and cross-comparisons. Furthermore, it is often neglected that alpha, gamma, and 78 beta diversities are often related in an additive or multiplicative way. Consequently, calculation of beta 79 diversity is not independent from alpha diversity and can lead to incorrect conclusions when comparing 80 beta values of regions with different alpha diversities (Jost, 2007). The framework of Hill numbers (Hill,81 1973) has been developed to avoid this problem. Hill numbers mathematically unify the broad array of 82 diversity concepts through incorporating relative abundance and species richness. As a result, all Hill 83 numbers are in the intuitive units of 'species', overcome many of the shortcoming of traditionally used 84 diversity metrics and different diversity levels calculated as Hill numbers can be directly compared to 85 each other. While the framework of Hill numbers is already relatively well-established in community 86 ecology (Chao, Chiu, & Jost, 2014;Chase et al., 2018;Jost, 2007), its advantages have only more recently 87 been suggested for DNA-based analysis of micro-organisms (Kang, Rodrigues, Ng, & Gentry, 2016) and 88 macro-organisms (Alberdi & Gilbert, 2019a).

89
In our study, we investigated how common practices of OTU stringency filtering are affecting diver-90 sity estimates based on Hill numbers derived from eDNA samples. We collected eDNA at many sites 91 across a river catchment, reflecting a spatially realistic scenario in riverine biomonitoring. At each site, 92 we collected three replicates of eDNA samples and used a metabarcoding approach by amplifying the 93 cytochrome c oxidase subunit I (COI) region to uncover a broad spectrum of eukaryotic diversity, con-94 taining both multiple currently used as well as potential future indicator groups. After bioinformatic 95 data preparation with established quality filtering and clustering criteria, we then studied how different 96 stringency treatments applied to the analysis of the eDNA data affect the interpretation and conclusions 97 with respect to biodiversity: The first set of strategies used were based on consistent presence of each op-98 erational taxonomic unit across filter replicates (Alberdi, Aizpurua, Gilbert, & Bohmann, 2018;Beentjes, 99 Speksnijder, Schilthuizen, Schaub, & van der Hoorn, 2018; Mächler et al., 2019). The other strategies 100 were based on a relative abundance cut-off value (Elbrecht & Leese, 2015;Macher et al., 2018;Taberlet, 101 Bonin, Zinger, & Coissac, 2018;Yamamoto et al., 2017) by comparing the data from two different MiSeq 102 runs, and the rarefaction of the data set to the smallest sequencing depth of a sampling site (Taberlet et 103 al., 2018). We investigated how these different selections of stringency treatments are affecting diversity 104 detected at the catchment and the site levels. We conducted all of our analyses in the Hill diversity 105 5 framework, which not only allowed comparison of these strategies, but also gives direct and intuitive 106 comparability of diversities at different levels, crucial to unify and generalize biomonitoring studies.

108
We used a data set specifically collected for the analysis of diversity patterns across a large riverine 109 ecosystem, as commonly surveyed in biomonitoring studies. We investigated how diversity estimates 110 derived by eDNA are affected by stringency decisions during analytical steps after the bioinformatic data 111 preparation. Below, we first describe the collection of the data set and second how we applied stringency 112 decisions during the data analysis. We used data on diversity of eukaryotes estimated by eDNA metabarcoding in a riverine network. To 115 do so, we sampled eDNA at 61 sites in a 740 km 2 catchment of the river Thur (Fig. 1A), northeastern 116 Switzerland, in June 2016. The streams at the selected sites range from first to seventh stream order and 117 cover an elevation gradient from 472 to 1241 m above sea level. All details of collection and preparation 118 of the herein presented data are extensively described in a previous study (Mächler et al., 2019), which 119 analysed the diversity of the taxonomic subgroups of may-, stone-, and caddisflies. Here we give only a 120 short overview of the general field and laboratory methods, necessary for the understanding of the data 121 and the experimental setup (Fig. 1B), and we then explain the analysis of the complete metabarcoding 122 data set covering all eukaryotic eDNA.

123
At each site, we filtered 250 mL of water on a single GF/F glass fiber filter (25 mm diameter, 0.7 124 µm pore size, Whatman International Ltd., Maidstone, U.K.). We replicated the filtration three times 125 (Fig. 1B). After filtration, we transferred the filters to independent Eppendorf tubes and stored them 126 immediately in the dark on ice. At the end of the sampling day, samples were moved to -20°C until 127 extraction. In total, we collected 183 filter samples over the whole sampling campaign. At the beginning 128 of each field day, we produced replicated filter controls (FC) consisting of 250 mL of UVC light treated 129 nanopore water, to estimate possible contamination of used material. Over the 11 sampling days, we 130 produced 33 filter controls. Additional to the eDNA samples, we measured environmental variables (pH,  The NGS data was demultiplexed with the MiSeq Reporter and the read quality was checked with 161 FastQC (Andrews, 2015). Next, raw reads were end-trimmed (usearch, version 10.0.240) and merged 162 (Flash, version 1.2.11). We then removed the primer sites (cutadapt, version 1.12) and quality-filtered 163 the data (prinseq-lite, version 0.20.4). We used UNOISE3 (usearch, version 10.0.240) to determine zero-164 radius OTUs (ZOTUs). ZOTUs are valid operational taxonomic units that are corrected for point errors 165 to retrieve reliable sequence clustering and are further filtered to remove chimeric sequences (Edgar, 2016). 166 To reduce sequence diversity, we conducted an additional clustering at 99% sequence identity. To ensure 167 an intact open reading frame, we checked the resulting ZOTUs for stop codons using the invertebrate 168 mitochondrial code. 169 We cleaned the data by following the description of Evans et al. (2017) due to indication of contami-170 nation (i.e., presence of ZOTUs in negative controls). We calculated the relative frequency of each ZOTU 171 appearing in the negative controls by dividing the total reads of an individual ZOTU in all negative con-172 trols by the total number of reads in the negative controls. We removed all ZOTUs that occurred with We conducted all our diversity analyses in the framework of Hill numbers, which mathematically unifies 181 diversity concepts based on relative abundances (Hill, 1973;Jost, 2007). As such, the sensitivity of Hill 182 numbers towards abundant taxa can be modulated with the parameter q, also known as the order of Hill 183 numbers. A Hill number with an order 0 (q = 0) is insensitive to abundance, and is analogue to the 184 classical species richness measure. A Hill number of order 1 (q = 1) takes the exact relative abundance 185 into account and is analogue to the exponential of Shannon diversity as alpha diversity measure. A Hill 186 8 number of order 2 (q = 2) gives abundant taxa more weight and is analogue to the inverse of the Simpson 187 index. The partitioning of biological diversities into the classical alpha, beta and gamma diversity can be 188 directly computed and expressed as Hill numbers, whereby alpha measures reflect the average diversity of 189 the sampling sites, gamma the overall diversity, and beta the dissimilarity between sites and the overall 190 diversity, for example of a catchment (Jost, 2007). All Hill numbers are calculated and interpreted as 191 units of 'species' (i.e., in our case as units of ZOTUs), can be directly compared across orders (Alberdi 192 & Gilbert, 2019a;Chao et al., 2014;Jost, 2007), and directly allow us to study differences between the 193 stringency treatments. 194 We then applied five different stringency treatments to identify how these affect diversity patterns in 195 the eDNA data, and thus may directly affect interpretation of biomonitoring and biodiversity-assessments.

196
The first three strategies are based on consistent appearance of ZOTUs across replicates (Alberdi et al.,197 2018; Beentjes et al., 2018;Mächler et al., 2019). We assessed how often each individual ZOTU was than PCR replicates (as in Alberdi et al., 2018). The fourth strategy was based on a relative abundance 203 cut-off values (Elbrecht & Leese, 2015;Macher et al., 2018;Taberlet et al., 2018;Yamamoto et al., 2017). 204 We plotted for each ZOTU its relative abundance in the first run against its relative abundance in the 205 second run (Fig. S1) and identified a valid threshold of 0.005%, as abundances between the two runs 206 were relatively similar above this value and it also corresponds to the recommended percentage to use for 207 studies without the positive control of a mock community (Bokulich et al., 2013). We subsequently refer 208 to this stringency strategy as the threshold treatment. As a fifth stringency treatment, we rarefied the 209 data set to the smallest sequencing depth of a sampling site (251,207 reads, all three replicates combined), 210 which we subsequently refer to as the rarefaction treatment. These five stringency treatments resulted 211 in five distinct data sets, that were analyzed subsequently in R (R Core Team, 2018, version 3.5.2) and We were interested in how the different stringency measures are affecting the partitioning of diversity 215 measures at the catchment level. We partitioned for alpha, gamma and beta diversity based on qualitative 216 (i.e., presence/absence, order 0) and quantitative (i.e., abundance based, from order 0.5-2, in steps of 217 0.5) Hill numbers for each of the five stringency treatments (function div part in the 'hilldiv' R package, 218 Alberdi & Gilbert, 2019b, version 1.5.0). As diversity measures tend to increase with increasing sampling 219 effort, ecologists implement rarefaction to correct for this bias (Gotelli, 2001;Gotelli & Colwell, 2011).

220
Classically, rarefaction is based on sample size, which can be represented by the number of collected 221 individuals (e.g., Simberloff, 1978), or the number of sites visited (e.g., Gotelli, 2001). Alternatively , Good 222 (1953, Good 222 ( , 2000 suggested a coverage-based approach, which standardizes samples based on their recovered 223 completeness (i.e., the recovered proportion of the total number of individuals and species of a community 224 represented in the sample). We calculated site-and coverage-based rarefaction and extrapolation curves 225 (Chao & Jost, 2012;Colwell et al., 2012) for the catchment diversity with the 'iNEXT' R package (Hsieh,226 Ma, & Chao, 2016, version 2.0.20). To do so, we utilized incidence data, which are based on relative 227 detection (presence/absence) over the whole catchment rather than relative abundances at the individual 228 sites (Alberdi & Gilbert, 2019a;Colwell & Coddington, 1994). sampling sites, we calculated Hill numbers for each of the five treatments using presence/absence data 232 (order 0) and abundance based data (order 1 and 2). To identify whether there were significant differences 233 in detected diversity between the stringency criteria or not, we performed the following tests: Within 234 each order, we tested whether the variance is equal between the stringency treatments by using a Bartlett 235 test. When significant, we performed a Kruskal-Wallis test to identify whether there was at least one 236 difference in means. If so, we then did a multiple mean comparison post-hoc test based on rank sums 237 (R package 'pgirmess', Giraudoux, 2018, version 1.6.9). We further investigated beta diversity between 238 the individual sites for the orders 0, 1, and 2 (function 'part div', Sørensen-type overlap from the R 239 package 'hilldiv', Alberdi & Gilbert, 2019b) and tested for differences between the mean of the stringency 240 treatments as explained above.

241
To identify how decisions and implementations of stringency criteria are affecting the diversity at 242 individual sites, we selected three different sampling sites (Gertenau, Niederuzwil, and Stocken) from the 243 catchment as proof-of-concept examples. These sites correspond to different stream orders (1, 4, and 7 244 respectively) and reflect distinct habitats, but are geographically still relatively close to each other, and 245 thus do not experience different large-scale environmental drivers. We investigated sample-size based 246 rarefaction/extrapolation curves with abundance data of the individual sampling sites.

248
The two Illumina MiSeq runs resulted in 17,495,861 and 21,026,108 raw sequences, respectively. After the 249 bioinformatic data preparation, 12,619,152 and 14,735,110 merged, primer-trimmed, and quality filtered 250 sequences remained for run 1 and run 2, respectively. In the combined data set, in total 26,544,120 counts 251 assigned to the samples remained for our analysis after the correction for ZOTUs in negative controls.

252
After the application of the five stringency treatments, the resulting five data sets differed mainly by 253 the number of ZOTUs. The biggest difference in numbers of ZOTUs occurred between the threshold 254 treatment and the additive stringency treatment, with the former containing only 12.4% of the ZOTUs 255 found in the latter (Table 1). Differences in counts were relatively moderate. We detected the biggest 256 difference between the additive and rarefied treatments, where the rarefied treatment had about 43% less 257 counts. The rarefaction curves were mostly saturating for sampling sites of the relaxed, strict, and rarefied  The partitioning of the five stringency-filtered data sets into Hill numbers of order 0 resulted in alpha and 262 gamma diversity (corresponding to richness) varying over nearly an order of magnitude (Fig. 2, Table   263 S1). The beta diversity estimated by Hill numbers quantifies how many times the diversity of the entire 264 catchment is richer than the average sampling site and varies greatly between the stringency treatments. 265 Beta values indicated that the difference between sites and the catchment data set was largest for the 266 strict data set, where the regional (gamma) diversity was nearly 12 times higher than the average local 267 (alpha) diversity, indicating large differences between the diversity at the individual sites. Contrary, the 268 richness for threshold treatment was only 2.39 times higher for the catchment than the sites, indicating 269 less differences between the richness of the individual sites.

271
Taking abundances into account (order 0.5 and higher), the partitioning values for all five data sets 272 became more similar when abundances were stronger weighted, indicating that all stringency data sets 273 were dominated by a few, very abundant ZOTUs (Fig. 2). All five stringency treatments ended up with 274 a beta diversity value of around 9 for order 2 ( Table S1). The beta value for the threshold data set was 275 increasing from order 0 to order 2, demonstrating that the ZOTUs of this data set were present in many 276 samples and dissimilarity between sites increased with abundance weighting.

277
Whether the accumulation of diversity with number of sampling sites was saturating or not depends 278 on the used order of Hill numbers (Fig. 3A). For the order q = 0, only the threshold data set was 279 saturating (even at a very low number of sampling sites), while the other stringency treatments did not 280 reach a saturation phase. With increasing order, the shape of the curves started to saturate with lower 281 number of samples and for the order q = 2 all stringency treatments were saturating.

282
The sample coverage of the catchment diversity was highest for the threshold treatment (99.9%), 283 followed by the additive (97.1%), rarefied (95.6%), relaxed (93.0%) and strict with a sample coverage of 284 90.3% (Fig. 3B). The coverage-based rarefaction curves turned towards a positive linear relationship with 285 increasing order of Hill numbers, revealing that with higher abundance weighting the coverage increases 286 only linearly with each additional sample. In biomonitoring studies, the diversity uncovered at the site level is important, for example, to define 289 management priorities between individual sites. When analysing the mean of the alpha diversity, no 290 significant difference were ever observed for the additive and rarefied stringency criteria, no matter how 291 strong abundances were weighted (Table 2, Fig 4). Additionally, these two stringency treatments were 292 always significantly different from the strict treatment. When comparing the mean of the recovered 293 diversities, we observed a dependency on the order of the Hill number (Table 2). The significant difference 294 of the relaxed and threshold treatment to others was depending on the Hill order. With q = 0, both 295 treatments were significantly different to the other treatments, but not each other. With q = 1, the 296 threshold was significantly different from the strict, but the relaxed stringency treatment was not different 297 from any other treatment. For q = 2, both treatments were indifferent from the other treatments. For 298 the beta diversity between sites we also observed different patterns for the various Hill orders. When 299 looking at q = 0 (i.e., Jaccard dissimilarity), all five stringency treatments had significantly different 300 means (Table S2, Fig. S3). With q = 1, three different groups appeared: additive, relaxed, and rarefied 301 were different from strict, as well as from the threshold treatment, but did not differ from each other, 302 and the strict treatment significantly differed from the threshold treatment. For the order q = 2 (i.e., the 303 Morisita-Horn index) there was no significant difference in means observed.

304
Sequence-based rarefaction curves behaved similarly between the three sites, but differed between 305 orders of Hill numbers (Fig. 5). For the order q = 0, the additive and rarefied stringency treatments did 306 not saturate, while the other three treatments did. However, for orders q > 0 the saturation for all five 307 stringency treatments occurred even at a low number of sequences. Surprisingly, the strict stringency 308 treatment resulted in highest alpha diversity for higher Hill orders q > 0 for the site of Niederuzwil ( presence-absence measures, and especially in the context of ecological functions and services provided by 321 an ecosystems, relative or absolute abundances become important. Thus, a coherent understanding of 322 diversity measures needs to consider also abundance-based estimates, and do so in ways that are com-323 parable across and within data sets. This is possible in the mathematical unifying framework of Hill 324 numbers when calculating diversity measures (Alberdi & Gilbert, 2019a;Chao et al., 2014;Jost, 2007).

325
In this study, we investigated how common practices of stringency filtering in the analysis of eDNA 326 metabarcoding data are affecting diversity estimates. The evaluation of different stringency treatments 327 within the Hill number framework showed that differences between treatments scale with the abundance 328 weighting of the ZOTUs. While we observed large differences between the individual stringency treat-329 ments for orders q = 0, the differences decreased or even disappeared with the stronger weighting of

337
The applied stringency measures already showed an effect in the overall number of ZOTUs. The 338 consistent presence in at least two (i.e., relaxed ) or three (i.e., strict) replicates reduced the numbers of 339 ZOTUs to half or even a third, respectively. Such a variation in detectability across replicates might be 340 due to stochasticity in the field (i.e., due to heterogeneous distribution of eDNA, Adrian-Kalchhauser & 341 Burkhardt-Holm, 2016; Macher & Leese, 2017), or in the lab through sub-sampling and stochasticity in 342 the PCR reaction. Especially errors or biases in PCR reactions can contribute to an inconsistency across 343 true replicates (Alberdi et al., 2018;Leray & Knowlton, 2017;Murray, Coghlan, & Bunce, 2015), which 344 cannot be eliminated during bioinformatic data preparation. The use of technical (e.g., PCR reactions, 345 Alberdi et al., 2018;Beentjes et al., 2018) or biological (e.g., Macher & Leese, 2017;Mächler et al., 2019) 346 replicates allows to apply such stringency treatments that are capable to infer validity of ZOTUs. Our 347 study highlights, consistent with other studies, the importance of taking multiple replicates, and to use 348 14 these replicates to get more coherent and robust conclusions on the diversity assessed (Alberdi et al.,349 2018; Leray & Knowlton, 2017).

350
The application of a customized abundance thresholds requires either the performance of at least two 351 separate sequencing runs or the implementation of a mock community (Bokulich et al., 2013), with the 352 latter being complicated for complex and broad taxonomic communities. As shown in our study, the two 353 sequencing runs on the same data set allow the identification of good estimates for a cut-off values of 354 minimal read numbers to be included and should be preferred over the selection of a random threshold. 355 We acknowledge that there are higher costs associated with additional sequencing runs, which may be 356 limiting in some monitoring programs. However, in the long run, monitoring programs are only valuable 357 when they provide reliable data, and thus these additional costs may be well-invested. Alternatively, 358 rarefaction can be applied without further investment in replication of samples or runs, but is often 359 heavily criticized (McMurdie & Holmes, 2014), as it tends to loose rare, but potentially essential, ZOTUs 360 that can reflect important aspects of communities (Ainsworth et al., 2015;Eren et al., 2015;Zhan et al., 361 2013).

362
Our analysis showed that differences, and thus conclusions on diversity, between the five stringency 363 treatments at the catchment level are large, given that the diversity estimates between treatments spanned 364 nearly an order of magnitude when considering richness only. Also, the differentiation between sites and 365 the overall catchment was limited for the threshold treatment for richness, indicating that remaining 366 ZOTUs of the threshold stringency treatment had a high prevalence in all the sampling sites, also reflected 367 in the fast saturation of the catchment diversity with accumulating number of samples. However, the 368 individual sites became more distinct with increasing abundance weighting (i.e., most ZOTUs are present 369 in all sampling sites but differ in abundances), an effect that was similar for all stringency treatments 370 for the order q ≥ 1. The observed effect might be a specific feature of our study, where we expected 371 little differences between the sites as they are connected and do not have substantial gradients in habitat 372 quality between sites. However, such a data set can reflect a biomonitoring study and it is important 373 to be aware that the application of cut-off values for abundances and the simultaneous use of richness 374 measurements can minimize differences between sites.

375
When comparing differences between the treatments at the site level, both alpha diversity and beta 376 dissimilarity estimates get more similar with increased weighting of abundances. Nevertheless, the mean 377 alpha diversity of the strict treatment was always significantly lower compared to the additive and rarefied 378 treatments, regardless how strong ZOTU abundances were weighted. But the situation can differ when 379 focusing on an individual site, as shown in one of the selected sampling location of Niederuzwil (Fig. 5). 380 We generally could observe that differences between the treatments disappeared with higher weighting 381 of abundances. This indicates that weighing abundances may be one way to make data sets more robust  (Beentjes et al., 2018;Buchner et al., 2019). In contrast, our study shows that the eDNA metabarcoding 387 results will differ to quite some extent depending on the used stringency treatment in the data analysis,

Figure 2
Diversity partitioning of the five stringency treated data sets for Hill numbers of different orders. The higher the order, the more weight is put on abundant ZOTUs; order q = 0 reflects presence-absence data. While gamma and alpha diversity behave similar for all the treatments, beta dissimilarity shows differences: For each stringency treatment, the differentiation between the individual sites (alpha) and the overall (gamma) increases with higher weighting of abundances. Hill numbers are presented on a logarithmic scale for better visualization. Exact numbers can be found in Table S1. The colors represent the different stringency treatments: additive (blue), relaxed (green), strict (red), threshold (orange), and rarefied (yellow).

Figure 3
Accumulation curves of Hill diversity based on number of samples (A) and sample coverage (B) over the whole catchment. Diversity for the Hill numbers are calculated for the order q = 0 (left panels), q = 1 (middle panels) and q = 2 (right panels). The solid lines indicate parameter space that is established by rarefaction of the measured data, while dashed lines indicate extrapolated values.

Figure 4
Hill numbers given for the five stringency measures (additive, relaxed, strict, threshold, and rarefied ) and the different orders of Hill numbers (indicated by the exponent of D given in the y-axis). Violin plots depict the distribution of values of the 56 sites. The inserted boxplots are indicating the 25%, 50% (bold) and 75% quantile, respectively. Stringency treatments with dissimilar letters (a, b, and c) are significantly different according to multiple mean comparison test after Kruskal-Wallis applied separately to the individual Hill number orders.

Figure 5
Rarefaction curves of Hill diversity based on number of retrieved sequences for three sampling sites: Gertenau (stream order 1, left panels), Niederuzwil (stream order 4, middle panels), and Stocken (stream order 7, right panels). Diversity for the Hill numbers are calculated for the order q = 0 (upper panels), q = 1 (middle panels) and q = 2 (lower panels). The solid lines indicate parameter space that is established by rarefaction of the measured data, while dashed lines correspond to extrapolated values.