# Morphological vs. DNA metabarcoding approaches for the evaluation of stream ecological status with benthic invertebrates: Testing different combinations of markers and strategies of data filtering

Albin Meyer, Frédéric Boyer, Alice Valentini, Aurélie Bonin, Gentile Francesco Ficetola, Jean-nicolas Beisel, Jonathan Bouquerel, Philippe Wagner, Coline Gaboriaud, Florian Leese, et al.

► To cite this version:

Albin Meyer, Frédéric Boyer, Alice Valentini, Aurélie Bonin, Gentile Francesco Ficetola, et al. Morphological vs. DNA metabarcoding approaches for the evaluation of stream ecological status with benthic invertebrates: Testing different combinations of markers and strategies of data filtering. Molecular Ecology, 2021, Environmental DNA for Biodiversity and Ecosystem Monitoring, 30 (13), pp.3203-3220. 10.1111/mec.15723. hal-03100929

HAL Id: hal-03100929

https://hal.univ-lorraine.fr/hal-03100929v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Morphological vs. DNA metabarcoding approaches for the evaluation of stream
ecological status with benthic invertebrates: testing different combinations of markers
and strategies of data filtering**

4

5

6 Albin MEYER (1), Frédéric BOYER (2), Alice VALENTINI (3), Aurélie BONIN (2,3,4), G.

7 Francesco FICETOLA (2,4), Jean-Nicolas BEISEL (5), Jonathan BOUQUEREL (1), Philippe

8 WAGNER (1), Coline GABORIAUD (3), Florian LEESE (6), Tony DEJEAN (3), Pierre

9 TABERLET (2,7), Philippe USSEGLIO-POLATERA\* (1)

10

11 \* Corresponding author's email address: philippe.usseglio-polatera@univ-lorraine.fr

12

13 (1) Université de Lorraine, CNRS, LIEC, F-57000 Metz, France

14 (2) Univ. Grenoble Alpes, CNRS, LECA, Laboratoire d'Ecologie Alpine, F-38000 Grenoble,

15 France

16 (3) SPYGEN, F-73370 Le Bourget du Lac, France

17 (4) Department of Environmental Science and Policy, Università degli Studi di Milano, Via

18 Celoria 26, Milano 20133, Italy

19 (5) Laboratoire Image Ville et Environnement (LIVE, UMR 7362), Strasbourg, France

20 (6) University of Duisburg-Essen, Aquatic Ecosystem Research, Universitaetsstr. 5, 45141

21 Essen
22 (7) UiT – The Arctic University of Norway, Tromsø Museum, Tromsø, Norway

23 Keywords: bioassessment, WFD, bulk, taxonomy, biotic index, metrics

24
25

26

27

28

29

30

31

32

33

34

35

36

37

Macroinvertebrate assemblages are the most common bioindicators used for stream
biomonitoring, yet the standard approach exhibits several time-consuming steps, including the
sorting and identification of organisms based on morphological criteria. In this study, we
examined if DNA metabarcoding could be used as an efficient molecular-based alternative to
the morphology-based monitoring of stream using macroinvertebrates. We compared results
achieved with the standard morphological identification of organisms sampled in 18 sites
located on 15 French wadeable streams to results obtained with the DNA metabarcoding
identification of sorted bulks of the same macroinvertebrate samples, using read numbers
(expressed as relative frequencies) as a proxy for abundances. We especially evaluated how
combining and filtering metabarcoding data obtained from three different markers (COI: BF1-
BR2, 18S: Euka02 and 16S: Inse01) could improve the efficiency of bioassessment.

38

39

40

41

42

43

44

45

46

One hundred and forty taxa were identified based on morphological criteria, and 127 were identified based on DNA metabarcoding of the three markers, with an overlap of 99 taxa. The threshold values used for sequence filtering based on the “best identity” criterion and the number of reads had an effect on the assessment efficiency of data obtained with each marker. Compared to single marker results, combining data from different markers allowed improving the match between biotic index values obtained with the bulk-DNA versus morphology-based approaches. Both approaches assigned the same ecological quality class to a majority (86%) of the site sampling events, highlighting both the efficiency of metabarcoding as a biomonitoring tool but also the need for further research to improve this efficiency.
47 **Introduction**

Human activities have major negative impacts on freshwater ecosystems with drastic consequences on their biodiversity at the global scale (Dudgeon et al., 2006, Vörösmarty et al., 2010). In this context, important directives have been implemented to assess the ecological status of freshwater systems (e.g. the Water Framework Directive, hereafter abbreviated “WFD”, in Europe [Directive 2000/60/EC]) or to protect their integrity (the Clean Water Act in the USA [Public Law 92-500]). There is therefore a major need for tools that allow large-scale, efficient monitoring of the ecological status of water bodies, with the ultimate objectives of identifying the underlying causes of the observed deterioration in water quality (or habitat suitability) and taking the appropriate measures to improve the ecological status of the monitored ecosystems. Such monitoring usually involves the survey of specific groups of organisms, i.e. Biological Quality Elements (BQEs) in the WFD: fishes, macroinvertebrates, macrophytes, phytoplankton and diatoms. Standard monitoring approaches are generally based on the morphological identification of sampled organisms using harmonized, intercalibrated protocols. Depending on the standards and the studied BQE, these approaches may exhibit several limitations: (i) they can be destructive/invasive, (ii) they can be resource-intensive, i.e. time-consuming and financially expensive, and (iii) they require taxonomic expertise for morphology-based identification, whereas such expertise is continually declining in many countries (Hutchings, 2017; Terlizzi, Bevilacqua, Fraschetti, & Boero, 2003).

DNA metabarcoding is an innovative molecular-based alternative for ecosystem monitoring. This approach consists first in extracting DNA from environment samples of sediment, soil, water, faeces, or directly from community bulks (Hering et al., 2018; Taberlet, Coissac, Pompanon, Brochmann, & Willerslev, 2012). DNA is then amplified using versatile molecular markers and sequenced through high-throughput sequencing. These sequences are compared to those found in reference databases, online and/or developed for a given project (e.g. Baird, Pascoe, Zhou, & Hajibabaei, 2011; Rimet et al., 2016), in order to obtain a
73

list of taxa potentially present in the monitored ecosystem or in the bulk sample of the surveyed community.

74

75

Compared to traditional morphology-based methods, the metabarcoding approach is usually considered as non-invasive (for instance, when directly extracting DNA from water samples; Dejean et al., 2012; Valentini, Pompanon, & Taberlet, 2009), and comparably inexpensive and rapid (Baird & Hajibabaei, 2012; Ji et al., 2013; Taberlet et al., 2012). It can also allow for better taxon identification than morphological expertise (Sweeney, Battle, Jackson, & Dapkey, 2011), and an overall better detection of all the species in aquatic ecosystems (Civade et al., 2016; Valentini et al., 2016). Moreover, metabarcoding also allows a sound estimate of the beta diversity (Ji et al., 2013; Serrana, Miyake, Gamboa, & Watanabe, 2019; Sweeney et al., 2011; Yu et al., 2012), and is a reliable source of information for policy making (Ji et al., 2013). Therefore, metabarcoding has been considered as a potential and credible alternative to morphology-based monitoring for both terrestrial and aquatic ecosystems (Baird & Hajibabaei, 2012; Elbrecht, Vamos, Meissner, Aroviita, & Leese, 2017; Shaw, Weyrich, & Cooper, 2017). In freshwaters, several studies have highlighted the biomonitoring potential of metabarcoding, which could efficiently discriminate streams according to their ecological quality (Gibson et al., 2015; Hajibabaei et al., 2011; Ji et al., 2013; Kuntke, de Jonge, Hesselsøe, & Lund Nielsen, 2020; Mächler et al., 2019; Serrana et al., 2019; Sweeney et al., 2011; Zizka, Geiger, & Leese, 2020). For instance, eDNA information on eukaryotic communities in bottom sediments has been strongly associated to land-use pressure types (Xie et al., 2017), and macroinvertebrate bulk data have been used to infer key gradients of stream condition, including dissolved oxygen, dissolved organic carbon, total nitrogen and conductivity (Emilson et al., 2017).

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

The main objective of this study was to test the ability of metabarcoding, applied to standardized bulk samples of benthic macroinvertebrates, to assess the ecological status of streams based on a large-scale biomonitoring program performed within the WFD context.
99 Using metabarcoding on bulk samples could allow to bypass the organism identification step, 100 which is time-consuming and a source of uncertainty due to operator misidentifications and 101 inter-operator identification variability (Metzeling, Chessman, Hardwick, & Wong, 2003). We 102 also investigated how to improve the bioassesment ability of metabarcoding by testing:

103 (i) combinations of different markers and primer sets, selected for their ability to identify a large 104 range of benthic macroinvertebrate groups. If individual markers can exhibit some bias in the 105 detection of taxa in a given sample (Elbrecht & Leese, 2015, 2017; Piñol, Mir, Gomez-Polo, & 106 Agustí, 2015), using multiple markers which efficiently amplify different taxonomic groups can 107 limit the global bias of detection. However, increasing the number of markers inflates the time 108 and cost of analyses (Clarke, Soubrier, Weyrich, & Cooper, 2014); and,

109 (ii) varying threshold values for metabarcoding data filtering. Namely, we investigated the 110 minimal threshold values for best identity match (the percentage of similarity between a 111 sequence and a barcode found in a reference database) and the number of reads for the 112 identified Molecular Operational Taxonomic Units (MOTUs). Classically, a single threshold 113 value is chosen for best identity matches, either one value per taxonomic level of interest (e.g. 114 species, genus) or a unique value for all the taxonomic levels. For arthropods, a unique 115 threshold value generally close to 97% is frequently found in the literature (e.g. 97% in Elbrecht 116 & Leese, 2017; Serrana et al., 2018; 2019; Yu et al., 2012; 97.5% in Carew, Kellar, Pettigrove, 117 & Hoffmann, 2018; 98% in Lobo, Shokralla, Costa, Hajibabaei, & Costa, 2017). Low 118 abundance (in terms of number of reads) MOTUs, e.g. exhibiting a relative abundance lower 119 than 0.003% (Elbrecht & Leese, 2017) or 0.005% (Bokulich et al., 2013; Carew et al., 2018) 120 in a given sample, are usually excluded from metabarcoding data, as they are considered 121 unreliable (e.g. Elbrecht & Leese, 2017; Elbrecht et al., 2017).

122
123 The selection of genetic markers (e.g. mitochondrial vs. nuclear, single vs. multicopy,
124 fast evolving vs. conserved, protein-coding vs. ribosomal), the respective primer sets as well
125 as the threshold values are important decisions to take in a metabarcoding experiment, but
126 for which no clear guidelines exist (but see Bokulich et al., 2013). Such decisions may have a
127 strong impact on the DNA-based bioassessment, which could be exacerbated by the high
128 phylogenetic diversity of benthic invertebrates in a single sample. For example, taxonomic
129 groups that are not as well amplified as other groups by a given primer set or marker could be
130 filtered out of the dataset due to their low numbers of reads. The amplification rate can indeed
131 vary greatly among the major taxonomic groups in benthic invertebrate assemblages, even
132 when accounting for biomass (Elbrecht et al., 2017).

133 We systematically studied the efficiency of the three selected markers and of all the
134 possible combinations of these markers, by comparing bulk-DNA and standard results through
135 the lens of the French Multimetric Invertebrate Index (I₂M₂; Mondy, Villeneuve, Archaimbault,
136 & Usseglio-Polatera, 2012). The I₂M₂ was designed as a WFD-compliant index for the
137 invertebrate-based ecological assessment of French wadeable streams. The standardized
138 taxonomic levels needed for calculating the I₂M₂ are mainly the genus level, except for Diptera,
139 Hirudinea and Turbellaria (family level) and Nematoda or Oligochaeta identified as such
140 (standard XP T90-388; AFNOR, 2010). We tested the usefulness of metabarcoding for
141 assessing and discriminating the ecological status of streams, based on a set of streams with
142 a wide range of ecological features.

143
144 **Material & Methods**

145

146 *Sampling sites and data acquisition*

147 Sampling was conducted in 18 sites on 15 streams belonging to five different stream
types defined according to stream order (Strahler, 1957) and French hydroecoregions
148 (Wasson, Chandesris, Pella, & Blanc, 2002; Wasson, Chandesris, Pella, Sauquet, & Mengin,
149 2006) (Fig. 1; see also Table S1 in the Supplementary Material). Streams were selected from
150 two national networks, surveying (i) reference sites (Réseau de Référence, RdR, about 400
151 sites) and (ii) the mean ecological quality of French streams via the long-term survey of a large
152 selection of sites (Réseau de Contrôle de Surveillance, RCS, about 1500 sites). These
153 surveys have allowed gathering a lot of information on chemical and hydro-morphological
154 pressures impairing water quality and habitats since 2007 (Larras et al., 2017; Mondy et al.,
155 2012). The selection of sites was based on three criteria: their geographic origin (two different
156 hydro-eco-regions x nine streams), pressure intensity and category (water quality degradation
157 or hydrological alteration), and stream type (Wasson et al., 2002). Pressure intensity ranged
158 from "very low" impairment corresponding to "Least Impaired River Reaches" (LIRRs;
159 following Dolédec & Statzner, 2008; Mondy et al., 2012; Statzner, Bady, Dolédec, & Schöll,
160 2005), to "moderate" or "strong" impairment, both corresponding to "Impaired River Reaches"
161 (IRRs). We selected IRRs which had been impaired by only one main pressure category
162 (water quality degradation or hydrological alteration) over the 2007-2012 period. Sites were a
163 *priori* grouped by triplet within a given stream type, including one LIRR and two IRRs, one IRR
164 exhibiting a moderate impairment and the other a stronger impairment level.
165

