Molecular identity crisis: environmental DNA metabarcoding meets traditional taxonomy—assessing biodiversity and freshwater mussel populations (Unionidae) in Alabama

Study area

We selected six sites on the Sipsey River (Fig. 1) where freshwater mussels had been sampled previously (Pierson, 1991; Haag & Warren, 2010; Atkinson & Forshay, 2022). In addition to sampling mussels with eDNA and traditional methods, we quantified environmental characteristics of the study area. We derived watershed areas for each sampling point using the Spatial Analyst Toolkit in ArcMap 10.6 (Environmental System Research Institute, Redlands, CA, USA) with a 30-m digital elevation model (DEM) from the National Elevation Dataset. We obtained land cover (30-m resolution) for the United States from the 2011 National Land Cover Database. We measured water temperature, pH, and conductivity with a YSI multiparameter probe. In addition, we filtered water samples through a Gelman A/E, GFF, 0.7-μm nominal pore size to determine water chemistry and placed them on ice. We determined NH₄-N and soluble reactive phosphorus (SRP) with a Lachat Quikchem 8000 flow-injection auto-analyzer using colorimetric methods (Lachat Instruments, Milwaukee, MN, USA). See Table 1 for a summary of study site characteristics.

Traditional surveys

We conducted quantitative mussel surveys in 60–80 m reaches between July-October 2016 across the six sampling sites. We first visually determined the extent of each mussel aggregation by snorkeling and then divided each total reach into 20 m sections (3–4 sections per site) in which three random numbers were selected to determine our transects. Transects were run perpendicular to stream flow within each 20 m section. Across each transect we excavated the sediment from a 0.25 m² quadrat to a depth of 15 cm every 1.5 m (6–20 quadrats per transect) retaining all the material dug into a mesh bag. To determine burial preference, prior to excavating the quadrat we lightly moved our hand over the sediment surface and recovered mussels that appeared at the surface and were not completely buried and placed them into a separate small bag that was placed within the quadrat’s mesh bag. Each bag was individually labeled, and all mussels were sorted through a series of sieves with species identification (Williams, Bogan & Garner, 2008) and lengths recorded for each individual. This sampling regime resulted in us sampling approximately 2.5% of the total area of the designated reach and we determined mussel species richness, density (indiv/m²), and estimated biomass using published length-dry mass relationships (see Atkinson et al., 2020) at each site.

We used fish community data collected by the Geological Survey of Alabama (GSA) and the Alabama Department of Environmental Management (ADEM) in 2011–2013 in the Sipsey River and major tributaries to evaluate fish species composition and abundance. Sampling combined seining and backpack electrofishing following the methodology of the Index of Biologic Integrity in the Hills and Coastal Terraces Ichthyoregion in Alabama (O’Neil & Shepard, 2011, Bearden et al., 2019). Briefly, surveyors performed a minimum of ten sampling efforts at each distinct type of habitat; riffle, pool, run or shoreline at each site. Freshwater mussel collection was conducted under US Fish and Wildlife Service permit #TE68616B-1 and Alabama Department of Conservation and Natural Resources permit #2016077745468680 to C.L. Atkinson.

Water sampling and eDNA extraction

We collected water samples for eDNA analysis at the six sampling sites on four occasions from August to November 2016 when we expected flow and turbidity to be low and when we hypothesized that eDNA would be most concentrated. We collected 4, 1-L samples of water at each site and sampled all sites in a day from downstream to upstream to avoid contamination. We repeated this sampling four times with each occasion separated by >2 weeks: August 3, August 17, October 12, and November 9, 2016. Conducting repeated visits separated by time allowed us to capture eDNA across a range of environmental conditions and through life history milestones (e.g., reproductive periods) of some of the target Unionidae community. This information is key to quantifying conditions that influence detection of mussel species from eDNA-based surveys (Schmidt et al., 2013).

Some sites were accessible on foot and others required a short canoe paddle. To prevent cross-site contamination, we entered sites from a downstream direction, and we decontaminated equipment (e.g., bottles, tweezers, waders) with a 50% bleach solution. With a gloved hand, each bottle was labeled, opened, and filled in the river at arm’s length from the shore, downstream but as close to the mussel quadrat surveys as possible, with the 1-L bottle mouth facing upstream in the center of the water column. Water depth at each site varied from 16 cm to 60 cm so the collection was anywhere from approximately 8 cm to 30 cm above the benthos. Water samples were kept in coolers on wet ice for transport back to the lab, where they were placed in a standard refrigerator until filtering.

