Reference databases, primer choice, and assay sensitivity for environmental metabarcoding: Lessons learnt from a re‐evaluation of an eDNA fish assessment in the Volga headwaters

Biodiversity monitoring via environmental DNA, particularly metabarcoding, is evolving into a powerful assessment tool for riverine systems. However, for metabarcoding to be fully integrated into standardized monitoring programmes, some current challenges concerning sampling design, laboratory workflow, and data analysis need to be overcome. Here, we review some of these major challenges and potential solutions. We further illustrate three potential pitfalls, namely the choice of suitable metabarcoding primers, the necessity of complete reference databases, and varying assay sensitivities, by a reappraisal of our‐own recently carried out metabarcoding study in the Volga headwaters. TaqMan qPCRs had detected catfish (Silurus glanis) and European eel (Anguilla anguilla), whereas metabarcoding had not, in the same samples. Furthermore, after extending the genetic reference database by 12 additional species and re‐analysing the metabarcoding data, we additionally detected the Siberian spiny loach (Cobitis sibirica) and Ukrainian brook lamprey (Eudontomyzon mariae) and reassigned the operational taxonomic units previously assigned to Misgurnus fossilis to Cobitis sibirica. In silico analysis of metabarcoding primer efficiencies revealed considerable variability among primer pairs and among target species, which could lead to strong primer bias and potential false‐negatives in metabarcoding studies if not properly compensated for. These results highlight some of the pitfalls of eDNA‐metabarcoding as a means of monitoring fish biodiversity in large rivers, which need to be considered in order to fully unleash the full potential of these approaches for freshwater biodiversity monitoring.

tial solutions. We further illustrate three potential pitfalls, namely the choice of suitable metabarcoding primers, the necessity of complete reference databases, and varying assay sensitivities, by a reappraisal of our-own recently carried out metabarcoding study in the Volga headwaters. TaqMan qPCRs had detected catfish (Silurus glanis) and European eel (Anguilla anguilla), whereas metabarcoding had not, in the same samples. Furthermore, after extending the genetic reference database by 12 additional species and re-analysing the metabarcoding data, we additionally detected the Siberian spiny loach (Cobitis sibirica) and Ukrainian brook lamprey (Eudontomyzon mariae) and reassigned the operational taxonomic units previously assigned to Misgurnus fossilis to Cobitis sibirica. In silico analysis of metabarcoding primer efficiencies revealed considerable variability among primer pairs and among target species, which could lead to strong primer bias and potential false-negatives in metabarcoding studies if not properly compensated for. These results highlight some of the pitfalls of eDNA-metabarcoding as a means of monitoring fish biodiversity in large rivers, which need to be considered in order to fully unleash the full potential of these approaches for freshwater biodiversity monitoring.

