Sampling sites located in Southern Norway (regions of Viken, Oslo, Vestfold/Telemark, Agder and Innlandet) were chosen randomly from sites with observations of amphibian species as reported in the Norwegian species database “Artsdatabanken”. Resampling for each month was done on the complete database, meaning that sites were sampled multiple times. As the density of sites with species observations is related to human presence, our sampling focused on anthropogenic impacted and potential Bd positive sites. Introduction of Bd is proposed to be related to human activities (Nielsen et al. 2019). During 2019, 110 bodies of water, including forest-, agricultural- and urban ponds, were sampled. To test for the optimal sampling season, we collected and analyzed samples from three different months in 2019: May or early June (37 samples), July (36 samples) and September (39 samples). These sampling times correspond to life history stages such as spawning, post-spawning and juveniles, respectively (See Suppl. material 1: Table S1). To compare the detection probabilities between the traditional and eDNA-based sampling, 91 bodies including forest, agricultural and urban ponds, and 10 swabs from various amphibian specimens were collected over a single season in 2020. Swabs were taken from frog and salamander skins at Hasseldalen (1 swab), Værmyr (3 swabs), Ringveien 52 (3 swabs) sites, while one swab from captured tadpoles was collected from each of Butjenna, Øyenkilen, and Sodra Skauen Gård waterbodies (See Suppl. material 2: Table S2). Swabs and tadpoles were preserved in 70% ethanol on site and stored at 4 °C until further analysis. Details about the survey design including geographical and temporal distribution of the samples are given in Suppl. material 3: Fig. S1.

Water samples were collected with a beaker attached to an expandable sampling pole (a device that can reach 4 meters). At least 5 samples (of 20–100 ml) per site were pooled together. We used cartridge filters (Sterivex filter units 0.22 μm pore diameter, Millipore, Merk, Darmstadt, Germany) and sterile 50-ml syringes to filter between 75 (minimum) and 1440 ml (maximum) of pooled water samples from each site (median volume of 480 ml). Sampling was performed by the same two persons throughout the study. A syringe - cartridge filter setup allowed for reproducible sampling without the need for heavy equipment and access to electricity. Filters were immediately frozen in dry shippers cooled with liquid nitrogen and then stored at -80 °C in the laboratory until further analysis. Chemical and physical parameters such as temperature, conductivity, pH, and turbidity were also measured on-site (See Suppl. materials 1 and 2).

All lab benches and equipment were cleaned by using 5% sodium hypochlorite and 70% ethanol in the field and laboratory. This included the three laboratories used for DNA extraction, prePCR or postPCR. Here pipettes and consumables were subjected to UV exposure for 30 min prior to usage. For water samples, DNA was extracted from the Sterivex filters using a Qiagen DNeasy PowerWater Sterivex Kit (Qiagen, Germany) including one unused filter as a negative control in one of the DNA extraction runs. The collected swabs were removed from ethanol and left to dry in clean Eppendorf tubes for 2 hours at 50 °C to evaporate residual ethanol. The DNA extraction was subsequently performed using the Qiagen Blood & Tissue DNA extraction kit. The quantity and quality (280/260 ratio) of the extracted DNA were assessed using NanoDrop.

Bd qPCR assays: We used two TaqMan assays amplifying different regions of the ribosomal RNA operon (ITS1 – 5.8S rRNA gene – ITS2). First the Boyle TaqMan assay (Boyle et al. 2004), targeting the Internal Transcribed Spacer 1 (ITS1) region where qPCR mixture contained 10 μL of 2× BioRad SsoAdvanced Universal Probes Supermix, 0.8 μM of the forward primer ITS1-3 Chytr (5′- CCTTGATATAATACAGTGTGCCATATGTC-3′) and the reverse primer 5.8S Chytr (5′- AGCCAAGAGATCCGTTGTCAAA-3′) and 0.24 μM of Chytr MGB2 (FAM-5′ TTCGGGACGACCC-3′-NFQ- MGB), 1 μL of internal positive contro (IPC, if not co-extracted with the sample), 1 μL IPC primer and probe mixture, 5 μL of environmental DNA sample and RNase/DNase free water to reach 25 μL as final volume. The second assay was based on the commercially available product Batrachochytrium dendrobatidis 5.8S rRNA Advanced Genesig Kit (Primer design Ltd, UK) which was carried out according to the manufacturer’s instructions. Both qPCR assays were conducted in a BioRad CFX96 thermocycler (USA) using the following program: 2 min at 95 °C, followed by 50 cycles of 10 s at 95 °C and 1 min at 60 °C. At least 8 serial standard dilutions were tested and validated before choosing five standard dilutions (1000, 100, 10, 1 and 0.1) to perform qPCR runs with environmental samples in triplicates. Bd -positive controls were added in the post-PCR room while the qPCR master mix was prepared in pre-PCR room.

