Composition and distribution of fish environmental DNA in an Adirondack watershed

Abstract

Background: Environmental DNA (eDNA) surveys are appealing options for monitoring aquatic biodiversity. While factors affecting eDNA persistence, capture and amplification have been heavily studied, watershed-scale surveys of fish communities and our confidence in such need further exploration.

Methods: We characterized fish eDNA compositions using rapid, low-volume filtering with replicate and control samples scaled for a single Illumina MiSeq flow cell, using the mitochondrial 12S ribosomal RNA locus for taxonomic profiling. Our goals were to determine: (1) spatiotemporal variation in eDNA abundance, (2) the filtrate needed to achieve strong sequencing libraries, (3) the taxonomic resolution of 12S ribosomal sequences in the study environment, (4) the portion of the expected fish community detectable by 12S sequencing, (5) biases in species recovery, (6) correlations between eDNA compositions and catch per unit effort (CPUE) and (7) the extent that eDNA profiles reflect major watershed features. Our bioinformatic approach included (1) estimation of sequencing error from unambiguous mappings and simulation of taxonomic assignment error under various mapping criteria; (2) binning of species based on inferred assignment error rather than by taxonomic rank; and (3) visualization of mismatch distributions to facilitate discovery of distinct haplotypes attributed to the same reference. Our approach was implemented within the St. Regis River, NY, USA, which supports tribal and recreational fisheries and has been a target of restoration activities. We used a large record of St. Regis-specific observations to validate our assignments.

Results: We found that 300 mL drawn through 25-mm cellulose nitrate filters yielded greater than 5 ng/µL DNA at most sites in summer, which was an approximate threshold for generating strong sequencing libraries in our hands. Using inferred sequence error rates, we binned 12S references for 110 species on a state checklist into 85 single-species bins and seven multispecies bins. Of 48 bins observed by capture survey in the St. Regis, we detected eDNA consistent with 40, with an additional four detections flagged as potential contaminants. Sixteen unobserved species detected by eDNA ranged from plausible to implausible based on distributional data, whereas six observed species had no 12S reference sequence. Summed log-ratio compositions of eDNA-detected taxa correlated with log(CPUE) (Pearson's R = 0.655, P < 0.001). Shifts in eDNA composition of several taxa and a genotypic shift in channel catfish (Ictalurus punctatus) coincided with the Hogansburg Dam, NY, USA. In summary, a simple filtering apparatus operated by field crews without prior expertise gave useful summaries of eDNA composition with minimal evidence of field contamination. 12S sequencing achieved useful taxonomic resolution despite the short marker length, and data exploration with standard bioinformatic tools clarified taxonomic uncertainty and sources of error.

Keywords: Barcode sequencing; Computational biology; Dam removal; Environmental DNA; Fisheries restoration; Metagenetics; Mitochondrial 12S ribosomal RNA; New York state.

Figure 1. Sampling locations within the St. Regis River watershed and DNA yield.

(A) Location of St. Regis sampling sites relative to the northeast United States. (B) Overview showing relative positions of 22 sampling sites from four contiguous regions (separated by dashed lines). The direction of flow northward to the St. Lawrence River is indicated by a blue arrow. The yellow star indicates the location of the USGS water gauge from which flow data was obtained. (C) Numbered site locations in the lower portions of the watershed, above and below the former Hogansburg Dam, NY, USA. (D) Numbered site locations in the remainder of the watershed. The number of samples exceeding 5 ng/uL, an approximate threshold for achieving strong sequencing libraries, is indicated for both (C) and (D) according to the legend. Note the scale of (D) is compressed relative to (C). (E) Detail showing increased sampling at sites 3 and 4 in the vicinity of the Hogansburg Dam, NY, USA (now removed), which is denoted by the yellow box. (F) Average DNA yield at each site, with sample standard error indicated by error bars.

Figure 2. Relation between eDNA concentration and library yield.

Library yield is shown on the primary axis and eDNA concentration and flow are shown on the secondary axis, with samples sorted by increasing library yield. Samples are grouped subjectively into three categories of library strength. The red horizontal line corresponds to an initial total DNA concentration of 5 ng/µL, as a reference. Samples greater than 2 ng/µL were diluted to that value prior to the 12S preamplification reaction.

Figure 3. Heat map of taxon abundance by sample.

Rows represent detected taxa and columns represent water samples with at least 1,500 12S sequence counts. Samples are arbitrarily sorted by increasing site number and then by sample date, and are unlabeled for image clarity. Color intensity in each cell is scaled by percentile from 0 (no color) to 100% (darkest color) and four potential contaminant species were removed (see Methods). Taxa are ordered by total prevalence above the threshold and then alphabetically. (A) Heat map based on raw counts. (B) Heat map based on scaled log-ratio compositions with a minimum taxon proportion of 0.1% imposed prior to transformation.

Figure 4. Taxon proportions in technical replicates correlate well overall but exhibit a strong dropout effect.

(A) Each point represents scaled log-ratio compositions of individual taxa in two replicates of a single biological sample. Points are pooled across the eight technical replicate pairs (of 12 total) that had at least 1,500 counts per library. (B) Histogram of scaled log-ratio compositions of taxa that were detected in one replicate but not the second. (C) Histogram of mean scaled log-ratio compositions of taxa that were detected in both replicates.

Figure 5. Outlier taxa with respect to average eDNA composition.

Scaled log-ratio compositions were averaged across all biological samples with at least 1,500 total counts in which the taxon was present at 0.1% or more. Moxostoma anisurum has notably higher average composition than the majority of taxa, whereas the CYPRINID3 bin has notably lower average composition.

Figure 6. Catch per unit effort (CPUE) correlates with eDNA composition.

CPUE was scaled to a total of 100% across all taxa and then log transformed. Scaled log-ratio compositions were summed across all samples in which the taxon was detected at 0.1% or greater

Figure 7. Among-sample similarity in eDNA composition.

Color scale represents pairwise values of Spearman’s rank correlation coefficient and the order of samples is based on clustering by Ward’s method. Technical replicate pairs are marked by matching colored boxes. Taxa with a prevalence of less than four samples at a minimum abundance of 0.1% (prior to transformation) were excluded.

Figure 8. Per-taxon changes in average eDNA composition above and below the Hogansburg Dam, NY, USA.

Averages are of five biological replicates at each site, and only taxa detected in at least two of the five sampling events at a single site were included. Asterisks indicate zero detections above the 0.1% threshold for that taxon at the corresponding site. Taxa are sorted by average log-ratio composition at site 3 (below the dam), in descending order.

Figure 9. Haplotype distribution for two groups of taxa show shifts coincident with the Hogansburg Dam, NY, USA.

Samples with greater than 1,500 total counts are shown ordered from downstream to upstream, by sampling date, with the dam location marked by a red line. For samples with technical replicates, only the replicate with the highest total counts is shown. Values are shown as proportion of counts for equivalence of scale across samples. (A) Shows four Etheostoma bins, including Etheostoma haplotype 1 (Etheostoma H1) which is of uncertain taxonomy but is closest by edit distance to E. nigrum (see text for details). (B) Shows aggregate values obtained for reference sequences of Ictalurus punctatus and for the novel haplotype attributed to that species (Ictalurus punctatus H1).