Primer design

The 16S rRNA gene has been suggested as a suitable DNA barcoding marker for amphibians (Vences et al. 2005). To design universal primers for amphibians’ metabarcoding analyses, a total of 1,034 sequences of the mitochondrial 16S rRNA region of amphibians from 58 (/75: 77.3%) families and 293 (/564: 52.0%) genera were downloaded from the National Center for Biotechnology Information database (NCBI, https://www.ncbi.nlm.nih.gov). These sequences were all amphibian sequences that were available from the NCBI, which included data from amphibians worldwide. To extract the taxonomic information of the 16S rRNA sequences from NCBI text (GenBank format), the programme Taxonparser was used (newly written for this study; available from https://github.com/RyosukeKakehashi/taxonparser). The 16S rRNA sequences were aligned using MAFFT (Katoh and Standley 2013) with the default parameters of the MAFFT server (https://mafft.cbrc.jp/alignment/software/). Based on the alignments, several candidate universal primer sets were designed from the conserved regions (Fig. 1; Table 1), using Primer3web (Untergasser et al. 2012) and visual observation.

In silico and in vitro evaluations of designed primer sets

To characterise each primer set, we performed in silico tests of the following three parameters: 1) universality of each priming site, 2) specificity of the priming site for target taxa and 3) resolution of the internal amplified regions. To examine the universality of each primer, the frequency of bases at each locus of the primers was visualised using the sequence logo from the aforementioned 1,034 sequences, which were the same as those used for the primer design. The specificity of each primer set was confirmed using in silico PCR, which was performed using the “search_pcr” command implemented in USEARCH v.10.0.240 (Edgar 2010), targeting the full-length mitochondrial DNA reference database of vertebrates. This database was created from the mtDNA data deposited in NCBI (see Supporting Information) and comprises 26,096 sequences (fish: 7861, amphibia: 959, reptile: 1170, bird: 2988, mammalian: 13,118). For each primer set, in silico PCR was performed with minimum and maximum amplicon sizes of 100 and 1,000, respectively and with varying maximum acceptable primer mismatches of 0, 1 and 2 (“maxdiff” option). The dissimilarity of the internal amplified regions of each primer set was calculated using the “dist.dna” function of the ape package in R ver. 3.6.3 (R Development Core Team 2019). The average length and dissimilarity of the internal amplified regions were then used to calculate the expected number of different bases amongst species in order to evaluate the taxonomic resolution of the primer sets. These numbers were compared amongst the primer sets by ANOVA and post-hoc Tukey-Kramer test.

In vitro tests were performed using tissue-derived DNA from 16 amphibians, including 13 species from the Anura clade, one from the Caudata clade and two from the Gymnophiona clade (Suppl. material 1: Table S1). The total DNA, used as PCR templates, was extracted from previous studies and newly obtained in the present study (Kurabayashi et al. 2006, 2008, Kurabayashi and Sumida 2009, 2013). Five newly designed and existing primer sets were used, which were chosen, based on the Tm value and presence of mixed bases (Table 1). The PCRs were conducted in a 10 µl solution containing 5 ng template DNA, 10× TaKaRa Ex Taq Buffer, 0.25 U TaKaRa Ex Taq DNA polymerase, 0.25 mM each dNTP and 0.5 µM forward and reverse primers (each). PCR conditions were as follows: initial denaturation at 94 °C for 30 s; 35 cycles of denaturation at 95 °C for 30 s, annealing at 52 °C for 30 s and extension at 72 °C for 30 s; and final extension at 70 °C for 5 min. Electrophoresis was used on the resultant PCR solutions (5 μl) on a 1.5% agarose gel to confirm amplification of the target DNA.

Field surveys