For obtaining the standard curves needed to quantify Bd in environmental samples with the Boyle TaqMan assay, we used a DNA extract from Bd spores provided by Professor Anssi Laurila’s research group at Uppsala University. This extract had an expected concentration of 100 genomes (or genome equivalent) per µL. The number of rRNA operons in each spore was not known, therefore concentration was given as genome equivalent. Samples collected in 2019 were analyzed with only the commercial assay while samples collected in 2020 were analyzed using both assays in triplicates including two negative controls, in 96 well plates. Samples were declared Bd positive if two or more of the three qPCR runs generated a positive Bd signal within the range of serial dilutions of the standard. If only one of the three qPCR runs was positive, an additional fourth qPCR run was conducted. Sanger sequencing was performed to confirm the amplification of the target fragment for samples with Ct values > 39.5. To compare our results with known occurrences of Bd, we included data from a previous Norwegian survey in 2017 (Taugbøl et al. 2021).

In addition, to further assess the sensitivity of this assay and estimate the target copy numbers (ITS1 copies), we constructed an extra standard curve using a synthetic gBlocks double-strain ITS1 fragment (“Bd_26-271”; this name refers to the positions in the reference sequence NR_119535; 246 bp) based on eight decimal serial dilutions from 0.39 to 3.9 ×10⁶ ITS1 copies per reaction. The limit of detection (LOD) was estimated using a discrete approach and defined by the highest dilution where at least a single replicate of the triplicated standard dilution of 5 serial concentrations showed amplification in the respective assay consistently throughout all qPCR experiments with at least 95% confidence (modified from Kylmus et al. 2019). To define the limit of quantification (LOQ), we identified the highest dilution where a R² > 0.98 and an efficiency between 85 and 115% was obtained when fitting a linear regression to log-transformed Bd quantity and Cq values of the dilution series for each individual experimental (qPCR) run.

Internal positive control: In order to determine the DNA extraction efficiency and the presence of inhibitors, a separate internal extraction control DNA, primers and probe mixture with different target region, IPC, was provided with the Batrachochytrium dendrobatidis 5.8S rRNA Advanced Genesig Primer design kit (United Kingdom). In 2019, the batch of samples collected in September were co-extracted with IPC by mixing 4 μL of IPC with Qiagen lysis buffer (as described in the manufacturer’s instructions) and tested in a qPCR run. In 2020, 30 randomly selected DNA samples were co-extracted with IPC while the remaining extracted DNA samples were spiked with the IPC. To spike DNA with the internal control, the internal control was diluted 1:20 and 1 μL was added to a subset of samples which did not undergo internal control co-extraction, prior to the qPCR. The internal control is detected through the VIC fluorophore channel. We assumed a 100% extraction efficiency if Cq value of 27 was observed with Cq values of 27±2 within the normal range. Samples with Cq values of 30 or higher were considered to have PCR inhibitors.

Mock communities are now widely used in metabarcoding as a positive control and to validate analysis procedures from library preparation to bioinformatics analyses. A mock community of five amphibian species was prepared from samples provided by the lab of Anssi Laurila (Uppsala University) and Nordens Ark, including DNA extracts of Rana arvalis and Bufo bufo, and swabs provided by Anssi Lab (preserved in 95% ethanol) of Bufotes variabilis, Epidalea calamita, and Pelophylax lessonae. DNA extraction from the swabs was conducted as described above. Quantification of the DNA concentration was performed using the PicoGreen assay (ThermoFisher, US). Approximately 10 ng of the extracted DNA from each species were mixed in equal amounts and diluted to prepare five dilutions: 1.0, 0.5, 0.1, 0.01 and 0.001 ng/μL.

