Parallel, targeted analysis of environmental samples via high-throughput quantitative PCR

When analyzing environmental samples for DNA from multiple taxa, researchers must usually decide between iterative analyses with single-taxon assays—which are reliable and sensitive, but also laborious to apply—and approaches such as metabarcoding that can simultaneously target multiple species, but which are less sensitive for detection across taxa. Here, we test an intermediate approach that allows efficient, parallel assessment of taxon-specific qPCR assays via high-throughput quantitative PCR (HT-qPCR). Based on an assessment of over 500 environmental samples, we found that sensitivity and specificity of our HT-qPCR approach were similar (con-cordance 0.900–1.000) to values achieved through single-species qPCR in six out of seven assays tested. Thus, HT-qPCR may provide analyses of similar quality as single-species qPCR analyses for environmental DNA, but at a lower cost per taxon. We see this approach as being a valuable addition to the eDNA sampling toolbox, particularly for situations where reliable inferences are needed for a defined suite of rare invasive or imperiled taxa.

qPCR despite failure to detect the species in eight replicate metabarcoding analyses. On the other end of the spectrum, approaches based on high-throughput sequencing (HTS) allow simultaneous detection of dozens to hundreds of taxa (e.g., eDNA metabarcoding; Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012). Although these sequencing-based methods have a broader taxonomic scope, the inference process is more complex, resulting in increased uncertainty about species presence or absence (Deiner et al., 2017). For example, in a typical metabarcoding protocol, the targeted genome region must be amplified across all target taxa (potentially confounded by amplification bias; Evans et al., 2016), correctly indexed and demultiplexed (Schnell, Bohmann, & Gilbert, 2015), bioinformatically processed to remove errors, and then assigned to taxa by comparison with a reference library (Murali, Bhargava, & Wright, 2018). The combined effects of primer bias and template competition may result in some templates not being sequenced. In addition, metabarcoding datasets are also susceptible to taxon identification errors due to inaccuracies in species DNA databases or cross-sample DNA contamination (Port et al., 2015).
An optimal solution for many applications would preserve the sensitivity, specificity, and simple interpretation of single-species approaches, but allow for efficient, parallel analysis for detection of multiple taxa. Some studies have sought to accomplish this goal through nested PCR, where DNA is initially enriched across a broad taxonomic group, followed by targeted PCR of those amplicons for individual taxa (Stoeckle, Das, & Charlop-Powers, 2018) or by simply multiplexing qPCR assays (Jo, Fukuoka, Uchida, Ashimaru, & Minamoto, 2020;Pilliod, Goldberg, Arkle, & Waits, 2013;Tsuji et al., 2018). Here, we test a related approach, known as high-throughput qPCR (HT-qPCR), with which large numbers (e.g., thousands) of individual reactions are run in parallel using microfluidic or assay-printed plates to assay samples simultaneously across different primer-probe sets. Because very small (<0.1 µl) reaction volumes are used, samples can be assayed many times and reagent costs are substantially reduced. However, the basic workflow and data output is the same as for single-species qPCR. Although HT-qPCR systems have been successfully applied to microbial eDNA (Shahraki, Heath, & Chaganti, 2019;Waseem et al., 2019), to our knowledge this is the first published application of the method to aqueous macrobial eDNA. Below, we validate the sensitivity and specificity of a suite of assays (Table 1) on an OpenArray (Life Technologies) HT-qPCR platform and then compare the outcomes of this platform with those from single-species qPCR across over 500 environmental samples.

| Samples
We assessed HT-qPCR in terms of specificity, sensitivity, and level of concordance with single-species qPCR analyses using a combination of DNA templates derived from tissue extractions, synthetic oligonucleotides, and environmental samples.

| Synthetic DNA
To test the sensitivity of the HT-qPCR platform, we generated a multiplexed dilution series, each containing 0.8-520.8 copies/µl of each assay template. We used a Qubit fluorometer to quantify synthetic gene fragments (gBlocks; Integrated DNA Technologies) and then used these to create a 5 × dilution series (10, 50, 250, 1,250, and 6,250 copies of each template per reaction). Each dilution level was assessed in triplicate.

