We collected water samples for eDNA isolation from Eagle Creek, a small first-order stream (Q = 0.1-324 L/s) in the Kalamazoo River watershed that flows approximately 12 km from its source in the Fort Custer Training Center (Michigan Army National Guard) to its confluence with the Kalamazoo River (42.3374N, 85.3377W) near the town of Augusta, Michigan. We sampled six locations in Eagle Creek and two sites in the Kalamazoo River (Figure 1). In the Kalamazoo river, we sampled one site upstream (Location 7) and one site downstream (Location 8) of the confluence of Eagle Creek to detect the new species that Eagle Creek contributed to the river and novel species carried by the Kalamazoo River.

Samples were collected on 29 April 2015, between 0940 and 1230 h. We sampled from the furthest downstream site (Location 8) to the most upstream site (Location 1) (Figure 1). At each sampling location, we collected three 250-ml replicate water samples to yield a total of 24 samples. To avoid contamination, we collected each surface water sample with an extendable pole with a dip bucket attached to the end. To maintain semi-sterile conditions, the pole and dip bucket were wiped with 10% bleach and rinsed with reverse-osmosis (RO) water prior to collecting each additional sample and between locations. For each sampling location, a 250-ml bottle filled with RO water in the laboratory prior to sampling, was transported alongside field sampling bottles to serve as a field control. One field control bottle was opened briefly at each location and resealed in the field. Sample bottles were tightly sealed, wiped with a 10% bleach solution and immediately placed in a cooler containing wet ice for transport to the laboratory. Before being opened and filtered, sample bottles were rinsed with deionised water as in our previous studies (Olds et al. 2016, Evans et al. 2017).

After collecting the water samples, we recorded the latitude and longitude coordinates for each sampling location using a handheld global positioning system (Garmin, Olathe, Kansas, USA). We then measured turbidity, dissolved oxygen (D.O.), pH, water temperature and flow velocity at each site using the cross-sectional approach (Gore and Banning 2017). Surface water turbidity was quantified using a portable turbidity meter (Hach, Loveland, Colorado, USA). Dissolved oxygen, pH and water temperature were quantified using a multi-parameter probe (Yellow Springs Instrument, Ohio, USA). The flow velocity at each site was measured at both the left and right edge of the wetted channel and the thalweg using a Marsh-McBirney Flo-mate flow meter and wading rod (Hach, Loveland, Colorado, USA). Discharge was calculated as mean velocity × width × depth for each sampling site. The distance between the sampling locations was calculated as the linear stream distance between the locations following the stream channel midpoint.

Figure 1.

Map of sampling locations along Eagle Creek. Location 1 is the most upstream sampling location in Eagle Creek. Location 6 is the most downstream sampling location in Eagle Creek. Locations 7 and 8 are in the Kalamazoo River. Location 7 is upstream of the confluence of Eagle Creek whereas Location 8 is downstream of the confluence.

In the laboratory within 6 hours of sample collection, water samples were vacuum-filtered on to 47 mm, 1.2-µm pore size polycarbonate membrane filters (EMD Millipore, Billerica, Massachusetts, USA). The filters were placed in 2.0 ml microcentrifuge tubes containing 700 µl of CTAB and stored at -20 °C as described by Evans et al. (2017). DNA was isolated following a modified Chloroform-Isoamyl alcohol (24:1, Amresco, Dallas, Texas, USA) extraction and an isopropanol precipitation (Olds et al. 2016; Renshaw et al. 2015). To remove potential PCR inhibitors, re-suspended DNA was treated with the OneStep™ PCR Inhibitor Removal Kit (Zymo Research, Irvine, California, USA).

To check for and quantify contamination during the experimental process, we included a series of field controls, mock communities and PCR negative controls (Suppl. material 1). Eight field controls (one for each sampling location) were extracted separately from the eDNA samples and mock community samples, amplified and visualised on agarose gel to detect potential contamination in the field. Four mock community samples were constructed independently and processed along with eDNA samples from DNA extraction to sequencing to detect any potential contamination during the experiment. To avoid the mock community species affecting the results, the mock community samples comprised 60 ng (10 ng per species) of tissue-derived DNA (measured with Qubit) from six Indo-Pacific marine fishes: Amphiprion ocellaris, Salarias fasciatus, Ecsenius bicolor, Centropyge bispinosa, Pseudanthias dispar and Macropharyngodon negrosensis. Mock communities 1 and 2 were processed along with samples from Locations 1-4. Mock communities 3 and 4 were processed along with samples from Locations 5-8. A no-template PCR negative control (sterile water) was used for each of the three gene regions amplified.

