Metabarcoding is (usually) more cost effective than seining or qPCR for detecting tidewater gobies and other estuarine fishes

Field sampling

Water samples were taken for eDNA at three sites known to contain tidewater gobies in the past: 16 samples Calleguas Creek (34 06′42″, 119 04′54″, Ventura County), 21 samples at Ormond Lagoon (34 8′23″, 119 11′20″, just west of Naval Base Ventura County; and 17 samples at Santa Clara River Mouth (34 4′8″, 119 15′53″, Ventura County). Site access was authorized by NBVC - Point Mugu agreement number Q2 DAR-Q N6923222MP001XY (Fig. 1).

Sampling supplies and analyses were from Jonah Ventures ® and included metabarcoding kits (sampling supplies, DNA extraction, PCR, sequencing of six pooled PCR replicates, and bioinformatics) and qPCR kits (sampling supplies, DNA extraction, and qPCR of three laboratory replicates). An additional one-time cost was required to validate a previously published cytochrome b tidewater goby assay (Schmelzle & Kinziger, 2016). At all sites, nearshore water samples were taken wearing latex gloves to reduce contamination with human DNA. Samples were then filtered through a 1-micron disc filter by pushing water through the filter with a 60cc Luer-lock syringe until clogging (mean sample volume: 174 cc +/- 0.141 S.D.). Filter capsules were purged of water before filling with preservative (tris-EDTA) then refrigerated until express shipped to Jonah Ventures ® for laboratory analysis.

To compare the efficacy of multi-species detection across methods, seines were taken at each sampling site for one of the estuaries (Calleguas Creek) that was sampled for eDNA (federal tidewater goby collection permit #PER0046428). Seine hauls were 2.4 m wide by 6.4 m distance on average in 0.6 m water depth. All captured fish were identified to species and released alive on site. Water temperature averaged 21.2 °C, conductivity was 32 mS (close to seawater), and dissolved oxygen was 8.3 mg/L.

Lab methods

As part of the eDNA sample cost, DNA was extracted from the DNeasy Blood & Tissue Kit (QIAGEN^©). Three qPCR replicates were run for the tidewater goby-specific cytB assay (forward primer: 5′-CCTCAATTCTCGTTCTACTAGTTGT-3′, reverse primer: 5′-GAGAATAAGTACGTCTGCTACTAGG-3′, probe: FAM-ACGTGCACTGACCTTCC GGCCTTTCTCC-MGB). Metabarcoding was done for the mitochondrial 12S ribosomal RNA (rRNA) gene which was PCR amplified from each genomic DNA sample using the MiFishUF (5′-GTCGGTAAAACTCGTGCCAGC-3′) and MiFishU-R (5′-CATAGTGGGGTATCTAATCCCAGTTTG-3′) primers with spacer regions. Amplicon size and PCR efficiency were visually inspected then cleaned by incubating. A second round of PCR was performed to complete the sequencing library construct. Final indexed amplicons from each sample were cleaned and normalized using SequalPrep Normalization Plates then pooled. Sample library pools were sent for sequencing on an Illumina NovaSeq 6000 (San Diego, CA, USA) at the Texas A&M Agrilife Genomics and Bioinformatics Sequencing Core facility (https://www.txgen.tamu.edu/) using the SP Reagent Kit v1.5 (500 cycles). Raw sequence data were then demultiplexed using pheniqs v2.10, primers were removed with Cutadapt v3.4, and read pairs were merged, denoised and chimeras were removed with VSEARCH v2.15.2. Amplicon sequence variants (ASVs) observed more than seven times were assigned taxonomy using a custom best-hits algorithm and a reference database that combined GenBank® and a Jonah Ventures® voucher sequence record, keeping any matching sequences with >90% agreement across top hits. Raw sequences were vouchered through the National Center for Biotechnology Information Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra; SRA# PRJNA924783).

Fish community assessment with metabarcoding

To estimate fish community composition, I determined a consensus taxonomy from the several hypothesized taxa per amplicon sequence variant (ASV). ASVs with <8 total reads were removed due to low confidence and, for a given ASV, only those hypothesized taxa with the maximum percent match per ASV were retained. Hypothesized taxa meeting these criteria were refined using geography and habitat information compiled from information on fish communities previously seined at Calleguas Creek from four sources (most of which were gray-literature reports), as well as those species reported within 1° latitude or longitude of the study site from the Global Biodiversity Information Facility (GBIF.org, 2024) or the list of invasive species from the Global Introduced Species Database (Poorter & Browne, 2005). Assigned taxa not on this combined plausible list were considered highly implausible and were reassigned to their closest plausible relative (e.g., Fundulus lima was re-assigned to the locally abundant, but unsequenced, Fundulus parvipinnis. The two are each other’s closest relatives, and so genetic similarity is expected.) In cases where there was more than one plausible taxon per ASV, the consensus was the lowest common taxon. Two consensus taxa not associated with estuarine habitats were considered valid assignments, but from an exogenous source (e.g., from bird feces).

Modeling cost with effort

I then created a cost model from the individual cost parameters, resulting in the following equation that considers estimated total cost H (sum of hour equivalents per location sampled) has a fixed cost plus a variable cost: ${\displaystyle fixed=\frac{P+G+A}{L}+B+WV+R}$ ${\displaystyle variable=k\left(v+n\left(tW+E+D+q\right)\right)}$

I parameterized the cost model based on >30 years sampling tidewater gobies and estuarine fishes (Table 1). Specifically, I considered that it takes time to maintain a seining permit (P) and build a plausible species list from background knowledge (B). It also costs money for field gear (G). Fixed costs per location sampled were assumed divisible by the number of locations visited (L) per sampling program. Field effort was multiplied by the number of field workers (W) that visit the location (with a trained fisheries biologist’s hourly rate equivalent to three field workers). Visiting a location (V) takes effort to prepare for as well as time and expenses for travel. At a location, one visits (k) sites. At each site, one takes (n) samples. Visiting a new site within a location takes time (v), as does taking a sample at each site (t) (seining takes a similar amount of time as collecting plus filtering a water sample). Seining can result in environmental impacts (E) related to incidental take and crushing burrows (this is the hardest cost to put into common units). At the lab, there may be costs to validating or developing an assay (A), depending on if the sequencer has validated the assay (no new cost), or if an assay has been published (moderate cost), or if a marker was not yet developed (large cost). Therefore, I investigate all three scenarios. There are also costs for sequencing each sample (q). For qPCR, there can be additional costs, per technical replicate (r), but these are now small enough ($5 USD) to be ignorable (I thus used the cost for three technical replicates per sample). Data entry and quality control (D) increases per site visit. The final cost was report writing (R), which, like data entry, is more labor intensive for multispecies than single species efforts. I expressed cost in the common currency of field-worker hours (assuming a 2023 technician’s hourly rate of $30 USD).

Modeling tidewater goby detection probability per sample by method

To estimate detection probability per sample for qPCR and seines, I used existing data for tidewater goby catch per unit effort (CPUE) and site-level detection obtained from 14 sites where tidewater gobies were present in a northern California survey (Schmelzle & Kinziger, 2016). To get conditional probabilities, I excluded sites where tidewater gobies were not detected. I also excluded the Virgin Creek and Pudding Creek sites, which only had three and five samples, respectively. I used the first three of Schmelzle & Kinziger’s (2015) six qPCR replicates to make the analysis comparable to typical studies. Detection per sample was coded as 0 (no detection) or 1 (detection), and estuary and sampling method (e.g., qPCR and seine haul) were factors.

To estimate joint detection probability per sample for qPCR and metabarcoding, I used observations from qPCR and metabarcoding at the J-Street Canal, Ormond Lagoon (34 8′23″, 119 11′20″) and the Santa Clara River Mouth (34 4′8″, 119 15′53″) so that qPCR and metabarcoding were factors.

For each methods contrast, a logistic regression was used to estimate the estuary-level coefficient associated with detecting a tidewater goby in a sample for a particular method where tidewater gobies were present. At sites where gobies are very easy to detect, any method would be suitable, whereas I was mostly interested in detecting gobies where they are rare. For that reason, rather than using mean detection to parameterize the cost-benefit model, I used the low-detection coefficient (and confidence intervals) from Schmelzle & Kinziger’s (2016) Ocean Ranch site (intercept (1.08) - 2.2 site effect) as a “rare-case” scenario. The Ocean Ranch coefficient was then transformed to a probability (using the plogis function in R) to parameterize the cost-benefit model. Detection per sample, γ, was thus used to calculate how the cumulative detection probability, d, increases with sampling effort (n*k). ${\displaystyle d=\left(1-{\left(1-\gamma \right)}^{nk}\right)}$

Then, to help guide method selection, I made three plots for each method: (1) detection probability (d) as a function of effort (n), (2) cost (H) as a function of effort (n), and (3) detection probability (d) as a function of cost (H). Because the costs of permitting and environmental effects associated with seining were uncertain, the cost plot was then reproduced for two additional assumptions about seining costs: no environmental costs, and no permitting or environmental costs. Also, given that fixed costs become relatively diminished as the number of locations sampled increases, I reproduced the cost plot for 10 locations (e.g., the effort it took to sample the southern tidewater goby populations) and again with 100 locations (the latter encompassing all potential tidewater goby populations).

Modeling species accumulation with effort

The number of species detected in a sample tends to increase with sampling effort to an asymptote. In addition, sampling effort can be applied at multiple scales. For instance, one can take many samples at a site and/or sample many sites. Diversity should increase with both efforts, but due to spatial heterogeneity, one might expect samples from the same site to be more similar than samples from different sites. I simulated the expected species accumulation with effort by taking 1,000 random orderings of the samples for seining and metabarcoding at Calleguas Creek (one or two replicates per site, and 1 to N sites), then calculating the average cumulative richness as a function of sample size, to generate a curve that started at the average richness per sample and leveled off at total cumulative richness. To make the collector’s curve more flexible for scenario modeling, I fit it using a multi-level model for species accumulation with effort using the Polce and Kunin method (Kunin et al., 2018), which can be represented as: ${\displaystyle S=\frac{\left(a+k^z{n}^c\right)}{\left(b+n^c\right)}}$ where S is the number of detected species, k is the number of sites (or area sampled), n is the number of samples per site, z increases with beta diversity among sites (or the species–area relationship independent of effort), c relates to the slope of the collector’s curve within a site, and a and b are constants.

To compare the efficacy of metabarcoding and seining for fish communities, I used non-linear least squares (with starting values = 1 for all Polce and Kunin parameters) to fit median log(S) as a function of n and k. Because metabarcoding had two samples per site, one of the samples was chosen at random to represent a site at each iteration. I then compared the predicted relationship between log(S) and k between metabarcoding and seining. Then, by linking effort to cost using the investment model, I estimated the relation between, log(S) and investment. Two metabarcoding strategies were evaluated. Double barcode refers to taking two samples per site (as done in my field samples), whereas single barcode projects a single sample per site (as commonly done by others). To compare among methods for species detection I plotted species accumulation against effort, cost against effort, and detection against cost. As for tidewater goby detection, I made separate plots to explore assumptions about the cost of seining and the cost efficiency of sampling multiple locations.

Statistics and data analysis, and model fitting were performed in R version 3.6.3 (R Core Team, 2020), cleaned using Tidyverse (Wickham et al., 2019), and visualized using ggplot2 (Wickham, Chang & Wickham, 2016). Original data, supplementary statistical tables and R code are available as a U.S. Geological Survey data release (Lafferty, 2024). The R code makes it possible to reproduce the analyses and figures and explore how different cost models or parameters might affect the results.

Tidewater goby assessment

Tidewater gobies were visibly present during sampling at Ormond Lagoon. All 18 qPCR samples at Ormond Lagoon and 15 of 20 qPCR Santa Clara River samples were positive for tidewater gobies, as were all 18 metabarcoding samples at Ormond Lagoon and 18 of 20 metabarcoding samples at Santa Clara River. Indeed, the logistic regression found no difference in detection between metabarcoding and qPCR at the Ormond Lagoon and the Santa Clara River Mouth sites (p = 0.99, Table S1). Thus, hereafter, I assumed detection by metabarcoding to be indistinguishable from qPCR. No tidewater gobies were detected by qPCR, metabarcoding or seining at the eight Calleguas Creek sites. For this reason, comparing the ability of seines and qPCR to sample tidewater goby used published data by Schmelzle & Kinziger (2016). Consistent with Schmelzle & Kinziger’s (2016) analysis of their data, detectability with qPCR was higher than with seining (p < 0.0001, Table S2). As described above, in the cost-benefit analysis, I used relatively low detection estimates from Ocean Ranch (qPCR and metabarcoding = 0.62, CI [0.42−0.78], vs. seining = 0.25, CI [0.12−0.43], Table S3). In general, there was more benefit to increasing sampling sites (k), than replicate samples within sites (n), so for the figures, n is set to one sample per site.

Cost effectiveness

Plotting detectability against sites sampled (Fig. 2, top panel) for Ocean Ranch showed that seining had less detectability than eDNA (especially at low density sites). However, eDNA and seining had nearly perfect detection for more than ten sites, even with relatively low detection probability. In other words, although eDNA had higher detection than seining, all methods achieved near-perfect detection under typical sampling efforts.

Costs increased with effort, but the different methods had different cost intercepts and slopes (Table 1, Fig. 2 middle panel). As an example, for a single sample at one location, metabarcoding cost 9 h and seining 42 h. Initial qPCR costs varied tremendously depending on whether an assay was validated or had to be developed. Indeed, having to develop assay from scratch made qPCR the least cost-effective method for a single location. Seining also had high fixed costs related to permitting and gear, and high environmental impacts per sample, though excluding permitting and environmental costs for seining did not make seining more cost effective than qPCR (Fig. 3). Metabarcoding was usually the most cost-effective method for detecting tidewater gobies due to its lower fixed costs than the other methods and its higher detection rate per sample than seining (Fig. 2 bottom panel). Exceptions occurred for cases where the sequencer had a specific tidewater goby marker in hand (true for future studies), or when many locations were visited (Fig. 4), in which case qPCR was slightly more cost-effective than metabarcoding.

Fish community assessment

Environmental DNA

Not counting one ASV (12 reads) vaguely assigned to Leuciscidae (e.g., G. orcuttii or L. armatus), metabarcoding detected DNA from 16 distinct fish species at Calleguas Creek (Table 2). Two of these were oceanic baitfish that seem most likely derived from bird feces. Of the remaining 14 estuarine fish species detected, 12 could be matched to a single species, and two were difficult to assign to the species rank. Specifically, although the Cypriniformes ASV was most likely Carassius auratus, Cyprinus was a plausible possibility too similar to exclude. Additionally, two gobies present in the system (Quietula y-cauda or Ilypnus gilberti) have not yet been sequenced and might be the source of the unidentified goby sequences. Metabarcoding detected to species all the fishes detected by seine except for Syngnathus leptorhynchus and, perhaps, Quietula y-cauda. Metabarcoding detected eight species not detected by seine. Moreover, it detected three previously unreported species for Calleguas Creek (Citharichthys stigmaeus, Mugil cephalus, and Pleuronichthys guttulatus). In addition to tidewater goby and the two unsequenced goby species, metabarcoding did not detect three previously reported freshwater fish species (Ameiurus nebulosus, Cottus asper, and Lepomis cyanellus). However, these species are detectable at other freshwater locations, indicating their DNA would have likely been detected had it occurred in the sample.

Species accumulation per sample by method

The Polce and Kunin species accumulation model fit to the qPCR and metabarcoding resampled accumulation curves made it possible to estimate parameters with error for incorporation into the methods comparison for species accumulation with effort (Table S4 for single sample, Table S5 for double sample).

Cost effectiveness

Unsurprisingly, taking two metabarcoding samples per site detected more species than taking a single sample, but not twice as many species (Fig. 5). And taking two samples per site was not twice as expensive as taking one. Thus, overall, there was little difference in the cost effectiveness of taking one or two metabarcoding samples per site with respect to detecting species. Nevertheless, given that taking two samples per site makes it possible to estimate within-site detection probabilities, taking two samples per site was the most informative strategy overall. Regardless, metabarcoding had higher species detection and lower costs than seining, and thus was much more cost effective. Assuming seining was less costly (Fig. 6) or several locations were visited (Fig. 7) made seining more cost effective, but not more cost effective than metabarcoding.