166 Macroinvertebrate community sampling occurred in autumn 2014 and spring 2015.
167 Thirty-six macroinvertebrate field sampling events were done following a standardized
168 protocol (French standard XP T90-333; AFNOR, 2009), commonly used in France in the
context of the WFD. This protocol advocates the sampling of macroinvertebrates in eight dominant habitats (i.e. with an individual share of at least 5% coverage at reach scale) and four marginal ones (i.e. with an individual share of less than 5% coverage) in three successive phases (Fig. 2A). Each habitat is characterized by its substrate type (among twelve categories) and its superficial current velocity range in front of the substrate (among four categories).

Sampled marginal substrates (phase A; samples 1-4) and the first four dominant ones (phase B, samples 5-8) are selected according to a decreasing gradient of "hosting capacity" (i.e. their ability to support a rich and diverse invertebrate assemblage; this gradient is defined by the norm XP T90-333). The four other samples performed on dominant substrates (phase C, samples 9-12) are made in proportion to their individual benthic coverages. A Surber sampler (net mesh size = 500 µm, opening area = 1/20 m²) was used to sample macroinvertebrates in each habitat. Samples were preserved in undenatured alcohol (ca. 70% final concentration), for up to one year. In laboratory, organisms were sorted, numbered and identified at the standardized taxonomic level (standard XP T90-388; AFNOR, 2010). Some individuals, difficult to identify at the required taxonomic level (e.g. early instars or organisms altered during the sampling process), were identified at the best taxonomic level possible (i.e. 29 individuals from eight taxa over 89157 individuals catch during the whole study). After taxonomic identification on morphological criteria, all the sorted organisms, i.e. the "bulk samples" (N = 432), were stored in 95° undenatured alcohol (for up to one year) and sent for metabarcoding.

## $I_2M_2$ index calculation

Faunal data (i.e. abundance per taxon x phase x site) were combined per site, as indicated in Fig. 2A, to calculate the values of the five individual metrics aggregated in the $I_2M_2$ index (Mondy et al., 2012). These five metrics are the total taxonomic richness, the Shannon-Weaver diversity index (Shannon & Weaver, 1963), the Average Score Per Taxon (ASPT;
Armitage, Moss, Wright, & Furse, 1983), and the proportions of ovoviviparous (Ovoviviparity) and polyvoltine (Polyvoltinism) organisms in the invertebrate assemblage. The I₂M₂ index was constructed and calibrated over 10 chemical and 7 hydromorphological pressure categories (see Mondy et al., 2012, for further details). One sub-index per pressure category is first calculated as the weighted mean of the values of the five metrics, with the weights equal to the ability of each metric to discriminate between LIRRs and IRRs for this pressure category (quantified by its "discrimination efficiency", DE; Ofenböck, Moog, Gerritsen, & Barbour, 2004). Then, the final I₂M₂ index value is calculated as the arithmetic mean of the 17 sub-index values (Mondy et al., 2012).

203

204 *Metabarcoding*

Bulk samples were homogenized using an IKA ULTRA-TURRAX Tube Drive control system with sterile 20 mL tubes and 10 steel beads (5 mm Ø) by grinding at 4,000 rpm for 15 min (IKA, Staufen im Breisgau, Germany). Complete samples were then incubated overnight at 56 °C in 5 mL of lysis buffer (Tris-HCl 0.1 M, EDTA 0.1 M, NaCl 0.01 M and N-lauroyl sarcosine 1%, pH 7.5–8.0). Extractions were then completed using the DNeasy Blood Tissue Kit (Qiagen GmbH, Hilden, Germany), according to the manufacturer's instructions. DNA extracts were recovered in a total volume of 300 µL. Two DNA extractions were performed per bulk sample. Negative extractions without samples were systematically performed to monitor possible contaminations.

Three primer pairs respectively corresponding to three different markers were used for each sample, Inse01 for a mitochondrial 16S rDNA region (Elbrecht et al., 2016; Taberlet, Bonin, Zinger, & Coissac, 2018), Euka02 for a nuclear 18S rDNA region (Guardiola et al., 2015) and BF1-BR2 for the Cytochrome c Oxidase I (COI) region (Elbrecht & Leese, 2017). DNA amplifications were performed in a final volume of 20 µL, using 2 µL of extract DNA as
219 template. The amplification mixture contained 10 µL of Applied Biosystems™ Master Mix
220 AmpliTaq Gold™ 360, 0.2 µg/µL of bovine serum albumin (BSA, Roche Diagnostic, Basel,
221 Switzerland) and 0.5 µM of each primer for COI and 16S primers or 0.2 µM for 18S primers.
222 Two PCR replicates were amplified for each DNA extraction and each primer pair, for a total
223 of four PCR replicates per bulk sample. The primers were 5'-labeled with an eight-nucleotide
224 tag unique to each replicate (with at least five differences between any pair of tags) allowing
225 the assignment of each sequence to the corresponding sample during sequence analysis. The
226 PCR mixture was denatured at 95°C for 10 min, followed by 35 cycles of 30 s at 95°C, 30 s at
227 52°C for COI and Inse01 or 45°C for Euka02 and 1 min at 72°C (1m 30s for COI), and followed
228 by a final elongation at 72°C for 7 min. Negative PCR controls (ultrapure water, with 12
229 replicates as well) were analyzed in parallel to the samples to monitor possible contaminations
230 during the PCR step.