From July to September 2019, field surveys were conducted at 160 agricultural ecosystem sites across 10 areas of Japan (Fig. 2; Suppl. material 1: Table S2). In each area, the study sites included various types of agricultural water environments, including streams, irrigation waterways, and ponds. These study areas covered a wide geographic region of Japan and included multiple local populations of some amphibian species. Most amphibians found in agricultural ecosystems are anurans and these areas harbour approximately half of the Japanese anuran species (19/43 species) (Ministry of Agriculture, Forestry and Fisheries 2009; https://www.maff.go.jp/j/nousin/keityo/tanbo/). In 122 of the 160 sites, both eDNA and physical surveys (capturing and visual surveys) were conducted. For the remaining 38 sites, only the eDNA survey was performed.

Water sampling was performed according to the Environmental DNA Sampling and Experiment Manual ver. 2.1 (The eDNA Society, 2019). We sampled 1 litre of surface water at each site and, after sampling, we added 1 ml benzalkonium chloride (final concentration = 0.1%) to prevent eDNA degradation (Yamanaka et al. 2017). To monitor potential contamination during the filtration and eDNA extraction processes, 1 litre of distilled water was used as the negative control for each sampling date.

We performed physical surveys of amphibian species at the 122 sampling sites (Suppl. material 1: Table S1) and compared the amphibian detectability between the physical and eDNA metabarcoding surveys. Water sampling for eDNA metabarcoding was performed prior to the physical surveys. Physical surveys were conducted with constant capture efforts; specifically, for each site, three people collected amphibians for 30 min in a 50 m section. The collections were performed using either a fixed net, hand net, net casting or boxed net. To ensure safety when conducting the surveys at the water’s edge, including physical surveys, water sampling and measuring water parameters, surveys were performed during daylight hours (07:00–16:00 h). The surveys were conducted in agricultural water environments, which included streams, irrigation waterways and ponds in each area.

Environmental DNA sample processing

Filtration and eDNA extraction were performed according to the Environmental DNA Sampling and Experiment Manual ver. 2.1 (The eDNA Society 2019). We filtered water samples using a glass fibre filter with nominal pore size of 0.7 μm (GF/F; GE Healthcare Life Science). However, two filters were used for processing samples with high turbidity to avoid the potential filter clogging, as 39 samples being processed this way. Filters were preserved individually (but paired filters from turbid samples were pooled together) in tubes at −30 °C. The eDNA on the filter was extracted using a Salivette (Sarstedt) and DNeasy Blood & Tissue Kit (QIAGEN Science, Hilden, Germany). Salivette was used to centrifuge the GF/F filter to elute the DNA solution from the filter. The eluate was then purified using a DNeasy Blood and Tissue kit. A final elution volume of 100 µl of DNA was obtained and then stored at –25 °C. To prevent cross-contamination, all equipment used in the water collection and filtration steps, including plastic bottles, filter funnels and tweezers, were decontaminated using > 0.1% sodium hypochlorite solution (The eDNA Society 2019).

To detect amphibian species from environmental samples, we amplified the partial 16S rRNA gene of amphibians and then performed high-throughput sequencing using a MiSeq platform (Illumina, San Diego, CA, USA). Amongst the primer sets, Amph16S, – consisting of 16S rRNA gene specific primers (Amph16S_1070_F + Amph16S_1340_R) – was used to amplify the gene fragments from the eDNA samples for the following reasons: 1) this primer combination successfully amplified the target gene fragment of all 16 species tested belonging to all three amphibian orders (Table 3) and 2) the amplified region of this primer set was estimated to have the highest resolution in species identification amongst the primer sets tested (see Results). The first-round PCR (1^st PCR) was performed in eight replicates with KOD-Plus-Neo polymerase (Toyobo, Osaka, Japan). Each PCR reaction (25 µl final volume) contained 300 nM primers, 2.5 µl 10× KOD buffer, 2.5 µl dNTPs (2 mM), 1.5 µl MgSO₄ (25 mM), 0.5 µl KOD-Plus-Neo polymerase, 1 µl DNA template and ultrapure water. The thermal cycle profile was 3 min at 95 °C; 40 cycles of 20 s at 98 °C, 15 s at 58 °C and 15 s at 72 °C; and 72 °C for 5 min. We used MiSeq-R1 tailed forward (5ʹ-ACACTCTTTC- CCTACACGACGCTCTTCCGATCTNNNNNN + gene-specific sequences (Amph16S_1070_F) – 3ʹ) and R2 tailed reverse (5ʹ-GTGACTGGAGTTCAGACGTGT- GCTCTTCCGATCTNNNNNN + gene-specific sequences (Amph16S_1340_R) – 3ʹ) primers. These primers included six random hexamers (N) to enhance cluster separation on the flow cells during initial base call calibrations on the MiSeq platform. Ultrapure water was used instead of eDNA in the eight reaction mixtures (non-template negative controls). After the eight technical replicates were pooled into a single tube, we removed unreacted reagents and primer dimers from the 1^st PCR products using the GeneRead Size Selection Kit (QIAGEN Science, Hilden, Germany) according to the Environmental DNA Sampling and Experiment Manual (Version 2.1; The eDNA Society 2019).