The 3 replicates from each sampling location (24 samples in total), 4 mock communities and a single PCR negative control were PCR amplified for 3 mitochondrial gene regions. The three PCR primer sets were used in independent reactions to amplify segments of the mitochondrial 16S rRNA gene (L1865/H2195, hereafter Ac16s), 12S rRNA gene (L909/H1155, hereafter Am12s) and Cytochrome B gene (L14912/H15149c, hereafter CytB) (Evans et al. 2016, 2017; Olds et al. 2016). Amplified gene fragments were prepared for Illumina sequencing following a two-step PCR-based approach as outlined in the Illumina 16S Metagenomic Sequencing Library Preparation guidelines (Illumina, Inc., San Diego, CA, USA).

Amplified DNA from the first PCR were electrophoresed through a 2% agarose gel, stained with ethidium bromide and then visualised on a UV light platform. Each amplified product was manually excised from the gels using single-use razor blades, cleaned with a QIAquick Gel Extraction Kit (Qiagen, Venlo, Netherlands) and eluted from spin columns with 30µl of Buffer EB (Qiagen, Venlo, Netherlands). Regardless of visual confirmation of amplification, for both the mock communities and PCR negative control, we excised a band from the agarose gel at the expected amplicon size. For each sample, the gel extracted PCR products from the three gene regions were combined in proportions determined from previous work with these markers to approximate equal read production across genes (Am12S = 2.5 ng; Ac16S = 3.75 ng; CytB = 18.75 ng). In a second PCR, the Illumina adaptor sequence was added, including individualised library barcodes per sample. The DNA concentration of each library was quantified via Qubit dsDNA HS Assay. Libraries were pooled in equal molar concentrations along with 25% PhiX DNA (v3, Illumina, San Diego, California, USA) and then 300-bp paired-end sequenced on an Illumina MiSeq in two MiSeq flow cells in the University of Notre Dame’s Genomics and Bioinformatics Core Facility (http://genomics.nd.edu/) using a MiSeq Reagent Kit v3 (600-cycle; Illumina, San Diego, California, USA).

Bioinformatics analysis followed Olds et al. (2016) (scripts available on GitHub: https://github.com/pfrender-laboratory/epps). Briefly, raw sequence reads were removed from subsequent analysis if the read matched MiSeq sequencing adaptors > 6-bp, read quality score was smaller than 20 with a sliding window = 10bp or read length was shorter than 50-bp using Trimmomatic version 0.32 (Bolger et al. 2014). Filtered paired-ended sequences were merged with USEARCH version 8.0.1623 (Edgar 2010). Identical merged reads were collapsed into one unique read using a custom Perl script (scripts available on GitHub: https://github.com/pfrender-laboratory/epps). The number of reads collapsed in unique reads was recorded. Unique reads were then clustered at 97% similarity with USEARCH version 8.0.1623 (Edgar 2010) to establish Operational Taxonomic Units (OTUs). The number of unique reads and their identical reads in each OTU was recorded and used as the abundance of the OTU. To accelerate the species assignment process, profile-hidden Markov models (HMMs) were used to remove sequences that are not of metazoan origin. A HMM model was constructed for each targeted gene region based on sequences of metazoan species from NCBI RefSeq Release 69 using HMMER version 3.1b1 (Wheeler and Eddy 2013). The centroid sequence of each OTU was mapped to the profile hidden Markov models to remove sequences that are not of metazoan origin. To establish the taxonomic affiliation of the OTUs, the centroid sequence of each OTU was assigned to species with SAP v1.9.4 (Munch et al. 2008) using NCBI NR database (version of April 2016) based on 97% sequence similarity to references (-minidentity 0.97). OTUs of mammal or bird origin were removed.

No observable bands (PCR products) were detected on agarose gels in the expected size ranges for the eight field negative controls. Therefore, the eight field negative controls were removed from further analysis. After species assignment, potential contamination was removed based on the threshold correction method (Hänfling et al. 2016; Evans et al. 2017). In brief, the relative read abundance was calculated for each species found in each genetic marker in the control samples (four mock communities and one PCR negative control). The relative read abundance then functioned as a minimum detection threshold for each marker and species. Species in the samples detected at a lower relative read abundance than this threshold were removed from subsequent analysis. The script for the threshold correction method can be found on Github (https://github.com/lyy005/threshold_correction_for_eDNA/). Mock communities 1 and 2 and the PCR negative controls were used to estimate the contamination thresholds for samples from Locations 1-4. Mock communities 3 and 4 and the PCR negative controls were used to set the contamination thresholds for samples from Locations 5-8.

To test whether the biodiversity patterns varied significantly between gene regions, an unfiltered and a filtered dataset were used for each region. A species (OTU) was considered detected if its DNA was found in any sample replicate or any marker from a site. This single occurrence criterion was used in the comparison of different gene regions since the “at least two markers” criterion in the moderate bioinformatics stringency (Evans et al. 2017) for species detection (see below) was not applicable for single gene regions. The unfiltered datasets used all the OTUs, including those not identified to species. In contrast, the filtered datasets included only the OTUs that could be identified to species-level. By this definition, inferences of species diversity will always be more conservative for the filtered than the unfiltered datasets because the filtered dataset is always a subset of the total OTUs detected.

Based on the sampling effort (three replicates for each site), the potential biodiversity in a relatively small environment and detection thresholds from previous studies, we utilised a moderate bioinformatics stringency (Evans et al. 2017) for species detection, namely a species was considered detected if its DNA was found in a least two replicate samples from a site or two markers from a single sample at a site.

All analyses of α- and β-diversity were conducted in R version 3.3.2 (R Core Team 2016). The R-scripts can be found in figshare (https://doi.org/10.6084/m9.figshare.5860851). Genetic data were treated as incidence data. Alpha-diversity was calculated as species richness (number of OTUs assigned to species level) at each sampling location. Beta-diversity was estimated as Jaccard distance using the “vegdist” command in the “vegan” package version 2.4-2 of R (Oksanen et al. 2013). PCoA plots were made based on the β-diversity matrix with “cmdscale” command in R to visualise the relationships amongst replicate samples and sites.

The correlation between β-diversity matrices and the longitudinal distance matrix was tested with a Mantel test using “mantel.rtest” command in the R package “ade4” version 1.7-6 (Dray et al. 2007). Since Location 7 in the Kalamazoo River was upstream of the Eagle Creek confluence, it was not included in the correlation analysis between biodiversity and linear stream distance.

To evaluate the biodiversity detected by each replicate, the β-diversity between replicates was also calculated with “vegdist” based on Jaccard distance. Differences in β-diversities within and between the sample locations were tested using the Multi-Response Permutation Procedure (MRPP) in the R package “vegan” version 2.4-2 (Oksanen et al. 2013) using Jaccard distance.

The physical properties, including discharge and longitudinal distance from the headwaters of each sampling location along Eagle Creek are summarised in Suppl. materials 2 and 3, respectively. In general, discharge increased with distance downstream as expected. The headwaters of Eagle Creek showed the lowest discharge (Location 1 Q = 0.10 L/s) whereas the most downstream site in Eagle Creek had much higher discharge (Location 6 Q = 324 L/s). The two locations (Location 7 and 8) in the Kalamazoo River had the highest discharge (Q = 898 L/s).

Illumina MiSeq sequencing generated 20.9 million raw reads. After filtering, merging and quality control, 2.1 million, 3.9 million and 3.0 million reads were used for downstream analysis for Ac16S, Am12S and CytB markers, across the total dataset (Table 1). After clustering, 54, 93 and 36 OTUs were detected with Ac16S, Am12S and Cytb markers, respectively. In these OTUs, 43, 79 and 36 OTUs were identified as vertebrates with HMM filtering. After the threshold correction, reads from mock community species were detected in the Am12S region of replicates 1a and 1b in Location 1 (Table 1). Therefore, Am12S marker of replicates 1a and 1b were removed from subsequent analysis. No other laboratory contamination was found. After the data filtering steps, 29, 33 and 15 OTUs were identified as fish species for Ac16S, Am12S and CytB markers, respectively (Table 1).

The identifications of fish species were refined based on their known geographic distributions. For example, the Ac16S and CytB markers were not diagnostic between Umbra limi and Umbra pygmaea. Since only U. limi is known from Michigan, all OTUs with an assignment to an Umbra spp. were assigned to Umbra limi. Similarly, the three gene regions were not diagnostic for the two closely related species, Ambloplites ariommus and Ambloplites rupestris. As only A. rupestri is found in Michigan, Ambloplites spp. was assigned to A. rupestri. The Am12S marker could also not distinguish between Catostomus commersonii and Catostomus catostomus. Since both are found in Michigan, the identification of these species was grouped at the generic level as Catostomus spp. Notropis topeka were removed as these species are not found in Michigan or nearby regions (NatureServe 2014).

Three gene regions were used for estimating α-diversity. When including all the OTUs (i.e. unfiltered dataset), the Am12S region detected the highest number of OTUs. When including only species-level assignments (i.e. filtered dataset), Ac16S detected the highest number of OTUs (Table 1). Only four species were detected by all three gene regions (Figure 2). The number of species co-detected by two markers was the highest for the pair Ac16S and Am12S (10 species). The pair Ac16S and CytB co-detected seven fish species, while Am12S and CytB had the lowest co-detection of species at four.

Figure 2.

Venn diagram of number of species detected by each genetic marker based on filtered datasets.

Different gene regions showed similar β-diversity patterns. For both the filtered and unfiltered data sets, β-diversity estimates, based on single gene regions and three markers, were positively correlated with each other. All estimates of β-diversity and longitudinal distance were positively correlated except between the Ac16S-filtered estimates and distance (all the r statistics and p-values are available in Suppl. material 4). The β-diversity estimates generated from the three unfiltered datasets (Ac16S unfiltered, Am12S unfiltered and CytB unfiltered) and three filtered datasets (Ac16S filtered, Am12S filtered and CytB filtered) were significantly correlated with each other. For each marker, the unfiltered datasets had a better correlation with the combined three marker dataset than with the corresponding single marker filtered dataset (Figure 3).

Figure 3.

Heat map of correlations (Mantel r) between β-diversities as determined by the different genetic markers as well as the physical distance between samples. “Three_markers” refers to the combination of three filtered datasets (Ac16S_filtered, Am12S_filtered and Cytb_filtered).

Based on the combined signal of all three genetic markers, 23 fish species were considered detected in the 6 sampling locations in Eagle Creek and 2 sampling locations in Kalamazoo River (Table 2). From upstream in Eagle Creek (Location 1) to downstream in the Kalamazoo River (Location 8), the number of fish species at a sampling site increased from 2 species to 11 species (Table 2). Location 1 had the lowest fish diversity (two species). Locations 2 and 3 had an intermediate fish diversity with seven and five species, respectively. Two species (Notemigonus crysoleucas and Perca flavescens) were only found in these two locations in Eagle Creek. Locations 4, 5 and 6 showed greater fish diversity with 8, 8 and 10 species detected, eight of which were shared between at least two of three locations. Furthermore, Erimyzon sucetta was only found in Location 5 and Lepomis gulosus was only found in Location 6. The two sampling locations in the Kalamazoo River (Locations 7 and 8) showed the highest number of species detected (9 and 11 species detected). Noturus flavus and Percina maculata were only detected at Location 8, while Cyprinella spiloptera was only detected at Location 7.

Based on PCoA (Figure 4), the most upstream Location 1 had a distinct fish assemblage from other locations. Locations 2 and 3 displayed similar biodiversity, as did Locations 5 and 6 and Locations 7 and 8. Fish diversity in Location 4 was intermediate between Locations 3 and 5. The longitudinal stream distance between sampling locations and β-diversity were positively correlated (Figure 5). In addition to distance, other environmental variables (temperature, pH, turbidity) significantly varied with β-diversity (Suppl. material 5). These variables (temperature, pH, D.O., discharge) are not independent of distance (Suppl. material 5) or other environmental variables. Based on the β-diversity measures amongst all the replicates, replicates from the same location were more similar than replicates between locations with a MRPP test (p-value < 0.01) (the PCoA plot and table with species and their read number in each replicate can be found in Suppl. materials 6 and 7).

Figure 4.

PCoA of the β-diversity of sampling locations with Jaccard distance based on all three genetic markers.

Figure 5.

Relationship between the distance between sampling locations and β-diversity based on the combined results of all three genetic markers.

In this study, we used eDNA metabarcoding to describe fish diversity along the course of a small stream. The local fish assemblage varied spatially from upstream to downstream, consistent with previous studies (Magalhaes et al. 2002; Vannote et al. 1980). In stream headwaters, fish species diversity tends to be low relative to downstream sites where habitat diversity increases and species tend to accumulate along the river continuum (Tornwall et al. 2015). The most upstream headwater site of Eagle Creek had a low-diversity and distinctive fish assemblage. The two locations in the much larger Kalamazoo River had the highest number of fish species, consistent with predictions of the River Continuum Concept (Vannote et al. 1980).

Fish species assemblages also showed absence, additions and turnovers (as measured using β-diversity) along the stream gradient. Samples collected from adjacent locations often showed similar species compositions and dissimilarity increased with distance between sites. The differences in species composition along the stream could be due to changes in habitat structure and availability, as well as changes in temperature, pH and dissolved oxygen along the stream course (Almeida and Cetra 2016). The two sampling locations (Locations 7 and 8) upstream and downstream of the confluence of Eagle Creek, respectively, showed different fish species composition. Eagle Creek likely contributed eDNA to the downstream site of the Kalamazoo River. For example, Etheostoma caeruleum and Micropterus salmoides were detected in Eagle Creek and at Location 8 in the Kalamazoo River where water from both systems had mixed but not at Location 7 upstream of the confluence. The absence of some fish species at Location 8 that were detected at Location 7 may be related to the distance between sites or dilution by incoming water from Eagle Creek.

We collected three replicate eDNA samples in each location to improve species inference following our practice (Olds et al. 2016; Evans et al. 2017) and as suggested by Goldberg et al. (2016). Based on a MRPP test and the PCoA analysis, replicates collected from the same location displayed similar biodiversity, although some variation was detected. We collected a modest number of replicates (3 replicates of 250 ml water) for each sampling location. Although more replicates could improve the detection of low-abundance species (Evans et al. 2017), the low variance amongst the replicates indicates that additional replicates may have had a minor effect on our results. We also found a positive correlation between β-diversity amongst sample sites and stream distance, indicating our site sampling scheme generated sufficient statistical power to infer the patterns of fish community similarity within Eagle Creek. Our previous studies suggested that use of a moderate bioinformatics stringency (“at least two markers or two samples” criterion) can minimise the effect of primer bias (Evans et al. 2016), minimise cross-sample contamination (Olds et al. 2016) and improve the reliability of species identification (Evans et al. 2017). In this study, we found that the overall biodiversity patterns (β-diversity estimates) inferred by the different genetic markers were similar. This finding is consistent with previous research that β-diversity estimates are less sensitive to the specific genetic markers employed than β-diversity estimates (Drummond et al. 2015). By comparing the filtered and unfiltered dataset of each genetic marker to the three-marker dataset, we found that including OTUs that did not have a species-level assignment (unfiltered datasets) added more information about community similarities amongst samples. In general, our results suggest that different genetic markers express the same β-diversity pattern along the longitudinal stream gradient in Eagle Creek.

Understanding the spatial distribution of eDNA in flowing environments is critical for accurate detection of species in stream ecosystems (Deiner et al. 2016; Jerde et al. 2016; Shogren et al. 2017). Environmental DNA transported from upstream may lead to false positive detection of species in downstream locations and it becomes difficult to infer species location relative to sampling location.

Based on previous studies, the average transport distance of eDNA can be highly variable, depending on stream discharge (Barnes and Turner 2016; Jerde et al. 2016; Shogren et al. 2016, 2017; Wilcox et al. 2016). For example, studies in low discharge (Q = 2-10 L/s) experimental streams found that eDNA transport distance was 50-1459 m (Jerde et al. 2016; Wilcox et al. 2016), but eDNA has been detected up to 10 km downstream in high discharge systems (Q > 3,000 L/s, Deiner and Altermatt 2014). In this study, we estimated β-diversity between sampling sites and showed that these estimates were correlated with stream distance. The implication of this result is that the downstream transport of eDNA in this low discharge stream is insufficient to completely homogenise the signal of fish diversity at individual sites. The low discharge in the headwaters of Eagle Creek likely limited the eDNA transport distance to less than the between site distances. For example, Locations 1 and 2 had no species overlap even though they were separated by only 829 m. These differences in species composition were probably not the result of low abundance and low detection probabilities since the most abundant species (based on sequence read number) in Location 1, Umbra limi, was not detected at Location 2. Given these results and the low discharges measured at Location 1 and 2 (0.1 L/s and 73.5 L/s respectively), it is likely that the DNA transport distance is less than 829 m in the headwaters of Eagle Creek. Sampling eDNA in low discharge streams is potentially a useful way for managers to measure diversity for more local scales. Higher-discharge downstream locations in Eagle Creek shared species including most of the high-abundance species (Table 2) so we cannot distinguish whether the detected species’ eDNA at any particular locality was derived from eDNA transport or the mobility of species (Pilliod et al. 2013).

Our goal was not to mechanistically describe why β-diversity changes in the river, but to test whether transport of eDNA obscures the signal of species turnover in lower discharge streams. The influence that environmental variables (turbidity, D.O., temperature etc.) have on the transport and detection of eDNA is not well understood, although, increasingly, studies are finding evidence that environmental variables can affect how long eDNA remains detectable in a system (Seymour et al. 2018; Shogren et al. 2017). Additional mechanistic tests of such variables are needed to understand when and how the environmental conditions in streams influence the detection and transport of eDNA (Deiner et al. 2017). By understanding the physical and ecological mechanisms of eDNA behaviour in lotic environments, we can better estimate the presence of species in rivers and streams using eDNA metabarcoding.