K E Y W O R D S
eDNA, false negatives, fish, inhibition, metabarcoding, primer bias, qPCR, real-time PCR 1 | eDNA-METABARCODING AS A TOOL FOR BIODIVERSITY ASSESSMENTS AND ITS ASSOCIATED CHALLENGES Biodiversity assessments are vital to every biomonitoring effort and aim to record all resident taxa of specific bioindicator groups in a target area, like a river catchment. A potentially powerful addition to conventional biodiversity surveys is DNA-based assessment. The utilization of environmental DNA (or eDNA; Taberlet, Coissac, Pompanon, Brochmann, & Willerslev, 2012) for biomonitoring has been promoted by numerous researchers, particularly for freshwater ecosystems (e.g., Hering et al., 2018;Senapati et al., 2018). Bioassessments utilizing eDNA have various positive assets that can compensate for the drawbacks and limitations associated with conventional sampling and species identification methods. These include non-invasiveness, better cost and time effectiveness, and increased taxonomic identification (Deiner et al., 2017). For overall biodiversity assessment of a target species community, eDNA metabarcoding is the ideal method. This approach uses universal primers to bind to the eDNA left by members of the target group in an environmental sample and amplifies a genetic barcode for each of these species in parallel during polymerase chain reaction (PCR). The amplicons then get sequenced using high-throughput sequencing, and species composition is reconstructed using the retrieved DNA sequences.
However, eDNA and metabarcoding approaches are still under development, and multiple issues need to be resolved before they can be fully integrated into standard biomonitoring systems. These issues concern the influence of various factors during eDNA sampling, laboratory workflow, and data analysis on detection probabilities. Even before sampling, researchers need to consider DNA persistence in the environment, which has been shown to be influenced by factors such as pH, temperature, UV radiation, flow, and chemical properties of bedrock (e.g., Barnes & Turner, 2016;Seymour et al., 2018;Stoeckle et al., 2017).
For eDNA sampling, storage, and extraction, a plurality of options exist, and the choice of method can significantly influence detection probabilities (e.g., Deiner, Walser, Maechler, & Altermatt, 2015;Spens et al., 2017), making it increasingly difficult to pick the optimal protocol without preliminary information on the study system. A crucial step in any metabarcoding analysis is the PCR. Ideally, PCR amplifies all target DNA at an equal and exponential rate. However, this is rarely true due to various reasons. First, eDNA samples frequently contain PCR inhibitors, hampering polymerase function and therefore increasing the risk of false negatives, particularly when working with species occurring at low abundance. Inhibitors might stem from abiotic or biological sources such as bedrock, algae, or organic input (McKee, Spear, & Pierson, 2015;Stoeckle et al., 2017). Inhibitor-removal kits have been shown to efficiently remove inhibitors but add costs. While inhibition influences PCR function in general, PCR bias favours particular template DNA strands from the mixed DNA template. PCR bias can stem from (a) PCR stochasticityrandom preference of individual DNA strands during PCR reaction (Kebschull & Zador, 2015); (b) preference of more abundant DNA templates, potentially masking the presence of low-abundant species (Bylemans, Gleeson, Duncan, Hardy, & Furlan, 2019); and (c) primer bias. Primer bias originates from sequence variation in the primer annealing sites among individual species, leading to higher efficiency of the PCR reaction for certain species in a mixed DNA template, while other species might get amplified less or not at all (Bylemans, Gleeson, Hardy, & Furlan, 2018;Elbrecht & Leese, 2015). While the first two canat least to some extentbe overcome by performing multiple PCR replicates (Ficetola et al., 2015), primer bias is best avoided by careful evaluation of the metabarcoding primers used, in order to avoid polymorphism in the DNA binding region. This goal might be complicated by additional primer criteria, such as small amplicon size, high taxonomic resolution of the amplified fragment, and target group specificity (Bylemans et al., 2018). Therefore, primer performance is best evaluated for the specific target community and geographic region under study (Elbrecht & Leese, 2017).
After data generation and raw data filtration, the retrieved DNA sequences are usually compared with a reference database of the expected species community in order to translate the obtained molecular operational taxonomic units into biological species for final data interpretation. However, such reference databases are still far from complete despite global initiatives like the international Barcode of Life (iBOL; Weigand et al., 2019). Additionally, the target fragment of iBOL, a part of the mitochondrial cytochrome oxidase subunit I (COI), has been shown to be less suitable for metabarcoding work for vertebrates (Deagle, Jarman, Coissac, Pompanon, & Taberlet, 2014). Therefore, other mitochondrial markers, particularly fragments of the mitochondrial 12S and 16S rRNA genes, have been the focus for metabarcoding primer design for this group (Evans et al., 2016;Miya et al., 2015). Databases for these gene fragments are less complete than for COI and have been shown to be heavily biased for certain geographic regions and taxonomic groups (Belle, Stoeckle, & Geist, 2019;Weigand et al., 2019), highlighting the need to expand databases for understudied areas and taxonomic groups.
An alternative to metabarcoding are species-specific approaches like real-time PCR (qPCR). Such assays are designed to specifically amplify only the target species from the bulk DNA in a sample. This approach does not require DNA sequencing and genetic reference databases. Instead, species presence is inferred via the amplification itself during the PCR step in real time. Whereas the overall informative value of qPCR is taxonomically narrow compared to metabarcoding, its application can circumvent some of metabarcoding's challenges. These include the elimination of primer bias, as the primers target one species only, easier detection of inhibition, as the amplification efficiency can be directly observed, and a higher potential to reliably correlate the results to species abundances (Doi et al., 2017;Lacoursière-Roussel, Côté, Leclerc, & Bernatchez, 2016). Additionally, such a targeted approach has been shown to be more sensitive than metabarcoding (Bylemans et al., 2019;Harper et al., 2018) and is more cost effective when the target taxa involve only a few species (Tsuji et al., 2018).
We tested for the influence of three potential problemsan incomplete reference database, PCR inhibition and primer biason the results of our own, recently carried out, metabarcoding study in the headwaters of the Volga River (Lecaudey, Schletterer, Kuzovlev, Hahn, & Weiss, 2019). Lecaudey et al. (2019), detected 23 fish species across the sampled stretches, but the number of detected species varied considerably between the three barcode markers used, namely 12S rRNA, 16S rRNA, and cytochrome B (CytB). The Ukrainian brook lamprey (Eudontomyzon mariae), blue bream (Ballerus ballerus), white-eye bream (Ballerus sapa), asp (Leuciscus aspius), Siberian spiny loach (Cobitis sibirica), wels catfish (Silurus glanis), and European eel (Anguilla anguilla) were not detected by any of the three markers used. The first five species lacked reference sequences for the 12S and 16S rRNA, which had an overall better performance than CytB for species detection. For the latter two species (S. glanis and A. anguilla), reference sequences were available and recent presence at the sampled stretches has been documented (Schletterer, 2006;Schletterer et al., 2018). Additionally, potential presence of inhibitors was observed in the laboratory during PCR steps. Therefore, in this re-evaluation, we aimed to (a) verify the absence of S. glanis and We performed TaqMan qPCRs for S. glanis and A. anguilla using previously published protocols (Jensen, Knudsen, Munk, Thomsen, & Møller, 2018;Roy, Belliveau, Mandrak, & Gagné, 2018). Sampling sites which either represent suitable habitat for target species or for which recent sightings were reported were identified (Figure 1). For each relevant site, the species-specific assay was applied to at least three randomly chosen sampling replicates and 10 qPCR replicates were carried out for each. An eDNA sample was recorded as positive when positive amplification was observed in at least two independent qPCR replicates. If a PCR replicate only yielded 1 positive amplification out of 10 amplifications trials, additional samples at that location were subjected to qPCR analysis until a clear result was obtained. We evaluated the potential presence of inhibitors by spiking positives controls with an aliquot of the eDNA extracts. Details on qPCR reaction set-up and inhibitor testing are given in Data S1.
2.1.2 | Expansion of genetic reference data and reanalysis of the metabarcoding data We produced reference sequences for the targeted 12S and 16S rRNA fragments for Lampetra planeri, E. mariae, B. ballerus, B. sapa, L. aspius, and C. sibirica. The third marker, CytB, was excluded for genetic reference data generation because of its generally poor performance in the previous metabarcoding analysis. Most
For this, we expanded the alignments of the reference databases at both the 5 0 and 3 0 ends, to include the primer binding regions of the respective primer sets. In total, 54, 56, and 43 species were included for this analysis for 12S and 16S rRNA and CytB, respectively, using a total of 153, 141, and 100 sequences (Data S1).
Primer efficiency analysis was carried out via the DECIPHER package v.2.10.2 (Wright, 2015) in R v.3.5.2 (R Core Team, 2018) with PCR amplification efficiency ranging from 0 (very low efficiency) to 1 (maximum efficiency possible) and calculated based on the hybridization and elongation efficiency of the primers to the target sites (settings: temp = 58 ºC for 12S and 16S rRNA, and 53 ºC for CytB, p = .0003, ions = 0.1558, TaqEfficiency = True, maxDistance = 0.4, maxGaps = 2 and processors = Null). We then calculated the lower mean efficiency (LM Eff ) value between the forward and reverse primer for each species. Using the data of the positive control (mock community), we tested for each marker separately for correlations between (a) LM Eff and the number of reads assigned to each species and (b) LM Eff , and whether the species was scored as 'detected' in the metabarcoding data, using both no threshold for detection (i.e., detection was scored as positive as soon as one read was assigned to that species) and applying the threshold levels used by Lecaudey et al. (2019). Correlations were calculated for tests of (a) using Pearson correlation coefficient, whereas for (b) point-biserial correlations were applied. All statistical testing was carried out in SPSS 24.0 (IBM Corp, 2016).