Paired-end sequencing on the Illumina MiSeq platform was based on two steps of PCR. The first step amplified the mitochondrial 12S rRNA gene of amphibian species using the following primers: forward batra-F 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNACACCGCCCGTCACCCT-3′ and reverse batra-R 5′-AGACGTGTGCTCTTCCGATCTNNNNNNGTAYACTTACCATGTTACGACTT-3′. Human blocking primer: 5′-TCACCCTCCTCAAGTATACTTCAAAGGCA-SPC3I-3′ was used to bind to human DNA and prevent its amplification (Valentini et al. 2016). This first PCR was performed in a final volume of 25 μL containing 5 μL of 5× Q5 reaction buffer, 0.2 μM “batra” primers, 0.2 mM dNTPs, 4 μM human blocking primer and 0.02 U/μL Q5 High-Fidelity DNA Polymerase (New England Biolabs) as well as 1 μL of positive control (mock community) or 5 μL of eDNA extract as a template. Amplification was performed in a Veritipro Thermal Cycler PCR Machine (Applied Biosystems, USA) under the following conditions: 30 s at 98 °C, followed by 35 cycles of 30 s at 98 °C, 30 s at 57 °C and 1 min at 72 °C and a final extension step at 72 °C for 7 min (Valentini et al. 2016). Each sample was amplified in triplicate and then pooled before the purification step. We also used two negative PCR products to monitor contamination and the mock community as a positive control of the amplification in each PCR run. 1% agarose gel electrophoresis was run to check the successful amplification of the target band (170 bp) before purifying the PCR products using AMPure XP magnetic beads (Beckman Coulter, US).

Twenty forward and reverse index primers developed based on Sinclair et al. (2015) were used for Illumina sequencing. This second PCR was performed in a final volume of 20 μl containing 4 μl of 5× Q5 reaction buffer, 0.2 mM dNTPs, a combination of the indexed forward and reverse primers (0.25 μM each), 0.02 U/μl Q5 High-Fidelity DNA Polymerase (New England Biolabs) and 2 μL of purified first step PCR product. Amplification conditions were 30 s at 98 °C, followed by 15 cycles of 10 s at 98 °C, 30 s at 66 °C, and 30 sec at 72 °C, then a final extension step at 72 °C for 2 min. The indexed target band was checked using 1% agarose gel electrophoresis. The second PCR product was purified using the same reagent as above, and quantified using the PicoGreen assay. The detailed DNA metabarcoding protocol for amphibians was published on protocols.io (dx.doi.org/10.17504/protocols.io.973h9qn). The first and second PCR mixtures were prepared in pre-PCR laboratory while the first purified PCR product was added in the post-PCR laboratory.

Equal amounts of DNA from each amplified sample were mixed to prepare the sequencing library. In addition to sequencing the five dilutions of the mock community, two or three environmental samples were spiked with one of the mock community species which was processed in the same sequence pool to confirm the presence and absence of the identified species. Furthermore, one random negative PCR product (entire sample) was included in each Illumina sequence run performed in 2019 and 2020. To detect all amphibian species in positive swab samples, the three swabs collected from positive Bd sites were sequenced while one swab was only sequenced from the other sites. The final pooled samples were visualized by gel electrophoresis to ensure that no extra unwanted fragments existed. The amplicon library was then sequenced on an Illumina MiSeq machine using the MiSeq v2 PE 150 Micro protocol, generating approximately 20 million paired-end reads of 150 bp length and exact matches to batra primers. Raw sequence data are available at NCBI through the SRA accession no. PRJNA821076 for 2019 and PRJNA822096 for 2020.