231 After PCR amplification, PCR products from the same marker were combined in equal
232 volumes and purified using the MinElute (Qiagen GmbH, Hilden, Germany) purification kit.
233 Purified amplicons were checked by high resolution capillary electrophoresis (QIAxcel System,
234 Qiagen GmbH, Hilden, Germany) and sent to Fasteris (Geneva, Switzerland) for library
235 preparation and sequencing. Libraries were prepared according to the PCR-free MetaFast
236 protocol (Taberlet et al., 2018; for further details, also see:
237 https://www.fasteris.com/dna/?q=content/metafast-protocol-amplicon-metagenomic-analysis),
238 which limits chimera formation. The Inse01 and Euka02 amplicons (three libraries each) were
239 sequenced on a HiSeq 2500 platform (Illumina, San Diego, CA, USA) with a paired-end
240 approach (2 × 125 bp), while the COI amplicons (two libraries) were sequenced on a MiSeq
241 platform (Illumina, San Diego, CA, USA) producing 2 × 250 bp paired-end reads.

242

243 Workflow for metabarcoding data
Sequences were processed using the OBITools software (Boyer et al., 2016). Each pair of raw reads was paired end merged with *illuminapairedend* to recover the full amplicon sequences. Pairs of reads that did merge with an alignment score above 40 (equivalent to align 10bp of maximal quality on both read ends) were discarded. For the Euka02 primer, as the read length did not allow to recover the full amplicon sequences for important taxa such as Gammaridae, Coleoptera and Trichoptera, pairs of reads whose ends could not be aligned (score < 40) were concatenated and kept separately for further processing. Recovered amplicon sequences were then assigned to their respective sample with *ngsfilter* and dereplicated to get MOTUs with *obiuniq*. MOTUs were then aligned against dedicated reference sequence databases (Ficetola et al., 2020) for each primer pair using *ecotag*. The Euka02 barcodes kept as concatenated sequences were processed separately and the alignment score was based on combining the alignments obtained for both ends of the barcode.

After the taxonomic assignment step, metabarcoding data were subjected to subsequent steps of preparation and filtering before the ultimate step of I₂M₂ index calculation (Fig. 2B). Discordant PCR replicates (i.e. that did not cluster when compared to other technical replicates of the same bulk sample) were identified using an iterative process. This iterative process was akin to minimizing the intra-sample distances (between PCR replicates of a given bulk sample) while maximizing the inter-samples distances. At each step of this “min-max” process, the algorithm identified which replicates were the most discordant, i.e. the replicates that were too distant from the other replicates of the same bulk sample, and it removed them before iterating. The assumption of this process is that PCR replicates from the same sample should be more similar to one another than to PCR replicates from other samples. To be more specific, at each iteration all the PCR replicates were projected on a 2D space using a Correspondence Analysis (implemented in the ade4 R package; Dray & Dufour, 2007) based on their square root transformed counts. Euclidian distances between all the PCR replicates
in this 2D space were then computed. Distances were partitioned depending on whether or
not they involved PCR replicates of the same sample or different samples. These two distance
distributions were then compared to pinpoint outlier distances for PCR replicates of the same
sample. PCR replicates responsible for these distances were removed and the remaining PCR
replicates were used again for a new iteration until no PCR replicates had to be removed.
Respectively, we removed 43.4, 7.5 and 9.7% of the replicates for the primer sets COI, Euka02
and Inse01. PCR replicates were thus available for 68.8, 100 and 97.9% of the bulk samples
amplified with COI, Euka02 and Inse01, respectively (mean numbers of available replicates
per bulk sample = 2.3, 3.7 and 3.6, respectively for COI, Euka02 and Inse01). The remaining
PCR replicates were then pooled together, by summing the numbers of reads per MOTU found
across all the replicates, for a given bulk sample.

Next, MOTUs were filtered based on the best identity percentage (i.e. the percentage
of similarity between the MOTU sequence and the closest one identified in the reference
database; Ficetola et al., 2020), then, after further sample pooling (Fig. 2A), they were filtered
based on their total number of reads (Fig. 2B). For the filtering step based on best identity
values, we tested several thresholds, between 80 to 100%. After this first filtering step, we
standardized the MOTUs data using the reference list of taxa taken into account for the
calculation of the I₂M₂ index value (norm XP T90-388; AFNOR, 2010). As a result, the reads
were either pooled by genus, sub-family or family level according to the taxa (or even at a
higher taxonomic level; e.g. for Oligochaeta), for the MOTUs which could be aggregated to
the taxonomic level requested by the standard. MOTUs identified at a taxonomic level too high
for being used in the bioevaluation process have been removed from the dataset (e.g. MOTUs
such as Metazoa, Neoptera or Holometabola). Figure S1, available in the Supplementary
Material, shows the remaining total number of reads per PCR replicate at this step of the
process. We also tested how a uniform sequencing/read depth would influence the
performance of the markers, by filtering out MOTUs with read abundance equal to or lower
296 than 0.003% per PCR replicate at this step of the process. Preliminary analyses had shown
297 that this added filtering step has only a negligible effect on the biomonitoring results compared
298 to the other filtering steps (namely for best identity and minimal number of reads; see the next
299 paragraph). Therefore, results presented in this study do not include this optional filtering step.

300 According to the standard used for the index calculation (XP T90-333; AFNOR, 2009;
301 Fig. 2A), the numbers of reads per MOTU were further pooled together within each group of
302 four samples corresponding to the three successive phases of the field sampling protocol (A,
303 B and C, respectively; Fig. 2A). A table summarizing the main identified MOTUs by each
304 marker is available in the Supplementary Material (Table S2). Taxa were then filtered
305 according to their total number of reads in each pooled sample for each marker. We went for
306 a naive and global approach, and systematically tested several threshold values of minimal
307 number of reads, ranging from 1 to 100 reads, five reads by five reads (i.e. 1, 5, 10, 15, up to
308 100).
309

309 As we tested all the possible combinations of threshold values for both filtering steps,
310 we obtained a total of 441 datasets for each marker, i.e. 21 threshold values for the best
311 identity percentage multiplied by 21 threshold values for the minimal number of reads. After
312 this step, read data were transformed either in relative frequencies (RF) or in
313 presence/absence (PA).

314

315 *Marker combinations*

316 For each mode of data expression (either RF or PA) we tested the individual marker
317 but also all the different combinations of two (2) or three (3) markers. In order to limit the
318 number of tested marker combinations over all the possible combinations of tested threshold
319 values, we first identified the threshold values that could best allow us to maximize the
correlation (i.e. exhibiting the highest adjusted-R²; see next section) between the results obtained with both the standard and the bulk-DNA approaches. Thus, we first selected the three "best" threshold values of best identity for each individual marker. Then, we selected the five best threshold values for the minimal number of reads, based on the datasets already filtered with the three best identity thresholds previously identified. As a result, we obtained 15 datasets (one per combination of best identity threshold x read minimal number threshold) for each individual marker, 225 datasets (15 x 15 combinations) for each pair of markers (e.g. [COI + Euka02]), and 3375 datasets (i.e. 15 x 15 x 15 combinations) for the combination of three markers ([COI + Euka02 + Inse01]).

When combining data from different markers, a taxon was considered as "present" if it was present at least in one of the two or three datasets included in the marker combination. For data coded in RF, for a given taxon, RF values were averaged over all the data provided by the two or three combined datasets.

## Statistical analyses

The values of the I₂M₂ index and its five individual metrics were calculated for all the available site sampling events, combinations of markers and selected thresholds, using relative frequencies of reads as a proxy for abundances (i.e. for all the metrics except total taxonomic richness and ASPT). For data expressed in presence/absence, the abundance was fixed as equal to one for each identified taxon. The values of the bulk-based index (B-I₂M₂) and its metrics were compared to the values provided by the standard approach ("Morphology-based" I₂M₂; M-I₂M₂ for the index). These values were expressed as Ecological Quality Ratios (EQR; range [0;1]). Reference values, known as "BEST" and "WORST" and needed for the expression of all the metric values as EQRs, were based on the revised I₂M₂ construction dataset (see next paragraph) for both the standard and the bulk-DNA approaches. We also
calculated the discrimination efficiency (DE; Ofenböck et al., 2004) of the I₂M₂ index on each dataset. Here, DE is the relative frequency of IRRs exhibiting I₂M₂ values lower than the first quartile of the distribution of the I₂M₂ values in the LIRRs.

The I₂M₂ was updated in 2016 thanks to a new and bigger available dataset (with 10,074 sampling events, versus only 4,132 sampling events in the original dataset used in Mondy et al., 2012). Based on this new dataset, more robust reference and DE values were (re)calculated, and combinations of phases used for individual metric calculation were revised in order to optimize the discrimination efficiency of metrics (see Fig. 2A). In this study we used this revised version of the I₂M₂. (ref?)

Linear regressions were calculated between the two sets of I₂M₂ index values (or individual metric values) obtained with the standard and bulk-DNA based approaches. For each regression we used the adjusted R² as a measure of the variance explained by the model (Crawley, 2007). The distributions of index values provided by both methods were statistically compared over all the site sampling events, with Wilcoxon signed rank tests for paired data. Friedman rank sum tests were used to identify whether or not the tested thresholds for best identity and minimal number of reads (after having selected the best three thresholds for best identity) led to significant differences in adjusted R² values over the whole range of tested values. All statistical analyses were done in R (version 3.6.3; R Core Team, 2020).
363 **Results**

364

365 *Taxonomic identification*

366 One hundred and forty morphotaxa were identified in the 36 site sampling events (Fig.

367 3). The bulk data filtered with the lowest threshold for best identity (80%), allowed to recover

368 75, 66 and 57 of the morphotaxa, for the markers COI, Euka02, and Inse01, respectively.

369 Twenty-eight morphotaxa were independently recovered by all three markers. Forty-one taxa

370 identified on morphological criteria were not recovered by at least one of the markers. Of these

371 forty-one taxa, six were not referenced as such in our marker-specific reference databases:

372 the phylum Nematoda, the family Rhagionidae, and the genera *Cyphon*, *Agriotypus*,

373 *Capnioneura* and *Lasiocephala*. Twenty-eight taxa were obtained by metabarcoding but not

374 based on morphological identification: 12 taxa for COI, 15 taxa for Euka02 and 9 taxa for

375 Inse01.

376

377 *Relative frequencies vs. presence/absence*

378 We first examined how the mode of data expression (RF or PA) affected the efficiency

379 of metabarcoding-based bioassessment. Whatever the mode of data expression, we generally

380 found significant and positive linear relationships between B-I₂M₂ and M-I₂M₂ values. For the

381 primer sets COI and Euka02, the relationship was significantly higher on RF data than on PA

382 data (Wilcoxon signed rank test for paired data; p < 0.001 for both primer sets; Fig. 4A). For

383 Inse01, adjusted R² were not significantly different between RF and PA data (Wilcoxon signed

384 rank test for paired data; p > 0.05). Nevertheless, the primer set Inse01 exhibited the highest

385 values of adjusted R² with PA data, with a maximum of 0.738. The maximum adjusted R² for

386 the primer sets COI and Euka02 were respectively 0.524 and 0.671, with RF data. Based on
these preliminary results, we chose to focus on data expressed in relative frequencies (RF) for the remaining analyses.

389

390 *Threshold selection for markers*

We used adjusted R² for examining how the selected threshold value for best identity modified the bioassessment efficiency of each marker. The relationship between the adjusted R² and the tested thresholds for the best identity greatly varied among markers (Friedman rank sum tests; p < 0.001 for each primer set; Fig. 5). For COI, the thresholds equal or lower than 84% of identity provided the best values of adjusted R² (medians > 0.50; Fig. 5A). Adjusted R² values decreased with increasing identity thresholds from 85 to 100%, albeit with a small increase of R² for identity comprised between 96 and 99%. For Euka02, the adjusted R² slowly increased with the threshold values until 96% of identity, with all the medians higher than 0.61 (Fig. 5B), while quickly decreased (median values close to 0.49) at very high values of identity (97-100%). For Inse01, the adjusted R² median values globally ranged between 0.58 and 0.62 for thresholds between 80 and 97% of identity. The highest R² values were observed with identity values of 98-100% (Fig. 5C). Based on these results, we selected the best identity threshold values of 80, 81 and 82% for COI; 94, 95 and 96% for Euka02; and 98, 99, and 100% for Inse01 (Fig. 5).

The relationships between the adjusted R² and the tested thresholds for the minimal number of reads exhibited different patterns of change according to the marker and the selected best identity threshold (Fig. 6). The adjusted R² varied significantly according to the tested thresholds for the minimal number of reads (Friedman rank sum tests; p < 0.001 for each primer set). For COI, adjusted R² was higher than 0.50 for thresholds equal to 5 and 10, and for thresholds higher than 50 reads, but was lower if the selected minimal number of reads was equal to 1 or comprised between 15 and 45 (Fig. 6A). For Euka02, adjusted R² increased
412 with increasing minimal number of reads until a plateau at about 65 reads (Fig. 6B). For both
413 markers, the relationships between adjusted R² and thresholds for the minimum number of
414 reads were similar for all the tested thresholds for best identity. The patterns of response of
415 Inse01 were similar for the 98 and 99% thresholds for best identity (Fig. 6C), but differed from
416 that obtained with a threshold of 100% of identity. The highest adjusted R² values were
417 observed for all the identity thresholds using thresholds comprised between 50 and 70 reads.
418 Based on these results, the threshold values 5, 60, 65, 70 and 100 for COI; 80, 85, 90, 95 and
419 100 for Euka02; and, 50, 55, 60, 65 and 70 for Inse01, were respectively selected (Fig. 6).

420

421 *Biomonitoring efficiency of the bulk-DNA approach*

422 After having identified, for each marker, the best mode of data expression (relative
423 frequencies) and the best ranges of identity threshold and minimal number of reads, we
424 analyzed the global bioassessment efficiency of each individual marker and each combination
425 of markers. We compared the values of the B-I₂M₂, its related metrics and the allocated
426 ecological quality classes with those provided by the standard approach. We also compared
427 the discrimination efficiency of the B-I₂M₂ calculated on each tested combination.

428

429 • *Biomonitoring efficiency of individual markers*

430 The marker Euka02 provided the best regression between M-I₂M₂ and B-I₂M₂ (Fig. 7A)
431 with a best adjusted R² of 0.671 for the combination of thresholds "96%:85" (for 96% of best
432 identity and a minimum of 85 reads, afterwards abbreviated min.r), among the 15 combinations
433 of selected thresholds for Euka02. The best adjusted R² were equal to 0.524 and 0.665 for
434 COI (80%:5 min.r) and Inse01 (100%:70 min.r), respectively and among the 15 combinations
435 of selected thresholds for each primer set. The I₂M₂ values differed significantly between
morphology-based versus bulk-DNA approaches for all three markers (Wilcoxon signed rank test for paired data; p-values < 0.05; Tab. 1). The ecological quality classes allocated to sites based on both approaches differed also in many cases: i.e. 5/36, 16/36 and 18/36 times for Inse01, Euka02 and COI respectively (Tab. 1). Both COI and Euka02 often allocated a worse ecological status to sites than did the morphology-based approach, while Inse01 tended to provide a good match between ecological classes assigned with both approaches (Tab. 1). The discrimination efficiency of the I₂M₂ index values provided by the bulk-DNA approach with COI (DE = 0.875; Tab. 1) was higher than that obtained with the I₂M₂ index values provided by the standard approach (DE<sub>standard</sub> = 0.833), while for the other markers the discrimination efficiency was poorer (DE = 0.75 for both markers).

The values of the ASPT and the relative frequency of polyvoltine organisms in the invertebrate assemblage (Polyvoltinism in Fig. 7) were correctly modelled from bulk-DNA data for all the markers. Adjusted R² ranged from 0.522 (Euka02) to 0.779 (Inse01) for ASPT, and between 0.457 (Inse01) and 0.768 (Euka02) for polyvoltinism. The relative frequency of ovoviviparous organisms (Ovoviviparity in Fig. 7) was correctly modelled by Euka02 and COI, with adjusted R² equal to 0.398 and 0.358, respectively, but not by Inse01. Richness was correctly modelled by COI (adj-R² = 0.520), but not by Euka02 (adj-R² < 0) and Inse01. Whatever the marker, the bulk-DNA approach did not correctly assess the Shannon-Weaver diversity, with adjusted R² ranging from -0.004 (COI) to 0.046 (Inse01).

455

456

• Biomonitoring efficiency of marker combinations

Combining bulk-DNA information provided by different markers improved the match between the standard and bulk-DNA based I₂M₂ values (Fig. 4B & 7B). The best adjusted R² between M-I₂M₂ and B-I₂M₂ increased from 0.671 (Euka02) to 0.717 when we combined all

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454
460 markers information (Fig. 4B & 7B). DE values of marker combinations ranged from 0.833 (for
461 [COI + Inse01]) to 0.958 (for [COI + Euka02] and [COI + Inse01 + Euka02]) (Tab. 1).

462 The increased quality of the model was readily explained by an increase in the adjusted
463 R² values for four out of the five individual metrics included in the I₂M₂ index, although some
464 of these increases remained modest (Fig. 7A & 7B). The adjusted R² of total richness and the
465 relative frequency of ovoviviparous organisms in the assemblage greatly increased from 0.520
466 (COI) to 0.634 [COI + Inse01], and from 0.398 (Euka02) to 0.610 [COI + Euka02], respectively.
467 In contrast, the adjusted R² of the Shannon-Weaver index and of the relative frequency of
468 polyvoltine organisms increased marginally from 0.046 (Inse01) to 0.064 [Inse01+ Euka02]
469 and from 0.768 (Euka02) to 0.770 [COI + Euka02], respectively. The adjusted R² of the ASPT
470 index decreased for combinations of primer sets from 0.779 (Inse01) to 0.710 [COI + Inse01].

471 Even if adjusted R² values increased, ecological quality classes allocated to tested
472 sites still differed for 9/36, 8/36, 5/36 and 9/36 sampling events respectively for [COI + Euka02],
473 [COI + Inse01], [Inse01 + Euka02] and for the combination of the three markers (Tab. 1).
474 Globally, the combinations of markers allocated a rather worse ecological status than the
475 standard approach (Tab. 1). Values differed significantly between M-I₂M₂ and B-I₂M₂ for [COI
476 + Euka02] and [COI + Inse01] (paired Wilcoxon's test; both p < 0.01; Tab. 1), but not for
477 [Inse01 + Euka02] and for the combinations of the three markers (both p > 0.05; Tab. 1).
478 **Discussion**

479

480 *Bulk-DNA and biomonitoring*

481 In this study, our objective was to evaluate the efficiency of bulk-DNA metabarcoding
482 of benthic macroinvertebrates to assess the ecological status of wadeable streams. We
483 therefore compared results based on morphological identification of organisms using a
484 standardized protocol (the I₂M₂; Mondy et al., 2012) to those provided by different
485 combinations of metabarcoding markers. We have confirmed the usefulness of bulk-DNA
486 metabarcoding for invertebrate-based stream bioassessment and have identified several
487 strategies to maximize the match between metabarcoding and standard approaches. In
488 addition, B-I₂M₂ and M-I₂M₂ values matched better for two out of the three markers when using
489 the relative number of reads (RF) for each MOTU, instead of their presence/absence only.

490 Aylagas, Borja, Muxika, and Rodríguez-Ezpeleta (2018) also observed better bioassessment
491 efficiency when using the number of reads instead of the presence/absence of MOTUs.

492 Nonetheless, results from Buchner et al. (2019) and Zizka, Geiger, and Leese (2020) have
493 suggested that macroinvertebrate presence/absence data can lead to similar bioassessement
494 results when compared to abundance-based data.

495 The good similarity observed between B-I₂M₂ and M-I₂M₂ could be explained by the
496 good correlations observed between the values of the individual metrics of the I₂M₂ obtained
497 with both bulk-DNA and standard approaches for four out of the five metrics. Several studies
498 confirmed that bulk data would readily allow to efficiently retrieve metric values initially
499 calculated on data achieved with a standard approach (Elbrecht et al., 2017; Emilson et al.,
500 2017; Gibson et al., 2015; Serrana et al., 2019).

501
503 Our results also confirmed that combining data from different markers allows improving
504 the appraisal of biodiversity based on bulk samples (Clarke et al., 2014; Hajibabaei, Spall,
505 Shokralla, & van Konynenburg, 2012). Overall, the marker COI (primer set BF1-BR2) exhibited
506 the best performance in terms of number of recovered taxa, but stream bioassessment based
507 on COI data alone was the least efficient. Combining metabarcoding data obtained with at
508 least two different markers has allowed us to increase both the number of taxa recovered with
509 metabarcoding and the quality of the bioassessment. Adding data obtained with a third marker
510 increased the number of recovered taxa, but not the quality of the bioassessment. We
511 therefore suggest that a minimum of two different markers should be used for the
512 biomonitoring of stream based on the bulk-DNA metabarcoding of benthic macroinvertebrates.

513 Moreover, our results have also highlighted that, depending on the type(s) of metrics
514 included in a biotic index, the choice of the markers should be quite logically also governed by
515 their ability to recover important indicator taxa (e.g. EPT) in order to minimize the risk of
516 missing such taxa. For instance, evaluating the ecological status of French streams with the
517 $I_2M_2$ index would need a combination of markers able to recover Gammaridae (and the
518 generally abundant genus *Gammarus*). This taxon is indeed highly contributing to the trait-
519 based metrics "ovoviviparity" and "polyvoltinism", and therefore its absence (when filtered out)
520 would explain why both the COI and Euka02 markers exhibited drops of their $R^2$ past certain
521 best identity thresholds (respectively at 85 and 98%). This observation also confirms the
522 importance of the threshold selection step (also addressed in the next section).

523
We have highlighted the importance of the selected thresholds for the best identity value and for the minimal (absolute or relative) number of reads when filtering the taxa which will be considered as present in the sample during the bioinformatics step of the metabarcoding approach. Slightly different thresholds can lead to far different results in terms of M-I₂M₂ and B-I₂M₂ correlations, at least when using individual markers, confirming the results observed by Tapolczai et al. (2019) for diatoms. Moreover, our results have shown that best threshold values were highly different among markers. Therefore, in order to maximize the robustness of metabarcoding-based biomonitoring indices, we suggest the preliminary selection of marker-specific thresholds, based – for instance – on correlation tests between metabarcoding data and abiotic information (e.g. metabarcoding-based index values calculated along a pressure gradient, Tapolczai et al., 2019) or biotic information (e.g. between metabarcoding- and morphology-based index values, this study).

537

### *Ecological quality class allocation*

The B-I₂M₂ values efficiently discriminate between least impaired and significantly impaired sites, but misclassifications (i.e. differences in the ecological quality class allocated by the metabarcoding-based and the morphology-based approaches) were observed. Comparing both methods, Dowle, Pochon, Banks, Shearer, and Wood (2016) indicated that such changes in quality class allocation were often due to changes in index values that were already close to an ecological quality class boundary. In estuarine and coastal sediment, Aylagas et al. (2018) observed changes in the quality class allocation for 14 of 18 stations, based on benthic macroinvertebrate communities. They also identified a lower rate of change for the stations allocated to the ecological quality classes exhibiting the widest range of values. It should be noted that in our study, the majority of the sampling events (22/36) were performed in sites of "high" ecological quality (based on their M-I₂M₂ value), which is also the quality class
550 exhibiting the largest extent. This "high" ecological status of many stream sites within the database may also have limited the number of misclassifications compared to other studies.
551
552 It may as well explain why worse evaluations of the ecological status were mainly observed
553 when classifying site sampling events based on the metabarcoding approach.

554 In addition, misclassifications were mainly due to discrepancies between the values of
555 the $I_2M_2$ metrics calculated with both approaches. Such discrepancies were due to: (i)
556 differences in the provided taxonomic lists; (ii) differences in the estimated relative
557 abundances of taxa (i.e. relative numbers of individuals vs. relative numbers of reads); and,
558 (iii) inadequacy of the "reference" values needed to calculate the EQR (calculations of
559 Ecological Quality Ratios, as advocated by the WFD, need to normalize metrics values in a 0-
560 1 range, using "reference" values). Potentially promising prospects for improvement will be
561 discussed in the following sections.

562

563 *On individual metrics*

564 Combining data obtained with at least two different markers has allowed to correctly
565 model the values of four out of the five metrics of the $I_2M_2$. Namely, we observed good results
566 for the ASPT, the total taxonomic richness, and the relative frequencies of the polyvoltine and
567 ovoviviparous organisms within the invertebrate assemblage, but not for the Shannon-Weaver
568 diversity index. Gibson et al. (2015) also observed a lack of significant, positive relationships
569 between data obtained with standard and bulk-DNA approaches, for Pielou's evenness and
570 Simpson's dominance indices. In contrast with the Shannon-Weaver index, the ASPT
571 calculation only needs information on the presence/absence of taxa identified at the family
572 level (Armitage et al., 1983). It probably explains the better robustness of this metric, and why
573 ASPT still exhibited good correlation between values calculated applying both approaches,
574 even if calculated on only 2/3 of the available data on the sampled invertebrate assemblage
575 (both the Shannon-Weaver index and the ASPT are calculated on eight sample units over 576 twelve, corresponding to only two of the three phases; see Fig. 2B and Mondy et al., 2012, for 577 further details).

578 As in previous studies (Hajibabaei et al., 2012; Serrana et al., 2019), morphotaxa not 579 recovered with metabarcoding were often “rare” taxa, i.e. taxa with less than 10 individuals 580 over the whole sampling design (observed for 29 out of the 41 missing taxa). Among the five 581 individual metrics combined in the I₂M₂, the Shannon-Weaver index and the relative 582 frequencies of polyvoltine and ovoviviparous invertebrates within the assemblage are the only 583 ones taking into account the taxon abundances. Missing the rarest taxa seems to have only a 584 limited impact on the calculated values of the trait-based metrics, but probably not on those of 585 the Shannon-Weaver diversity index.

586 Metabarcoding also identified taxa absent in the morphotaxa list. Their overall 587 frequency was moderate (20% of the global faunal list, including both identified morphotaxa 588 and DNA-based taxa) when compared, for instance, to the results provided by Serrana et al. 589 (2019), who had observed nine taxa identified only with metabarcoding for a total of 20 taxa. 590 Such taxa could especially lead to bias in values of metrics closely depending on taxonomic 591 richness; i.e. total richness and Shannon-Weaver diversity in the I₂M₂. Nevertheless, in our 592 study, a non-negligible part of these taxa corresponded to taxa identified at heterogeneous 593 taxonomic levels among samples; i.e. at the family level in some samples (e.g. Ephemerellidae) 594 and the genus level (e.g. *Ephemerella*) in other samples from the same site. This bias could 595 be partially addressed by standardizing the taxonomic list at the site scale, instead of at the 596 sample scale as it is currently done in the data preparation step of the I₂M₂ calculation.

597

598 *On reference values and discrimination efficiency*
For best managing the potential bias in taxonomic identification and abundance estimation when using the bulk-DNA approach (Hering et al., 2018), stream-type reference values, "BEST" and "WORST", could have been specifically redefined for the calculation of the B-I₂M₂ index. During preliminary analyses we tested such an approach, but preliminary results showed that it was not actually pertinent: performances were similar, but we actually observed more misclassifications in ecological quality class allocation when using metabarcoding-specific "BEST" and "WORST" reference values. "BEST" and "WORST" values used routinely for the calculation of the M-I₂M₂ were indeed defined per stream type (for "BEST" values) and on the available national dataset (values currently in use were defined on a dataset containing more than 10,000 site sampling events), whereas "BEST" and "WORST" values could only have been defined on our modest bulk dataset. Therefore "BEST" and "WORST" values defined for metabarcoding data would have been less robust than those used for the M-I₂M₂, and would have represented an additional source of discrepancy between B-I₂M₂ and M-I₂M₂ approaches, both in terms of correlation and ecological status assessment.

Similarly, the I₂M₂ index is calculated as the arithmetic average of 17 sub-indices (one per pressure category), each corresponding to the mean of the individual metrics weighed by their respective DE for the corresponding category of pressure (Mondy et al., 2012). The discrimination efficiency of individual metrics has been calculated at the national scale. Here, the sampling design was too small to have the possibility of calculating robust values of DE specifically allocated to the calculation of bulk-based I₂M₂ values. Defining such DE values in the future should also improve the relationship between M-I₂M₂ and B-I₂M₂.

## Perspectives and conclusion

This study supports the bulk metabarcoding approach as a promising method for stream biomonitoring based on bulk-DNA from benthic macroinvertebrate samples (Aylagas,
Borja, Rodríguez-Ezpeleta, & Consuegra, 2014; Beentjes et al., 2019; Carew et al. 2018; Carew, Pettigrove, Metzeling, & Hoffmann, 2013; Elbrecht et al., 2017; Emilson et al., 2017; Gibson et al., 2014, 2015; Hajibabaei et al., 2011, 2019; Kuntke et al., 2020; Serrana et al., 2019; Zizka et al., 2020). However, more work is needed before implementing bulk-metabarcoding in the routine monitoring of streams. Indeed, we have highlighted a series of biases leading to the reclassification of several site sampling events in terms of ecological quality class. Solutions exist for reducing these biases, and we have focused on several of them in the previous sections. For instance, a reference DNA barcoding database, including 578 different taxa (62% identified at the species level) has been specifically built for this study (Ficetola et al., 2020), in order to work with a database as complete as possible, as recommended by several authors (Aylagas et al., 2014; Hering et al., 2018). Several studies (Gibson et al., 2015; Ji et al., 2013; Mächler et al., 2019; Serrana et al., 2019; Sweeney et al., 2011) have also suggested that new indices could be constructed specifically on bulk-DNA information for stream bioassessment based on macroinvertebrate assemblages, as already been done for benthic diatoms (e.g. Vasselon, Domaizon, Rimet, Kahlert, & Bouchez, 2017). On one hand, even if retrieved MOTUs are not assigned to taxa, MOTUs can still be used to efficiently discriminate between impaired and reference situations (e.g. Emilson et al., 2017). However, this approach would need a huge preliminary sampling phase in order to construct a new index based on a database that would include the majority of the MOTUs that could be recovered, for instance, for any stream found at the national scale. On the other hand, when molecular methods will prove to be mature enough, they could be used to obtain species-level taxonomic lists. Based on these lists, the observed and reference values of taxonomy-based and trait-based metrics could be refined to further improve the discrimination efficiency of DNA-based indices. Other promising sources of DNA for stream bioassessment are the ethanol used for sample/invertebrate preservation (e.g. Hajibabaei et al., 2012; Martins et al., 2020; Zizka, Leese, Peinert, & Geiger, 2019) or even DNA directly extracted from unsorted samples (Pereira-da-Conceicao et al., 2019). Both sources would allow to avoid the time-
651 consuming steps of invertebrate sorting and counting. Supervised machine learning is also
652 another promising approach, as it could be used to directly model a standard biotic index value
653 based on metabarcoding data (Cordier, Lanzén, Apothéloz-Perret-Gentil, Stoeck, &
654 Pawlowski, 2019; Frühe et al., 2020).

655 Last, one major challenge for an optimal match between standard and bulk-DNA
656 results in a stream bioassessment context is the optimization of taxon abundance estimations.
657 Taxon abundances are taken into account in many invertebrate-based bioassessment
658 methods (e.g. in all the intercalibrated European methods; Bennet et al., 2011). In this study,
659 we did not investigate how to better correlate the relative abundances of morphotaxa in
660 samples with information provided by the relative numbers of sequence reads. However,
661 several studies have shown that the numbers of reads could be correlated with taxon
662 abundances or biomasses, albeit frequently with a poor fit (Carew et al., 2013; Deagle,
663 Thomas, Shaffer, Trites, & Jarman, 2013; Dowle et al., 2016; Elbrecht & Leese, 2015, 2017;
664 Serrana et al., 2019), strengthening the idea that better estimating the relative abundances of
665 taxa based on their relative numbers of reads in a bulk sample could improve the modelled
666 values of a biotic index.

667
668 Acknowledgments

669 This study was supported by a grant from Agence Nationale de la Recherche (aquaDNA;

670 ANR-13-ECOT-0002-01).

671 We thank SPYGEN staff for its help for DNA extraction.

672 We thank Vasco Elbrecht for his insightful comments. Thanks also go to Jan Pawlowski,

673 Christian Moritz and three anonymous reviewers for improving a previous version of the paper.

674
675 References

676

677 AFNOR. (2009). Qualité de l'eau – Prélèvement des macro-invertébrés aquatiques en rivières peu profondes. XP T90-333.
678

679

680 AFNOR. (2010). Qualité de l'eau – Traitement au laboratoire d'échantillons contenant des macro-invertébrés de cours d'eau. XP T90-388.
681

682

683 Armitage, P., Moss, D., Wright, J., & Furse, M. (1983). The performance of a new biological water quality score system based on macroinvertebrates over a wide range of unpolluted running-water sites. Water Research, 17(3), 333-347. https://doi.org/10.1016/0043-1354(83)90188-4
686

687

688 Aylagas, E., Borja, Á., Muxika, I., & Rodríguez-Ezpeleta, N. (2018). Adapting metabarcoding-based benthic biomonitoring into routine marine ecological status assessment networks. *Ecological Indicators*, *95*, 194-202. https://doi.org/10.1016/j.ecolind.2018.07.044

691

692 Aylagas, E., Borja, Á., & Rodríguez-Ezpeleta, N. (2014). Environmental Status Assessment Using DNA Metabarcoding: Towards a Genetics Based Marine Biotic Index (gAMBI). *PLoS ONE*, *9*(3), e90529. https://doi.org/10.1371/journal.pone.0090529

695

696 Baird, D. J., & Hajibabaei, M. (2012). Biomonitoring 2.0: A new paradigm in ecosystem assessment made possible by next-generation DNA sequencing. *Molecular Ecology*, *21*(8), 2039-2044. https://doi.org/10.1111/j.1365-294X.2012.05519.x

699

700 Baird, D. J., Pascoe, T. J., Zhou, X., & Hajibabaei, M. (2011). Building freshwater macroinvertebrate DNA-barcode libraries from reference collection material: Formalin preservation vs specimen age. *Journal of the North American Benthological Society*, *30*(1), 125-130. https://doi.org/10.1899/10-013.1

704

705 Beentjes, K. K., Speksnijder, A. G. C. L., Schilthuizen, M., Hoogeveen, M., Pastoor, R., & van der Hoorn, B. B. (2019). Increased performance of DNA metabarcoding of macroinvertebrates by taxonomic sorting. *PLOS ONE*, *14*(12), e0226527. https://doi.org/10.1371/journal.pone.0226527

709

710 Bennett, C., Owen, R., Birk, S., Buffagni, A., Erba, S., Mengin, N., ... Wasson, J.-G. (2011). Bringing European river quality into line: An exercise to intercalibrate macro-invertebrate classification methods. *Hydrobiologia*, *667*(1), 31-48. https://doi.org/10.1007/s10750-011-0635-2

714

715 Birk, S., Bonne, W., Borja, A., Brucet, S., Courrat, A., Poikane, S., ... Hering, D. (2012). Three hundred ways to assess Europe's surface waters: An almost complete overview of biological methods to implement the Water Framework Directive. *Ecological Indicators*, *18*, 31-41. https://doi.org/10.1016/j.ecolind.2011.10.009

719

720 Bokulich, N. A., Subramanian, S., Faith, J. J., Gevers, D., Gordon, J. I., Knight, R., ... Caporaso, J. G. (2013). Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. *Nature Methods*, *10*(1), 57-59. https://doi.org/10.1038/nmeth.2276

723
724 Boyer, F., Mercier, C., Bonin, A., Le Bras, Y., Taberlet, P., & Coissac, E. (2016). obitools : A unix-inspired software package for DNA metabarcoding. *Molecular Ecology Resources*, **16**(1), 176-182. https://doi.org/10.1111/1755-0998.12428

727 Buchner, D., Beermann, A. J., Laini, A., Rolauffs, P., Vitecek, S., Hering, D., & Leese, F. (2019). Analysis of 13,312 benthic invertebrate samples from German streams reveals minor deviations in ecological status class between abundance and presence/absence data. *PLOS ONE*, **14**(12), e0226547. https://doi.org/10.1371/journal.pone.0226547

732 Carew, M. E., Kellar, C. R., Pettigrove, V. J., & Hoffmann, A. A. (2018). Can high-throughput sequencing detect macroinvertebrate diversity for routine monitoring of an urban river? *Ecological Indicators*, **85**, 440-450. https://doi.org/10.1016/j.ecolind.2017.11.002

736 Carew, M. E., Pettigrove, V. J., Metzeling, L., & Hoffmann, A. A. (2013). Environmental monitoring using next generation sequencing : Rapid identification of macroinvertebrate bioindicator species. *Frontiers in Zoology*, **10**(1), 45. https://doi.org/10.1186/1742-9994-10-45

740 Civade, R., Dejean, T., Valentini, A., Roset, N., Raymond, J.-C., Bonin, A., ... Pont, D. (2016). Spatial Representativeness of Environmental DNA Metabarcoding Signal for Fish Biodiversity Assessment in a Natural Freshwater System. *PLOS ONE*, **11**(6), e0157366. https://doi.org/10.1371/journal.pone.0157366

745 Clarke, L. J., Soubrier, J., Weyrich, L. S., & Cooper, A. (2014). Environmental metabarcodes for insects : *In silico* PCR reveals potential for taxonomic bias. *Molecular Ecology Resources*, **14**(6), 1160-1170. https://doi.org/10.1111/1755-0998.12265

749 Crawley, M. J. (2007). *The R Book*. John Wiley & Sons, Ltd. https://doi.org/10.1002/9780470515075

752 Cordier, T., Lanzén, A., Apothéloz-Perret-Gentil, L., Stoeck, T., & Pawlowski, J. (2019). *Embracing Environmental Genomics and Machine Learning for Routine Biomonitoring. Trends in Microbiology*, **27**(5), 387-397. https://doi.org/10.1016/j.tim.2018.10.012

756 Deagle, B. E., Thomas, A. C., Shaffer, A. K., Trites, A. W., & Jarman, S. N. (2013). Quantifying sequence proportions in a DNA-based diet study using Ion Torrent amplicon sequencing : Which counts count? *Molecular Ecology Resources*, **13**(4), 620-633. https://doi.org/10.1111/1755-0998.12103

761 Dejean, T., Valentini, A., Miquel, C., Taberlet, P., Bellemain, E., & Miaud, C. (2012). Improved detection of an alien invasive species through environmental DNA barcoding : The example of the American bullfrog *Lithobates catesbeianus: Alien invasive species detection using eDNA*. *Journal of Applied Ecology*, **49**(4), 953-959. https://doi.org/10.1111/j.1365-2664.2012.02171.x

767 Dolédec, S., & Statzner, B. (2008). Invertebrate traits for the biomonitoring of large European rivers : An assessment of specific types of human impact. *Freshwater Biology*, **53**(3), 617-634. https://doi.org/10.1111/j.1365-2427.2007.01924.x

771 Dowle, E. J., Pochon, X., Banks, J., Shearer, K., & Wood, S. A. (2016). Targeted gene enrichment and high-throughput sequencing for environmental biomonitoring : A case study using freshwater macroinvertebrates. *Molecular Ecology Resources*, **16**(5), 1240-1254. https://doi.org/10.1111/1755-0998.12488

776
777 Dray, S., & Dufour, A.-B. (2007). The ade4 Package: Implementing the Duality Diagram for Ecologists. *Journal of Statistical Software*, 22(4). https://doi.org/10.18637/jss.v022.i04

779

780 Dudgeon, D., Arthington, A. H., Gessner, M. O., Kawabata, Z.-I., Knowler, D. J., Lévêque, C., ... Sullivan, C. A. (2006). Freshwater biodiversity: Importance, threats, status and conservation challenges. *Biological Reviews*, 81(02), 163.
783 https://doi.org/10.1017/S1464793105006950

784

785 Elbrecht, V., & Leese, F. (2015). Can DNA-Based Ecosystem Assessments Quantify Species Abundance? Testing Primer Bias and Biomass—Sequence Relationships with an Innovative Metabarcoding Protocol. *PLOS ONE*, 10(7), e0130324.
788 https://doi.org/10.1371/journal.pone.0130324

789

790 Elbrecht, V., & Leese, F. (2017). *Validation and development of COI metabarcoding primers for freshwater macroinvertebrate bioassessment* [Preprint]. PeerJ Preprints. https://doi.org/10.7287/peerj.preprints.2044v5

793

794 Elbrecht, V., Taberlet, P., Dejean, T., Valentini, A., Usseglio-Polatera, P., Beisel, J.-N., Coissac, E., Boyer, F., & Leese, F. (2016). Testing the potential of a ribosomal 16S marker for DNA metabarcoding of insects. *PeerJ*, 4, e1966. https://doi.org/10.7717/peerj.1966

797

798 Elbrecht, V., Vamos, E. E., Meissner, K., Aroviita, J., & Leese, F. (2017). Assessing strengths and weaknesses of DNA metabarcoding- based macroinvertebrate identification for routine stream monitoring. *Methods in Ecology and Evolution*, 8(10), 1265-1275.
801 https://doi.org/10.1111/2041-210X.12789

802

803 Emilson, C. E., Thompson, D. G., Venier, L. A., Porter, T. M., Swystun, T., Chartrand, D., ... Hajibabaei, M. (2017). DNA metabarcoding and morphological macroinvertebrate metrics reveal the same changes in boreal watersheds across an environmental gradient. *Scientific Reports*, 7(1), 12777. https://doi.org/10.1038/s41598-017-13157-x

807

808 Ficetola, G. F., Boyer, F., Valentini, A., Bonin, A., Meyer, A., Dejean, T., Gaboriaud, C., Usseglio- Polatera, P., & Taberlet, P. (2020). Comparison of markers for the monitoring of freshwater benthic biodiversity through DNA metabarcoding. *Molecular Ecology*. https://doi.org/10.1111/mec.15632

812

813 Frühe, L., Cordier, T., Dully, V., Breiner, H., Lentendu, G., Pawlowski, J., Martins, C., Wilding, T. A., & Stoeck, T. (2020). Supervised machine learning is superior to indicator value inference in monitoring the environmental impacts of salmon aquaculture using eDNA metabarcodes. *Molecular Ecology*, mec.15434. https://doi.org/10.1111/mec.15434

817

818 Gibson, J. F., Shokralla, S., Curry, C., Baird, D. J., Monk, W. A., King, I., & Hajibabaei, M. (2015). Large-Scale Biomonitoring of Remote and Threatened Ecosystems via High-Throughput Sequencing. *PLOS ONE*, 10(10), e0138432.
821 https://doi.org/10.1371/journal.pone.0138432

822

823 Gibson, J., Shokralla, S., Porter, T. M., King, I., van Konynenburg, S., Janzen, D. H., ... Hajibabaei, M. (2014). Simultaneous assessment of the macrobiome and microbiome in a bulk sample of tropical arthropods through DNA metasystematics. *Proceedings of the National Academy of Sciences*, 111(22), 8007-8012. https://doi.org/10.1073/pnas.1406468111

827
828 Guardiola, M., Uriz, M. J., Taberlet, P., Coissac, E., Wangensteen, O. S., & Turon, X. (2015).

829 Deep-Sea, Deep-Sequencing: Metabarcoding Extracellular DNA from Sediments of Marine

830 Canyons. PLOS ONE, 10(10), e0139633. https://doi.org/10.1371/journal.pone.0139633

831

832 Hajibabaei, M., Porter, T. M., Robinson, C. V., Baird, D. J., Shokralla, S., & Wright, M. T. G.

833 (2019). Watered-down biodiversity? A comparison of metabarcoding results from DNA

834 extracted from matched water and bulk tissue biomonitoring samples. PLOS ONE, 14(12),

835 e0225409. https://doi.org/10.1371/journal.pone.0225409

836

837 Hajibabaei, M., Shokralla, S., Zhou, X., Singer, G. A. C., & Baird, D. J. (2011). Environmental

838 Barcoding: A Next-Generation Sequencing Approach for Biomonitoring Applications Using

839 River Benthos. PLoS ONE, 6(4), e17497. https://doi.org/10.1371/journal.pone.0017497

840

841 Hajibabaei, M., Spall, J. L., Shokralla, S., & van Konynenburg, S. (2012). Assessing

842 biodiversity of a freshwater benthic macroinvertebrate community through non-destructive

843 environmental barcoding of DNA from preservative ethanol. BMC Ecology, 12(1), 28.

844 https://doi.org/10.1186/1472-6785-12-28

845

846 Hering, D., Borja, A., Jones, J. I., Pont, D., Boets, P., Bouchez, A., ... Kelly, M. (2018).

847 Implementation options for DNA-based identification into ecological status assessment under

848 the European Water Framework Directive. Water Research, 138, 192-205.

849 https://doi.org/10.1016/j.watres.2018.03.003

850

851 Hutchings, P. (2019). An advocate for taxonomic research in Australia. Pacific Conservation

852 Biology, 25(1), 34. https://doi.org/10.1071/PC17033

853

854 Ji, Y., Ashton, L., Pedley, S. M., Edwards, D. P., Tang, Y., Nakamura, A., ... Yu, D. W. (2013).

855 Reliable, verifiable and efficient monitoring of biodiversity via metabarcoding. Ecology Letters,

856 16(10), 1245-1257. https://doi.org/10.1111/ele.12162

857

858 Kuntke, F., de Jonge, N., Hesselsøe, M., & Lund Nielsen, J. (2020). Stream water quality

859 assessment by metabarcoding of invertebrates. Ecological Indicators, 111, 105982.

860 https://doi.org/10.1016/j.ecolind.2019.105982

861

862 Larras, F., Coulaud, R., Gautreau, E., Billoir, E., Rosebery, J., & Usseglio-Polatera, P. (2017).

863 Assessing anthropogenic pressures on streams: A random forest approach based on benthic

864 diatom communities. Science of The Total Environment, 586, 1101-1112.

865 https://doi.org/10.1016/j.scitotenv.2017.02.096

866

867 Lobo, J., Shokralla, S., Costa, M. H., Hajibabaei, M., & Costa, F. O. (2017). DNA

868 metabarcoding for high-throughput monitoring of estuarine macrobenthic communities.

869 Scientific Reports, 7(1), 15618. https://doi.org/10.1038/s41598-017-15823-6

870

871 Mächler, E., Little, C. J., Wüthrich, R., Alther, R., Fronhofer, E. A., Gounand, I., ... Altermatt,

872 F. (2019). Assessing different components of diversity across a river network using eDNA.

873 Environmental DNA, 1(3), 290-301. https://doi.org/10.1002/edn3.33

874

875 Martins, F. M. S., Porto, M., Feio, M. J., Eger, B., Bonin, A., Serra, S. R. Q., Taberlet, P., &

876 Beja, P. (2020). Modelling technical and biological biases in macroinvertebrate community

877 assessment from bulk preservative using multiple metabarcoding markers. Molecular Ecology.

878 https://doi.org/10.1111/mec.15620

879
880 Metzeling, L., Chessman, B., Hardwick, R., & Wong, V. (2003). Rapid assessment of rivers using macroinvertebrates: The role of experience, and comparisons with quantitative methods.
*Hydrobiologia*, 510(1-3), 39-52. https://doi.org/10.1023/B:HYDR.0000008500.34301.a0

881
882
883
884 Mondy, C. P., Villeneuve, B., Archaimbault, V., & Usseglio-Polatera, P. (2012). A new macroinvertebrate-based multimetric index (I2M2) to evaluate ecological quality of French wadeable streams fulfilling the WFD demands: A taxonomical and trait approach. *Ecological Indicators*, 18, 452-467. https://doi.org/10.1016/j.ecolind.2011.12.013

888
889 Ofenböck, T., Moog, O., Gerritsen, J., & Barbour, M. (2004). A Stressor Specific Multimetric Approach for Monitoring Running Waters in Austria Using Benthic Macro-Invertebrates. In D. Hering, P. F. M. Verdonschot, O. Moog, & L. Sandin (Éds.), *Integrated Assessment of Running Waters in Europe* (p. 251-268). Springer Netherlands. https://doi.org/10.1007/978-94-007-0993-5 15

894
895 Pereira-da-Conceicoa, L., Elbrecht, V., Hall, A., Briscoe, A., Barber-James, H., & Price, B. (2019). *Metabarcoding unsorted kick-samples facilitates macroinvertebrate-based biomonitoring with increased taxonomic resolution, while outperforming environmental DNA* [Preprint]. https://doi.org/10.1101/792333

899
900 Piñol, J., Mir, G., Gomez-Polo, P., & Agustí, N. (2015). Universal and blocking primer mismatches limit the use of high-throughput DNA sequencing for the quantitative metabarcoding of arthropods. *Molecular Ecology Resources*, 15(4), 819-830. https://doi.org/10.1111/1755-0998.12355

904
905 R Core Team. (2020). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

907
908 Rimet, F., Chaumeil, P., Keck, F., Kermarrec, L., Vasselon, V., Kahlert, M., ... Bouchez, A. (2016). R-Syst::diatom: An open-access and curated barcode database for diatoms and freshwater monitoring. *Database*, 2016, baw016. https://doi.org/10.1093/database/baw016

911
912 Serrana, J. M., Miyake, Y., Gamboa, M., & Watanabe, K. (2019). Comparison of DNA metabarcoding and morphological identification for stream macroinvertebrate biodiversity assessment and monitoring. *Ecological Indicators*, 101, 963-972. https://doi.org/10.1016/j.ecolind.2019.02.008

916
917 Serrana, J. M., Yaegashi, S., Kondoh, S., Li, B., Robinson, C. T., & Watanabe, K. (2018). Ecological influence of sediment bypass tunnels on macroinvertebrates in dam-fragmented rivers by DNA metabarcoding. *Scientific Reports*, 8(1), 10185. https://doi.org/10.1038/s41598-018-28624-2

921
922 Shannon, C. E., & Weaver, W. (1963). The mathematical theory of communication. Univ. Illinois Press, Urbana, USA.

924
925 Shaw, J. L. A., Weyrich, L., & Cooper, A. (2017). Using environmental (e)DNA sequencing for aquatic biodiversity surveys: A beginner's guide. *Marine and Freshwater Research*, 68(1), 20. https://doi.org/10.1071/MF15361

928
929 Statzner, B., Bady, P., Dolédec, S., & Schöll, F. (2005). Invertebrate traits for the biomonitoring of large European rivers: An initial assessment of trait patterns in least impacted river reaches. *Freshwater Biology*, 50(12), 2136-2161. https://doi.org/10.1111/j.1365-2427.2005.01447.x

932
933 Strahler, A. N. (1957). Quantitative analysis of watershed geomorphology. *Transactions, 934 American Geophysical Union, 38*(6), 913. https://doi.org/10.1029/TR038i006p00913

935

936 Sweeney, B. W., Battle, J. M., Jackson, J. K., & Dapkey, T. (2011). Can DNA barcodes of 937 stream macroinvertebrates improve descriptions of community structure and water quality? 938 *Journal of the North American Benthological Society*, 30(1), 195-216. 939 https://doi.org/10.1899/10-016.1

940

941 Taberlet, P., Bonin, A., Zinger, L., & Coissac, E. (2018). *Environmental DNA: For biodiversity 942 research and monitoring* (First edition). Oxford University Press.

943

944 Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C., & Willerslev, E. (2012). Towards 945 next-generation biodiversity assessment using DNA metabarcoding. *Molecular Ecology*, 946 21(8), 2045-2050. https://doi.org/10.1111/j.1365-294X.2012.05470.x

947

948 Tapolczai, K., Vasselon, V., Bouchez, A., Stenger-Kovács, C., Padisák, J., & Rimet, F. (2019). 949 The impact of OTU sequence similarity threshold on diatom-based bioassessment: A case 950 study of the rivers of Mayotte (France, Indian Ocean). *Ecology and Evolution*, 9(1), 166-179. 951 https://doi.org/10.1002/ece3.4701

952

953 Terlizzi, A., Bevilacqua, S., Fraschetti, S., & Boero, F. (2003). Taxonomic sufficiency and the 954 increasing insufficiency of taxonomic expertise. *Marine Pollution Bulletin*, 46(5), 556-561. 955 https://doi.org/10.1016/S0025-326X(03)00066-3

956

957 Valentini, A., Pompanon, F., & Taberlet, P. (2009). DNA barcoding for ecologists. *Trends in 958 Ecology & Evolution*, 24(2), 110-117. https://doi.org/10.1016/j.tree.2008.09.011

959

960 Valentini, A., Taberlet, P., Miaud, C., Civade, R., Herder, J., Thomsen, P. F., ... Dejean, T. 961 (2016). Next-generation monitoring of aquatic biodiversity using environmental DNA 962 metabarcoding. *Molecular Ecology*, 25(4), 929-942. https://doi.org/10.1111/mec.13428

963

964 Vasselon, V., Domaizon, I., Rimet, F., Kahler, M., & Bouchez, A. (2017). Application of high- 965 throughput sequencing (HTS) metabarcoding to diatom biomonitoring: Do DNA extraction 966 methods matter? *Freshwater Science*, 36(1), 162-177. https://doi.org/10.1086/690649

967

968 Vörösmarty, C. J., McIntyre, P. B., Gessner, M. O., Dudgeon, D., Prusevich, A., Green, P., ... 969 Davies, P. M. (2010). Global threats to human water security and river biodiversity. *Nature, 970 467*(7315), 555-561. https://doi.org/10.1038/nature09440

971

972 Wasson, J.-G., Chandesris, A., Pella, H., & Blanc, L. (2002). *Les hydro-écorégions de France 973 métropolitaine — Approche régionale de la typologie des eaux courantes et éléments pour la 974 définition des peuplements de référence d'invertébrés*.

975

976 Wasson, J.-G., Chandesris, A., Pella, H., Sauquet, E., & Mengin, N. (2006). *Appui scientifique 977 à la mise en œuvre de la Directive Européenne Cadre sur l'Éau. Typologie des cours d'eau 978 de France métropolitaine*.

979

980 Xie, Y., Wang, J., Yang, J., Giesy, J. P., Yu, H., & Zhang, X. (2017). Environmental DNA 981 metabarcoding reveals primary chemical contaminants in freshwater sediments from different 982 land-use types. *Chemosphere*, 172, 201-209.

983 https://doi.org/10.1016/j.chemosphere.2016.12.117
984
985 Yu, D. W., Ji, Y., Emerson, B. C., Wang, X., Ye, C., Yang, C., & Ding, Z. (2012). Biodiversity

986 soup: Metabarcoding of arthropods for rapid biodiversity assessment and biomonitoring:

987 Biodiversity soup. *Methods in Ecology and Evolution*, 3(4), 613-623.

988 https://doi.org/10.1111/j.2041-210X.2012.00198.x

989

990 Zizka, V. M. A., Geiger, M. F., & Leese, F. (2020). DNA metabarcoding of stream invertebrates

991 reveals spatio-temporal variation but consistent status class assessments in a natural and

992 urban river. *Ecological Indicators*, 115, 106383. https://doi.org/10.1016/j.ecolind.2020.106383

993

994 Zizka, V. M. A., Leese, F., Peinert, B., & Geiger, M. F. (2019). DNA metabarcoding from

995 sample fixative as a quick and voucher-preserving biodiversity assessment method. *Genome*,

996 62(3), 122-136. https://doi.org/10.1139/gen-2018-0048

997
998 **Data Accessibility**

999 All datasets will be available at figshare upon manuscript acceptance:
1000 https://doi.org/10.6084/m9 figshare.13110899 (bulk-DNA data) and
1001 https://doi.org/10.6084/m9 figshare.13110692 (faunal data)

1002 **Author Contributions**

1003 Designed the research: Pierre Taberlet, Tony Dejean, Gentile Francesco Ficetola, Frédéric Boyer, Alice Valentini, Aurélie Bonin, Albin Meyer, Philippe Usseglio-Polatera, Jean-Nicolas Beisel, Florian Leese.

1006 Performed research: Alice Valentini, Aurélie Bonin, Albin Meyer, Coline Gaboriaud, Pierre Taberlet, Tony Dejean, Gentile Francesco Ficetola, Frédéric Boyer, Philippe Usseglio-Polatera, Philippe Wagner, Jonathan Bouquerel, Jean-Nicolas Beisel.

1009 Analyzed data: Albin Meyer, Philippe Usseglio-Polatera, Gentile Francesco Ficetola, Frédéric Boyer, Aurélie Bonin.

1011 Wrote the paper: Albin Meyer wrote the first draft of the paper, with subsequent contribution of all the authors.

1013
1014 Tables

1015

1016 Tab.1 Total abundance and number of reads for all the taxa identified with the morphology-

1017 and/or bulk-based approaches. Bulk data were filtered with a homogeneous best identity

1018 threshold of 0.80.

<table><thead><tr><th>Higher rank</th><th>Family</th><th>Taxon</th><th>Abundance</th><th>Reads (COI)</th><th>Reads (EUKA02)</th><th>Reads (INSE01)</th></tr></thead><tbody><tr><td>Oligochaeta</td><td></td><td></td><td>6971</td><td>87902</td><td>4023314</td><td>2091067</td></tr><tr><td rowspan="3">Hirudinea</td><td>Erpobdellidae</td><td></td><td>158</td><td>4169</td><td>313919</td><td>1567</td></tr><tr><td>Glossiphoniidae</td><td></td><td>92</td><td>4274</td><td>27341</td><td>0</td></tr><tr><td>Piscicolidae</td><td></td><td>6</td><td>181</td><td>4380</td><td>0</td></tr><tr><td>Hydracarina</td><td></td><td></td><td>86</td><td>32</td><td>102</td><td>0</td></tr><tr><td rowspan="6">Crustacea</td><td>Gammaridae</td><td>Echinogammarus</td><td>0</td><td>0</td><td>985</td><td>0</td></tr><tr><td>Gammaridae</td><td>Gammarus</td><td>20318</td><td>209364</td><td>943643</td><td>0</td></tr><tr><td>Astacidae</td><td>Pacifastacus</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Cambaridae</td><td>Procambarus</td><td>0</td><td>65</td><td>0</td><td>0</td></tr><tr><td>Asellidae</td><td></td><td>70</td><td>943</td><td>0</td><td>1975</td></tr><tr><td>Chrysomelidae</td><td></td><td>0</td><td>0</td><td>0</td><td>1039</td></tr><tr><td rowspan="14">Coleoptera</td><td>Curculionidae</td><td></td><td>27</td><td>0</td><td>44087</td><td>15211</td></tr><tr><td>Dytiscidae</td><td>Colymbetinae</td><td>14</td><td>0</td><td>100</td><td>0</td></tr><tr><td>Dytiscidae</td><td></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Dytiscidae</td><td>Hydroporinae</td><td>127</td><td>5041</td><td>79</td><td>0</td></tr><tr><td>Dytiscidae</td><td>Laccophilinae</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Elmidae</td><td></td><td>0</td><td>0</td><td>516</td><td>0</td></tr><tr><td>Elmidae</td><td>Elmis</td><td>1171</td><td>5488</td><td>25160</td><td>227850</td></tr><tr><td>Elmidae</td><td>Esolus</td><td>963</td><td>2911</td><td>23131</td><td>68410</td></tr><tr><td>Elmidae</td><td>Limnius</td><td>1309</td><td>23073</td><td>142941</td><td>1017298</td></tr><tr><td>Elmidae</td><td>Macronychus</td><td>2</td><td>0</td><td>0</td><td>128</td></tr><tr><td>Elmidae</td><td>Oulimnius</td><td>581</td><td>2792</td><td>1315</td><td>18411</td></tr><tr><td>Elmidae</td><td>Riolus</td><td>716</td><td>12770</td><td>22435</td><td>96015</td></tr><tr><td>Gyrinidae</td><td>Orectochilus</td><td>31</td><td>128</td><td>0</td><td>1144</td></tr><tr><td>Haliplidae</td><td>Haliplus</td><td>8</td><td>142</td><td>0</td><td>0</td></tr><tr><td>Hydraenidae</td><td>Hydraena</td><td>336</td><td>829</td><td>61998</td><td>16424</td></tr><tr><td>Hydrophilidae</td><td></td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Hydrophilidae</td><td>Hydrophilinae</td><td>0</td><td>0</td><td>383</td><td>0</td></tr><tr><td>Hydrophilidae</td><td>Sphaeridiinae</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Scirtidae</td><td>Cyphon</td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Scirtidae</td><td>Helodes</td><td>70</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Scirtidae</td><td>Hydrocyphon</td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td rowspan="14">Diptera</td><td>Anthomyidae</td><td></td><td>3</td><td>539</td><td>0</td><td>0</td></tr><tr><td>Athericidae</td><td></td><td>624</td><td>20069</td><td>0</td><td>2895</td></tr><tr><td>Blephariceridae</td><td></td><td>16</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ceratopogonidae</td><td></td><td>296</td><td>3989</td><td>0</td><td>0</td></tr><tr><td>Chaoboridae</td><td></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Chironomidae</td><td></td><td>18992</td><td>14658</td><td>0</td><td>18513</td></tr><tr><td>Culicidae</td><td></td><td>0</td><td>1450</td><td>0</td><td>582</td></tr><tr><td>Dixidae</td><td></td><td>5</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Empididae</td><td></td><td>218</td><td>414</td><td>0</td><td>0</td></tr><tr><td>Ephydridae</td><td></td><td>4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Limoniidae</td><td></td><td>508</td><td>8357</td><td>375057</td><td>428638</td></tr><tr><td>Psychodidae</td><td></td><td>147</td><td>534</td><td>0</td><td>0</td></tr></tbody></table>
1021 Tab.1 (continued)

<table><thead><tr><th>Higher rank</th><th>Family</th><th>Taxon</th><th>Abundance</th><th>Reads (COI)</th><th>Reads (EUKA02)</th><th>Reads (INSE01)</th></tr></thead><tbody><tr><td rowspan="7">Diptera</td><td>Ptychopteridae</td><td></td><td>0</td><td>13</td><td>0</td><td>0</td></tr><tr><td>Rhagionidae</td><td></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Simuliidae</td><td></td><td>4287</td><td>27203</td><td>54558</td><td>780508</td></tr><tr><td>Stratiomyidae</td><td></td><td>9</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Tabanidae</td><td></td><td>8</td><td>14036</td><td>0</td><td>0</td></tr><tr><td>Tipulidae</td><td></td><td>14</td><td>5924</td><td>0</td><td>18095</td></tr><tr><td>Ephemeroptera</td><td>Baetidae</td><td></td><td>0</td><td>0</td><td>0</td><td>21088</td></tr><tr><td rowspan="14">Ephemeroptera</td><td>Baetidae</td><td>Baetis</td><td>7063</td><td>282764</td><td>3766695</td><td>3664057</td></tr><tr><td>Baetidae</td><td>Cloeon</td><td>4</td><td>187</td><td>0</td><td>0</td></tr><tr><td>Caenidae</td><td>Brachycercus</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Caenidae</td><td>Caenis</td><td>330</td><td>0</td><td>25517</td><td>529</td></tr><tr><td>Ephemerellidae</td><td>Ephem器ella</td><td>2757</td><td>33602</td><td>204825</td><td>31769</td></tr><tr><td>Ephemerellidae</td><td></td><td>0</td><td>57357</td><td>686</td><td>0</td></tr><tr><td>Ephemerellidae</td><td>Torleya</td><td>99</td><td>0</td><td>576400</td><td>0</td></tr><tr><td>Ephemeridae</td><td>Ephemera</td><td>53</td><td>25640</td><td>57405</td><td>566199</td></tr><tr><td>Heptageniidae</td><td>Ecdyonurus</td><td>89</td><td>63253</td><td>489674</td><td>1250160</td></tr><tr><td>Heptageniidae</td><td>Electrogena</td><td>0</td><td>0</td><td>0</td><td>649</td></tr><tr><td>Heptageniidae</td><td>Epeorus</td><td>300</td><td>5092</td><td>264696</td><td>83753</td></tr><tr><td>Heptageniidae</td><td>Heptagenia</td><td>1</td><td>2</td><td>681</td><td>259</td></tr><tr><td>Heptageniidae</td><td>Rhithrogena</td><td>780</td><td>75947</td><td>1149041</td><td>2130344</td></tr><tr><td>Leptophlebiidae</td><td>Habroleptoides</td><td>708</td><td>22839</td><td>0</td><td>741183</td></tr><tr><td>Leptophlebiidae</td><td>Habrophlebia</td><td>100</td><td>1280</td><td>291009</td><td>115316</td></tr><tr><td>Leptophlebiidae</td><td>Leptophlebia</td><td>0</td><td>0</td><td>171</td><td>0</td></tr><tr><td>Leptophlebiidae</td><td></td><td>0</td><td>0</td><td>15370</td><td>0</td></tr><tr><td>Leptophlebiidae</td><td>Paraleptophlebia</td><td>6</td><td>0</td><td>223154</td><td>21646</td></tr><tr><td>Potamanthidae</td><td>Potamanthus</td><td>0</td><td>0</td><td>21</td><td>10</td></tr><tr><td>Siphlonuridae</td><td>Siphlonurus</td><td>0</td><td>0</td><td>936</td><td>11</td></tr><tr><td rowspan="5">Hemiptera</td><td>Aphelocheiridae</td><td>Aphelocheirus</td><td>49</td><td>2492</td><td>34159</td><td>129833</td></tr><tr><td>Corixidae</td><td>Corixinae</td><td>1</td><td>0</td><td>692</td><td>1336</td></tr><tr><td>Corixidae</td><td>Micronecta</td><td>119</td><td>0</td><td>559</td><td>0</td></tr><tr><td>Mesoveliidae</td><td>Mesovelia</td><td>5</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Agriotypidae</td><td>Agriotypus</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Hymenoptera</td><td>Crambidae</td><td></td><td>1</td><td>0</td><td>1917</td><td>114</td></tr><tr><td>Lepidoptera</td><td>Sialidae</td><td>Sialis</td><td>8</td><td>467</td><td>0</td><td>0</td></tr><tr><td rowspan="10">Megaloptera</td><td>Aeshnidae</td><td>Aeshna</td><td>0</td><td>2316</td><td>0</td><td>0</td></tr><tr><td>Aeshnidae</td><td>Anax</td><td>2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Calopterygidae</td><td>Calopteryx</td><td>24</td><td>270</td><td>0</td><td>66203</td></tr><tr><td>Coenagrionidae</td><td></td><td>2</td><td>0</td><td>345</td><td>0</td></tr><tr><td>Cordulegasteridae</td><td>Cordulegaster</td><td>16</td><td>283</td><td>0</td><td>49980</td></tr><tr><td>Gomphidae</td><td>Gomphus</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Gomphidae</td><td>Onychogomphus</td><td>1</td><td>982</td><td>0</td><td>0</td></tr><tr><td>Platycnemididae</td><td>Platycnemis</td><td>1</td><td>804</td><td>0</td><td>825</td></tr></tbody></table>
<table><caption>1023 Tab.1 (continued)</caption><thead><tr><th>Higher rank</th><th>Family</th><th>Taxon</th><th>Abundance</th><th>Reads (COI)</th><th>Reads (EUKA02)</th><th>Reads (INSE01)</th></tr></thead><tbody><tr><td rowspan="20">Plecoptera</td><td>Capniidae</td><td>Capnia</td><td>2</td><td>0</td><td>24</td><td>0</td></tr><tr><td>Capniidae</td><td></td><td>0</td><td>0</td><td>0</td><td>32</td></tr><tr><td>Capniidae</td><td>Capnioneura</td><td>14</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Chloroperlidae</td><td>Chloroperla</td><td>27</td><td>278</td><td>0</td><td>26272</td></tr><tr><td>Chloroperlidae</td><td></td><td>0</td><td>28</td><td>0</td><td>0</td></tr><tr><td>Chloroperlidae</td><td>Siphonoperla</td><td>559</td><td>7380</td><td>1711584</td><td>436685</td></tr><tr><td>Leuctridae</td><td>Leuctra</td><td>1563</td><td>2704</td><td>779840</td><td>1232887</td></tr><tr><td>Nemouridae</td><td>Amphinemura</td><td>1018</td><td>4517</td><td>0</td><td>362313</td></tr><tr><td>Nemouridae</td><td>Nemoura</td><td>278</td><td>355</td><td>0</td><td>2335</td></tr><tr><td>Nemouridae</td><td></td><td>0</td><td>76</td><td>941894</td><td>7</td></tr><tr><td>Nemouridae</td><td>Nemurella</td><td>2</td><td>0</td><td>0</td><td>10378</td></tr><tr><td>Nemouridae</td><td>Protonemura</td><td>3469</td><td>129089</td><td>0</td><td>2148653</td></tr><tr><td>Perlidae</td><td>Dinocras</td><td>13</td><td>3938</td><td>31865</td><td>27300</td></tr><tr><td>Perlidae</td><td>Perla</td><td>36</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Perlodidae</td><td>Dictyogenus</td><td>5</td><td>0</td><td>0</td><td>32162</td></tr><tr><td>Perlodidae</td><td>Diura</td><td>2</td><td>0</td><td>132202</td><td>225296</td></tr><tr><td>Perlodidae</td><td>Isoperla</td><td>175</td><td>20271</td><td>0</td><td>89548</td></tr><tr><td>Perlodidae</td><td>Perlodes</td><td>3</td><td>11896</td><td>0</td><td>78653</td></tr><tr><td>Perlodidae</td><td></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td rowspan="10">Trichoptera</td><td>Taeniopterygidae</td><td>Brachyptera</td><td>4</td><td>0</td><td>1845</td><td>4901</td></tr><tr><td>Taeniopterygidae</td><td>Rhabdiopteryx</td><td>0</td><td>0</td><td>712</td><td>0</td></tr><tr><td>Taeniopterygidae</td><td>Taeniopteryx</td><td>62</td><td>148</td><td>3123</td><td>0</td></tr><tr><td>Beraeidae</td><td>Beraeodes</td><td>0</td><td>0</td><td>13651</td><td>0</td></tr><tr><td>Beraeidae</td><td>Ernodes</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Brachycentridae</td><td></td><td>0</td><td>4465</td><td>0</td><td>0</td></tr><tr><td>Brachycentridae</td><td>Brachycentrus</td><td>8</td><td>0</td><td>210</td><td>0</td></tr><tr><td>Brachycentridae</td><td>Micrasema</td><td>670</td><td>8496</td><td>17315</td><td>0</td></tr><tr><td>Glossosomatidae</td><td>Agapetus</td><td>581</td><td>0</td><td>1230</td><td>0</td></tr><tr><td>Glossosomatidae</td><td>Glossosoma</td><td>157</td><td>3385</td><td>68511</td><td>0</td></tr><tr><td>Goeridae</td><td>Goera</td><td>5</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Goeridae</td><td></td><td>1</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Goeridae</td><td>Lithax</td><td>60</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Goeridae</td><td>Silo</td><td>92</td><td>1890</td><td>26334</td><td>398</td></tr><tr><td>Hydropsychidae</td><td>Cheumatopsyche</td><td>79</td><td>1247</td><td>1127</td><td>0</td></tr><tr><td>Hydropsychidae</td><td>Hydropsyche</td><td>3120</td><td>91288</td><td>7170</td><td>378</td></tr><tr><td>Hydroptilidae</td><td>Hydroptila</td><td>208</td><td>435</td><td>0</td><td>0</td></tr><tr><td>Hydroptilidae</td><td>Orthotrichia</td><td>0</td><td>1474</td><td>0</td><td>0</td></tr><tr><td>Hydroptilidae</td><td>Ptilocolepus</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Hydroptilidae</td><td>Stactobia</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lepidostomatidae</td><td>Crunoecia</td><td>6</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lepidostomatidae</td><td>Lasiocephala</td><td>4</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lepidostomatidae</td><td>Lepidostoma</td><td>316</td><td>1463</td><td>0</td><td>0</td></tr></tbody></table>
<table><caption>1025 Tab.1 (continued)</caption><thead><tr><th>Higher rank</th><th>Family</th><th>Taxon</th><th>Abundance</th><th>Reads (COI)</th><th>Reads (EUKA02)</th><th>Reads (INSE01)</th></tr></thead><tbody><tr><td rowspan="15">Trichoptera</td><td>Leptoceridae</td><td>Adicella</td><td>20</td><td>326</td><td>269</td><td>0</td></tr><tr><td>Leptoceridae</td><td>Athripsodes</td><td>7</td><td>0</td><td>0</td><td>122</td></tr><tr><td>Leptoceridae</td><td>Mystacides</td><td>40</td><td>258</td><td>187</td><td>0</td></tr><tr><td>Leptoceridae</td><td>Oecetis</td><td>4</td><td>178</td><td>304</td><td>0</td></tr><tr><td>Limnephilidae</td><td>Drusinae</td><td>80</td><td>4998</td><td>25279</td><td>0</td></tr><tr><td>Limnephilidae</td><td></td><td>10</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Limnephilidae</td><td>Limnephilinae</td><td>1133</td><td>66589</td><td>100892</td><td>186822</td></tr><tr><td>Odontoceridae</td><td>Odontocerum</td><td>253</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Philopotamidae</td><td>Philopotamus</td><td>72</td><td>970</td><td>0</td><td>0</td></tr><tr><td>Philopotamidae</td><td>Wormaldia</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Polycentropodidae</td><td>Holocentropus</td><td>40</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Polycentropodidae</td><td>Plectrocnemia</td><td>105</td><td>0</td><td>0</td><td>8023</td></tr><tr><td>Polycentropodidae</td><td>Polycentropus</td><td>137</td><td>171</td><td>2654</td><td>0</td></tr><tr><td>Psychomyiidae</td><td>Type</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Psychomyiidae</td><td>Psychomyia</td><td>530</td><td>1041</td><td>12772</td><td>119483</td></tr><tr><td>Psychomyiidae</td><td>Tinodes</td><td>70</td><td>0</td><td>4177</td><td>173</td></tr><tr><td>Rhyacophilidae</td><td>Rhyacophila</td><td>677</td><td>8312</td><td>108836</td><td>910479</td></tr><tr><td>Sericostomatidae</td><td>Sericostoma</td><td>821</td><td>51100</td><td>0</td><td>546599</td></tr><tr><td>Sericostomatidae</td><td></td><td>0</td><td>0</td><td>87001</td><td>0</td></tr><tr><td>Hydrozoa</td><td></td><td></td><td>0</td><td>2</td><td>0</td><td>0</td></tr><tr><td rowspan="2">Bivalvia</td><td>Sphaeriidae</td><td>Pisidium</td><td>444</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Sphaeriidae</td><td>Sphaerium</td><td>0</td><td>0</td><td>323460</td><td>2810</td></tr><tr><td rowspan="15">Gastropoda</td><td>Acroloxidae</td><td>Acroloxus</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Ancylidae</td><td>Ancylus</td><td>342</td><td>5113</td><td>699087</td><td>425</td></tr><tr><td>Bithyniidae</td><td>Bithynia</td><td>6</td><td>114</td><td>4477</td><td>0</td></tr><tr><td>Hydrobiidae</td><td>Potamopyrgus</td><td>145</td><td>961</td><td>273236</td><td>0</td></tr><tr><td>Lymnaeidae</td><td>Galba</td><td>3</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lymnaeidae</td><td>Lymnaea</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Lymnaeidae</td><td></td><td>1</td><td>0</td><td>14</td><td>0</td></tr><tr><td>Lymnaeidae</td><td>Radix</td><td>76</td><td>11975</td><td>484006</td><td>0</td></tr><tr><td>Neritidae</td><td>Theodoxus</td><td>0</td><td>82</td><td>0</td><td>0</td></tr><tr><td>Physidae</td><td>Physa stricto-sensu</td><td>17</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Physidae</td><td>Physella</td><td>10</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Physidae</td><td></td><td>0</td><td>530</td><td>223269</td><td>0</td></tr><tr><td>Planorbidae</td><td></td><td>10</td><td>0</td><td>960</td><td>0</td></tr><tr><td>Valvatidae</td><td>Valvata</td><td>7</td><td>0</td><td>209</td><td>0</td></tr><tr><td>Nemathelmintha</td><td></td><td></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td rowspan="2">Turbellaria</td><td>Dugesiidae</td><td></td><td>296</td><td>221</td><td>270076</td><td>0</td></tr><tr><td>Planariidae</td><td></td><td>527</td><td>1847</td><td>1671815</td><td>0</td></tr><tr><td>Porifera</td><td>Spongillidae</td><td></td><td>0</td><td>0</td><td>108</td><td>0</td></tr></tbody></table>

<caption>1026</caption>
Tab. 2. Summary of the results corresponding to the best match between morphology-based and bulk-DNA based approaches, obtained with each individual primer set and each combination of primer sets. Wilcoxon tests and adjusted R² (from linear regressions) are applied/calculated between the values of the I₂M₂ index obtained with morphology-based versus bulk-DNA based approaches. The thresholds for best identity and for minimum number of reads are respectively provided in brackets for each primer set or combination of primer sets. DE = Discrimination Efficiency.

<table><thead><tr><th rowspan="2">Bulk data</th><th rowspan="2">Adjusted R² DE</th><th rowspan="2">Wilcoxon tests (p-value)</th><th colspan="6">Reclassifications</th><th rowspan="2">Total</th></tr><tr><th colspan="4">Worse predicted state</th><th colspan="2">No changes</th><th colspan="2">Better predicted state</th></tr><tr><th></th><th></th><th></th><th>-4</th><th>-3</th><th>-2</th><th>-1</th><th>0</th><th>+1</th><th>+2</th><th></th></tr></thead><tbody><tr><td>COI (83:0)</td><td>0.396 0.833</td><td>0.119</td><td>1</td><td></td><td></td><td>4</td><td>27</td><td>4</td><td></td><td>9/36</td></tr><tr><td>EUKA02 (96:85)</td><td>0.544 0.751</td><td>0.001</td><td></td><td></td><td>3</td><td>11</td><td>19</td><td>2</td><td>1</td><td>17/36</td></tr><tr><td>INSE01 (100:10)</td><td>0.611 0.875</td><td>0.010</td><td></td><td>1</td><td></td><td>7</td><td>26</td><td>2</td><td></td><td>10/36</td></tr><tr><td>COI (82:5) + EUKA02 (95:85)</td><td>0.648 0.875</td><td>0.002</td><td></td><td>1</td><td>1</td><td>7</td><td>25</td><td>1</td><td>1</td><td>11/36</td></tr><tr><td>COI (84:10) + INSE01 (100:10)</td><td>0.740 0.875</td><td>0.003</td><td></td><td></td><td>1</td><td>9</td><td>23</td><td>2</td><td></td><td>12/36</td></tr><tr><td>INSE01 (100:5) + EUKA02 (96:85)</td><td>0.697 0.875</td><td>0.000</td><td></td><td></td><td>2</td><td>9</td><td>23</td><td>2</td><td></td><td>13/36</td></tr><tr><td>COI (82:5) + EUKA02 (95:85) + INSE01 (100:5)</td><td>0.698 0.958</td><td>0.002</td><td></td><td></td><td>1</td><td>8</td><td>25</td><td>2</td><td></td><td>11/36</td></tr></tbody></table>

1034
1035
<figure><img src="image_2.png" /></figure>

Fig. 1: Location of the sampled sites. Delineations indicate the main French hydro-eco-regions (Wasson et al. 2002). The number in the code name of each sampled site gives the information about its group membership (from 1 to 6), whereas the letter gives information about its status, either a LIRR ("Reference"; "R") or a IRR, exhibiting either a Moderate "M" pressure intensity or a Strong "S" pressure intensity (see first section of Material & Methods for further details). All three sites found in a group exhibit the same stream type. Further information about the sites are available in the Supplementary Material (Table S1).
<figure><img src="image_3.png" /></figure>

1047

1048 Fig. 2: Workflow. Dashed double-arrowed lines indicate comparisons of the results obtained

1049 with both the morphological-based approach and the bulk DNA one.

1050

1051
<figure><img src="image_4.png" /><figcaption>Fig. 3: Boxplots describing the distributions of adjusted R² values derived from linear regressions between I₂M₂ index values calculated on relative abundances (morphology-based approach) versus presence/absence (PA) or relative frequencies of reads (RF) for each primer set. N = 21 thresholds for best identity × 21 thresholds for minimal number of reads = 441 values per couple of “primer × data type”. Statistical differences between PA vs RF, for a given primer, were investigated with Wilcoxon tests for paired data, with *** if p < 0.001. Each boxplot</figcaption></figure>
1060 represents minimum/Q25/median/Q75/maximum values, respectively. Outliers (open circles)

1061 are outside the 1.5 interquartile range of the corresponding adjusted R² value distribution.

1062
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>
Fig. 4: Boxplots of the values of adjusted R² of linear regression between values of I₂M₂ based on morphology-based data and I₂M₂ calculated on bulk DNA for each primer set (A = COI, B = EUKA02, C = INSE01) and for different thresholds for best identity. Reads were expressed as relative frequencies. N = 21 for each boxplot. Each boxplot represents minimum/Q25/median/Q75/maximum values, respectively. Outliers (open circles) are outside the 1.5 interquartile range of the corresponding adjusted R² value distribution. In grey, selected thresholds for the following analyses.
<figure><img src="image_8.png" /></figure>

1072

1073

1074

1075

1076

1077

1078

<caption>Fig. 5: Values of the adjusted R² of linear regressions between I₂M₂ calculated on morphology-based data and I₂M₂ values calculated on bulk DNA-based data for each primer set (reads are expressed in relative frequencies) according to different values of threshold for the minimal number of reads, and calculated for the three previously selected thresholds for best identity (cf. Fig. 4). Arrows indicate the five selected thresholds for minimal number of reads for the following analyses.</caption>
<figure><img src="image_9.png" /><figcaption>Fig. 6: Linear regressions (solid lines) between the values of the I₂M₂ index (and its associated metrics) calculated from the morphology-based versus bulk DNA-based data, for each primer set and combinations of primer sets. Results presented here are for the best combination (= maximum adjusted R²) in terms of thresholds for best identity and minimum number of reads within the range of tested values. All metric values are expressed in Ecological Quality Ratios (EQR). Dashed lines (in right graphs) are lines with a slope equals to 1. N = 36 (18 streams x 2 years). Top-left values are adjusted R².</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 7: Boxplots of the values of adjusted R² corresponding to linear regressions between I₂M₂ values calculated on morphology-based and bulk DNA-based data, for each combination of primer sets. Read numbers were first transformed in relative frequencies within data provided by each primer set separately, before being summed for each primer set combination and once again transformed in relative frequencies. N = 225 for each pair of primer sets, and N = 3375 for the combination of the three primer sets.</figcaption></figure>
1095 Supplementary material

1096

1097 Tab. S1. Detailed location of the sampled sites. Stream types are according to Wasson et al.

1098 (2002).<sup>a</sup> this site was not considered as a reference for the definition of bulk-DNA based

1099 "BEST" values, as it was an acidic stream (pH = 6.4). Pressure intensity: "M" = moderate, "S"

1100 strong.

1101

<table><thead><tr><th>Code</th><th>Stream</th><th>Site</th><th>Longitude (WGS84)</th><th>Latitude (WGS84)</th><th>Stream type</th><th>Pressure category</th></tr></thead><tbody><tr><td>1M</td><td>GELON</td><td>VILLARD</td><td>6.1215540</td><td>45.4770249 TP2</td><td></td><td>Hydrological alteration (M)</td></tr><tr><td>1R</td><td>CROP</td><td>VILLARD-BONNOT</td><td>5.9749229</td><td>45.2224414 TP2</td><td></td><td>Reference</td></tr><tr><td>1S</td><td>GELON</td><td>VERNEIL</td><td>6.1417991</td><td>45.4646720 TP2</td><td></td><td>Hydrological alteration (S)</td></tr><tr><td>2M</td><td>BOURBRE</td><td>ST-CLAIR</td><td>5.4574619</td><td>45.5892812 P5</td><td></td><td>Organic matter (M)</td></tr><tr><td>2R</td><td>SEMINE</td><td>BELLEYDOUX</td><td>5.7885934</td><td>46.2489425 P5</td><td></td><td>Reference</td></tr><tr><td>2S</td><td>BOURBRE</td><td>CESSIEU</td><td>5.3987061</td><td>45.5647726 P5</td><td></td><td>Organic matter (S)</td></tr><tr><td>3M</td><td>OIGNIN</td><td>BRION</td><td>5.5419970</td><td>46.1744351 M5</td><td></td><td>Organic matter (M)</td></tr><tr><td>3R</td><td>CHERAN</td><td>CARLET</td><td>6.1829477</td><td>45.6528586 M5</td><td></td><td>Reference</td></tr><tr><td>3S</td><td>OIGNIN</td><td>MATAFELON</td><td>5.5525989</td><td>46.2655298 M5</td><td></td><td>Hydrological alteration (S)</td></tr><tr><td>4M</td><td>RAVINE</td><td>RAON L'ETAPE</td><td>6.9292740</td><td>48.4087970 TP4</td><td></td><td>Acidification (M)</td></tr><tr><td>4R</td><td>PLAINE</td><td>RAON SUR PLAINE</td><td>7.0612097</td><td>48.4986129 P4</td><td></td><td>Reference</td></tr><tr><td>4S</td><td>GENTIL SAPIN</td><td>SENONES</td><td>7.0400186</td><td>48.4536782 TP4</td><td></td><td>Acidification (S)</td></tr><tr><td>5M</td><td>VEZOUZE</td><td>CIREY</td><td>6.9373024</td><td>48.5896737 P4</td><td></td><td>Organic matter (M)</td></tr><tr><td>5R</td><td>SAVOUREUSE</td><td>LEPUIX</td><td>6.8261808</td><td>47.7677172 P4</td><td></td><td>Reference</td></tr><tr><td>5S</td><td>VEZOUZE</td><td>BLAMONT</td><td>6.8384378</td><td>48.5884879 P4</td><td></td><td>Organic matter (S)</td></tr><tr><td>6M</td><td>CLEURIE</td><td>LA FORGE</td><td>6.7042109</td><td>48.0559488 P4</td><td></td><td>Toxic/Chemical (M)</td></tr><tr><td>6R</td><td>MENAURUPT</td><td>SAPOIS</td><td>6.7868727</td><td>48.0290126 P4</td><td></td><td>Reference<sup>a</sup></td></tr><tr><td>6S</td><td>CLEURIE</td><td>LE THOLY</td><td>6.7252253</td><td>48.0891638 P4</td><td></td><td>Toxic/Chemical (S)</td></tr></tbody></table>

1102

1103