| TaqMan qPCRs
Throughout both TaqMan qPCR assays, no amplification of negative controls was observed. A total of 16 eDNA samples stemming from the five relevant sampling sites were typed with the S. glanis TaqMan assay. Four samples showed positive amplification in two or more out of the 10 PCR replicates ( Table 1). Three of these samples stem from Mel'nikovo at the Tvertsa River, whereas one is from the Volga River, at Rzhev. For A. anguilla, a total of 18 eDNA samples from four locations were tested (Table 2) testing in duplicate, suggesting inhibition. After repeating these five samples with five additional PCR replicates, only one sample-assay combination (sample 2 at Volga-Staritsa for the A. anguilla assay) still yielded a ΔC T ≥ 2 (2.58).

| Extension of genetic reference data and reanalysis of the metabarcoding data
We successfully sequenced target fragments of the 12S and 16S rRNA of two L. planeri, one E. mariae, two B. sapa, one B. ballerus, two L. aspius, and one C. sibirica (GenBank accession numbers: Data S1), providing the first reference sequences for the target fragments in these species. The extended genetic databases for metabarcoding consisted of a total of 235 reference sequences from 68 species for 12S rRNA, 193 sequences from 63 species for 16S rRNA, and 326 sequences from 66 species for CytB. After metabarcoding data re-analysis, we retrieved significant blast hits on two additional species, namely C. sibirica and E. mariae. The latter was detected in the Moksha River from the 12S rRNA dataset. The operational taxonomic units (OTUs) assigned to E. mariae had all been assigned to L. planeri in Lecaudey et al. (2019). C. sibirica was detected at two sampling sites in the Volga (Staritsa and Rhzev) from both 12S and 16S rRNA data.
The OTUs assigned to C. sibirica had been assigned to C. taenia or the family Cobitidae in Lecaudey et al. (2019). Additionally, all OTUs assigned to the Misgurnus genus in Lecaudey et al. (2019) were now assigned to the newly added C. taenia sequences, eliminating M. fossilis from the final species list.
Considering the additions and one removal of species from the detected species list, the overall number of detected species in the studied free-flowing section of the headwaters of the Volga River now amounts to 26 (Table 3).