| Rattlesnake Creek
To test concordance of the HT-qPCR platform with conventional qPCR analyses, we assessed 48 environmental samples collected from the Rattlesnake Creek basin in Western Montana, USA ( Figure 1). Bull trout, brook trout, brown trout, rainbow trout, and westslope cutthroat trout are known to occupy this basin, and western pearlshell mussel are thought to be absent.

| Bull Trout Inventory
To test concordance of the HT-qPCR platform with conventional qPCR analysis over a larger dataset, we also assessed 476 samples collected from streams across Western Montana, USA, as part of the Rangewide Bull Trout eDNA Project (Young et al., 2017;Figure 1). All environmental samples in this study were collected by filtering 5 L of stream water through a 1.5-µm pore size glass microfiber filter (GE Health Sciences) in the field, which was then stored in silica desiccant, as described in . We then extracted DNA from one half of each filter using a modified DNEasy Blood and tissue DNA extraction protocol (Carim, Dysthe, Young, McKelvey, & Schwartz, 2016). Extracted DNA was then stored at −20°C for future analysis. All sample handling and extraction were done in dedicated low-DNA laboratory spaces following stringent quality control protocols.
Assays for the HT-qPCR chip were synthesized as described in their original publications, except for the westslope cutthroat trout and western pearlshell mussel assays. Because the HT-qPCR chips used were unable to accommodate mixed bases within primers, a degenerate base within the reverse primer sequence for each assay was converted to a single base (loci indicated in Table 1). These primer modifications could reduce sensitivity of the assays for some haplotypes, but single base mismatches internal to a primer typically have a modest impact on amplification efficiency (Wright et al., 2014).

| qPCR
To test concordance of the HT-qPCR platform with single-species qPCR, most of the environmental samples were assessed with both approaches. All 48 environmental samples from Rattlesnake Creek were analyzed with taxon-specific qPCR assays using singlespecies qPCR for detection of bull trout (mitochondrial marker), brook trout, brown trout, rainbow trout, and westslope cutthroat trout following laboratory protocols specific to each assay Dysthe, Rodgers, et al., 2018;Wilcox, Zarn, et al., 2018). Of the 476 Bull Trout Inventory samples, 465 (97.7%) had been previously analyzed for the presence of bull trout DNA using either the mitochondrial or nuclear bull trout marker. Note: Target species, original citation, locus gene, and oligonucleotide sequences (F = forward primer, R = reverse primer, P = hydrolysis probe) for each. All hydrolysis probes include a nonfluorescent quencher (NFQ), minor-grove binding (MGB) moiety, and a FAM fluorophore. All assays are designed for a 60°C annealing temperature. Underlined base converted from mixed base "R" in original assay.

| High-throughput quantitative PCR
We used targeted pre-amplification to enrich samples prior to HT-qPCR analysis. This was necessary to retain sensitivity because we expected eDNA to be present at low concentrations (<3 copies/ µl) and the reaction chambers for HT-qPCR are small (<0.1 µl). This based on their position on the OpenArray. As recommended by the manufacturer, we avoided these wells when we designed the plate.
The genotyping plates require two allele-specific hydrolysis probes; we used an identical probe sequence as the original assays for both the FAM-and VIC-labeled probes (same sequence in both probes). Plate F I G U R E 2 Study design showing indicating environmental Samples, Molecular analyses, and Concordance tests across detection platforms. Environmental samples for this study were collected from Rattle Snake Creek samples (n = 48, blue) and the Bull Trout Inventory samples (n = 476, orange). Samples from Rattlesnake Creek were analyzed on a single-species qPCR platform for the presence of bull trout (mitochondrial marker), rainbow trout, brown trout, brook trout, and westslope cutthroat trout DNA. Samples from the Bull Trout Inventory were analyzed on a single-species qPCR platform for the presence of bull trout DNA only (mitochondrial or nuclear marker; 465/476 samples). All samples were analyzed on a HT-qPCR platform for the presence of bull trout (mitochondrial and nuclear markers), rainbow trout, brown trout, brook trout, westslope cutthroat trout, and western pearlshell mussel DNA. Using these data, we conducted concordance tests of detection (+) and nondetection (−) on the single-species qPCR (y-axis) and HT-qPCR (x-axis) platforms. Confusion matrices indicate results for (a) bull trout mitochondrial, (b) bull trout nuclear, (c) brown trout, brook trout (d), and rainbow trout (e) DNA in Rattlesnake Creek samples and bull trout mitochondrial (f) and bull trout nuclear (g) DNA in Bull Trout Inventory samples. Westslope cutthroat trout analysis results from the HT-qPCR platform were discarded due to poor amplification. One plate of HT-qPCR platform results for bull trout mitochondrial DNA (81 samples) was also discarded due to amplification in some negative controls

