High-throughput environmental DNA analysis informs a biological assessment of an urban stream.

There is growing interest in the use of DNA barcoding and metabarcoding approaches to aid biological assessments and monitoring of waterbodies. While biodiversity measured by morphology and by DNA often has been found correlated, few studies have compared DNA data to established measures of impairment such as multimetric pollution tolerance indices used by many bioassessment programs. We incorporated environmental DNA (eDNA) metabarcoding of seston into a rigorous watershed-scale biological assessment of an urban stream to examine the extent to which eDNA richness and diversity patterns were correlated with multimetric indices and ecological impairment status designations. We also evaluated different filtering approaches and taxonomic classifications to identify best practices for environmental assessments. Seston eDNA revealed a wide diversity of eukaryotic taxa but was dominated by diatoms (36%). Differentiation among sites in alpha and beta diversity was greater when operational taxonomic units (OTUs) were classified taxonomically, but coarse resolution taxonomy (kingdom) was more informative than finer resolution taxonomy (family, genus). Correlations of DNA richness and diversity with multimetric indices for fish and macroinvertebrates were generally weak, possibly because Metazoa were not highly represented in our DNA dataset. Nonetheless, sites could be differentiated based on ecological impairment status, with more impaired sites having lower eDNA diversity as measured by the Shannon index, but higher taxonomic richness. Significant environmental drivers of community structure, as inferred from constrained ordination analyses, differed among kingdoms within the eDNA dataset, as well as from fish and macrobenthos, suggesting that eDNA provides novel environmental information. These results suggest that even a simple seston eDNA filtering protocol can provide biodiversity information of value to stream bioassessment programs. The approach bears further investigation as a potentially useful rapid assessment protocol to supplement more intensive field sampling efforts.

Introduction

Biological assessments of aquatic communities, or bioassessments, are approaches to evaluating the condition of waterbodies by surveying the composition, distribution and/or abundance of resident biota, often in conjunction with measures of physical habitat (US EPA, 2011). Over many decades, they have proved effective for assessing the level of anthropogenic disturbance in aquatic ecosystems, as well as for measuring progress of restoration efforts (Karr, 1981; Yoder and Rankin, 1995; Simon, 2002; Carlson et al., 2018). Various assemblages may be evaluated, but surveys in lotic systems have focused on the fish, benthic macroinvertebrate and, to a lesser extent, periphyton (biofilm) and amphibian communities (Ruaro and Gubiani, 2013). Benthic macroinvertebrates, including many insect larvae, are particularly sensitive to anthropogenic disturbances and form the basis of many stream bioassessment programs (Resh et al., 1995; Buss et al., 2014).

While bioassessments have become a critical component of water quality monitoring in the US and elsewhere, bioassessments have important limitations. Sampling for biota in aquatic systems can be labor intensive and specimen identifications are often highly dependent on trained taxonomy experts. As a result, many groups, including many microeukaryotes that may be functionally important such as microfungi (Gessner and Chauvet, 1994) are typically not assessed for lack of resources or taxonomic expertise. Other important groups, such as aquatic insect larvae, require careful microscopic examinations for identifications that set limits on realistic sample sizes and add to processing time and costs (Pfrender et al., 2010; Stein et al., 2014). These resource limitations necessitate compromises in the number of organisms identified and counted per sample, taxonomic precision of identifications, as well as where and how often samples are collected. Finally, inter-laboratory comparison studies have shown that taxonomic disagreements among experts can be significant (Stribling et al., 2008; Haase et al., 2010). Each of these factors impacts the utility of bioassessment results for adaptive resource management.

DNA-based specimen identification has been proposed as an approach that could increase throughput, precision, and taxonomic breadth of bioassessments (Pfrender et al., 2010). Indeed, several studies have shown that DNA-based approaches can enhance taxonomic identifications of larval invertebrates (Pilgrim et al. 2011, Hajibabaei et al. 2012; Carew et al., 2013, Gibson et al., 2015; Elbrecht et al. 2017), possibly at comparable cost (Stein et al., 2013). Nonetheless, DNA-based approaches also have limitations and uncertainties that hinder greater adoption of these methods. Standard “DNA barcoding” approaches that identify single specimens based on individual DNA sequencing reactions (e.g., Pilgrim et al., 2011) require as much time and effort to sample and sort organisms as traditional methods do, providing little time savings at higher cost (Cameron et al., 2006). On the other hand, ”DNA metabarcoding” (Taberlet et al., 2012) approaches which attempt to characterize whole communities in bulk samples are hindered by primer bias and competition in the PCR, causing some species to be absent from the measured species pool (false negatives) and estimates of relative abundances of others to be biased relative to biomass (Elbrecht and Leese, 2015). Reliance on databases to match DNA sequences to taxonomic names that are often sparse (Kvist, 2013; Trebitz et al., 2017; Leese et al., 2018) and sometimes inaccurate (Collins and Cruikshank, 2012) also cause significant concerns for interpretation of biodiversity patterns and environmental condition from DNA.

To date, few studies have directly compared results of DNA-based biodiversity analyses with assessments of aquatic ecosystem condition obtained from established bioassessment protocols. Such studies are critically important to the acceptance of these new methods in established biomonitoring programs. We sought to determine whether a simple protocol for filtering seston DNA (environmental DNA or “eDNA” represented in suspended particulate matter) followed by high throughput 18 s rRNA gene sequencing could provide data of value to stream bioassessment programs. This approach allows a broad survey of diversity and is different than previous studies that evaluated DNA of individual or bulked macroinvertebrate specimens which had been netted and then manually separated from other biota and organic matter (Stein et al., 2014; Elbrecht et al. 2017). It is unreasonable to think that this simple protocol could duplicate or replace the biodiversity information obtained by traditional bioassessment programs. For one thing, because eDNA methods lack the rich history of ecological research that led to development of sensitive multimetric indices and empirical models, they currently must rely on simplistic richness and diversity statistics that provide only a coarse assessment of biological condition. Nonetheless, eDNA methods could enhance efforts that are anchored in rigorous morphological specimen identification and enumeration. For example, DNA information can provide a snapshot of biodiversity information for the many taxa that are not amenable to morphological identification. If DNA metabarcoding provides information that correlates to traditional bioassessment metrics used to estimate environmental condition, it could be used as a rapid method to provide further spatial context by expanding to additional areas, providing more complete geographic coverage. Similarly, this approach could provide a rapid-turnaround assessment that can be deployed between traditional assessments or in target areas where environmental problems are suspected and may need quick management action.