Amphibian 12S rRNA gene sequences were obtained by searching NCBI and BOLD using species names of the Norwegian amphibian species and ‘12S’ as search criteria. Next, a homology search using Basic-Local-Alignment-Search-Tool (BLAST) was used to extract rRNA gene sequences of high (90%) homology. This two-step process led to a 12S rRNA gene sequence database of mostly amphibian sequences. To clean up the database, sequences were aligned and the “batra” primers were matched. This resulted in a database of 7 sequences including all species that have been reported in the wild in Norway.

Traditional observational data of amphibian diversity was obtained from the Norwegian public biodiversity databank “Artsdatabanken” (https://artskart.artsdatabanken.no/). Data for amphibians was downloaded on 2020-05-11 using “Amphibia” as a query resulting in 4086 species observations in the region overlapping with our samples. Species observations from the last 20 years were kept and then grouped if GPS coordinates overlapped as defined by resolution to the second decimal corresponding to a place up to 1.1 km².

Raw sequences were first processed with cutadapt v1.18 to remove PCR primers, and then analyzed with the R package dada2 v1.14.1 (Callahan et al. 2016) for denoising and sequence-pair assembly. After manual inspection of quality score plots, forward and reverse reads of the amplicon sequencing run were trimmed to 50 bp length. This assured an almost complete overlap of the forward and reverse reads as the amplicons without primers were 51–54 bp long. Additional quality-filtering steps removed any sequences with unassigned base pairs and reads with a single Phred score below 20. After reads were dereplicated, forward and reverse error models were created in dada2 with a subset of the sequences (10⁸ nucleotides). Amphibian 12S rRNA gene amplicons were assembled by merging the read pairs. Chimeras were removed using ‘removeBimeraDenovo’ from the dada2 package, which resulted in the final Amplicon Sequencing Variant (ASV) table for performing taxonomic assignments. Taxonomy (to species level) was assigned using blastn and an in-house 12S rRNA gene database. Sequences were assigned to specific species using cutoff criteria of identity > 96% (which corresponds to approximately two mismatches) and alignment length > 45 bp.

Heatmaps were constructed using the ggplot2 R package, plotting the number of samples uniquely detected by either eDNA or traditional observation of species deposited to the Norwegian species database “Artsdatabanken” (ABD method), as well as those samples where the species was detected using both methods. Correspondence of species observations between eDNA and those recorded in Artsdatabanken (using GPS data rounded to the second decimal) were evaluated by correlation analysis, fisher exact test and Pearson’s Chi-square test using R (version 3.6.2). In order to trace the outcome of the sequence analysis we used the sequenced amphibian mock community.

Generalized linear (glm; base R) and additive (gam; “mgcv” package in R) models were used to explore detection rates of amphibian species in relation to time of sampling and pH. We also performed co-occurrence analysis using the R package “cooccur”. The latter did not reveal any significant relationships, probably due to sparse Bd detection.

We tracked the distribution of the invasive chytrid fungus Bd in the sampled water bodies. By using two TaqMan assays we were able to detect Bd in six water bodies and could confirm its presence in most cases by repeated sampling (Table 1). Still, our results indicate that spring is the season of highest detection probability for Bd, as samples taken during other seasons were all negative.

We used two TaqMan assays to amplify different regions of the rRNA operon (ITS1 and 5.8S rRNA gene). qPCR efficiency of the Bd 5.8S rRNA Genesig Advanced Kit and the Boyle TaqMan assay was on average 98.82% and 99.34%, respectively. The limit of quantification (LOQ) of the assays depending on the run varied between 100 and 10 gene copies for the commercial Genesig kit and 1 and 0.1 genome equivalents for the Boyle TaqMan assay (Fig. 1; Suppl. material 3: Table S3). The limit of detection (LOD), as assessed from at least three independent runs, varied between 10, 1 and 0.1 gene copies for the commercial Genesig kit and 1, 0.1 and 0.01 genome equivalents for the Boyle TaqMan assay (see Suppl. material 3: Table S3). Thus it can be speculated that the Boyle TaqMan assay is slightly more sensitive than the Genesig kit assay. The higher sensitivity of the Boyle TaqMan assay was also confirmed using dilutions of the gBlocks fragment Bd_27-271 as standards, as the lowest concentration detected (in 2 of 3 triplicates) corresponds to 0.39 ITS1 copies (Fig. 1B), while the LOQ (based on a single qPCR experiment) seems to be two orders of magnitude higher (39 ITS1 copies). It is well known that rRNA operon copy number per genome (or spore) is highly variable between strains of the same fungal species, as reported by Longo et al. (2013) for Bd genomes which contain from 10 to 144 ITS1 copies.

Figure 1. 

Examples of standard curves used to determine the limits of detection (LOD) and quantification (LOQ) for the Boyle TaqMan assay using spore DNA (providing genome equivalents) (A) and the synthetic gBlocks fragment Bd_27-271 (providing ITS copy numbers) (B) as standards. Crosses represent standard dilutions above the LOD but below the LOQ while circles represent standard dilutions used for determining the corresponding standard curves and their R² and efficiency values (detailed on the panels). For the spore DNA curve (A), detailed statistics on multiple runs are provided in Suppl. material 3: Table S3.

Both qPCR assays generated positive Bd signals from water samples of three locations: Østre Støkken, Ringveien 52, and Odden, and from one swab from the location Værmyr (Table 1). Two additional locations, Butjenna and Rokkevannet, generated positive Bd signals when using the Boyle TaqMan assay. However, as a high Cq value (as in the case of Butjenna and Rokkevannet) can represent unspecific amplification and thus false positive Bd results, both PCR products were purified and Sanger-sequenced for validation. The resulting sequences were used to query the NCBI nt-database using blastn, which resulted in 100% matches with sequences from Bd, and thus the presence of Bd was concluded.

We used internal positive controls (IPC) consisting of a fully synthetic DNA, primers and a TaqMan probe, which are included in qPCR reactions or co-extracted with the DNA extract to detect PCR inhibition and false negative results due to inhibitory substances in the environmental samples (Hoorfar et al. 2004; Phillips 2004). In 2019 and 2020, five and six samples, respectively, showed delayed Cq values above 27±2 for the IPC. This variability in Cq results strongly suggests the presence of inhibitory compounds in a minor fraction of samples, which needs to be taken into account when interpreting negative results (Suppl. material 3: Fig. S2). The IPC was detected in all co-extracted samples indicating the efficiency of the extraction protocol used. Three of these samples generated a normal IPC-qPCR signal after dilution, which is one way to overcome the effect of inhibitors. It should be mentioned that only half of these samples were successfully sequenced and one amphibian species was identified in each of them, thus metabarcoding results from these samples are likely biased as well. By reporting IPC-qPCR signal, potential inhibition can be observed providing a good indicator of false-negative PCR results. A way to solve some false-negative detection might be by sample dilution or optimizing qPCR conditions (Kamal et al. 2017; Lance and Guan 2019) or the use of an inhibitor removal kit (McKee et al. 2015).

Amphibian diversity was assessed using a metabarcoding approach based on amplification of the mitochondrial 12S rRNA gene of amphibian species and subsequent sequencing. Positive controls such as the mock community including DNA from R. arvalis, B. bufo, B. variabilis, E. calamita, and P. lessonae showed amplification down to 10^-2 ng of template DNA and resulted in sequences for most mock species (Fig. 2). The exception was B. variabilis that could not be detected when the targeted DNA was below 10^-1 ng of template DNA. Most of the species were detected in low concentration indicating that the metabarcoding approach is still sensitive to very low amounts of genomic DNA.

Figure 2. 

Sequencing results from the dilution series of an amphibian mock community in 2019 (A) and 2020 (B). Both plots showing the five added species in the original, 5× and 10× diluted mock samples, while Bufotes variabilis disappeared in the 100× diluted sample.

Next, we evaluated species discrimination of native Norwegian amphibian species by the DNA metabarcoding approach. The comparison of the available reference sequences by pairwise alignment revealed sufficient divergence amongst native Norwegian species for annotation at species level. The closest species were R. arvalis and R. temporaria with a minimum of two nucleotide differences in the amplified 12S rRNA gene region while all other species differed in at least four positions. Species discrimination may, however, become an issue with this metabarcoding approach in areas of high amphibian diversity including multiple closely related species, such as in the tropics.

Analyses of the 110 and 91 samples from 2019 and 2020, respectively, resulted in all of the samples being well amplified as reflected by the amplified amplicon size (approx. 50 bp, expected to the target region). No amplification was obtained from negative control filters, thus allowing us to exclude the possibility of cross-contamination among filters and confirming the clean DNA extraction and purification steps. Moreover, there was no amphibian species detected in any negative PCR control of the metabarcoding after the end-point (no resulting sequences), even if a band appeared in gel electrophoresis due to the primer dimer formation, confirming a clean PCR setup.

A stringent raw sequence data quality filtering resulted in a loss of approximately 60% of the reads from an average of 15,736 raw paired end reads per sample (range from 830 to 87,543). From an average of 6,181 (range from 587 to 29,720) quality-filtered paired end reads per sample, denoising and merging resulted in a mean number of 5,676 high quality sequences (range from 567 to 29,136). Amphibian sequences were only detected in 68 of 110 samples collected in 2019, and in 69 of 91 samples from 2020. Unspecific amplification of non-targeted species could be observed in almost all amplified samples. In 2019 and 2020 the mean percentages of sequences per sample assigned to amphibian taxa were 9% (range 0 to 96%) and 30% (range 0 to 100%), respectively, corresponding to 947 (range 0 to 26,327) and 1,862 (range 0 to 12,667) amphibian reads per sample. Closer inspection of the non-amphibian sequences showed various taxonomic assignments from bacteria to fish. This indicates unspecific amplification in the PCRs and emphasizes the need for a high sequencing depth to detect amphibian species when using the batra primers. Still, as we obtained on average more than 5,000 paired reads per sample, false negatives as a result of sequence library preparation could be minimized.

Looking at the 2019 data, amphibian detection was significantly lower in September compared to May/June and July (estimate = -1.635, z-value = -3.710, p-value = 0.0021 as given by generalized linear models). This was corroborated by the total number of detected amphibian occurrences among all sites (Table 2) and at six sites that were sampled during the three sampling occasions in 2019. In autumn, amphibians could not be detected at any of the six sites while in May/June and July there were multiple species detections. Thus it can be concluded that sampling during spawning and post-spawning results in higher likelihood of detection when compared to autumn sampling. In a previous study, a comparison of the detection of the pool frog (P. lessonae) using both eDNA and traditional methods revealed that traditional methods gave a higher rate of observation in June, whereas eDNA gave at least as many or more observations during other months (Eiler et al. 2018).

Five amphibian species B. bufo, L. vulgaris, R. arvalis, R. temporaria, and T. cristatus are prevalent in Norwegian waterbodies. This limited number of species was observed by our eDNA approach with more than one amphibian species being detected in several samples. A study on the influence of pH on amphibian diversity in Norway reported that R. temporaria was detected in all pH levels but there is a decrease in the frequency of B. bufo and T. vulgaris (Dolmen 1988), while R. arvalis and T. cristatus seemed to be sensitive to pH changes. Using our eDNA data and generalized additive models revealed a significant decrease in both R. arvalis (Estimate = -1.388, z-value = -2.958, p-value = 0.0031) and T. cristatus (Estimate = -3.566, z-value = -3.059, p-value = 0.0023), stronger than for L. vulgaris (Estimate = -0.979, z-value = -2.186, p-value = 0.0288) and R. temporaria (Estimate = -0.789, z-value = -2.013, p-value = 0.0441). Furthermore, our results confirm previous conclusions that acidic areas are lower in amphibian diversity in the inland/highland region of Southern Norway as given by generalized additive models revealing a significant decrease in overall species observations with decreasing pH (Estimate = -3.350, z-value = -3.644, p-value = 0.0003).

Five amphibian species were detected by both traditional observation methods (ABD) and eDNA (Fig. 3). A much larger number of sites has been covered by ABD when compared to eDNA methods. Still, there were 133 sites where overlapping data was available for comparative analysis. An all data heatmap is shown in Fig. 3. Fisher exact and Pearson’s Chi-squared tests indicate a significant association between eDNA and ABD in detecting four of the amphibian species while no significant association could be observed for R. temporaria (Table 3). To conclude, correspondence of species observations between data from a national species reporting ABD-based program and eDNA metabarcoding was significant but weak. This may reflect both an incomplete overlap in species detection and differences in the number of total detections of individual species, and the different timeframes (20 years vs 2 years).

Figure 3. 

Comparison of amphibian species detection between our metabarcoding results (2 years of data) and species inventories in the Norwegian species database “Artsdatabanken” (20 years of data). ND – number of sites where the respective species was not detected by any of either method; shared – the number of sites where the respective species was detected by both methods; unique eDNA – the number of sites where the respective species was detected only by eDNA; unique ABD - the number of sites where the respective species was found only in the database. Detailed statistics on correspondence between methods are given in Table 3.

Screening amphibian occurrence by field techniques usually requires a suite of tools, which introduce several sources of errors. One of the errors is the high variability of amphibian detection in different seasons throughout the year (Heyer et al. 1993) and another is error in the sampling technique. For example, different types of amphibian traps affect detection probabilities (Dodd 2010; Petitot et al. 2014). This has been proposed by previous observations such as in neotropical streams in Brazil where both ABD and eDNA metabarcoding methods were used (Sasso et al. 2017). Our results confirm previous observation that amphibians should be monitored in late spring/early summer in Scandinavia when using eDNA methods (Eiler et al. 2018). Thus, national monitoring programs for amphibians based on eDNA are best performed from late May to early July.

The high dispersion and low detection probability represents a challenge for monitoring programs in risk assessment, as large areas including high numbers of systems need to be monitored. There is also a challenge to integrate the monitoring results in management strategies to protect amphibian diversity. To map large areas, such as across Norway, a random sampling strategy should be set up. In our case for the risk assessment of Bd, we focused on systems impacted by humans since the main pathway of introduction of Bd is related to human activities (VKM report). Since sampling of a vast amount of systems is necessary, efficient sampling techniques with low likelihood of contamination need to be applied. We chose a low-tech approach using a sterile beaker attached to an expandable (easy to clean) sampling pole, sterile syringes and cartridge filters. This kept the time spent per site to a minimum (20–45 minutes) including the assessment of water parameters with online sondes.

For sample preservation, we recommend flash freezing of filters in liquid nitrogen on-site using a dry shipper and transfer to -80 °C freezers in the lab as this allows the long-term preservation of all kinds of environmental and tissue samples. This has been the “gold standard” for downstream DNA and RNA analysis for at least three decades (Seutin et al. 1991; Reiss et al1995; Kilpatrick 2002; Anchordoquy and Molina 2007; Frampton and Sam 2008; Wong et al. 2012). For DNA extraction, we recommend commercial kits that are optimized for filter cartridges and remove potential PCR inhibitors while still keeping DNA yield high, i.e. the DNeasy PowerWater Sterivex kit. This is an important initial step of the analysis workflow that needs to be kept constant as extraction protocols determine to a large extent the outcome of metabarcoding studies (Majaneva et al. 2018). For library preparation in the case of metabarcoding we strongly argue for a two-step PCR. There are several reasons: 1) Short adapters with undefined nucleotides between target primer and sequencing primer minimize biased PCR amplification when compared to a 1-step (fusion) PCR adding adapters and index, 2) a two-step PCR minimizes PCR cycles in the first step limiting PCR primer bias as cycle number can be kept low (also in the case of using a ligation based library preparation) (Thompson et al. 2002). 3) Indexing by PCR, which adds adapter sequence, minimizes index jumping when compared to ligation based library preparation. 4) Compared to a ligation-based protocol the two-step PCR protocol is cheaper.