| Negative controls
The HT-qPCR plate for analysis of Rattlesnake Creek included 24 notemplate controls, and each of the six plates used for analysis of Bull Trout Inventory samples included 12 no-template controls for a total of 96 negative controls across the study (14% of samples analyzed).

| Concordance tests
We used the irr package in R (Gamer, Lemon, & Singh, 2019; R Core Development Team, 2018) to calculate Cohen's kappa (Cohen, 1960) between HT-qPCR and single-species qPCR in terms of detection/nondetection for each environmental sample. We also calculated raw concordance (proportion of detection/nondetection results in agreement) between the two methodologies. For the 48 Rattlesnake Creek samples, we assessed concordance for bull trout, brook trout, brown trout, and rainbow trout (westslope cutthroat trout dropped from the analysis due to poor assay performance on the OpenArray system, see below). For the 465 Rangewide Bull Trout Inventory samples with paired single-species qPCR data, we assessed concordance of the two bull trout markers (mitochondrial and nuclear) with those for single-species qPCR (both markers used to generate single-species qPCR records). An outline of samples and analyses used in these concordance tests is provided in Figure 2.

| RE SULTS
The westslope cutthroat trout assay performed poorly on the OpenArray system. Most amplification curves, even among the positive controls, lacked a clear linear amplification phase. We disregarded data from this assay in subsequent analyses.

| Sensitivity and specificity
All assays except that for westslope cutthroat trout provided unambiguous amplification curves, and no tissue-derived DNA samples from nontarget species resulted cross-amplification. All assays other than that for westslope cutthroat trout resulted in 100% amplification across all 15 standard curve dilutions down to 0.8 copies/µl template (10 copies/reaction).

| Rattlesnake Creek samples
Results for the 48 environmental samples from the Rattlesnake Creek basin using single-species qPCR and the HT-qPCR system were highly concordant. When comparing bull trout, brook trout, brown trout, and rainbow trout assay across platforms, Cohen's kappa ranged from .832 to 1.000 and was always significant (p < .01). Raw concordance ranged from 0.957 to 1.000 for the same assays ( Figure 1, Table 2). The spatial patterns of species detection for bull trout, brook trout, rainbow trout, and brown trout were also consistent with a previous understanding of fish distributions within the basin (Montana Fish, Wildlife, & Parks, 2019). Western pearlshell mussel are not known to inhabit Rattlesnake Creek, and all analysis results on the HT-qPCR system were negative. None of the 24 negative controls exhibited amplification.

| Bull Trout Inventory samples
In Note: In all cases, Cohen's kappa (κ) was highly significant (p < 0.01), and agreement (raw concordance) between methods was > 0.9. HT-qPCR amplification curves for westslope cutthroat trout were ambiguous. Western pearlshell are not known to occupy the Rattlesnake Creek basin, and all OpenArray analysis results were negative.