Here, we compare biodiversity information obtained from eDNA in the water column of an urban stream to bioassessment data and impairment status designations obtained from a morphological-taxonomy based survey of fish and macroinvertebrates. Our goal was to determine whether the information collected from this simple eDNA filtering protocol is similar to and can complement information obtained from a bioassessment based on morphology, and whether it could be used to differentiate sites by level of biological impairment inferred from this traditional bioassessment. Because we ultimately want to design an efficient, informative, and transferable assay, we evaluated several parameters that could affect biodiversity estimates, including filter pore size and month of sampling. We also evaluated the importance of taxonomic resolution, or using taxonomy at all, in correlating DNA-based biodiversity metrics to traditional bioassessment endpoints, as precise taxonomic information for DNA sequences is unavailable for many groups and could limit the utility of this approach.

Materials and methods

Study Location

The Mill Creek watershed encompasses 42,994 ha within southwestern Ohio, USA (Fig. 1). The mainstem of Mill Creek flows south ~45 km from headwaters in rural-suburban plots through industrialized sections of the city of Cincinnati until it empties into the Ohio River. It is a low gradient, low volume, but flashy system, dropping 107 m over its length, with a median daily flow of 0.54m3s−1. Historically, Mill Creek has been designated highly impaired for aquatic life use, recreation (primary contact) and fish consumption due to extensive pollution from municipal point sources, storm water runoff, sewer overflows, and legacy contaminated sediments (Ohio E.P.A, 1994). Habitat is impaired in many areas and portions of the waterway are extensively modified, including concrete channel hardening of a portion of the mainstem and a barrier dam near the Ohio River for flood control. Over the last two decades, the city and watershed action groups have initiated actions to restore the stream by reducing chemical pollutant loadings, improving habitat, and removing toxic materials adjacent to the stream. Recently, there is evidence of significant success in many parts of the watershed, as documented by a return of many native fish and macroinvertebrate species (MBI, 2011, 2016). Nonetheless, many areas of the watershed remain impaired due to siltation, hydrological modifications, nutrients, chlorides, and polycyclic aromatic hydrocarbons (PAHs). In its current state, Mill Creek provides a strong gradient of physical habitat, water quality, and biological condition (see Supplementary Table 1) that makes it very useful for evaluating bioassessment metrics, DNA-based biodiversity measures and environmental correlates.

Fig. 1.

Location of sample sites in the Mill Creek watershed, Cincinnati, OH, USA. Shaded areas in the inset map represent areas of high urban intensity.

Sample collection

We sampled 26 sites within the watershed, including 21 sites on the mainstem of Mill Creek and five on two of the larger tributaries (East Fork and West Fork). Two 1-L water samples were collected once per site in early August and then again in mid-September 2016 by carefully submersing sterile bottles and capping them under water. One of the samples was filtered onsite using a portable vacuum pump attached to a 250-ml disposable filter funnel (Thermo Scientific) fitted with a 3-μm pore size cellulose nitrate (CN) filter (MilliporeSigma); the filter was then transferred to a 2-ml microcentrifuge tube containing 95% ethanol before being stored at −20 °C. The other water sample was transported on ice to the EPA laboratory in Cincinnati OH, where it was stored overnight at 4 °C before filtering with a similar filter funnel but fitted with an 0.4-μm pore size polycarbonate (PC) filter, then preserved in 95% ethanol and stored at −20 °C. The 0.4-μm PC filters were prone to clog before the full 1-L sample was filtered; in these cases, a second filter was used to increase the volume filtered. Still, for 25 samples, the total volume filtered was < 1 L, ranging between 500 and 950 ml. Each day that a sample was collected, molecular biology grade water was poured into two 1-L bottles at one of the sites visited and processed along with other samples as field blanks.

Water column chemistry, sediment chemistry, and fish and macroinvertebrate sampling are described elsewhere (MBI, 2016). Briefly, water column chemistry was characterized by grab samples collected 6 times at each site between June 16-October 15, 2016. Sediment chemistry was determined from samples of the upper 4 in. of bottom material, collected once per site in October of that year. Fish were sampled by electroshocking once from late July to early August and again in early to mid-October; fish species, composition, abundance, and condition were then used to calculate an Index of Biotic Integrity (IBI; Karr et al., 1986) for each site. Benthic macroinvertebrates were sampled using modified multiple-plate artificial substrate (Hester-Dendy) samplers for quantitative samples and dip net sampling for qualitative samples following Ohio EPA methodology (Ohio E.P.A., 2015); these data were then used to calculate an Invertebrate Community Index (ICI; DeShon, 1995) for each site. Artificial substrates were colonized over 6 weeks from mid-July to mid-September 2016 before removal, disassembly, fixation of biological specimens and subsequent taxonomic identification at MBI’s laboratory. Sites that did not meet minimum stream flow requirements for Hester-Dendy samplers were assigned a qualitative ICI based on dip net samples. For purposes of analysis, we assigned quantitative values for the qualitative ICI scores as the midpoint of the applicable quantitative range for that score: “Very Good” = 43; “Good” = 35; “Marginally Good” = 27. IBI and ICI scores are used as numerical biocriteria and each site was assessed as fully, partially, or not attaining its designated aquatic life use (MBI, 2016). Aquatic life use designations are tiered in Ohio, ranging from “limited resource” to “exceptional” for warm water habitat. In Mill Creek, all sites were previously classified as Warm Water Habitat (WWH) or Modified Warm Water Habitat-Channelized (MWH-C) by Ohio EPA. Thus, the combination of aquatic life use designation and attainment status within each aquatic life use describes an environmental degradation gradient.

DNA sequencing