We filtered all four samples collected at each site and date with Cellulose Acetate filters but varied the pore size such that two samples were filtered with 3 µm and two with 5 µm pore size for each collection. Methods to filter water samples in the field were not feasible in our study area because of the high levels of turbidity that characterize many Gulf Coast river systems. Instead, we developed a vacuum pump manifold system set to 20 PSI in a laboratory setting to filter water samples within 24 h of collection to minimize degradation of DNA. Numerous precautions were incorporated into laboratory practices to minimize risk of contamination between samples, including using new gloves between samples, and cleaning equipment and the workstation with bleach prior to filtering and between samples. A negative control water sample was also run through the filtering process and included on the array. After vacuum filtration, each filter was carefully rolled with a gloved hand and placed in a labeled and sealed vial and frozen prior to shipment on dry ice to the United States Department of Agriculture (U.S.D.A) Forest Service Pacific Northwest Research Station in Corvallis, Oregon for DNA extraction and subsequent analysis. We conducted pilot tests on other nitrocellulose filter pore sizes of 0.45 µm (Hauck et al., 2019) and 1.2 µm but found that both were too fine to allow filtering of water from our study area to be completed within 24 h. Even with the vacuum pump, the turbid water samples clogged the filters and stopped the filtering process. DNA was extracted and processed as described in Hauck et al. (2019). Briefly, the eDNA was extracted using the inhibitor removing Power Water DNA extraction kit (Qiagen, Hilden, Germany) per manufacturer’s instructions and the samples were cleaned and concentrated post-extraction with the Zymoclean© Large Fragment DNA Recovery Kit (Zymo Research, Irvine, CA, USA).

Microfluidic PCR amplification and massively parallel sequencing

Microfluidic high throughput PCR amplification of the eDNA was performed on the Fluidigm 48.48 Access Array™ (Fluidigm, San Francisco, CA, USA; Fluidigm Corporation, 2016) as previously described in Hauck et al. (2019). Specifically, the Access Array uses integrated fluidic circuits and a 4-primer amplicon tagging scheme in which target-specific primer pairs amplify 48 different targets in combination with sample-specific barcoded primer pairs in 48 different samples allowing for the simultaneous amplification of barcoded targets in each of the 2,304 individual reaction chambers on each array. Many of the taxon-general universal metabarcoding primers we used in Hauck et al. (2019) were designed to amplify diverse classes of organisms (e.g., ray-finned fishes (Teleostei) and amphibians (Batrachia), Chondrostei fish, mussels (Bivalvia), and insects), by targeting the 12S rDNA region used by Valentini et al. (2016), as well as 16S, 18S, and internal transcribed spacer (ITS) rDNA. Due to their universal nature, they were used again in this study. Additional primers were designed to target organisms found in the Gulf of Mexico coastal rivers. The new primers were developed as described in the previous study (Hauck et al., 2019), and targeted mostly freshwater mussels, crayfishes, warm water fish, and turtles using mitochondrial genes commonly used for eDNA metabarcoding purposes: cytochrome B (CytB), 12S and 16S ribosomal genes, Cytochrome Oxidase Subunit I (COI), NADH dehydrogenase Subunit 1 (ND1), subunit 2 (ND2), subunit 4 (ND4), and Subunit 5 (ND5). See Table S1 for all primer sequences and the associated species and gene targets used in this study. We determined previously through test PCR reactions and preliminary Fluidigm experiments that the optimal Fluidigm submission concentration for eDNA is 12–15 ng/µl; unfortunately, many of these samples did not contain enough DNA to be submitted at the optimal concentration. As a result, we submitted the DNA samples in the following categories: 12, 6, 3, 3 > 1, or <1 ng/µl (see Table 2 for final sample submission details).

The negative control consisted of samples N1R3 (plate 1) and N4R3 (plate 2), which are ultrapure water carried through the sample filtration process. Our positive controls followed the methods in our previous study (Appendix 1, Hauck et al., 2019) with the goal of 2 × 10⁴ molecules per amplicon per 31.35 nl volume (the size of a Fluidigm Access Array reaction, after accounting for primer and mastermix volume). Target species DNA was derived from organisms that were collected by field sampling, through other researchers, or from acquiring tissues from museum specimens. All samples were sent to the USDA Forest Service Pacific Northwest Research Station in Corvallis, Oregon for DNA extraction. Unfortunately, the resulting DNA was often degraded, particularly when derived from the museum samples which had been stored in formaldehyde for a significant amount of time. Amplification was still attempted, and a basic standard was composed of amplicons for those few species whose DNA we were able to generate amplicons from (see Table S2 for a list of positive control species). A modified positive was also added to one plate; it consisted of the basic standard spiked with the same copy number of amplicons generated for the metabarcoding primers used in both the previous study and this one.