| D ISCUSS I ON
We found that HT-qPCR worked with six of seven assays developed for single-species qPCR, without further optimization. These six assays generated specific and sensitive species detection from environmental samples and were highly concordant with single-species qPCR analyses across hundreds of samples and with known spatial patterns of species' occupancies in a well-characterized stream.
This work demonstrates that HT-qPCR is an efficient, multi-taxa approach for eDNA sampling that retains the sensitivity and specificity found in single-species qPCR analyses.
Metabarcoding is currently the method of choice for multi-species detection using eDNA sampling. However, when compared with qPCR assays for detection of specific species, metabarcoding tends to be less sensitive and accurate. In metabarcoding, biased species detection is a persistent problem (Kelly, Shelton, & Gallego, 2019), despite extensive efforts to minimize bias introduced during template enrichment (Bylemans et al., 2019). Here, we describe a targeted, parallel method, which takes a fundamentally different approach from that of metabarcoding. HT-qPCR approaches rely on targeted species detection and achieve rates of detection and specificity comparable to those of single-species qPCR. However, because HT-qPCR performs analyses in parallel for several to dozens of templates, it is much more cost-effective than single-species qPCR when questions of species presence involve more than a handful of species.
The levels of concordance between single-species qPCR and the HT-qPCR platform were high, especially given the low concentration copies/reaction for the mitochondrial and nuclear marker, respectively, with 76% of samples having <10 mitochondrial copies/reaction. When these 42 positive samples were run twice with the same mitochondrial marker, raw concordance was 0.857 and concordance with the nuclear marker was 0.923 (Dysthe, Rodgers, et al., 2018).
For the hundreds of samples drawn from the same survey for this analysis, raw concordance was comparable (0.900-0.988), leading us to conclude that there was no difference in either sensitivity or specificity between single-species qPCR and HT-qPCR.
Although the cost of analysis per sample is greater for HT-qPCR when compared to single-species qPCR, the cost per taxon can be substantially lower, depending on the number of species of interest.
In terms of laboratory time and reagents, the pre-amplification step is roughly equivalent to a single-species qPCR analysis. The custom OpenArray genotyping chip, preloaded with assays, is approximately $700 USD at the time of writing. Given these costs and our analytical approach, we estimate that it costs as much to run the described panel (6-7 taxa) as it would to run three separate, single-species qPCR analyses per sample. We also note that the combination of assays and samples is flexible. For example, an OpenArray chip could include assays for 26 different taxa against 32 samples (three technical replicates each, not counting control samples). Other HT-qPCR platforms such as the Fluidigm BioMark or Takara Bio SmartChip (formerly WaferGen), which do not require preprinting of hydrolysis assays into a suite of manufactured chips, could be even more flexible for assay optimization and small-scale analyses. Further, because pre-amplification for the entire suite of taxa requires only as much template as for a single analysis with single-species qPCR, this approach does not deplete DNA extracts as rapidly, which is important if the intent is to archive samples for future analyses (Dysthe, Rodgers, et al., 2018).
The cause of poor performance of the westslope cutthroat trout assay is unknown. We tested for an impact of the pre-amplification step by running pre-amplified, diluted samples on single-species qPCR for this assay and did not identify any issues (i.e., qPCRs all resulted in normal amplification curve morphology; data not shown).
One contributing factor could be the inability to accommodate degenerate bases in the primers for this assay on the OpenArray genotyping platform. Of note, OpenArray plates are also offered in a "gene expression" formulation that may provide optimal conditions for eDNA sampling applications (Grigorenko et al., 2014), including incorporation of degenerate bases in primers and probes. However, the gene expression chips are roughly double the cost of the genotyping plates that we tested (>$10,000 USD per order of 10 chips).
Other platforms such as the Takara Bio SmartChip also have the flexibility to include degenerate bases in assay primers and probes (Takara Bio, 2018).
Start-up costs associated with an assay-based approach could be limiting for some multi-taxa questions. If taxon-specific assays are not already available, this HT-qPCR approach has a high initial cost relative to other multi-taxa methods such as metabarcoding and capture enrichment ; development of each taxon-specific qPCR assays requires extensive and labor-intensive testing and validation. However, this development burden is rapidly being reduced as taxon-focused projects have produced an abundance of qPCR assays in the published literature. Thus, the library of validated qPCR assays for eDNA-currently on the order of several hundred-is rapidly growing (Tsuji et al., 2018).
The pre-amplification step, which is necessarily to provide sufficient template abundance for downstream analysis, results in intrinsically greater laboratory-derived contamination for HT-qPCR than for single-species qPCR, which can typically be done without any handling of PCR products. Like other multi-taxa methods that require handling of PCR products, additional precautions are necessarily to avoid contamination (e.g., dedicated UV hoods for postamplification handling, as in this study). We observed false inference errors in seven negative controls, or 1.6% of control analyses, of which 71% were isolated to a single marker on a single OpenArray plate. Where reliable inferences are critical, repeated analyses from samples collected over time and multi-laboratory validation of results (Sepulveda et al., 2020) may further improve confidence.
Pre-amplification could also present a problem if oligonucleotide interactions caused the formation of products (e.g., primer dimers, chimeric sequences) that either reduced sensitivity or generated false-positive results. In this study, we used low primer concentrations and minimal pre-amplification cycles to reduce this risk. Although assessing and mitigating all potential oligonucleotide interactions would be difficult (e.g., 182 pairwise primer combinations were present in our pre-amplification solutions), we did not observe any problematic results. In previous studies of targeted pre-amplification, researchers have found that assay performance and consistency actually improve with increasing number of pre-amplification targets because primer interactions are diluted across a greater number of pairwise combinations (Andersson et al., 2015).
Multiplex qPCR, which has been used in several eDNA sampling studies, also sometimes seeks to combine multiple, taxon-specific assays for different taxa into a single qPCR solution in order to save on time and reagent costs. This approach can also be sensitive and cost-effective, and does not require specialized equipment, as is the case for HT-qPCR. However, the number of parallel analyses possible in multiplex qPCR will typically be limited by the number of available fluorescent channels (typically 3-6). Multiplexing conditions are likely also much more stringent than required for low-level pre-amplification prior to HT-qPCR analysis. In multiplex qPCR, all assays, including probes, are co-amplified for the entirety of the analysis, whereas in HT-qPCR, after an initial pre-amplification with primers only, reactions are partitioned by individual assay.
Finally, an important limitation of all assay-based eDNA approaches, relative to sequencing, is that inferences about habitat occupancy can only be made about species that are targeted. This contrasts with metabarcoding, for which generic primers can be designed to cover a large suite of potential taxa (e.g., all "teleosts" or all "anurans"). Similarly, new haplotypes of known taxa can be problematic for assay-based approaches, resulting in failure to detect a species, whereas in HTS approaches the new haplotype could be both detected and characterized for future studies. Thus, HT-qPCR currently represents a useful tool that can significantly enhance the power and efficiency of eDNA surveys that rely on conventional qPCR approaches, may be optimal when species that are only distantly related are of simultaneous interest, and can potentially serve as a complimentary tool to more community-wide approaches based on metabarcoding (Harper et al., 2019;Thomsen et al., 2012).