Filters from water samples in this study were processed along with filters from 85 other aquatic samples from the region to help build a consistent regional database. Total genomic DNA was extracted using a modification of the Sambrook et al. (1989) phenol-chloroform-isoamyl alcohol (PCI) method. Filters in 2- ml microcentrifuge tubes were incubated overnight with agitation at 56 °C in 900 μl of lysis buffer (10 mM NaCl, 0.5% SDS, 100 mM EDTA, 100 mM Tris 8.0) supplemented with 400 mg proteinase K (Fisher Scientific). The tubes were heated to 91 °C for 10 min to inactivate the proteinase then centrifuged at 15,000g for 1 min to pull the filters to the bottom of the tubes. Following addition of 900 μl PCI (25:24:1; Fisher Scientific), the tubes were vortexed for 5 s to create an emulsion and centrifuged at 15,000g for 5 min. Aliquots of 700 μl from the aqueous layer were transferred to new 2 ml microcentrifuge tubes containing 700 μl chloroform-isoamyl alcohol (24:1), vortexed for 5 s, then again centrifuged at 15,000g for 5 min. Aliquots of 500 μl were transferred to fresh 2-ml microcentrifuge tubes along with 20 μl of 5M NaCl and 1.25 ml chilled ethanol, then inverted to mix and allowed to precipitate overnight at −20 °C. Samples were then centrifuged at 15,000g at 4 °C for 10 min to pellet the DNA before decanting the liquid, washing the pellet with 1 ml 70% ethanol, decanting again and vacuum drying the DNA pellet at 45 °C for 15 min or until no visible liquid remained. DNA samples were rehydrated in 200 μl TE (10 mM Tris, 1 mM EDTA) and concentrations estimated using Picogreen on a BioTek Microplate Reader before 1:100 dilution for PCR.

We used a dual-index PCR strategy to provide template for Illumina MiSeq sequencing (Kozich et al., 2013). In the first round of PCR, we used primers targeting the V4 region of the 18S rRNA gene (SSU_F04 GCTTGTCTCAAAGATTAAGCC and SSU_R22 GCCTGCTGCCTTCCTTGGA; Creer et al. 2010) and incorporating 5′ adapter sequences for a second round of PCR by dual indexing (5′-ACACTGACGACATGGTTCTACA-3′ and 5′TACGGTAGCAGAGACTTGGTCT-3′, respectively). The first round PCR contained 2 μl DNA template, 2 μl 10X PCR buffer, 0.6 μl MgCl2 (25 mM), 2 μl dNTPs (10 mM), 0.5 μl each 18S primer with 5′ adapter (10 mM), 4 μl 1X BSA and 0.1 μl Taq polymerase (5U/μl; Qiagen) and 9.9 μl ultrapure water. PCR was performed under cycling conditions consisting of an initial 2.5-min denaturing step at 94 °C, followed by 32 cycles of 30 s at 94 °C, 1 min at 50 °C, and 1 min at 72 °C and a final elongation step of 10 min at 72 °C. Agarose gel electrophoresis was used to confirm the absence of visible amplification products in all field blanks and to confirm the correct size of amplification products in test samples. PCR products were purified with the Qiagen Qiaquick 96 PCR Purification kit and quantified with Picogreen on a BioTek Microplate reader before normalizing to l0 ng/μl with Qiagen EB buffer. The second round of PCRs used index primers homologous to amplicon adapters (1 μl each), 1 μl of PCR amplicon and 17 μl Accuprime Master Mix (Thermofisher Scientific). Index PCRs were performed as follows: 95 °C for 3 min followed by 8 cycles of 95 °C for 30 s, 55 °C for 30 s; and 72 °C for 30 s; followed by a final extension step of 72 °C for 5 min. Index PCR products were then purified using the AM-Pure XP kit (Beckman Coulter Life Sciences), quantified with picogreen as before and normalized to the lowest nanomolar concentration in Qiagen EB buffer. Normalized Index PCR plates were then pooled into a single sample by combining 3 μl from each well into a 1.5 ml microcentrifuge tube. Amplicons were then sequenced using a 2 × 300 600-cycle Illumina MiSeq sequencing kit according to manufacturer’s protocols. Raw sequence data have been deposited in the NCBI Short Read Archive under Bioproject PRJNA542144, Biosample accessions SAMN11613569-SAMN11613667.

Read filtering and taxon assignment

We used a Usearch (v9.2; Edgar 2010) bioinformatic pipeline to process demultiplexed reads for 184 of 314 samples that came off a single MiSeq DNA sequencing run (other samples were not relevant to this study). Approximately 7.2 M paired reads were merged and primers removed with Cutadapt v 1.14 (Martin 2011). Full-length sequences that were smaller than 341 bp or possessed one or more expected errors based on cumulative Phred quality scores were filtered out, leaving 5.1 M sequences. Remaining sequences were dereplicated and unique sequences (1.8 M) were identified; those with < 4 observations in the total data set were removed as possible sequencing artifacts. The remaining sequences were then screened to remove chimeras and clustered into 2854 Operational Taxonomic Units (OTUs) of 97% similarity or more, after which all quality-filtered reads were mapped to these OTUs. Following removal of sequences from the 85 regional samples that were not relevant to this study (see Section 2.3), 3.3 million sequences were mapped to 2006 remaining OTUs distributed among 99 samples. On average, 33 thousand reads per sample were mapped to OTUs, with a minimum of 23,711 and a maximum 47,246 reads per sample.

Several approaches were employed to obtain taxonomic classifications for OTUs. First, OTUs were compared to the SILVA ribosomal RNA reference database (SSU Ref NR 128, September 2016) using SINA Online (Pruesse et al. 2012) to identify the lowest common ancestor within an identity threshold set at 0.8. Classifications based on the EMBL database using SINA Online were also retrieved. Two other approaches made use of a reference database of NCBI sequences composed of sequences with 80% or better identity to OTUs over their full sequence length which were obtained using Blast + (Camacho et al., 2009). NCBI sequences obtained (downloaded Nov 14, 2017) were trimmed to the length of aligned OTUs and their associated taxonomy information was downloaded using Entrez efetch (Kans 2017). This database was then queried using the USEARCH sintax and RDP v1.9.1 (Wang et al. 2007) taxonomic classifiers.

To create a consensus taxonomy, the lowest taxonomic rank was retrieved from each of the four classifications (SINA-SILVA, SINA-EMBL, sintax-NCBI, RDP-NCBI) for each OTU and analyzed with the classification function in the R package taxize (version 0.9.4, Chamberlain and Szöcs, 2013), using NCBI’s taxonomy database. The function provided a full rank taxonomy using a standardized set of taxonomic ranks, allowing evaluation of taxonomic precision of different classifiers as well as conflicts among them. If there were no conflicts between classifications, the most precise (lowest level) taxonomy was assigned the consensus taxonomy for the OTU, unless the classification was to genus or lower and the identity score was < 0.95, in which case a family level resolution was used. In cases where the classifiers gave conflicting results, taxonomies were generally moved to the lowest rank in which there was no conflict. Exceptions were made if one of the classifiers had an identity score of 1.0, in which case it was used, or when one of the classifiers had an identity score ≥0.95 while the conflicting classifier had < 0.90 identity. For comparisons of different taxonomic ranks, taxa that do not have taxonomic names at a rank were treated as distinct from other taxa at that rank if they were distinct at lower ranks. This had the possible effect of over-splitting some taxa at higher ranks but was judged to be more acceptable than lumping distinct taxa together at these ranks simply because an accepted name does not exist. For purposes of comparing metrics for different taxonomic ranks, a suffix-associated rank (e.g., supra-, infra-, sub-) was substituted for a main rank if it was not available. If a kingdom was not identified by the classification function, a higher-level rank (e.g., Opisthokonta or Stramenopiles) was substituted, if available.