We carried out a second-round PCR (2^nd PCR) using the purified products from the 1^st PCR as templates. The 2^nd PCR was performed using P5-i5-R1 (5ʹ – AATGATACGGCGACCACCGAGATCTACAXXXXXXXXACACTCTTTCCCTACACGACGCTCTTCCGATCT – 3ʹ) and P7-i7-R2 (5′ – CAAGCAGAAGACGGCATACGAGATXXXXXXXXGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT – 3′) primers to add MiSeq adapter sequences and 8-bp index sequences to both ends of the amplicons. The octo-X segments represent dual-index sequences. The 2^nd PCR was conducted in a 12 μl solution containing 6.0 μl of 2× KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, WA, USA), 3.6 pmol each primer and 1 μl 1^st PCR product. The PCR profile was as follows: 3 min at 95 °C; 15 cycles of 20 s at 98 °C, 15 s at 72 °C and 72 °C for 5 min.

The 2^nd PCR products were pooled in equal volumes into a single 1.5-ml tube. Target amplicons were selected by electrophoresis using E-Gel SizeSelect 2% (Thermo Fisher Scientific, Waltham, MA, USA) with the E-Gel Precast Agarose Electrophoresis System (Thermo Fisher Scientific, Waltham, MA, USA). DNA concentration was measured using real-time PCR assays (QuantStudio3; Thermo Fisher Scientific, Waltham, MA, USA) and the library size distribution was confirmed using TapeStation 4200 (Agilent, Tokyo, Japan). The concentration of the DNA library was adjusted to 4 nM. Finally, the library was sequenced using an Illumina MiSeq v.2 Reagent kit for 2× 250 bp PE (Illumina, San Diego, CA, USA).