| CON CLUS ION
For many biological inventories in aquatic systems, robust estimates of habitat occupancy of several to dozens of taxa are the goal. The development of single-species qPCR analyses has yielded highly accurate estimates of species occupancy, but the method is costprohibitive when large numbers of taxa are of interest (e.g., suites of invasive species or aquatic communities containing a number of native species of conservation concern; Harper et al., 2019; . Here, we use extensive validation across hundreds of samples to demonstrate that HT-qPCR can produce inferences that are highly sensitive and concordant with single-species qPCR, facilitating design of efficient multi-taxa panels.

ACK N OWLED G M ENTS
We thank Kristy Deiner and Joseph Craine for suggesting pre-amplification prior to HT-qPCR analysis and for productive conversations about the approach early on in this project. We also thank Xin Guan,

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
TMW, KSM, MKY, MKS, and RFL conceived and designed the study.
TMW, MKY, KSM, MKS, and CE contributed to data acquisition, analysis, and interpretation. All authors contributed to writing the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw data from the Rattlesnake Creek analyses and broader bull trout assessments along with the associated R scripts for analysis are available on GitHub [https://github.com/taylo rmwil cox/HT-qPCR].
All other species data from the Rangewide Bull Trout eDNA Project samples will be incorporated into the Aquatic eDNA Atlas Database (https://www.fs.fed.us/rm/boise /AWAE/proje cts/the-aquat ic-eD-NAt las-proje ct.html).