Biodiversity analyses

Counts of each OTU in each sample were divided by sample totals to provide relativized OTU abundances, which were the basis of all biodiversity analyses. The Shannon index and taxa richness (rarefied to the minimum sample size) were calculated with the functions diversity and rarefy of the R package vegan (version 2.5–2, Oksanen et al. 2018), respectively. These alpha diversity estimates were normalized using Tukey transformations with the transformTukey function of the R package rcompanion (version 2.0.0) prior to evaluation by mixed model analysis. Mixed models were constructed with the R packages lme4 (version 1.1–17) and lmerTest (version 3.0–1), treating filter type and sample month as fixed effects and sample sites as random effects. Conditional and marginal coefficients of determination for these models were calculated with the function r.squaredGLMM in the R package MuMIn (version 1.42.1). Additional models included impairment status (designated aquatic life use and biocriteria attainment) as a fixed effect, with sample sites nested within. For this analysis, a single site classified as partially meeting the MWH-C designation was grouped with other sites that fully meet the MWH-C designation. Impairment status was then ordinally classified as MWH-C, MWH-nonattainment, MWH- partial attainment, and MWH-full attainment. Spearman rank correlations were used to compare the Shannon index and taxa richness estimates for each site based on molecular analysis with biodiversity metrics based on morphological analyses. For these analyses, the conditional mode of the alpha diversity metric for each site was extracted from the mixed model described above using the function ranef of the lme4 package, providing a single estimate for each site after accounting for fixed effects of filter type and date.

Assessments of beta diversity patterns were performed with the R package vegan. We Hellinger-transformed (square root of relative frequency in sample) abundances for unclassified OTUs and for OTUs classified at each taxonomic rank to reduce the importance of rare taxa on analyses. Transformed abundances were converted to pairwise matrices of Bray-Curtis dissimilarities for unclassified OTUs and for each taxonomic rank before a PERMANOVA was conducted using the adonis function to examine effects of filter type, month of sampling, environmental impairment status, and site differences on sample differentiation. Multivariate homogeneity of group variances was evaluated using the function betadisper. In addition, the Bray-Curtis matrix was evaluated graphically with nonmetric multidimensional scaling (NMDS) using the metaMDS function. We then performed redundancy analysis (RDA) on Hellinger-transformed abundances for each site in vegan to evaluate the explanatory power of abiotic variables on the distributions of taxonomic groups. To compare metabarcoding data to traditional bioassessment methods, we also performed RDA on Hellinger-transformed abundances of sampled macroinvertebrates and fish. Water chemistry, sediment chemistry, and spatial variables were included in the RDA to determine which factors best explained diversity patterns. This was followed by a partial RDA to determine whether chemistry effects were significant when spatial effects (drainage area, stream gradient) were controlled. To quantify and compare responses of taxonomic groups to different environmental dimensions, we performed separate RDAs for macroinvertebrates, fish, and the most common kingdoms identified by metabarcoding, using forward model selection to identify the linear combinations of abiotic variables that were useful in explaining variation in measured biodiversity patterns for each taxonomic group. Blanchet et al.’s (2008) double stopping criterion was used to minimize model overfit in developing these parsimonious RDA models.

Results

Following sequencing and bioinformatic analysis, all 2006 OTUs retrieved were classified as Eukaryota, 1476 of the OTUs were assigned to 38 kingdoms, 1339 of the OTUS were assigned to 106 phyla, 1288 of the OTUs were assigned to 162 classes, 1091 OTUs were assigned to 269 orders, 900 OTUs were assigned to 365 families, and 748 OTUs were assigned to 489 genera or species (Table 1). Of the 3.3 million sequences mapped to OTUs, 95.6% could be classified to phylum level and 69.0% to either genus or species. Stramenopiles (heterokonts) accounted for slightly more than half (52.4%) of all sequences, while Metazoa (animals), Viridiplantae (plants), and Fungi accounted for 12.6%, 12.4% and 6.6%, respectively (Table 1). Diatoms (Bacillariophyta) accounted for the majority of stramenopiles and for 36.0% of all sequences. Green algae (Chlorophyta) dominated among plants and represented 12.1% of all sequences. Among the metazoans, arthropods, nematodes, gastrotrichs, and bryozoans were the most dominant phyla at 3.6%, 3.2%, 2.3%, and 1.5% of all sequences, but Platyhelminthes, Mollusca, Annelida, Rotifera, Porifera, Cnidaria, Nemertea, Chordata, and Entoprocta sequences were also observed. The relative abundance of stramenopiles tended to increase moving downstream (Supplementary Fig. 1). Metazoans were most abundant in one tributary (East Fork), especially downstream of a wastewater treatment plant on the tributary, and were relatively rare farther downstream on Mill Creek.

Table 1

Relative abundance (%) of kingdoms averaged over sites (± standard error), counts of total OTUs in each kingdom, and counts of different taxa identified for each taxonomic rank. Taxonomic names are derived from the NCBI taxonomy database. Taxa names that begin with “_” are not actual kingdoms but lower ranks for which a kingdom was not identified. For this analysis, they were treated as separate kingdoms (see text).