DNA samples and primers were sent to the Roy J. Carver Biotechnology Center at the University of Illinois for amplification on their Fluidigm instrument and sequencing on their MiSeq instrument using 2X250 V2 chemistry. Amplification was done using a 2-step Fluidigm cycling protocol and a modified template-specific annealing temperature to 58 °C. Bovine serum albumin (BSA) was added at 0.2 μg/μl to alleviate PCR inhibition from contaminants in environmental DNA (Romanowski, Lorenz & Wackernagel, 1993; Widmer, Seidler & Watrud, 1996) and to mitigate inhibitor-driven bias (Valentini et al., 2016). Samples were barcoded with 10 basepair (bp) indices during amplification and then pooled for gel cleanup and size selection prior to sequencing. All raw sequence reads from this dataset are accessible to the public, they have been uploaded into the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) and can be accessed by searching for BioProject Accession Number PRJNA902311.

Bioinformatic analysis

Sequences were returned from the Roy J. Carver Biotechnology Center demultiplexed by index and by primer using a custom pipeline they designed to sort Fluidigm data. For this run, they allowed one mismatch in the index, and two mismatches in the primer. Due to anomalies found in the read two data during joining of the paired-end sequences, only ‘read 1’ sequences were used in subsequent analyses. Trimming of the ‘read 1’ sequences was done in primer groups using Trimmomatic v. 0.33 (Bolger, Lohse & Usadel, 2014). Using amplicon specific information, the 5’ end of the read was trimmed the length of the primer plus five nucleotides and the overall length of the read was cut off at 200 bp to avoid analyzing the 3’ end of the sequence where quality scoring of read one sequencing began to drop off. The taxa assigning pipeline that we used for this dataset is KMA version 1.2.23, which uses a kmer-based conclave scoring method (Clausen, Aarestrup & Lund, 2018). Once run through KMA, the data was then piped through CC-Metagen version 1.2.2 (Marcelino et al., 2020) for the final assignment of taxonomy to the read. The CC-Metagen pipeline was originally developed to assess gut microbiome and fungal populations using singular, short reads with high conservation within species. Because our data consists of many different gene targets of varying conservation, length, and taxonomic origin, we needed to modify the thresholds to better capture species richness. We found the following settings best fit this particular data set: class = 79.0, order = 80.0, family = 84.0, genus = 91.4, and species = 94.4. Any reads showing up in the negative control may be the result of index hopping, sequencing error in the indices or contamination. We used the negative control to calculate a 95th percentile minimum read threshold. Any reads at or below this threshold were removed from the dataset. Counts above the threshold were considered “positive hits” for that taxon in that sample × primer combination. A Krona chart for hierarchical viewing of all counts above the threshold for all sample sites was generated using a krona excel template (Ondov, Bergman & Phillippy, 2011).

Phylogenetic analysis

Based on the classification results of our eDNA data we decided to investigate the reads from unexpected species detections. Species assignments of Fusonaia cerina and F. flava and Pleurobema decisum and P. chatanoogaense, and Elliptio sp. eDNA reads from our dataset were investigated through phylogenetic analysis at the mitochondrial ND1 gene. ND1 sequences in our dataset were derived from multiple primers targeting different, and sometimes overlapping, fragments of ND1 ranging from 99 bp to 220 bp. Our reads reflect consensus sequences pulled from KMA alignments, and they are named by sample and preliminary KMA alignment classification. All available NCBI reference sequences for P. decisum (four), P. chatanoogaense (two) and F. cerina (seven) were downloaded for inclusion in alignments. F. flava had a very large number of sequences in NCBI, so seven representative sequences were selected with a preference for vouchered specimens. In addition, to provide an outgroup and guide tree branching, up to six representative sequences were pulled from NCBI for the following taxonomic groups within Ambleminae: Quadrula (outgroup), Obovaria, Lampsilis and Reginaia.