Bioinformatic workflows for national monitoring programs should denoise reads based on quality scores and optimize for each sequencing run, for example as implemented in the dada2 algorithm (Callahan et al. 2016) instead of clustering to improve species resolution. Using such methods allows discrimination of species with two nucleotide mismatches. Reference databases should be developed specifically for the national monitoring project taken from national species inventories including potential invasive species. This was rather simple in our case, with only seven amphibian species being present in Norway, but might be more challenging in species-rich taxonomic groups or locations. Species resolutions of the metabarcoding assay can be tested in silico using for example database-against-database homology searches and lowest common ancestor analysis (Huson et al. 2007).

In the case of qPCR, we recommend the use of IPC to identify samples that exhibit PCR inhibition. These samples can be flagged as problematic and run again after dilution which can lead to positive detection of targets. There is also a need to decrease LOD and LOQ, keeping amplification efficiency between 90–110%. The lowest number of replicate PCR reactions should be three, but this can be increase if allowed by budget constraints (Piggott et al. 2016). In addition, high Ct samples should always be evaluated by sequencing to detect false positives.

Besides these recommendations there are still many obstacles standing in the way of the implementation of eDNA-based biodiversity monitoring. International coordination is paramount to ensure comparability of data in time and space and to avoid unnecessary duplication of efforts and resulting delays in wide application. Awareness and knowledge of the new methods among practitioners is essential for implementation and to provide platforms for critical discussion on when, where and how the eDNA-based methods should be applied. Standardization efforts need to advance and cover all steps from sampling design to sequence data interpretation. Additionally, in order to allow for comprehensive mapping of genetic diversity, the number of published genomes (at least mitochondrial genomes) should cover all known species. Modelling and analysis tools should be applied and developed in tandem to facilitate efficient decision making.