Kingdom	Relative Abundance	OTUs	Phyla	Classes	Orders	Families	Genera-Species
Alveolata∣Ciliophora	3.34 ± 0.32%	113	2	6	20	36	55
Alveolata∣Perkinsea	0.03 ± 0.00%	4	1	1	1	1	0
Alveolata (Other)	0.97 ± 0.15%	88	7	7	14	23	25
Amoebozoa∣Discosea	0.00 ± 0.00%	2	2	2	2	2	2
Amoebozoa∣Discosea∣Flabellinia	0.00 ± 0.00%	7	1	1	1	1	3
Amoebozoa∣Discosea∣Flabellinia∣Dactylopodida	0.09 ± 0.03%	16	2	2	2	2	8
Amoebozoa∣Discosea∣Stygamoebida	0.00 ± 0.00%	1	1	1	1	1	1
Amoebozoa∣Gracilipodida	0.02 ± 0.00%	3	1	1	1	1	2
Amoebozoa∣Mycetozoa	0.00 ± 0.00%	1	1	1	1	1	1
Amoebozoa∣Mycetozoa∣Myxogastria	0.23 ± 0.08%	5	1	1	1	1	1
Amoebozoa∣Tubulinea	0.00 ± 0.00%	4	2	2	2	3	3
Amoebozoa∣Tubulinea∣Euamoebida	0.00 ± 0.00%	6	1	1	1	1	4
Amoebozoa∣Tubulinea∣Euamoebida∣Tubulinida	0.02 ± 0.00%	14	2	2	2	2	5
Amoebozoa (Other)	0.01 ± 0.00%	12	4	4	4	4	5
Apusozoa	0.12 ± 0.03%	17	4	4	4	6	8
Centroheliozoa	0.13 ± 0.02%	13	3	3	3	3	7
Haptophyceae	0.03 ± 0.00%	4	3	3	3	3	3
Jakobida	0.00 ± 0.00%	1	1	1	1	1	1
Metazoa	13.3 ± 3.76%	165	13	38	64	83	93
Opisthokonta∣Nucleariidae and Fonticula group	0.03 ± 0.00%	7	1	1	1	1	4
Opisthokonta (Other)	0.81 ± 0.11%	75	4	4	7	9	16
Rhizaria∣Cercozoa∣Imbricatea∣Plasmodiophorida	0.00 ± 0.00%	5	1	1	1	1	4
Rhizaria∣Cercozoa∣Imbricatea∣Silicofilosea	0.23 ± 0.11%	11	2	2	2	5	8
Rhizaria∣Cercozoa∣Imbricatea∣Spongomonadida	0.00 ± 0.00%	2	1	1	1	1	1
Rhizaria∣Cercozoa∣Thecofilosea	0.87 ± 0.35%	20	6	6	6	7	8
Rhizaria∣Cercozoa (Other)	2.01 ± 0.28%	111	7	7	7	11	21
Rhodophyta	0.01 ± 0.00%	3	3	3	3	3	3
Stramenopiles∣PX clade	1.65 ± 0.55%	2	1	1	2	2	2
Stramenopiles (Other)	49.88 ± 3.82%	308	15	18	35	53	88
Fungi	6.91 ± 0.92%	311	9	23	42	51	50
Viridiplantae	12.15 ± 1.32%	111	2	15	34	47	55
_Cryptophyta	3.38 ± 0.63%	19	1	1	3	5	10
_Glaucocystophyceae	0.04 ± 0.01%	1	1	1	1	1	1
_Heterolobosea	0.02 ± 0.00%	3	1	1	2	2	2
_Katablepharidophyta	0.20 ± 0.03%	5	1	1	2	2	4
_Dimorpha	0.00 ± 0.00%	1	1	1	1	1	1
_Paratrimastix	0.00 ± 0.00%	3	1	1	1	2	2
_Collodictyonidae	0.04 ± 0.00%	2	1	1	1	1	1
Eukaryota (Other)	3.29 ± 0.32%	530	0	0	0	0	0
Total	100.00%	2006	106	162	269	365	489

Alpha diversity

Filter type strongly impacted both the Shannon index and taxonomic richness, with the 3-μm CN filter demonstrating higher diversity than the 0.4-μm PC filter for all taxonomic ranks (Supplementary Table 2). Both the Shannon index and richness tended to be higher in August compared to July, but differences were not significant for some taxonomic ranks. Within kingdoms, only diversity of Fungi and Viridiplantae was not affected by filter type, while only Metazoa and Viridiplantae were not affected by sample month (Supplementary Table 3).

An important objective of our study was to evaluate how important taxonomic resolution of DNA identifications is for discriminating sites based on biological condition. To investigate this, we estimated variance due to sites for the Shannon index and richness at each taxonomic rank. Greater variance among sites with greater taxonomic resolution would infer better ability to discriminate site-to-site differences. Fig. 2 shows these variances as a proportion of total variance for different levels of taxonomic resolution (see also Supplementary Table 2), after accounting for effects of filter type and month of sampling. Results for neither index suggest greater variance among sites with increasing taxonomic resolution. In fact, for taxonomic richness, site variance was only significantly different from zero for kingdom-level richness and for OTU richness. Variance among sites was larger for the Shannon index and was significantly greater than zero at each taxonomic rank. However, as was the case for richness, variance among sites was proportionally highest at the kingdom level and site variance for OTUs was higher than for family and genus-species taxonomic ranks. Among different kingdoms, OTUs for Fungi had particularly high site-to-site variation, with more than 60% of variance for the Shannon index and richness attributable to site differences (Supplementary Table 3).

Fig. 2.

Proportion of total variance attributed to sample sites for Shannon index and Richness for each taxonomic rank. OTUs classified to genus or species were combined. Estimates that are significantly greater than zero are marked with asterisks (*: p < 0.05, **: p < 0.01, ***: p < 0.001).

Spearman rank correlations among Shannon index values and among richness values for different taxonomic ranks were positive and significant (Fig. 3). Few correlations with traditional bioassessment metrics were significant, the major exceptions being correlations between Shannon index and ICI, which were positive and relatively high for higher taxonomic ranks, as well as correlations between the number of EPT species (insects in the orders Ephemeroptera, Plecoptera, and Trichoptera) and the Shannon index, which were also positive (Fig. 3, Supplementary Fig. 2). Negative correlations were observed between ICI and richness at kingdom and phylum levels.

Fig. 3.

Spearman correlation coefficients between various morphological bioassessment metrics and DNA based alpha diversity measures for different taxonomic ranks. The first 9 rows represent morphological metrics while the rest are DNA metrics (SI: Shannon Index, RI: Richness). Correlation coefficients with significance greater than p ≤ 0.05 are color coded relative to the scale bar.

We added impairment status to the mixed model as a fixed effect, with random sites nested within, to see how DNA-based alpha diversity measures corresponded. For this analysis, we focused only on the kingdom rank and on taxonomically unclassified OTUs, as previous results suggested they had larger site-to-site variation for taxonomic richness than other ranks. Alpha diversity of kingdoms as measured by the Shannon index was smaller for sites classified as biologically impaired (MWH-C) than for sites classified as in better biological condition (fully or partially meeting WWH biocriteria; p < 0.05). Differences in Shannon diversity based on OTUs were not significant, however. In contrast, taxonomic richness was significantly higher for the more impaired MWH-C sites (p < 0.05), and this was true both for analyses based on kingdoms and for OTUs (Fig. 4).