| Testing for primer bias
The observed PCR efficiencies of primers showed overall high mean efficiency values, but individual efficiencies of certain primer-species combinations revealed significant drops (Table 4, Figure 2), potentially leading to strong primer bias. A total of seven species-primer combinations revealed a primer efficiency of zero: one species for the 12S rRNA reverse primer, one species for the 16 s rRNA forward and reverse primers, respectively, and four species for the CytB reverse primer.
Among the three markers, the 12S rRNA primer pair had the highest mean primer efficiencies for both forward and reverse primers. All but five species analysed (9.3%) showed almost optimal annealing conditions (LM Eff > 0.9), for both 12S rRNA primers. Concerning the 16S rRNA primers, 10 species yielded an LM Eff < 0.9 (17.9%). CytB had the lowest overall mean efficiency for both forward and reverse primers and the highest variability among species analysed. In total, 22 species (51.1%) revealed an LM Eff < 0.9 for CytB.
Overall, the mean primer pair efficiencies correlated well with the number of species detected by them in the original study (Table 4).
Among the six species of the mock community, LM Eff varied between 0 and 0.999 (Data S1). The 16S rRNA primers had the lowest mean LM Eff across all six species compared to 12S rRNA and CytB primers (0.518, 0.831, and 0.702, respectively). This went along with these primers being able to recover only three out of the six species of the mock community when no detection threshold was applied, whereas 12S rRNA and CytB both recovered five of the six species.
However, when the respective detection threshold was applied, the combined, the correlation between LM Eff and read count (p onetailed = .04; r = 0.42, df = 16) was significant, as was correlation between LM Eff and detection with and without a threshold (p = .009, r = 0.6, df = 16 and p < .0001, r = 0.83, df = 16, respectively).