All sequence alignments and trees were created using Geneious® software version 10.2.6 created by Biomatters (available from https://www.geneious.com). Alignments were conducted using MUSCLE with a maximum number of 100 iterations with first iteration kmer4_6, and terminal gaps penalties were turned off to accommodate multiple sequence lengths in the alignment. Prior to conducting an alignment our consensus sequences were screened for residual primer or adapter sequence and trimmed if necessary. Maximum likelihood trees were drawn from the alignments using RAxML 8.2.11 with the nucleotide model GTR CAT I, which includes an estimate of proportion of invariable sites. The tree algorithm used was ‘Rapid Bootstrapping’ and ‘search for best-scoring ML tree’ with 1,000 random replicates to generate the final tree and bootstrap values. The tree was rooted with the Quadrula sequences. As unionid mussels of the Quadrulini tribe, not the Pleurobimini tribe we were investigating, they are a valid outgroup for tree rooting.

Data analysis

We used summary statistics and graphics to discuss observed differences in species detection between our two sampling methods. This was necessary due to the dissimilar structures of the data sets themselves, the relatively small sample size of our data, and differences in the timing of data collection between the traditional samples and the eDNA samples. Traditional survey data was summarized to the site and species level. Our six sites used different numbers of quadrats during sampling (range 103 at Site 3 to 206 at Site 6) as a result of the differences in stream width and mussel aggregation length; therefore the total number of individual mussels encountered as well as total biomass was normalized to a per m² value. We summarized eDNA results by site when comparing the eDNA detections directly to the traditional survey results. We also summarized the eDNA data by sample, replicate and filter size to evaluate the total number of unique taxonomic identities detected and the number of reads for each sample; and did comparisons on the amount of DNA submitted (as a categorical variable) to evaluate the impact of DNA quantity on the total number of reads and unique taxonomic identities returned.

Mussel diversity and density from traditional quadrat surveys

Across our six sites, we sampled 779 quadrats and measured 6,195 individual unionid mussels in 2016, observing 17–29 species per site (Fig. 2). Unionid community densities (including all species) ranged 15.1–38.8 individuals per m² and community biomass ranged 6.6–15.9 g per m². The most abundant species encountered were Pleurobema decisum, Cyclonaias asperata, and Fusconaia cerina, respectively (Fig. 2).

eDNA yield and metabarcoding results

The yields from eDNA extractions were rather low, and despite having larger pore sizes, the 5 µm filters consistently yielded more DNA than the 3 µm filters. Due to the poor yields, Site 3 is missing one of the 3 µm filter replicates for each collection time and is underrepresented in the overall eDNA pool. Site 1 is also missing one of the 3 µm filter replicates for the fourth collection point but is fully represented at the other timepoints. In addition, many samples were submitted below the optimal concentration of 12–15 ng/µl (Table 2). This likely led to increased missingness in our dataset and will require trouble shooting for future eDNA studies in turbid water systems.

Amplification of the standard positive control failed, and limited amplifications occurred in the modified positive control. This limits our ability to define false negative results. The only primer pair that produced reads in the negative control was a highly efficient pair targeting Phytophthora, a genus of plant pathogens. The calculated minimum threshold was five reads. To avoid false positives any primer x taxon x sample combination that had five or fewer total reads was dropped from the data set. After processing, 1.76 million reads were classified and included in the data set (1.67 million reads with controls removed). The KMA/CC-Metagen pipeline for classifying organisms yielded very few unexpected results. Of those classified reads, only 57 reads were identified as non-native fish in the genus Oncorhynchus, and 714 reads classified to salamanders found only in the Pacific Northwest. The rest of the reads were classified as organisms, or closely related organisms, that we would expect to find in a Southeastern US river. The source of the salamander reads likely resulted from barcode swapping or sorting errors with another array that was run concurrently. Since there were no reads from these species in the negative control, contamination was unlikely. We removed the salamander reads from the data set.

We detected 126 different species with this array, including taxa representing unionid mussels, fish, Phytophthora, fungi, amphibians, arthropods (insects and crustaceans), cyanophora, oomycota, apicomplexia, and bacteria. Additional detections were only taxonomically classified to family (21), genus (33) or order (2). We did not detect any turtles despite targeting them directly with six primer pairs and known occurrences of numerous species of freshwater turtles in the watershed. Not surprisingly, when all detections are summed, 99.6% of the reads are composed of six taxonomic groups. Of the 48 primer pairs that we deployed on the array, 18 targeted mussels and 15 targeted fish, which are the two taxonomic groups that we detected with the greatest number of reads (Fig. 3). The single Phytophthora primer pair captured the next highest number of reads. All read classifications summarized by site can be viewed in an interactive Krona Graph (File S1), see Fig. 4 to view a static image of the Krona Graph displaying a summary of all classifications for all sites.

We did not observe strong patterns for DNA yield by site. There was a high variability between samples and overall number of reads that made it through the classification pipeline. This included numerous samples submitted at high concentration returning a low number of reads, and low concentration samples returning a relatively high number of reads. This may be a reflection of the overall post-extraction quality of the eDNA sample. We observed that samples filtered using the 5 µm filter often resulted in more DNA for submission and tended to have higher classified read returns (Fig. 5). The number of classified reads per sample ranged from 0 to 96,996 with an average of 18,442; and the number of unique taxonomic identities per sample ranged from 0 to 31 with an average of 5.25. The average number of unique taxonomic identities coming from the samples submitted at the optimal concentration of 12 ng/µl was 6.58, which was higher than the average overall (Fig. 5). Samples submitted at less than 1 ng/µl underperformed on average.

Comparison between traditional surveys and eDNA detections

Fish

We had 402,493 reads that classified as actinopteri with 26 species and nine genus level detections representing eight different families of ray-finned fishes. Of those reads, approximately 70% originated from universal primers that targeted the highly conserved 12S or 16S regions. These universal primers led to some classification issues where most of those reads were classified to family or genus but not to species. Most of the detections came from the following fish families: Cyprinidae, Centrarchidae, Ictaluridae, and Percidae. We loosely compared some of our fish detection results from eDNA to traditional surveys conducted 3–7 years prior by the Alabama Geological Survey (GSA), 2011–2013 (Bearden et al., 2019). The GSA survey occurred at two sites at or near to our study area. Our Site 6 corresponds with GSA Site SR14, separated only by a bridge. GSA Site Turkey Ford was approximately 4 km upstream from Site 1. Traditional surveys identified the Rock Darter Etheostoma rupestre in the Sipsey River whereas eDNA made a detection at the genus level for Etheostoma. Yellow Perch Perca flavescens, Bigmouth Buffalo Ictiobus cyprinellus, and two madtoms (Notorus munitus and N. stifmosus) were detected with eDNA, but not captured earlier in the traditional surveys (Fig. 6).

Mussels

In general, the detection of mussel species with eDNA approximated the assemblages recorded by traditional sampling methods with higher eDNA detections for mussel species with larger biomass and greater densities as observed in visual quadrat sampling (Fig. 7). We note, however, that no species had complete congruence across all sites. The species detected at the most sites using eDNA were Pleurobema decisum (five out of six sites), Fusconaia cerina (four out of six sites), Obovaria unicolor (four out of six sites), and Cyclonaias asperata (four out of six sites). Similarly, the most abundant species observed in the field surveys were P. decisum, C. asperata, and F. cerina, respectively (Fig. 2). Comparing our traditional samples to eDNA samples, we had 30 unique site/species combinations where the same species were identified by both sampling methods at the same site. Site 4 had the most overlapping species detections with eight species detected by eDNA which is 35% of the 23 observed species (Elliptio arca, Elliptio crassidens, Fusconaia cerina, Lampsilis straminea, Obovaria unicolor, Pleurobema decisum, Cyclonaias asperata and Truncilla donaciformes). Site 3 followed with detection of seven overlapping species, 44% of the 16 field observed species (Lampsilis ornata, Ligumia recta, Megalonaias nervosa, Obovaria unicolor, Pleurobema decisum, Cyclonaias asperata and Villosa lienosa). We detected six overlapping species at sites 1 (32% of the 19 field observations) and 2 (35% of the 17 field observations) (Elliptio crassidens, Fusconaia cerina, Lampsilis ornata, Lampsilis straminea, Leptodea fragilis and Pleurobema decisum) and (Cyclonais asperata, Fusonaia cerina, Ligumia recta, Obovaria jacksoniana, Pleurobema decisum and Tritogonia verrucosa) respectively. The sites with the poorest detections were Site 5 and Site 6, which only detected three of the expected species each (Fusconaia cerina, Obovaria unicolor and Pleurobema decisum) and (Cyclonaias asperata, Lampsilis ornata and Obovaria unicolor) respectively. These detections encompass only 16% of the 19 observed species at Site 5, and 11% of the 27 species observed at Site 6 and are substantially lower than the detection rates from the other four sites.

During the traditional survey, we recorded the location of mussels as they were observed (surface vs. buried). We hypothesized that a higher number of mussels on the surface would yield greater DNA signal in the sample, but the results are inconsistent and indicate that factors other than density at the surface are influencing species detection (Fig. 8). For example, Leptodea fragilis at Site 1 had 100% of the mussels at the surface but only 26 reads counted; however, the average biomass for this mussel at the same site is exceptionally small at 0.58 mg/m². The species with the highest eDNA detections overall, Pleurobema decisum, is often buried with an average of ~30% exposed at the surface; however, this species also has the largest average biomass and average density overall. This indicates that biomass and density may influence detection more than substrate location.

Many eDNA detections are classified to the taxonomic level of family (Unionidae) or genus (e.g., Venustaconcha, Fusconaia). It is possible that some of the missing detections are located at higher taxonomic levels and that insufficient information was captured or available in the database to bioinformatically assign these reads to species. We also noticed that our eDNA classifications returned taxa that may have changed nomenclature per Williams et al. (2017) or are not known to exist in the Sipsey River or the Mobile Basin and were not detected by traditional sampling. This was particularly evident in the Pleurobemini tribe containing Fusconaia and Pleurobema taxonomic groups. For example, within Fusconaia, F. cerina was the only species of that genus observed with traditional sampling; however, eDNA detected other Fusconaia species with multiple markers (ND1, COI). Of particular note, F. flava was detected at all sites with eDNA, but not quadrat sampling. Within Pleurobema, the most common species observed was P. decisum; however, eDNA data reflected a high number of P. chattanoogaense reads at four sites with multiple markers although not mentioned in the traditional survey. While it is possible that traditional methods missed some mussel species or misidentified them, it is more likely that our eDNA classifications generated the discrepancies between quadrat observations and eDNA detections.

Other species detections that were surprising included high numbers of sequences at Site 3 and Site 4 that most closely matched the genus Venustaconcha. To the best of our knowledge, no species in that genus has been recorded in the Sipsey. Since the match was not at a close enough% identification to be called at the species level those reads will remain undescribed. We also detected Lamspsilis satura and Elliptio ictarina, with moderately high numbers of reads at Site 4, and neither species was observed in the field. These may be misclassified detections of other species in that genus that were observed in the field.

Phylogenetic analysis of pleurobema, fusconaia and eliption species

We further examined sequences classified to the Fusconaia, Pleurobema, and Elliptio taxonomic groups from our ND1 eDNA pool to evaluate congruence between eDNA-based detections and field identification. By adding reference sequences from NCBI to our eDNA sequences, we constructed a maximum likelihood tree that resolves the two prominent Fusconaia species in a clade and the two prominent Pleurobema species in a clade (Fig. 9).

Our eDNA sequences for P. decisum and P. chattanoogaense resolved in clades with reference sequences from both species interleaved together with 93% bootstrap support. The resolution of sequences from these two species into two clades indicates that we captured two haplogroups that are represented in both species. The F. flava and F. cerina sequences also resolved into a single clade, with the reference sequences of both species interleaved together, although with lower bootstrap support (71%). The topology of Fusonaia sequences reflect three potential haplogroups, but as with the Pleurobema group, these are not well supported by bootstrap values. The mean intraspecific distances for the Pleurobema and Fusonaia species ranged from 1.15–1.77%; while the interspecific distance between P. decisum and P. chattanoogaense was 1.84% and between F. cerina and F. flava was 1.27% (Table 3).

We also generated alignments and maximum likelihood trees using consensus sequences identified as Elliptio sp. from our data. There was insufficient resolution at the locus we amplified to identify those reads to the level of species. Our amplified sequences did not group with E. arca or E. crassidens (both of which was observed in the field), or with any other closely related Elliptio species with sequences available in the NCBI reference database. Instead, they formed a separate clade closest to E. arctata which has not been seen in the Sipsey for many years. This indicates that the Sipsey River taxon represented by our eDNA sequences is not represented by a sequence in NCBI for an identification match. The failure to match E. crassidens or E. arca from NCBI also indicates that (a) eDNA primarily sampled another Elliptio species besides E. crassidens or E. arca, (b) E. crassidens or E. arca has greater haplotype diversity than is represented by NCBI, or (c) that there is an Elliptio sp. we detected in the Sipsey River that is a different taxon than those observed in the field study.