Fig. 4.

Least square means of alpha diversity statistics for different designated uses (MWH-C: Modified Warmwater Habitat - Channelized, WWH- Warmwater Habitat) and impairment status (nonattainment, partial or full attainment of biocriteria standards). One site designated MWH-C that was classified as partially attaining the biocriteria was combined with others that were fully meeting the biocriteria in this analysis. Bars that have different letters are significantly different (p < 0.05).

Beta diversity

Compositional patterns in OTU relative abundances assessed by PERMANOVA were similar to alpha diversity patterns (Supplementary Table 4). Filter type and sample month contributed to differentiation among samples but had smaller impacts on sample variation than impairment status (9.5%) and sample sites (36%). Applying taxonomy to OTUs increased model informativeness, as measured by the size of residuals (Supplementary Table 5). Taxonomy at the genus through phylum level provided relatively similar model informativeness, mostly by increasing the amount of differentiation due to sites, while differentiation at the kingdom level was most informative, with filter, month, impairment status and sample sites explaining 68% of differences, and impairment status alone explaining 14.1%. NMDS plots demonstrated substantial overlap in ordination space for groups classified by impairment status (Supplementary Fig. 3). Differences in multivariate group dispersions at least partially explained the impairment status differences in PERMANOVA results (p < 0.001 for all taxonomic ranks), with less dispersion among the more impaired MWH-C sites and more dispersion among WWH sites, especially for those that fully or partially met biocriteria.

The RDA model using kingdom-level taxonomy was able to explain much of the site to site differentiation (adjusted R2 = 0.70), with the first 2 axes being significant (p < 0.015). Ordinations using other taxonomic ranks were slightly less informative (0.58 ≤ R2 ≤ 0.62), while unclassified OTUs were least informative (R2 = 0.52). In comparison, RDA models for macroinvertebrates and fish collected from Mill Creek were substantially less predictive, with adjusted R2 s of 0.40 and 0.29, respectively. Fig. 5A shows the ordination of environmental variables with kingdoms. On the first axis, several sediment metals and PAH are associated with drainage size (e.g., lower watershed) while nutrient-related factors are associated with stream gradient (upper watershed). The second axis indicated dispersion between dissolved oxygen, sediment iron, and lead. Only those kingdoms with large projections (more dispersion) are shown in the figure; most were intermediate, indicating little relation to the environmental factors projected, and not shown for clarity. A subsequent partial RDA that controlled for geographic drainage and gradient factors was still significant (p = 0.011) and generally separated sediment chemistry from water chemistry (Fig. 5B). Viridiplantae and Cryptophyta were associated with sediment metals and PAHs while Rhizaria∣Cercozoa were more closely associated with nutrients and other water column chemistry. Stramenopiles, Fungi, and Metazoa were less strongly associated with chemistry data. Similar ordinations for macroinvertebrates and fish are shown in Fig. 5C and D, respectively. Note that quantitative macroinvertebrate samples could only be collected at 21 of the 26 sample sites, and therefore this is based on more limited data. The most pollution tolerant organisms (Oligochaeta, Cricotopus bicinctus) segregated with sediment metals in this ordination while mildly intolerant organisms (Ceratopshche morosa, Chimarra obscura) segregated with several water chemistry parameters. Fish identified in Fig. 5D are generally considered moderately pollution tolerant.

Fig. 5.

Redundancy analysis (RDA) ordination illustrating relationships among key spatial, water column, and sediment variables and sampled biological communities. A) Ordination for major kingdoms identified by metabarcoding; B) partial RDA for kingdoms identified by metabarcoding, controlling for effects of drainage and gradient; C) ordination for sampled macroinvertebrates; D) ordination for sampled fish. For clarity, taxa near the centroids that did not show strong relationships with environmental variables are not shown. Black open circles represent sample sites.

As expected, drainage area and gradient accounted for most of the variation in the parsimonious RDA model of relative abundances of kingdoms, but water column zinc and nitrate concentrations were also significant (Table 2). These variables also were commonly explanatory in RDA models of OTU relative abundances within individual kingdoms. Other significant variables were restricted to single kingdoms: water column lead was explanatory for Stramenopiles, sediment PAHs were significant for Rhizaria∣Cercozoa, sediment copper was significant for Metazoa and sediment iron was significant for Viridiplantae. In comparison, sediment copper and sediment iron, as well as conductivity were explanatory for macroinvertebrates, while only drainage area was explanatory for fish.

Table 2

Coefficients of determination for full RDA models and parsimonious RDA models after forward selection of significant environmental variables. Presented are results for analyses of relative abundances of kingdoms, OTUs within individual kingdoms, macroinvertebrates collected from artificial substrates, and electroshocked fish. Chemistry results refer to water column chemistry unless otherwise stated.

	Adjusted Full Model R2	Selected Variables	Adjusted Cumulative R2	F	P
Kingdoms	0.704	Gradient	0.467	22.86	0.001
		Zinc	0.564	6.35	0.001
		Drainage	0.658	7.35	0.001
		Nitrate	0.705	4.53	0.001
Metazoan OTUs	0.352	Nitrate	0.156	5.61	0.001
		TotalPhosphorus	0.272	4.83	0.001
		SedimentCopper	0.327	2.88	0.002
Fungal OTUs	0.451	DrainageArea	0.157	5.65	0.001
		Sulfate	0.243	3.75	0.001
		Gradient	0.311	3.25	0.001
		Nitrate	0.363	2.80	0.001
		Zinc	0.390	1.91	0.011
Rhizaria∣Cercozoan OTUs	0.598	TotalPhosphorus	0.234	8.63	0.001
		DrainageArea	0.359	5.68	0.002
		DO	0.421	3.45	0.007
		SedimentPAH	0.474	3.23	0.003
		Gradient	0.511	2.58	0.009
		Nitrate	0.556	3.02	0.004
Stramenopile OTUs	0.558	Gradient	0.251	9.36	0.001
		Nitrate	0.353	4.79	0.001
		DrainageArea	0.436	4.39	0.001
		Zinc	0.513	4.49	0.001
		Lead	0.541	2.27	0.014