Raw sequencing reads were converted to FASTQ format using Illumina bcl2fastq2 v.2.17 software allowing zero mismatches. To perform the species identification from the MiSeq output, FASTQ data were processed using the pipeline of the metabarcoding programme package Claident version 0.2.2017.05.22 (Tanabe and Toju 2013; downloaded from https://www.claident.org/). We demultiplexed the data and removed any reads with low-quality index sequences (i.e. Phred score < 30) using the clsplitseq command with the option minqualtag = 30. Paired-end reads were then merged using the clconcatpair command with the option maxnmismatch = 20 and minovllen = 20. We then trimmed low-quality tails until the Phred scores of the last base were 30 or higher. We removed low-quality sequences (Phred score < 30) using the clfilterseq command with the options minqual = 30 and maxplowqual = 0.1. In addition, apparently noisy and chimeric reads were removed using the clcleanseqv command with the options primarymaxnmismatch = 0, secondarymaxnmismatch = 1 and pnoisycluster = 0.5. After quality control, we classified the reads with > 99% sequence similarity into a molecular operational taxonomic unit (MOTU) using the clclassseqv command with the option of minident = 0.99. We obtained denoising sequence data and identified each MOTU as a species using online BLASTN (https://blast.ncbi.nlm.nih.gov/Blast.cgi), based on a 97% homology criterion.

MOTUs with less than 10 reads per sample were discarded because of potential contamination. The remaining sequence reads assigned to amphibians were vetted, based on habitat and species assignments were finalised. For all samples, the read depth was sufficiently large to saturate the number of amphibian species detected (checked using the “rarecurve” function in the vegan package version 2.5–4 (Oksanen et al. 2019)). Therefore, rarefaction was not performed.

To compare the monitoring results between eDNA metabarcoding and physical surveys, the following analyses were performed using the vegan package version 2.5–4 and lme4 package version 1.1–21 (Bates et al. 2015) in R ver. 3.6.3 (R Development Core Team 2019). We used the VennDiagram package version 1.6.20 and ggplot2 package version 3.2.0.9000 to draw figures. The site-species read matrix was converted to the presence/absence of each species and the data of 122 sites, where both eDNA survey and physical surveys were carried out, were used for comparison (Suppl. material 1: Table S2). First, we compared the species composition of the monitoring methods. Non-metric multidimensional scaling (NMDS) with 10,000 permutations depicted the similarity of species composition amongst sites and methods measured using the Jaccard Similarity Coefficient and PERMANOVA, with 10,000 permutations, was conducted using the “adonis” function in vegan. In addition, to compare the beta dispersion of the methods, PERMDISP analysis was performed with the “betadisper” function in vegan. Second, we compared the number of amphibian species detected between the monitoring methods. We fitted a generalised linear mixed model with a Poisson distribution, using the function “glmer”. In this model, the number of species was used as the response variable. The monitoring methods (i.e. eDNA metabarcoding or physical surveys) were set as explanatory variables and the survey area ID was set as a random effect to consider differences in the number of species inhabiting the survey areas. Finally, we visualised the detected community composition between the two methods in each area by drawing a Venn diagram.

We tested five sets of universal primer candidates for the amphibian 16S rRNA genes (Table 1). The sets of “MiAmphiS” and “Amphi16S” had the shortest and longest amplicon sizes, respectively. The universality of each primer set was shown by the sequence logo (Suppl. material 1: Fig. S1). Both forward and reverse primers of “MiAmphiL” and “MiAmphiS” had no degenerated bases in both forward and reverse primers and high universality at the 3ʹ-end of both primers; however, there were several mismatches in the region between the 5ʹ-end and the centre of both primers, which was located on a variable region amongst species. The primers of “Modified 16Sar” and “16Sar_mod2” contained some degenerated bases, resulting in high universality. “Amph16S” had mixed bases in the forward primer and no mixed base in the reverse primer and also showed high universality.

Regarding the power of amplification, we performed in vitro amplification tests for all five primer sets using DNA templates extracted from the 16 amphibian species belonging to all three amphibian orders and all primer sets amplified the target fragments tested. We performed in silico PCR using vertebrate 16S rRNA gene data reported to date (7861 fishes, 959 amphibians, 1170 reptiles, 2988 birds and 13,118 mammals). The specificity of each primer set, examined by in silico PCR, is shown in Suppl. material 1: Table S3. When we did not allow any mismatch between the primers and target 16S rRNA gene sequences (maxdiff = 0), Amph16S showed the highest amplification rate (925/959 sequences; 96.4%). However, this primer set also had the highest amplification rate for non-amphibian taxa. MiAmphiL and MiAmphiS had the lowest amplification rates for non-amphibian taxa (3.1% and 5.4%, respectively), but they had the lowest amplification rates for amphibians. Even when we allowed two mismatches (maxdiff = 2), Amph16S showed the highest amplification rate for amphibians (98.7%). Under maxdiff = 2, the amplification rates for reptiles and birds were clearly lower for MiAmphiL and MiAmphiS (8.2% and 28.9% for reptiles and 0% and 4.9% for birds, respectively) than for Amph16S (90.3% and 99.4%, respectively), whereas the amplification rates for fishes and mammals did not differ significantly amongst these primer sets.

In the context of taxonomic resolution, Amph16S had the highest expected number of bases that differ amongst species (mean: 92.90, range: 68.77–163.23) (Fig. 3, Suppl. material 1: Fig. S2), whereas MiAmphiS had the lowest value (mean: 23.00, range: 17.17–34.53). The means of 16S_mod2, MiAmphiL and Modified 16Sar were 43.61, 43.50 and 44.89, respectively.

Amph16S had the highest universality for PCR amplification and taxonomic resolution in eDNA metabarcoding. Therefore, we regarded this primer set as the most useful and was applied in the subsequent field surveys.

In the actual metabarcoding analysis, using eDNA with Amph16S for the 160 field sites, we detected a total of 15 Anuran and one Caudate species from 122 water samples (see Discussion and Suppl. material 1: Table S4). None of the negative controls detected any of the species. Seven Anuran and one Caudate species were detected in the physical surveys (Suppl. material 1: Table S5). Combining the eDNA metabarcoding and physical surveys, 16 amphibian species were detected (Table 3). Amongst them, seven species were detected by both eDNA metabarcoding and physical surveys. Eight species (Buergeria buergeri, Bufo japonicus, Bufo torrenticola, Pelophylax sp., Rana ornativentris, R. tagoi tagoi, Zhangixalus arboreus and Z. schlegelii) were detected only in the eDNA metabarcoding and only one species (Cynops pyrrhogaster) was detected by physical surveys. Seven sequences from multiple geographic areas were assigned to the single anuran species Glandirana rugosa (Suppl. material 1: Fig. S3), although they showed a > 3% nucleotide divergence. This is consistent with a previous report of many local populations with high genetic divergence in this species (Sekiya et al. 2013).

Community composition differed significantly between the eDNA metabarcoding and physical surveys (PERMANOVA: p < 0 .001; Fig. 4; Table 4). The subsequent PERMDISP analysis further indicated that the amphibian community composition was significantly heterogeneous between the monitoring methods (p < 0.001; Table 5) and the number of species detected by eDNA metabarcoding was higher than that of physical surveys (GLMM: p < 0.001; Table 6). The amphibian community composition detected by the eDNA metabarcoding method encompassed that of the physical survey in all areas excluding area J (Fig. 5; Table 3).