| TaqMan qPCRs
TaqMan qPCRs revealed unambiguous detection of S. glanis and A. anguilla at two and one sampling sites, respectively, which were missed with the metabarcoding data in Lecaudey et al. (2019). In Lecaudey et al. (2019), one OTU had been assigned to S. glanis at Mel'nikovo using 12S rRNA, but the assigned read count (6) had fallen below the set threshold for false-positive elimination and therefore was not counted as positive detection. These results are concordant with previous studies that showed a higher sensitivity of targeted qPCR compared to metabarcoding (Bylemans et al., 2019;Harper et al., 2018), indicating a better performance of targeted approaches particularly for species of low abundance. However, a recently developed TaqMan qPCR assay did not detect sterlet (Acipenser ruthenus) in the same system, despite recent reports of its presence at individual sites , potentially indicating the border of sensitivity of such assays for a species occurring at very low densities. The sensitivity of eDNA assays will always be influenced by the filtered water volume, storage and extraction methods, and the number of PCR replicates performed, as well as both metabarcoding thresholds or qPCR detection thresholds (as also discussed by Goldberg et al., 2016;Harper et al., 2019). Standardizations across different systems will be hard to achieve, as environmental conditions affecting the sensitivity of eDNA assays will very much differ. However, good scientific practice, such as meticulous reporting of the sampling and laboratory protocols and on experimentally assessed limits of detection or quantification, can nevertheless deliver comprehensible and interpretable results.
Overall, we did not detect strong signals of inhibition during our TaqMan qPCR assay. This is presumably due to both a spin-column treatment and the use of an inhibitor-tolerant enzyme. Therefore, we recommend using inhibitor-tolerant enzymes as well as the implementation of inhibition controls.

| Extension of reference databases and reanalysis of metabarcoding data
The extension of the genetic reference database by Sanger sequencing and using published sequences (12 species in total) resulted in the detection of two additional species, namely C. sibirica and E. mariae.
C. sibirica was detected in the Volga River, and it is also known to occur in the Tudovka River but this could not be verified in the present study. E. mariae yielded detection signals in the Moksha River, where the occurrence of E. mariae has recently been reported near the sampling site (Artaev & Ruchin, 2017).
In addition to the newly detected taxa, we observed reassignment of OTUs based on the extended reference database, demonstrating importance of exhaustive reference sequences, particularly when a 'best match' approach like BLAST is used for taxonomic assignment.
Not only does the lack of individual reference sequences result in potential false negatives, but they may also lead to erroneous detection of closely related sister-taxa.

| Primer bias and primer efficiencies of the metabarcoding primers
The observed primer efficiencies varied considerably among (a) the three metabarcoding primer pairs used, (b) among the target species, and (c) in some cases between forward and reverse primers, all leading to potential bias and associated problems in metabarcoding analyses. each researcher should carry out critical primer efficiency analysis on the target taxa before carrying out a metabarcoding study.

| CONCLUSIONS
The methodology and protocols around eDNA and metabarcoding surveys are still under development, and therefore it is imperative to F I G U R E 2 Mean primer efficiencies for the three metabarcoding primer sets analysed (12S rRNA, 16S rRNA, and CytB) for each species. Shown are the mean values for each species separated by forward (green triangle pointing upwards) and reverse (blue triangle pointing downwards) primers. Error bars represent standard deviations within species [Colour figure can be viewed at wileyonlinelibrary.com] recognize the pitfalls and challenges to overcome before fully integrating them into standard biomonitoring efforts. Here we reevaluated the results of our own metabarcoding assessment in the headwaters of the Volga River, as a case study. We showed that single-species TaqMan qPCR showed higher sensitivity than metabarcoding for selected rare species and highlighted the importance of complete reference databases to receive a full picture of the present species diversity. Finally, we showed that primer bias can affect the outcomes of metabarcoding analyses leading to potential false negatives. Thus, primer bias should be properly evaluated on each reference data set and primer pair.