Alveolate OTUs	0.418	TotalPhosphorus	0.145	5.23	0.001
		Gradient	0.225	3.48	0.002
		DrainageArea	0.299	3.42	0.001
		Nitrate	0.371	3.51	0.001
		Zinc	0.409	2.36	0.007
Viridiplantae OTUs	0.423	Gradient	0.226	8.28	0.001
		DrainageArea	0.306	3.80	0.001
		Zinc	0.376	3.58	0.001
		SedimentIron	0.416	2.515	0.01
Macroinvertebrates	0.399	SedimentCopper	0.125	3.862	0.001
		SedimentIron	0.191	2.558	0.005
		Conductivity	0.259	2.639	0.002
		Gradient	0.314	2.355	0.009
Fish	0.289	DrainageArea	0.155	5.407	0.001

Discussion

DNA barcoding and metabarcoding approaches have long been offered as ways to increase the accuracy and precision of aquatic biological assessments, as well as to provide a path to lower costs (Pfrender et al., 2010; Hajibabaei et al., 2011; Stein et al., 2014). There is now overwhelming evidence that sampled DNA can reveal biodiversity data relevant to these assessments (Shelton, 2016; Deiner et al., 2017). Yet, to date, few studies have attempted to directly compare biodiversity metrics derived from sampled DNA to endpoints that established bioassessment programs use to describe levels of environmental impairment, such as multimetric fish and benthic invertebrate indices. Information on the complementarity, efficiency, redundancy, and novelty of DNA indicators will be required before existing bioassessment programs are likely to fully incorporate these new methods. Those studies that have evaluated biotic indices have either evaluated DNA barcodes of individually picked macroinvertebrate specimens (Stein et al., 2014) or compared bulked specimens obtained from sieved or netted macroinvertebrate samples (Aylagas et al., 2016; Elbrecht et al., 2017). These studies found good correlations between morphological and DNA-based taxonomic identifications but still require substantial manual field sampling and sample processing effort.

This study sought to determine whether a simple seston DNA (eDNA) sampling protocol, requiring minimal field labor and no presorting of specimens, could provide information relevant to a stream bioassessment. The state of Ohio uses tiered aquatic life use designations in which Modified Warmwater Habitat represents a more degraded condition than Warmwater Habitat (Yoder and Rankin, 1995). Standard alpha diversity and beta diversity metrics for eDNA clearly differentiated these levels of impairment within an urban stream. Further, multivariate analyses of seston community structure pointed to different environmental drivers than for traditional bioassessment endpoints (fish and benthic macroinvertebrates), suggesting that it provides nonredundant information. Of course, our protocol was not entirely simple, as the laboratory and bioinformatic procedures were complex, but these are procedures that can be easily standardized and, to a large extent, automated.

The depth and breadth of biodiversity that can be screened using DNA approaches is largely determined by the DNA barcode primers deployed. Studies that aim to differentiate taxa within the macroinvertebrate and fish communities often target the mitochondrial gene COI, as it has been shown to have high species-level discriminatory ability (Hebert et al., 2003). The hypervariable V4 region of the small subunit of the 18S rRNA gene that we targeted provides greater taxonomic breadth, but at a cost of reduced depth of information on metazoans (Drummond et al., 2015). There is not yet consensus on which region of the 18S rRNA gene is most suitable for standardized DNA barcode development, although the V2, V4 and V9 regions appear most informative (Hadziavdic et al., 2014). Choosing among primers that target these regions involves trade-offs in terms of the universality of the primers (primer bias), taxonomic depth available for assessment, and suitability to different sequencing platforms (amplicon length). The primer set we used has been deployed frequently in aquatic biodiversity surveys and we have found it to perform well in a recent lake plankton survey (Banerji et al., 2018). Nonetheless, other 18S rRNA primer sets may provide better representation of seston DNA diversity; research into optimized primers continues (e.g., Hugerth et al., 2014; Bradley et al., 2016) and empirical comparisons will be useful in future biodiversity surveys.

We identified relatively few fish and invertebrates, and relatively few metazoans in general, with our approach. Instead, our seston DNA dataset was dominated by stramenopiles, particularly diatoms, which may partially explain why alpha diversity of unclassified OTUs was not strongly correlated with alpha diversity and associated metrics for benthic invertebrates and fish. It is not clear that use of DNA barcode primers that target mitochondrial loci would have enhanced assessment of macrobenthos since our sampling was in the water column. Nonetheless, stream seston DNA integrates biodiversity information over a broad geographic area (Deiner et al., 2016) and is likely to reflect, at least partially, organisms that reside in the water column. As such, it presents a valuable contrast to assessment data for macrobenthos. The greater taxonomic breadth of 18S rRNA compared to mitochondrial genes also provides potential advantages since the component biota reflect a greater range of niche diversity, and therefore may be responsive to a wider array of environmental stresses.

The relative paucity of taxonomic information in DNA databases for non-metazoan eukaryotic kingdoms is a potential drawback to their use in bioassessments. High taxonomic resolution has been found useful to enhancing ecological signal (Lenat and Resh, 2001; Waite et al., 2004; Feio et al., 2006, Bennett et al., 2014) and most bioassessment programs strive to balance costs of taxonomic precision with environmental discriminating power. Knowing this, we invested effort in using multiple classification approaches and databases to taxonomically classify OTUs to the highest resolution possible. Interestingly, we did not find great value in high-resolution taxonomic classifications for this 18S rRNA gene dataset. Generally, sites could be differentiated nearly as well with unclassified OTUs as with taxonomic classifications. This suggests utility for development of “taxonomy-free” eDNA biomonitoring approaches (e.g., Apothéloz-Perret-Gentil et al., 2017) that calculate ecological values directly from OTU occurrences or relative abundances, and therefore avoid data limitations due to poorly resolved taxonomies.

For taxonomically classified OTUs, biodiversity measures tended to be at least as good at differentiating site-to-site variation when based on relative abundances of kingdoms as when based on lower taxonomic ranks, including genus level classifications. Just why kingdoms assorted by environmental differences as well or better than component phyla and other taxonomic ranks is unclear. While only 69% of OTUs could be classified to genus or species, this still represented 748 taxa available for analysis, which is quite large. It is possible that taxonomic classifications at lower taxonomic ranks (e.g., family, genus, species) were more error-prone due to the slow evolutionary rate, and thus coarser taxonomic resolution, of the 18 s rRNA gene. However, this is unlikely to explain why, for example, kingdom-level classifications were better able to discriminate site-to-site differences in richness than phylum, class, or order-level classifications, as taxonomic resolution for these ranks should be well within the capabilities of this barcode gene. It may be that the strong environmental gradient in this urban watershed created a strong kingdom sorting mechanism. There is some evidence to suggest that strong disturbance gradients decrease the importance of high taxonomic resolution for uncovering ecological signal (Bowman and Bailey, 1997; Olsgard et al., 1998; Marshall et al., 2006). The pattern we uncovered in beta diversity analyses in which communities residing in the most impaired sites were much more similar to each other than were communities residing in less impaired sites is consistent with this interpretation. Clearly, additional studies are needed to judge the generality of these results. If our finding that high taxonomic resolution of OTUs is not necessary can be replicated in other studies, the implications for bioassessment would be favorable, as it would suggest value in use of a seston DNA approach even without further enhancement of genetic taxonomy databases.