Our study shows that eDNA surveys provide detailed insights into the distribution of amphibians and an associated pathogen, thus allowing the assessment of risks associated with the invasive chytrid fungus Bd. Reproducible detection of Bd over multiple years strongly suggests that this pathogen is established at multiple sites in Southern Norway. Bd positive samples were obtained throughout an area of at least 1000 km² south of Oslo with Bd detection highly dispersed within the area, which is comparable to similar climatic regions in Sweden and the UK (Rosquist 2020).

Multiple methods have been proposed to assess risks for biodiversity including Species Distribution Models (SDMs) combined with species specific biotic indices (Rödder et al. 2009), expert opinion (Non-native Risk Assessment scheme; Roy et al. 2014) and an integrative analytical approach combining ecological traits and anthropogenic risk factors (Munstermann et al. 2021). Bd is known to be a particularly temperature and moisture dependent species (Lips et al. 2006, Woodhams et al. 2008). The in vitro growth optimum range of Bd is at 17–25 °C, whereas temperatures higher than 29 °C, freezing and desiccation are lethal (Piotrowski et al. 2004). These findings are supported by observations in the field (Kriger et al. 2007). A previous study on the geographic extent of Bd’s climatic niche based on SDM suggested the presence but low risk of Bd for anuran amphibian species in Norway (Rödder et al. 2009). In another recent study it was argued that transmission from an environmental Bd reservoir and climate change could increase the ability of Bd to invade new amphibian populations and increase their extinction risk. Still, a recent study concludes that Bd-induced extinction dynamics were far more sensitive to host resistance and tolerance than to Bd transmission (Wilber et al. 2017). Infection-tolerant host species seem to dominate in Norway, as suggested by previous European studies (Smith 2014). Hence, the overall likelihood of a negative impact of Bd on indigenous amphibian species has been suggested (with moderate confidence) to be minimal, while for an introduced frog species (P. lessonae) there was a moderate risk identified (Nielsen et al. 2019. Considering previous outbreaks of chytridiomycosis in Europe (Bosch and Martínez-Solano 2006) and extensive studies on Bd in similar climate regions, as well as its low prevalence in our study coupled with co-occurrences of various indigenous species over multiple years, it can be assumed that chytridomycosis outbreaks in Norway are unlikely. Still, an annual follow-up of Bd presence at least for the positive sites and nearby areas may be advisable, to keep an eye on the spread of this fungal pathogen.