At the kingdom level, environmental quality, as determined by designated aquatic life use (MWH-C vs WWH), was easily distinguished by both the Shannon index and taxonomic richness. Although not significant, trends in both measures also were consistent with assessments of biocriteria attainment based on fish and macroinvertebrate multimetric indices (IBI, ICI). While the Shannon index for kingdoms increased with environmental quality, kingdom richness decreased with environmental quality. A similar response has been reported for benthic estuarine communities, in which the most polluted estuaries assessed by metabarcoding demonstrated the most richness (Chariton et al., 2015). The response may reflect a greater affinity for eutrophic or high sediment habitats of many microeukaryote kingdoms compared to the fish and benthic macroinvertebrates typically used to measure impairment. Of course, correspondence between our DNA-based indices and measures of biotic condition does not prove that they are useful (Hubert, 1971). Bioassessment programs rely either on biological indices built on key metrics of the biological assemblage that respond predictably to disturbance (e.g., IBI, Karr et al., 1986) or use empirical models to compare observed taxa to those that should be expected (O/E) in the absence of disturbance (Wright, 2000). Further development and validation of seston DNA as a bioassessment tool will require development of similar environmental response profiles based on metabarcoding data (Pawlowski et al., 2018). Possible approaches to achieve this include (i) identification and analysis of sensitive indicator taxa, analogous to EPT taxa in traditional bioassessments; (ii) development of “taxonomy-free” approaches based on ecological correlations of OTUs (e.g., Apothéloz-Perret-Gentil et al., 2017); (iii) or use of phylogenetic relationships between OTUs and well-characterized indicator taxa to develop new biotic indices (Keck et al. 2018).

Beta diversity patterns can highlight environmental response profiles that may be valuable to future bioassessment metrics and models. Ordination analyses indicated that different taxonomic groups assorted along different environmental gradients with, for example, plants (green algae) and cryptophytes assorting with many sediment metals and cercozoans assorting with many water column parameters. This differentiation may be used to help determine ecological condition, as evidenced by the explanatory power of impairment status in the PERMANOVA model. Clearly, more work needs to be done to understand the causality of these relationships before they are implemented in bioassessment programs as, for example, it is counterintuitive that autotrophs did not align more closely with nutrient gradients. Further, the variables we selected may be proxies for other, unmeasured environmental drivers. The most parsimonious RDA model for OTUs within kingdoms suggests that biodiversity within each kingdom was associated with different sediment and water column factors (nitrogen, phosphorus, copper, sulfate, zinc, lead, iron, PAHs), and that these were different from responses of traditional macroinvertebrate and fish endpoints, although this requires more exploration with larger datasets to be definitive. The sorting of different OTUs within Rhizaria∣Cercozoa with phosphorus was particularly strong and merits further research as a potential indicator of nutrient impact.

An important operational decision is whether to collect water samples in the field then transport them to the laboratory for sterile filtration and downstream processing or to filter the samples in the field. Filtration in the field simplifies transport and avoids concerns about biodiversity changes that may occur in the sample between sampling and filtration (Goldberg et al., 2016), but can greatly lengthen time at each site, limiting the number of sites that can be visited in a day. For turbid or eutrophic sites, 0.4-μm filters clog and must be periodically replaced, which can significantly increase time on site. Use of a larger pore size significantly reduces filtration time, to < 5 min in our experience for 3-μm filters, without clogging. The tradeoff, of course, is that smaller cells and biological materials may not be retained. Modeling of different pores sizes suggests there should be slightly more eDNA capture from an 0.4-μm filter with our nominal water volumes (Turner et al., 2014); however, clogging and restriction on the number of filters per sample (maximum of two for our protocol) reduced actual volumes filtered by 0.4-μm PC filters, contributing to lower DNA recoveries overall. We found higher biodiversity with the 3-μm CN filters but did find that one moderately abundant genus in the class Xanthophyceae was absent from 3-μm filtered samples. In general, we feel that field filtering with the 3-μm CN filters performed better than lab filtering with 0.4-μm PC filters, but either method should be acceptable if applied consistently for all samples.

Our study was limited to a single, highly disturbed watershed, and our results must be interpreted in that context. More studies are needed that evaluate DNA-based approaches in relation to bioassessment criteria and these studies need to evaluate different biotic community selection techniques (e.g., sieves, filters, centrifugation) and environmental matrices (e.g., sediment, biofilm), as well as different ecological resources (e.g., waterbody types, ecoregions). Spatial and biological relationships between seston DNA, benthic DNA, and traditional bioassessment metrics, in particular, merit further study. Studies are also needed to assess how best to transition these methods into established programs to optimize information content relative to resource and program management costs.

We believe that eDNA metabarcoding approaches can be complementary to traditional assemblage-focused (e.g., fish, macroinvertebrates, diatoms) bioassessments and would be best deployed in conjunction with those efforts, not in lieu of them. Use of the 18S rRNA gene sequencing analysis, in particular, provides a method to increase the taxonomic breadth of biodiversity surveys greatly with minimal additional field sampling. However, multiple barcode loci can and likely should be screened from the same filters, as combined screening of mitochondrial genes such as COI for Metazoa, the 16S rRNA gene for prokaryotes, and other genes (rbcL, ITS) likely will provide greater environmental informativeness (Stat et al., 2017). To the extent that relationships between traditional bioassessment data and eDNA are well-characterized, seston eDNA methods could be deployed as a supplementary rapid assessment approach to provide greater geographic and temporal context within a program that is anchored by more intensive traditional morphology-based bioassessments. This approach may provide a mechanism for more adaptive management, particularly at smaller scales of watersheds or municipalities, as target sites could be screened and assessed rapidly, allowing more rapid and effective management interventions.