Introduction
Over the past 20 years, environmental DNA (eDNA) metabarcoding with high-throughput sequencing (HTS) has become a powerful tool for capturing the diversity of free-living species (Hupało et al. Reference Hupało, Majaneva, Victoria, Lucas, Marquina, Lijtmaer and and Deiner2021; Meeus et al. Reference Meeus, Silva-Rocha, Adriaens, Brown, Chartosia, Claramunt-López and and Groom2023; Pollock et al. Reference Pollock, Howe, Irizarry, Lorusso, Kruger, Himmler and Struwe2015; Ruch et al. Reference Ruch, Karns, McMurray, Moore-Palm, Murphy, Namestnik and Roth2010). However, this approach has only recently emerged in the field of parasitology, with few assessments of individual parasite groups and even fewer characterizing biodiversity across diverse parasite communities (Hupało et al. Reference Hupało, Blasco-Costa, Trujillo-González, Leese, Smit and Sures2025). There are many benefits—both demonstrated and anticipated—from parasite community eDNA metabarcoding. Environmental sampling allows discovery and estimation of parasite diversity without having to collect, euthanize, and dissect hosts. The resulting sequences provide a characterization of diverse communities without the need for taxonomic experts for each parasite group, though this benefit relies on comparisons with reference libraries populated by sequences generated by traditional parasite surveys and descriptions. eDNA metabarcoding data can help taxonomists and ecologists generate a more comprehensive diversity assessment, especially for species with low prevalence and abundance that are seldom detected during necropsies. Furthermore, the relative ease of eDNA surveys enables resampling over time and thus, facilitates monitoring and early detection of changes in parasite communities, including emergence and losses that may occur in ecosystems either naturally (e.g., after hurricanes) or due to anthropogenic impact (e.g., oil spills, dam emplacement/removal).
The lag in the development of eDNA metabarcoding workflows for parasites can be attributed mainly to the fact that ‘parasites’ are a functional, not a monophyletic group, with high genetic heterogeneity. Therefore, there are no universal primers to detect parasite biodiversity that encompasses taxa from Chromista to Animalia (Blasco-Costa et al. Reference Blasco-Costa, Cutmore, Miller and Nolan2016; Poulin et al. Reference Poulin, Hay and Jorge2019). However, metabarcoding primers have been developed for specific parasite groups—mainly arthropods, helminths, fungi, microsporidians, myxozoans, and protists—targeting different genetic markers, including nuclear small subunit (18S) and internal transcribed spacer region (ITS) ribosomal RNA (rRNA) and mitochondrial 12S and 16S rRNA and cytochrome c oxidase I (COI) genes (reviewed in Hupało et al. Reference Hupało, Blasco-Costa, Trujillo-González, Leese, Smit and Sures2025), but these have not been applied widely. Primers used to barcode protists, nematodes, myxozoans, and microsporidians commonly target the fourth variable region (V4) of the 18S rRNA gene (e.g., Aivelo et al. Reference Aivelo, Harris, Cadle and Wright2018; Doliwa et al. Reference Doliwa, Grabner, Sures and Dunthorn2023; Lisnerová et al. Reference Lisnerová, Holzer, Blabolil and Fiala2023). Reference taxonomic data for the 18S rRNA gene is relatively rich, as this gene has been used for parasite identification since the beginning of the molecular era. In the case of helminths, primers that target the mitochondrial 12S and 16S rRNA and COI gene regions have been designed using an ecosystem-wide approach rather than focusing on a specific parasite taxon (Chan et al. Reference Chan, Saralamba, Saralamba, Ruangsittichai and Thaenkham2022a, Reference Buchnerb; Thomas et al. Reference Thomas, Milotic, Vaux and Poulin2022). Additional primers are being developed continuously, yet challenges remain, including the need to optimize metabarcoding workflows—from sample collection to data analysis—hence the value of field testing these when possible.
There is a fundamental need to understand the limitations and account for biases inherent in eDNA sampling methods. Template DNA in environmental samples is present in highly variable concentrations, with seasonal and daily fluctuations depending on the parasite’s ecology, the sample matrix (e.g., water, sediment), and habitat (reviewed in Hupało et al. Reference Hupało, Blasco-Costa, Trujillo-González, Leese, Smit and Sures2025). Therefore, the timing of sample collection would affect the probability of detection, especially for low-abundance parasites. Most aquatic parasites have free-living stages that may represent high DNA concentrations in the water column, either as whole organisms or extracellular DNA, whereas parasites without buoyant or motile stages might be more efficiently found in sediments as eggs (e.g., helminths) (Bass et al. Reference Bass, Christison, Stentiford, Cook and Hartikainen2023) and resistant spores (e.g., microsporidians) (Lisnerová et al. Reference Lisnerová, Holzer, Blabolil and Fiala2023). Moreover, the sample matrices can differ in their capacity to preserve DNA. For instance, bacteria and sediment-dwelling organisms can destroy or sequester forms deposited in the sediment (Barnes et al. Reference Barnes, Turner, Jerde, Renshaw, Chadderton and Lodge2014). In the water column, DNA degradation depends on the environmental conditions including temperature, salinity, UV exposure, and presence of suspended solids that can adsorb DNA and/or reduce the negative effects of UV exposure (Barnes et al. Reference Barnes, Turner, Jerde, Renshaw, Chadderton and Lodge2014; Strickler et al. Reference Strickler, Fremier and Goldberg2015). Methods also may be biased toward abundant taxa and leave rare species undetected when: sampling effort is too low, there is preferential amplification of certain parasite taxa, and sequencing depth is insufficient. Further, there are likely biases for certain life stages based on their physical properties (e.g., tough spores vs. soft tissues during DNA extraction). Lastly, reference DNA sequences are lacking for many parasite species—and likely entire clades—largely because parasites themselves are still neglected among the scientific community, and perhaps more significantly, because of the high phylogenetic diversity that characterizes parasite groups, which complicates the use of single genetic markers across manifold taxa.
Despite its limitations and biases, the eDNA approach can complement traditional methods, which remain necessary for formal species description and molecular characterization, which, in turn, allows for the population of reference databases. It is thus imperative to develop and apply metabarcoding tools to more accurately capture parasite biodiversity. Accordingly, we tested the efficacy of parasite eDNA metabarcoding as part of a ParasiteBlitz, a BioBlitz targeting parasites (de Buron et al. Reference de Buron, Hill-Spanik, Atkinson, Vanhove, Kmentová, Georgieva, Díaz-Morales, Kendrick, Roumillat and Rothman2025) at a locality that had not been previously surveyed for parasites. Given the limitations and context dependencies listed above, we used multiple primers targeting different parasite groups, considering 4 habitats, 3 sampling methods, and 2 matrices. To our knowledge, this is the first such combined approach, as published studies have been performed typically in isolation, focusing on certain parasite groups or ecosystems. Herein, we examined the performance of primers designed to target 5 major groups: Microsporidia, Myxozoa, Platyhelminthes, Nematoda, and Eukaryota/Protista, at 4 adjacent and connected aquatic habitats: a wetland, a freshwater pond, a brackish impoundment, and a tidal creek. In each habitat, we considered 2 sample matrices: sediment and water (filtered actively, or passively through the deployment of filters).
Materials and methods
Habitats and sampling sites
A detailed description of the ParasiteBlitz sampling habitats is given in de Buron et al. (Reference de Buron, Hill-Spanik, Atkinson, Vanhove, Kmentová, Georgieva, Díaz-Morales, Kendrick, Roumillat and Rothman2025). Briefly, we sampled 4 contiguous aquatic habitats at Stono Preserve, Charleston, SC, USA (32°44′06″N 80°10′48″W): an ephemeral wetland adjacent to a freshwater pond, and a brackish former rice impoundment that drains into a tidal creek of the Stono River (Figure 1), hereafter referred to as the ‘wetland’, ‘pond’, ‘impoundment’, and ‘tidal creek’, respectively. Environmental sampling occurred in 2023, on the mornings of April 23 (impoundment), 25 (tidal creek during flood tide), 27 (wetland), and 30 (pond), prior to disturbances generated from the collection of hosts for the traditional part of the ParasiteBlitz. Three sites within each habitat were sampled, which represented a range of ecological conditions (e.g., shallow and deep water; Figure 1). Salinity and temperature were recorded using a YSI Professional Plus (Yellow Springs, OH, USA) at 1 or 2 sites within each habitat, and water was collected for measurement of pH immediately upon return to the laboratory (Table 1).
Map of Stono Preserve, Charleston, South Carolina, USA showing sampling sites (triangles) across the 4 habitats: ephemeral wetland, freshwater pond, saltwater impoundment, and tidal creek in the Spartina grass mudflats of the Stono River.
Physicochemical parameters of Stono Preserve during the April 2023 ParasiteBlitz. Parameters were taken at two sites within the impoundment and pond habitats, shown as values separated by commas
Active eDNA sampling
Water samples were taken at each site (500 mL per site, with a total of 1.5 L per habitat). Water was collected into new Nalgene bottles by gloved hand at the surface. Samples were placed on ice during in situ filtration using a battery-powered peristaltic pump (Vampire Sampler, Bürkle GmbH) and self-preserving eDNA Filter Packs (Smith Root, PES membrane, 1.2 μm; Thomas et al. Reference Thomas, Nguyen, Howard and Goldberg2019). Flow-through was collected in a volumetric flask to measure the volume of water filtered. Sites varied in the amount of suspended particulates; therefore, some samples required more than 1 filter: once clogged, we replaced the filter and continued until the 500 mL target volume for each sample was reached. The wetland sites were the exception, with only ~250 mL per site filtered (750 mL total), as this ephemeral water body was much smaller and shallower in comparison to the others. As negative control, 500 mL of ultrapure water was filtered in situ after filtration of samples.
Passive eDNA sampling
eDNA was sampled passively at each site using filter assemblies suspended in the water column for 8 h, proximal to where active sampling occurred. We used a modified method based on Bessey et al. (Reference Bessey, Jarman, Simpson, Miller, Stewart, Keesing and Berry2021): each passive filter assembly consisted of a cellulose acetate filter membrane (0.45 μm, 47 mm diameter) centre-punched with a sterile scalpel, then threaded on nylon monofilament (Figure 2). Site depth was measured, and zip ties were placed on either side of the filter to keep it at ~ half the depth of the water. One end of the monofilament was tied to a lead weight to keep the assembly submerged, and the other was tied to a metal rod driven into the substrate to ensure the assembly would stay in place (Figure 2). Immediately upon retrieval, membranes were preserved in 100% EtOH in a 50-mL nuclease-free Falcon tube. As a negative control, a membrane was placed in situ into a nuclease-free Falcon tube that contained ultrapure water, retrieved, exposed to air for ~20 s, then preserved in 100% EtOH in a different tube.
Passive water sampling. (a), Image of assembly just prior to deployment; (b), Schematic of assembly while deployed; (c), Image of assembly while deployed.
Sediment sampling
Four sediment cores were collected haphazardly at each site. Surface sediment (top 3–5 mm) was collected using corers consisting of sterile 10-cc syringes with the Luer end removed; sediment was then expelled into 15-mL nuclease-free Falcon tubes, which were placed on ice until storage at –80°C.
DNA extraction
DNA extraction was performed using a modified version of Buchner (Reference Buchner2022a) for lysis, Buchner (Reference Buchner2022b) for extraction, and Buchner (Reference Buchner2023) for extract clean-up. DNA extraction from active and passive filters was performed using a guanidine-based approach under a UV-sterilized hood. Before proceeding with the extraction, a passive filter was placed in an individual sterile Petri dish, covered with a lid, and allowed to air-dry for 1 h. Each active filter was removed from the casing and transferred to a sterile Petri dish. The filters were torn into 1–3 cm2 pieces with sterile forceps in the dish, then transferred to a 2-mL twist-top tube filled with 900 μL TNES buffer (50 M TRIS, 400 M NaCl, 20 M EDTA), 100 μL proteinase K (10 mg mL–1 Pro-K 7, BioScience Catalog #RP100B), and ~30 1-mm and 10 2-mm Zirconia beads (BioSpec Products). Technical negative controls were used by integrating 8 controls at random positions among samples. The samples were then bead-beaten for 2 min at 2400 rpm, incubated at 55°C for 1 h in a shaker at 1400 rpm, and stored at –20°C until further extraction. For extraction, the lysate was cleared at 11,000 g at 20°C for 3 min, and 400 μL of lysate was mixed with 800 μL of guanidine hydrochloride (GuHCl) binding buffer (3 M GuHCl, 10 mM Bis-Tris, 90% (v/v) EtOH) in an EconoSpin® 96 Well DNA Binding Plate. The samples were passed through the membrane using a vacuum manifold for 2 min, followed by two washing steps of 600 μL wash buffer (10 mM Tris, 80% (v/v) EtOH) each. The membranes were dried by running the vacuum for 10 min. DNA was eluted by adding 50 μL 10 mM Tris to the membrane and vacuum was applied for 1 min. Eluted DNA was purified using a 96 Channel Portable Electronic Pipette (Integra) to mix 50 μL DNA with 100 μL clean-up bead solution (0.2 mg mL–1 of Sera-Mag SpeedBeads carboxylate modified particles dissolved in PEG-NaCl buffer; for details see Buchner [Reference Buchner2023]), incubated at room temperature for 5 min, then the beads separated by magnet for 3 min between each step, and the supernatant aspirated and discarded. The beads were washed in 100 μL wash buffer twice and incubated at room temperature for 30 s each time. After the second wash, buffer was removed and the beads dried for 5 min at room temperature, then DNA eluted with 50 μL elution buffer, and stored at –20°C. Sediment samples were homogenized via vortexing and/or mixing using a sterile spatula. DNA was extracted from each replicate using a Qiagen PowerSoil Pro Kit (Valencia, CA, USA) following the manufacturer’s instructions and eluted in 100 μL buffer C6.
Extracted DNA from every sample was amplified with a universal eukaryote PCR assay (details below) to ensure that the DNA was suitable for amplification (e.g., not degraded) and test for the presence of PCR inhibitors. Replicates of each sample type (active/passive/sediment) from each habitat were then pooled prior to quantification and PCR. Pools (n = 12) were quantified using a Nanodrop 2000 spectrophotometer (Thermo Fisher, Waltham, MA, USA). Two pools each were made using the sediment and water DNA negative controls, respectively, and were subjected to library preparation and sequencing in the same way as the samples.
Library preparation and high-throughput amplicon sequencing
We followed a modified version of the Illumina 16S metagenomic library preparation protocol (Illumina, 2013). Five parasite groups were targeted for DNA metabarcoding and HTS, Microsporidia, Myxozoa, Nematoda, Platyhelminthes, and Eukaryota/Protista, as primers had been developed previously for these groups (Table 2). Illumina overhang adapters (San Diego, CA, USA) were included on all primers (forward: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3′, reverse: 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3′) for use in the subsequent indexing PCR, except for the primers used in the first PCR of nested approaches. All PCRs were done in triplicate 25-μL reactions, which contained 1x TaKaRa Ex Taq buffer (Clontech, Mountain View, CA, USA), 0.2 mM dNTPs, 0.4 mg mL−1 bovine serum albumin, 0.2–0.5 μM each primer, 1 U TaKaRa Ex Taq DNA polymerase (Clonetech), and 1–5 μL template DNA. All PCRs contained a positive and negative control. PCRs for protists (Pagenkopp Lohan et al. Reference Pagenkopp Lohan, Fleischer, Carney, Holzer and Ruiz2016), nematodes (Weigand et al. Reference Weigand, Kremers and Grabner2016), and microsporidians (Aivelo et al. Reference Aivelo, Harris, Cadle and Wright2018) targeted a portion of the 18S rRNA gene. The myxozoan PCR followed the nested approach of Lisnerová et al. (Reference Lisnerová, Holzer, Blabolil and Fiala2023), but for the sake of efficiency, we multiplexed the second PCR. The first PCR used universal eukaryote 18S rRNA gene primers, and the second used a pool of multiple forward and reverse primers (Table 2). Triplicate reactions from the first PCR were pooled and then used as template in the second PCR, which was also done in triplicate.
Primers tested and used in this study, including their sequence, gene region amplified and sequenced, expected product sizes in base pairs (bp), and reference(s)
18S = the small subunit ribosomal RNA (rRNA), 18S V4 = variable region 4 of the 18S rRNA gene, ITS2 = second internal transcribed spacer region of the 18S rRNA gene, COI = mitochondrial cytochrome c oxidase subunit I gene, 12S = mitochondrial 12S rRNA gene, 18S V9 = variable region 9 of the 18S rRNA gene. Primers in grey were tested but resulting amplicons (if any) were not used in the final data analysis.
* Only pooled sediment and actively filtered water samples were included in this library (i.e., no passively sampled water).
† Used only in the first PCR of nested approaches and thus did not include Illumina overhang adapters.
‡ Sediment samples only were used to generate a Platyhelminthes library using the standard PCR approach.
§ When paired with COI-ASmit2.
For platyhelminths, previously developed assays targeting portions of the mitochondrial COI (Thomas et al. Reference Thomas, Milotic, Vaux and Poulin2022) and 12S rRNA genes (Chan et al. Reference Chan, Saralamba, Saralamba, Ruangsittichai and Thaenkham2022a, Reference Chan, Saralamba, Saralamba, Ruangsittichai, Chaisiri, Limpanont, Charoennitiwat and Thaenkhamb) (Table 2) did not yield amplicons of the expected size or did not amplify at all. Therefore, we tested a nested approach (Vanhove et al. Reference Vanhove, Pariselle, Van Steenberge, Raeymaekers, Hablützel, Gillardin and and Huyse2015) not used previously in metabarcoding studies. The first-round primers were JB3 (= COI-ASmit1) + Cox1_schist-3′(Bowles et al. Reference Bowles, Blair and McManus1995; Littlewood et al. Reference Littlewood, Rohde and Clough1997; Lockyer et al. Reference Lockyer, Olson, Ostergaard, Rollinson, Johnston, Attwood and and Littlewood2003; Table 2), with cycling: 95°C for 5 min, then 40 cycles at 94°C for 1 min, annealing at 44°C for 1 min, and extension at 72°C for 1 min, with final extension at 72°C for 7 min. Triplicate reactions were pooled and then used as template in the second-round PCR using primers JB3 (= COI-ASmit1) + JB4.5 (= COI-ASmit2) (Bowles et al. Reference Bowles, Blair and McManus1995; Littlewood et al. Reference Littlewood, Rohde and Clough1997; Table 2), also done in triplicate, with cycling: 95°C for 5 min, then 40 cycles at 94°C for 1 min, annealing at 50°C for 1 min, and extension at 72°C for 1 min, with final extension at 72°C for 7 min. We also tested a non-nested approach on the sediment DNA only using the JB3 (= COI-ASmit1) + JB4.5 (= COI-ASmit2) primers (Bowles et al. Reference Bowles, Blair and McManus1995; Littlewood et al. Reference Littlewood, Rohde and Clough1997; Table 2), which was also done in triplicate, with cycling as above except with 35 cycles.
Each triplicate reaction was pooled (normalized based on band intensity as visualized on a 1% agarose gel) and then dual-indexed with Illumina Nextera indices. All PCR negatives were pooled and indexed, resulting in 7 indexed negative controls (field water control from each habitat (4), DNA extraction controls (2), PCR negative control (1)). Each 25-μL reaction contained 1x KAPA HiFi HotStart ReadyMix (Kapa Biosystems, Wilmington, MA, USA), 0.4 μM of each adapter, and 1 μL of pooled PCR product. Initial denaturation at 95°C for 5 min was followed by 6 cycles of denaturation at 98°C for 20 s, annealing at 60°C for 45 s, and extension at 72°C for 45 s, and then by a final extension at 72°C for 5 min. Indexed products were purified using 1.8x Agencourt AMPureXP Beads (Beckman Coulter, Brea, CA, USA) or extracted from an agarose gel using a QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). All cleaned and indexed products were quantified using a High Sensitivity dsDNA kit on a Qubit 4.0 (Thermo Fisher Scientific) and were pooled in equimolar amounts. The 7 negative controls were also added to the pool (1 μL each).
The amplicon size of the equimolar pool was assessed using a Bioanalyzer TapeStation (Agilent Technologies, Santa Clara, CA, USA) and was quantified using a KAPA library quantification kit (Kapa Biosystems, Inc.) prior to sequencing on a MiSeq using the 2 × 300 bp kit (Illumina) at the University of Rhode Island, Rhode Island Idea Network of Biomedical Research Excellence (URI RI-INBRE) Centralized Research Core Facility at Rhode Island Genomics Sequencing Center (Kingston, RI, USA). Cluster density was low on the first run, so the library was run again and sequence data from both runs were used in subsequent analyses. Raw sequences and associated metadata were deposited as Project #PRJNA1333354 into the Sequence Read Archive database (Leinonen et al. Reference Leinonen, Sugawara and Shumway2010).
Bioinformatics
Sequences were demultiplexed by the sequencing facility. Primers and adapters were removed using Cutadapt (V3.7, Martin Reference Martin2011). The remainder of preprocessing and data analyses were performed using R (V4.1.2, R Core Team 2021) and RStudio (V2021.09.2, RStudio Team 2021). We used DADA2 (V1.18.0) to resolve amplicon sequence variants (ASVs) because of its advantages over operational taxonomic unit (OTU) clustering methods, which include improved resolution and reusability, allowing for comparison across studies (Callahan et al. Reference Callahan, McMurdie, Rosen, Han, Johnson and Holmes2016; Callahan et al. Reference Callahan, McMurdie and Holmes2017). DADA2 was used for quality filtering (fastqPairedFilter maxEE = 2–3, minLen = 175 bp, truncLen R1 = 250 bp, truncLen R2 = 175–200 bp, depending on primers used), ASV inference, PhiX and chimera removal, and merging of sequence reads. Sequences from each of the 2 MiSeq runs were preprocessed separately, using the same methods, then merged prior to removal of chimeras. To obtain a second estimate of parasite richness, DECIPHER (V2.22.0, Wright Reference Wright2015) was used to align, generate distance matrices, and cluster ASVs into OTUs given 3 similarity thresholds (97–99%) using UPGMA (unweighted pair group method with arithmetic mean) hierarchical clustering, but all subsequent analyses used ASVs.
ASVs were used as queries against the NCBI nt database (downloaded on 14 November 2024) using BLAST (V2.16.0), retrieving a maximum of 10 hits per query (-max_target_seqs 10 -max_hsps 1) along with the NCBI taxonomic identifiers and scientific names (-outfmt ‘7 std qcovhsp staxids sscinames’) (Camacho et al. Reference Camacho, Coulouris, Avagyan, Ning, Papadopoulos, Bealer and Madden2009). Subject sequences were then categorized with empirically determined bitscore cutoffs: 700 = high; 600 = medium; 400 = low; none = other; BLAST results were filtered according to those cutoffs. For each ASV, the number of hits above a given threshold was counted and used to calculate a confidence score (0–1, with a maximum score of 1 meaning 10 out of 10 hits). ASVs that had hits above a given threshold were then removed from the pool, and a lower threshold was examined until no thresholds were remaining. All 18S ASVs received species level identification if high bitscore, confidence >0.1, and resulting BLAST percent similarity (% ID) was >99%; genus if high or medium bitscore, confidence ≥0.1, and % ID >97.5%, family if low or medium bitscore with confidence ≥0.1, and % ID >97.5%; and class if bitscore was ‘other’. ASV taxonomic classifications were also double-checked for accuracy using the lowest common ancestor (LCA) assignment using the taxonkit (V0.18.0) LCA command (Shen & Ren Reference Shen and Ren2021). Decontam (Davis et al. Reference Davis, Proctor, Holmes, Relman and Callahan2018) was used to identify any likely contaminants based on prevalence of ASVs found in the negative controls.
All taxonomic assignments were then examined manually for parasites. Of those assignments, we filtered to keep the following parasite taxa: ‘Apicomplexa’, ‘Microsporidia’, ‘Perkinsozoa’, ‘Monogenea’, ‘Trematoda’, ‘Myxozoa’, ‘Poecilostomatoida’, ‘Mermithidae’, ‘Loxothylacus’, ‘Gyrinicola’ using phyloseq (V1.38.0, McMurdie & Holmes Reference McMurdie and Holmes2013). We did not include myco- or phytoparasites, or parasites of plants such as Cryptomycota, Endomyxa, and Oomycota, given the lack of information regarding these groups and our lack of expertise therewith.
We then compared each parasite ASV to those in GenBank (accessed between 27 August and 5 September 2025) using BLASTn (Altschul et al. Reference Altschul, Gish, Miller, Myers and Lipman1990). For the 18S (and sole arthropod COI) ASVs, we decreased the level of the original taxonomic assignment only in the following cases: species level was assigned if the resulting % IDs were ≥99% to sequence(s) associated with a peer-reviewed publication from only 1 species, genus if ≥98% to sequence(s) from only 1 genus, family if ≥95% to sequences from only 1 family, and class if ≥90% to sequences from only 1 class, all with ≥90% query coverage. If these criteria were not met, the assignment was increased to the higher taxonomic level relative to that assignment (e.g., if the original assignment was to genus level, but the manual BLAST result revealed that the ASV was ≥ 98% similar to sequences from more than 1 genus from the same family, then the taxonomic assignment was changed to family level). Otherwise, the assignment remained the same. For identification and delimitation of platyhelminth COI parasite sequences, alignments with sequences retrieved from GenBank for each taxonomic group were made using MAFFT (V7, Katoh & Standley Reference Katoh and Standley2013) as an online execution under default setting parameters, and taxonomic classifications were based on genetic distances (uncorrected p-distance) calculated in MEGA (V11, Tamura et al. Reference Tamura, Stecher and Kumar2021). ASVs were considered different species when dissimilarity was >3% and different genera when dissimilarity was >15%. For parasitic brackish/marine taxa, we modified the taxonomic classifications resulting from NCBI to follow those of the World Register of Marine Species (WoRMS 2025). For others (e.g., terrestrial, freshwater), we followed the most recent peer-reviewed literature.
Phyloseq (McMurdie & Holmes Reference McMurdie and Holmes2013) was used to merge sequences resulting from all libraries based on method*habitat, resulting in 12 samples (3 methods × 4 habitats). iNEXT (V3.0.1, Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Ma, Colwell and Ellison2014, Hsieh et al. Reference Hsieh, Ma and Chao2016) was used to plot rarefaction and extrapolation curves of species richness estimates and calculate richness estimates (q = 0, 1,000 bootstrap replicates, confidence = 95%, sampling depth = minimum number of reads x 2) based on the observed number of ASVs.
Bray-Curtis dissimilarities were calculated on proportionally transformed ASV tables in phyloseq (McMurdie & Holmes Reference McMurdie and Holmes2013), and hierarchical clustering analysis using the resulting dissimilarity matrix and the average linkage method was done in R (V4.1.2, R Core Team 2021). Permutational analysis of variance (PERMANOVA) was also calculated (using 999 permutations) to examine if parasite communities differed significantly among habitats and sampling methods using vegan (V2.6-2, Oksanen et al. Reference Oksanen, Simpson, Blanchet, Kindt, Legendre, Minchin and Weedon2022). Heatmap plots were generated using NeatMap (Rajaram & Oono Reference Rajaram and Oono2010) in phyloseq (ordination = principal coordinates analysis, distance = Bray-Curtis; McMurdie & Holmes Reference McMurdie and Holmes2013) to visualize differences in parasite ASV abundances across sampling methods for each parasite group (excluding Arthropoda and Nematoda because so few ASVs were detected for these 2 taxa) and from each habitat. Because some samples did not include any ASV of a particular parasite group, subsets of the data were plotted; for Microsporidia and Myzozoa: sediment and active water samples from all habitats were used; for Myxozoa: only samples from the impoundment; and for Platyhelminthes: all samples from only the impoundment and pond. All plots and graphs were generated using ggplot2 (Wickham Reference Wickham2016) and VENN (http://bioinformatics.psb.ugent.be/webtools/Venn/), and Inkscape (www.inkscape.org) was used to customize figures.
Results
Sequencing
DNA yields from the passive water samples were low (undetectable to 2.4 ng μL–1) compared to the active water (3.2–17.8 ng μL–1) and sediment samples (104.6–316.9 ng μL–1). Nevertheless, amplifications were successful for all sample types for all assays, except for the microsporidian assay on DNA from passive filters. For sequencing, we pooled 60 sample libraries (all 3 sampling methods × 4 habitats × 4 metabarcoding assays (Eukaryota/Protista, Myxozoa, Nematoda, Platyhelminthes nested) + 2 sampling methods (sediment and active water) × 4 habitats × 1 assay (Microsporidia) + 1 sampling method (sediment) × 4 habitats × 1 assay (Platyhelminthes)) and 7 control libraries.
Both sequencing runs yielded 5,831,732 sequences grouped into 12,941 ASVs after QA/QC, with the number of sequences per library ranging from 2 to 284,149 (mean = 100,981) for 54 of 60 sample libraries (5/12 myxozoan and 1/12 nematode libraries produced no sequence; Supplementary Tables 1 and 2). After filtering to keep only parasite taxa and removing 1 contaminant ASV that was found in both a negative control and a sample, 36/56 successfully sequenced libraries yielded parasite sequences (Supplementary Tables 1 and 2). The Platyhelminthes PCR assay that we tested only on 4 sediment samples yielded 782,685 sequences, but only 1 ASV (90 sequences) was identified as a parasite (a nematode, Gyrinicola Yamaguti, 1938) with the remaining ASVs identified as Oomycota. Therefore, we removed all sequences resulting from the Platyhelminthes (not nested) PCR assay, leaving 9,018 ASVs (5,046,029 sequences), which corresponded to 5,023–7,264 OTUs (97–99% similarity threshold).
Primer performance: sensitivity and specificity
Overall, the microsporidian primers had the highest specificity: 50.0% of all ASVs and 91.2% of sequences were identified as either microsporidian or myzozoan parasites (Table 3), and of those, 99.5% (733/737) were ASVs of Microsporidia and 0.5% were ASVs of Myzozoa (Figure 3). The remaining 50.0% of ASVs (9.8% of sequences) that could not be categorized as a parasite were mostly unclassifiable (i.e., NA (208 ASVs) or unclassified Eukaryota (122 ASVs)).
Percent and number of parasite and target parasite amplicon sequence variants (ASVs) and sequence reads along with the total number of ASVs and sequence reads resulting from each PCR assay and the number of libraries that yielded parasite ASVs (of 12 total samples for each except for the Microsporidia assay where only 8 samples amplified successfully (see Supplementary Table 1))
18S V4 = variable region 4 of the small subunit ribosomal RNA gene, COI = mitochondrial cytochrome c oxidase subunit I gene.
The other assays yielded the remaining 30% of parasite ASVs with variable success and specificity (Table 3). The myxozoan assay yielded the lowest number of sequences (64,122) and only yielded 3 parasite ASVs (2 myxozoans and 1 nematode) and unexpectedly, more myxozoan ASVs were recovered by other assays; the ‘non-parasite’ ASVs were classified as free-living arthropods or nematodes, or as unclassifiable. For the nematode assay, only 0.8% of ASVs were parasitic nematodes, but the primers also amplified ASVs from 4 other parasite groups (Table 3, Figure 3), plus arthropods and ‘unclassified’ ASVs that could not be confidently categorized as parasites. The Platyhelminthes nested COI assay yielded 78 parasite ASVs, of which 77 were of Platyhelminthes and 1 of Arthropoda; most ASVs resulting from this assay were unclassifiable (580), but there were also ASVs of Arthropoda (9) and Bacillariophyta (14) that were not categorized as parasitic. Lastly, the assay used to detect protists and other eukaryotic parasites yielded 191 parasite ASVs from 3 groups: Myzozoa, Myxozoa, and Nematoda (Figure 3). Most ASVs from the universal eukaryote assay were unclassifiable (2,338), and the rest could not confidently be categorized as parasitic: Bacillariophyta (417), Chlorophyta (344), Chitridiomycota (233), Ciliophora (277), Oomycota (129), and unclassified Fungi (221). Some ASVs had the same taxonomic identification, but were derived from different libraries, and thus only differed due to length or gene region. We kept all myxozoan ASVs as there were so few (36) across 3 different assays. However, we estimated that these likely only reflect 16–24 different species based on OTU clustering analysis (Supplementary Table 3). We did remove 20 other ASVs that were possibly redundant among the other 18S rDNA and COI assays.
Taxa detected
Our final dataset, therefore, contained 1,036 parasite ASVs (1,361,606 sequences) corresponding to 541–712 OTUs (Supplementary Table 3) from 6 parasite groups: Arthropoda (1 ASV), Microsporidia (733 ASVs), Myxozoa (36 ASVs), Myzozoa (176 ASVs), Nematoda (7 ASVs), and Platyhelminthes (76 ASVs of Trematoda, 7 of Monogenea) (Table 4). ASVs of Microsporidia made up 70% of all parasite ASVs, which were mainly detected via the microsporidian assay (Figure 3). We were able to identify ASVs to 9 species (11 ASVs), 28 genera (131 ASVs), 39 families (228 ASVs), 16 orders (303 ASVs), and 11 classes (420 ASVs; Table 4). Of note, no ASVs of Acanthocephala, Cestoda, Euhirudinea, or Nematomorpha were detected.
Taxonomic identification of all parasite amplicon sequence variants (ASVs) in our final dataset. Total number of ASVs per identification is in parentheses
NA = taxonomic level could not be assigned. Classifications of brackish and marine taxa based on those found in World Register of Marine Species (WoRMS 2025); freshwater and terrestrial taxa classifications were based on most recent peer-reviewed literature.
The sole arthropod ASV was 100% similar with 100% coverage to a sequence of the sacculinid Loxothylacus panopaei (Gissler, 1884) from the Harris mud crab Rhithropanopeus harrisii (Gould) collected from the Neuse River in North Carolina, USA (GenBank accession number AY265372, Glenner et al. Reference Glenner, Lützen and Takahashi2003) and was only detected in the tidal creek via active water filtration using the Platyhelminthes COI assay. Identified microsporidian ASVs were from 3 orders, 10 families, 15 genera, and 7 species: Agmasoma penaei (Sprague, 1950) detected in the sediment from the tidal creek and impoundment (2 ASVs); Enteropsectra longa Zhang and Félix, Reference Zhang, Sachse, Prevost, Luallen, Troemel and Félix2016 found in pond sediment; Euplotespora binucleata Fokin, Di Giuseppe, Erra, and Dini, Reference Fokin, di Giuseppe, Erra and Dini2008 detected in the pond via active water filtration (1 ASV); Mrazekia macrocyclopis Issi, Tokarev, Voronin, Seliverstova, Pavlova, and Dolgikh, Reference Issi, Tokarev, Voronin, Seliverstova, Pavlova and Dolgikh2010 detected in sediment and active water samples from the tidal creek, impoundment, and pond (1 ASV); Nosema ceranae Fries, Feng, da Silva, Slemenda, and Pieniazek, Reference Fries, Feng, da Silva, Slemenda and Pieniazek1996 detected in impoundment sediment (1 ASV); Paranosema grylli (Sokolova, Selezniov, Dolgikh & Issi, Reference Sokolova, Selezniov, Dolgik and and Issi1994) detected in wetland sediment (1 ASV), and Perezia nelsoni (Sprague, 1950) from impoundment sediment (2 ASVs). Two ASVs were 99.7–100% similar with 100% query coverage to a sequence of Agmasoma penaei from white shrimp Penaeus setiferus (L.) off the coast of Louisiana, USA in the Gulf of Mexico (KF549987, Sokolova et al. Reference Sokolova, Pelin, Hawke and Corradi2015) and in the Charleston Harbor watershed, SC, USA (OL467311, OL467312, Zuidema et al. Reference Zuidema, de Buron, Kingsley-Smith, Hill-Spanik, Fanani and Kendrick2023); 1 ASV was 99.7% similar with 92% coverage to a sequence of Enteropsectra longa from a unidentified free-living nematode species of Oscheius Andrassy in Iceland (KX360142, Zhang et al. Reference Zhang, Sachse, Prevost, Luallen, Troemel and Félix2016); 1 ASV was 100% similar with 100% coverage to a sequence of Euplotospora binucleata from a free-living ciliate Euplotes woodruffi Gaw collected from a brackish lagoon in northeastern Italy (DQ675604, Fokin et al. Reference Fokin, di Giuseppe, Erra and Dini2008); 1 ASV was 99.4% similar with 100% coverage to a sequence of Mrazekia macrocyclopis from copepod Macrocyclops albidus albidus (Jurine) collected from a freshwater pond in St. Petersburg, Russia (FJ914315, Issi et al. Reference Issi, Tokarev, Voronin, Seliverstova, Pavlova and Dolgikh2010); 1 ASV was 100% similar with 100% coverage to several sequences of N. ceranae (= Vairimorpha ceranae in GenBank, but see Bartolomé et al. (Reference Bartolomé, Higes, Martín Hernández, Chen, Evans and Huang2024) for debate regarding the validity of genus synonymy) from bees in Europe (e.g., GU131064, Sagastume et al. Reference Sagastume, Del Águila, Martín-Hernández, Higes and Henriques-Gil2011; KC680620, Roudel et al. Reference Roudel, Aufauvre, Corbara, Delbac and Blot2013) and Asia (e.g., JN872261, Li et al. Reference Li, Chen, Wu, Peng, An, Schmid-Hempel and Schmid-Hempel2012; LC510236, Takashima et al. Reference Takashima, Ohari and Itagaki2021); 1 ASV was 99.8% similar with 100% coverage to a sequence of Paranosema grylli collected from experimentally-infected cricket Gryllus bimaculatus De Geer (AY305325, Sokolova et al. Reference Sokolova, Dolgikh, Morzhina, Nassonova, Issi, Terry and and Vossbrinck2003); and 2 ASVs were 99.2% and 99.7% similar with 96% and 91% coverage, respectively, to a sequence of Perezia nelsoni from white shrimp Penaeus setiferus collected from the Gulf of Mexico in Louisiana, USA (KX856426, Sokolova & Hawke Reference Sokolova and Hawke2016). Myxozoan ASVs comprised 2 orders, 3 families, 2 genera, and 1 species: 1 ASV detected in the tidal creek and impoundment via active water sampling, that was 100% similar with 100% coverage to a sequence of Ellipsomyxa adlardi Whipps & Font, Reference Whipps and Font2013 from the naked goby Gobiosoma bosc (Lacepède) collected in Lake Pontchartrain, Louisiana, USA (JX443488, Whipps & Font Reference Whipps and Font2013). Myzozoan ASVs comprised 2 subphyla, 4 classes, 5 orders, 12 families, and 4 genera; nematode ASVs: 2 classes, 2 orders, 2 families, 1 genus; and platyhelminth ASVs: 1 subphylum, 2 classes, 4 orders, 9 families, 3 genera (Table 4). Collectively, the parasite-targeted assays also recovered ASVs of potential host taxa including annelids (e.g., Diopatra Audouin & Milne Edwards, Capitella Blainville), arthropods (e.g., copepods such as Macrocyclops albidus albidus, insects, and arachnids), echinoderms (e.g., Asteroidea de Blainville, Holothuroidea de Blainville), bivalve and gastropod mollusks (e.g., Crassostrea virginica (Gmelin) and Stylommatophora A. Schmidt), sponges (e.g., Mycale Gray), and vertebrates including fishes (e.g., Pagrus Cuvier).
Sampling methods
Sediment and actively filtered water eDNA samples were the most successful in generating parasite ASVs and sequences: 15/20 and 14/20 libraries yielded 596,921 (36.3% of the total) and 495,310 (24.2%) parasite sequences, respectively (Table 5A). The actively filtered water samples yielded ASVs of all 6 groups compared to only 4 from sediment (no parasitic arthropod or nematode), despite detection of ~3× more ASVs in the sediment (837 vs. 250) (Table 5A, Figure 4). Passively sampled water yielded the fewest parasite ASVs (27) and sequences (296,375 or 19.8% of the total), partially due to the microsporidian assay not amplifying these samples, as mentioned above, but also because 2/16 libraries did not yield any usable sequences and 7/16 libraries did not yield parasite sequences, i.e., only 7/16 libraries yielded parasite sequences (Table 5A, Supplemental Tables 1 and 2). ASVs of only 3 groups were detected from the passively filtered water samples: Platyhelminthes, Myxozoa, and Myzozoa (Figure 4). There were only 2 ASVs (ASV 1 and 1109) shared among all 3 sampling methods (Figure 4), both identified as the haploporid digenean Saccocoelioides Szidat, 1954. There were 6 ASVs shared among actively and passively sampled water, 1 was identified as myxozoan Ellipsomyxa Køie, 2003, 2 as Saccocoelioides, 2 as unclassified species from Class Perkinsea, and 1 as an apicomplexan in Order Eucoccidiorida. Each method also yielded unique ASVs (Figure 4). The 19 unique ASVs resulting from the passive water sampling were identified as Saccocoelioides (5), gyrodactylid monogeneans (5), and apicomplexans (9), which included a species of Eimeriidae (1) and unclassified species of Eucoccidiorida (6), Perkinsea (1), and Conoidasida (1). When comparing only sediment and active water sampling, sediment captured more parasite ASVs compared to active water for Microsporidia (per sample means = 202 and 30 ASVs, respectively) and Myzozoa (per sample means = 28 and 12, respectively) across all habitats (Supplementary Figure 1). The opposite was true for Platyhelminthes detection in the impoundment and pond (the only habitats for which we could make comparisons), where more ASVs were detected in active water than in sediment (per sample means = 31 and 4 ASVs, respectively) and for Myxozoa in the impoundment (total ASVs = 44 and 7, respectively; Supplementary Figure 1).
Percent and number of parasite amplicon sequence variants (ASVs) and sequence reads along with the total number of ASVs and sequence reads; the number of libraries that yielded parasite ASVs of the total of samples that amplified successfully (see Supplementary Table 1) per substrate/method (A) and habitat (B)
Proportion of parasite amplicon sequence variants (ASVs) resulting from each sampling method/substrate used during the April 2023 ParasiteBlitz at Stono Preserve, South Carolina, USA, and the number of unique and shared ASVs.
Habitats
A similar number of ASVs was found in each habitat (Figure 5), with each having a similar number of unique ASVs, despite uneven sequencing and a varying number of libraries yielding parasite sequences across habitats (Table 5B, Figure 6). Richness estimates for each habitat varied according to sampling matrix, as sediment had higher richness than water samples (Figure 6). In the sediment, estimated richness was highest in the wetland (268.13 ± 2.59 ASVs) and impoundment (265.21 ± 3.02), followed by the tidal creek (243.56 ± 4.39) and pond (208.74 ± 4.23; sampling depth = 195,660 sequences; Figure 6). In active water samples, estimated richness was highest in the impoundment (92.21 ± 0.89 ASVs), followed by the wetland (61.25 ± 2.82 ASVs), pond (51.74 ± 0.56), and lastly the tidal creek (48.95 ± 0.45; sampling depth = 16,626 sequences; Figure 6). Sequencing depth was very low for 2 passive water samples (tidal creek = 11 reads, pond = 99 reads), so we did not calculate richness estimates for this sampling method. ASVs of microsporidians were dominant in all habitats, followed by myzozoans (Figure 5). There were only 5 ASVs shared among all habitats (Figure 5), all of which were microsporidians (3 identified as unclassified species, 1 assigned to family level only (Nosematidae Labbé, 1899), and 1 to genus level as Nosema Nägeli, 1857). No parasitic nematode ASV was detected in the brackish habitats (tidal creek and impoundment), and no myxozoan ASV was recovered from freshwater (pond and wetland). There were more ASVs shared between the tidal creek and impoundment and between the pond and wetland than among the more dissimilar habitats (Figure 5). This relationship also was reflected in the hierarchical clustering analysis where parasite communities from the sediment and active water from each habitat formed a distinct cluster (PERMANOVA P = 0.001 for habitat, P = 0.686 for sampling method; Figure 7). The community from the passive water sample from the impoundment also clustered with those from the other sample types from this habitat, but the other passive water sample communities did not group with their respective habitats (Figure 7).
Proportion of parasite amplicon sequence variants (ASVs) detected in each habitat during the 2023 ParasiteBlitz at Stono Preserve, South Carolina, USA, and the number of unique and shared ASVs.
Rarefaction (solid lines) and extrapolation (dashed lines) curves of estimated species richness based on the number of observed amplicon sequence variants for each method*habitat with 1,000 bootstrap replicates generated using iNEXT (Chao et al. Reference Chao, Gotelli, Hsieh, Sander, Ma, Colwell and Ellison2014, Hsieh et al. Reference Hsieh, Ma and Chao2016). The shaded area represents the 95% confidence intervals around the estimates.
Hierarchical clustering based on Bray-Curtis dissimilarities among parasite communities (relative abundance of amplicon sequence variants) found in each method*habitat using the average linkage method.
Discussion
The methods employed in this study proved successful for the detection of parasites across different habitats and complementary to traditional approaches for parasite discovery (see other papers in this Special Collection). Overall, we identified ~1,000 parasite ASVs corresponding to ~600 parasite OTUs out of 1.36 million parasite sequences (and 5.8 million total sequences). The exact number of parasite species likely lies somewhere between these two estimates, as our ASVs might include sequences of the same parasites (due to possible sequencing error, for instance) and thus may represent an overestimation of parasite richness when compared to the OTUs. However, we consider our estimate of parasite richness to be rather conservative because we only included taxa that we were certain to be parasitic; in other words, we were strictly selective and excluded ASVs corresponding to taxa not definitely known to be parasitic due to lack of biological data and those for which we lacked the necessary expertise, such as plant parasites.
Microsporidians were the most diverse parasite group detected both overall and in each habitat, but we detected fewer taxa (733 ASVs or 392–517 OTUs) compared to previous metabarcoding studies using water samples: Chauvet et al. (Reference Chauvet, Debroas, Moné, Dubuffet and Lepère2022) detected 1,472 OTUs from 2.9 million reads across 15 marine and freshwater sites; Chauvet et al. (Reference Chauvet, Monjot and Lepère2023) detected 1,853 OTUs from 2.9 million reads in samples collected from a freshwater lake over 1 year, and Dubuffett et al. (Reference Dubuffet, Chauvet, Moné, Debroas and Lepère2021) detected 1,531 OTUs from 323,904 reads from 2 freshwater lakes over 3 months. In addition to sampling across larger spatial or temporal scales, these authors filtered a significantly larger volume of water (2.5–22 L per sample) and filtered water successively to concentrate spores, culminating in a smaller final pore size (0.8 μm), which could account for at least part of the higher richness estimates compared to our study.
Only 1 of the 7 microsporidian species we detected has been previously reported from SC: Agmasoma penaei in its type host, white shrimp Penaeus (Litopenaeus) setiferus (see Zuidema et al. Reference Zuidema, de Buron, Kingsley-Smith, Hill-Spanik, Fanani and Kendrick2023). Perezia (Nosema) nelsoni has also been reported in the white shrimp in the Gulf of Mexico (Sokolova & Hawke Reference Sokolova and Hawke2016) and in its type host, the brown shrimp Farfantepenaeus aztecus (Ives) from Ossabaw Sound, Georgia, USA (Overstreet, Reference Overstreet1973), a neighbouring state to SC, but not reported from either shrimp species in SC. Nosema ceranae was first described in the Asian honeybee Apis cerana Fabricius from China (Fries et al. Reference Fries, Feng, da Silva, Slemenda and Pieniazek1996) but since has been reported as an emerging parasite of honeybees worldwide (MacInnis et al. Reference MacInnis, Keddie and Pernal2022). It has been present in the United States since 1995 (Chen et al. Reference Chen, Evans, Smith and Pettis2008), where it is now common in the domestic honeybee Apis mellifera L. (Chen et al. Reference Chen, Evans, Smith and Pettis2008; Williams et al. Reference Williams, Shafer, Rogers, Shutler and Stewart2008), and while we could find no report of its presence in SC, honeybees are known to be infected by this microsporidian in Virginia, also a southeastern US state (Traver & Fell Reference Traver and Fell2011). Paranosema grylli has not been reported outside of Russia where it was described from laboratory crickets (Sokolova et al. Reference Sokolova, Selezniov, Dolgik and and Issi1994). However, its type host, the black cricket Gryllus bimaculatus, is an introduced species considered a threatening pest in the USA, where it is known to be sold to the public in pet stores (Weissman et al. Reference Weissman, Gray, Pham and Tijssen2012). While no feral population is yet reported from the USA, the species could become invasive (Weissman & Gray Reference Weissman and Gray2019). Significantly, most US cricket breeders are located in the Southeast (Weissman et al. Reference Weissman, Gray, Pham and Tijssen2012), so detection of Parasonema grylli at Stono Preserve could be an indication of a parasite introduced, possibly along with its host in the wild, or the parasite may have spilled over to native species of Gryllus that are common to our area (B. Scholtens, College of Charleston, pers. comm.). The detection of 3 other microsporidian species that have only ever been reported once may reflect potential introductions as well: Euplotespora binucleata was described from ciliate Euplotes woodruffi from a brackish lagoon in Italy (Fokin et al. Reference Fokin, di Giuseppe, Erra and Dini2008), Mrazekia macrocyclopis was described from copepod Macrocyclops albidus in a freshwater pond in Russia (Issi et al. Reference Issi, Tokarev, Voronin, Seliverstova, Pavlova and Dolgikh2010), and Enteropsectra longa was described from a soil nematode Oscheius sp. in Iceland (Zhang et al. Reference Zhang, Sachse, Prevost, Luallen, Troemel and Félix2016). It is possible these parasites were introduced to the SC region, likely via anthropogenic means given the large geographic distance. Alternatively, the 18S rRNA gene region can be too conserved to differentiate species in some cases (Bojko et al. Reference Bojko, Reinke, Stentiford, Williams, Rogers and Bass2022) and the ASVs assigned to these species thus could be from closely related but different species.
Also because of the use of the 18S rRNA gene region, our results are likely an underestimate of microsporidian richness at Stono Preserve. Bojko et al. (Reference Bojko, Reinke, Stentiford, Williams, Rogers and Bass2022) recommended the use of the partial/complete 28S rRNA and the ITS region of rRNA genes as additional markers; however, adding these markers herein was over our resource budget. Microsporidia is a large and diverse group of ~1,400 known species that infect metazoans and protozoans (Murareanu et al. Reference Murareanu, Sukhdeo, Qu, Jiang and Reinke2021), and there are likely thousands of unknown or uncharacterized taxa (Bojko et al. Reference Bojko, Reinke, Stentiford, Williams, Rogers and Bass2022). In our case, 594 of 733 microsporidian ASVs could not be assigned to a lower taxonomic level, which illustrates the evident gap in microsporidian reference sequences that are needed to make molecular identifications.
One of the 2 Platyhelminthes assays generated the second-highest level of specific taxon recovery: 9.8% of ASVs and 29.4% of sequences recovered from the COI nested protocol (Vanhove et al. Reference Vanhove, Pariselle, Van Steenberge, Raeymaekers, Hablützel, Gillardin and and Huyse2015) were Platyhelminthes. Unexpectedly, this assay detected the only arthropod parasite ASV in our samples. It was identified as the invasive rhizocephalan Loxothylacus panopaei, which has been reported previously from flatback mud crab Eurypanopeus depressus (Smith) in SC (O’Shaughnessy et al. Reference O’Shaughnessy, Freshwater and Burge2014). The non-nested Platyhelminthes assay using the same second round of primers as the nested assay failed, which we suspect was due to low target DNA concentrations and not an issue with the nested approach, given the enrichment of the target in the first round of PCR. Adapting existing primers that are already used for taxonomic identification of parasites can be a starting point for metabarcoding primer development, but the obvious problem is the limitation in sequence length of current short-read HTS technologies—the maximum length is 600 bp including adapters and indexes, severely limiting the number of informative sites especially for already fairly conserved gene regions like 28S. Although long-read sequencing (e.g., Pacific Biosciences and Oxford Nanopore) has improved in accuracy, its use in metagenomics is still in development (Marx Reference Marx2023). In the meantime, further development of assays targeting the COI, and potentially the ITS2 region of the rRNA gene, is necessary. The length of the ITS2 region is variable but less than 600 bp, and because it has been widely applied in the identification and differentiation of parasitic platyhelminths (e.g., Bowles et al. Reference Bowles, Blair and McManus1995; León-Règagnon et al. Reference León-Règagnon, Brooks and Pérez-Ponce de León1999; Morgan & Blair Reference Morgan and Blair1995), reference sequences are available in GenBank, though currently there are fewer than COI (roughly 12,000 vs. 22,000 sequences as of 11 September 2025).
The no-to-low sequence yields for the myxozoan assay were undoubtedly related to the non-optimal PCR design rather than the diversity of these parasites at Stono Preserve based on our fish and annelid examinations (see other papers in this Special Collection). We attempted to use existing myxozoan primers (Lisnerová et al. Reference Lisnerová, Holzer, Blabolil and Fiala2023) in a novel multiplex PCR; however, the assay yielded low-quality results, which we attribute to the possible formation of cross-primer dimers that interfered with amplification. We knew this was a risk, but the alternative, which was to perform 3 nested PCRs in triplicate—9 total PCRs per sample (Lisnerová et al. Reference Lisnerová, Holzer, Blabolil and Fiala2023)—was again impractical given our financial and time constraints and did not align with our goal of efficiency in the framework of a ParasiteBlitz. Despite this, we found between 16 and 36 putative myxozoan species out of 22,855 reads, encompassing 2 orders, 5 families, 4 genera (Ellipsomyxa, Myxolobus Bütschli, 1882, Parvicapsula Shulman, 1953, and Kudoa Meglitsch, 1947), and 1 species (Ellipsomyxa adlardi). For comparison, Lisnerová et al. (Reference Lisnerová, Holzer, Blabolil and Fiala2023) detected 17 myxosporean OTUs from 260,720 reads in a single sediment sample from the Douro River estuary (Portugal). To our knowledge, Ellipsomyxa adlardi has been reported only from its type locality on the Gulf of Mexico coast of Louisiana, USA and only from its type host, the naked goby Gobiosoma bosc (see Whipps & Font Reference Whipps and Font2013). Although not encountered during the traditional part of this ParasiteBlitz (de Buron et al. Reference de Buron, Hill-Spanik, Atkinson, Vanhove, Kmentová, Georgieva, Díaz-Morales, Kendrick, Roumillat and Rothman2025), the naked goby is commonly found in the Charleston Harbor watershed, and given the multiple examples of overlap in organismal diversity between the Gulf of Mexico and the Southern Atlantic Bight, including parasites (e.g., Euzet & de Buron Reference Euzet and de Buron2010; Rosas-Valdez et al. Reference Rosas-Valdez, Morrone, Pinacho-Pinacho, Domínguez-Domínguez and García-Varela2020; Ruhnke et al. Reference Ruhnke, Daniel and Jensen2020), presence of this myxozoan at Stono Preserve is not surprising. No myxozoans were detected in the freshwater habitats, which likely reflects both the aforementioned problems with the multiplex PCR and the pond and wetland having very low fish diversity—only 1 fish species, Gambusia holbrooki Girard, was encountered in the traditional part of the survey (de Buron et al. Reference de Buron, Hill-Spanik, Atkinson, Vanhove, Kmentová, Georgieva, Díaz-Morales, Kendrick, Roumillat and Rothman2025), and no myxozoans were observed in this potential host.
The primers used in the nematode assay lacked specificity but still yielded some non-target parasite sequences, allowing for the detection of myxozoans (which supplemented results of the myxozoan-targeted assay). However, the lack of specificity resulted in low recovery of nematode sequences and even fewer parasitic nematode sequences, with 6 ASVS identified as Gyrinicola. Another nematode ASV (a mermithid) was detected using the general eukaryote assay (Bråte et al. Reference Bråte, Logares, Berney, Ree, Klaveness, Jakobsen and Shalchian-Tabrizi2010). Unfortunately, the second nematode assay targeting the ITS2 region of the rRNA gene (Davey et al. Reference Davey, Utaaker and Fossoy2021) failed to amplify any targets. While this could mean there are few parasitic nematodes at Stono Preserve, it is more likely that our methods were not optimal for recovering either the DNA from these worms or their various life stages. Most studies characterizing nematode communities via eDNA and those evaluating methods for nematode eDNA metabarcoding from faeces, soil, and sediment start by isolating and concentrating the nematodes prior to DNA extraction (e.g., Avramenko et al. Reference Avramenko, Redman, Lewis, Yazwinski, Wasmuth and Gilleard2015, Davey et al. Reference Davey, Utaaker and Fossoy2021, Kawanobe et al. Reference Kawanobe, Toyota and Ritz2021, Porazinska et al. Reference Porazinska, Giblin-Davis, Faller, Farmerie, Kanzaki, Morris, Powers, Tucker, Sung and Thomas2009); this has been shown to enhance DNA recovery and reduce non-target amplification (Davey et al. Reference Davey, Utaaker and Fossoy2021). Additionally, it is notoriously difficult to extract DNA from nematodes due to their multilayered cuticle (Dawkins & Spencer Reference Dawkins and Spencer1989; Seesao et al. Reference Seesao, Audebert, Verrez-Bagnis, Merlin, Jérôme, Viscogliosi and and Gay2014), and the DNA extraction step is crucial for detection of nematode DNA and PCR efficiency (Högberg et al. Reference Högberg, Baltrušis, Enweji and Höglund2022). For holistic surveys, like BioBlitzes, prior isolation of nematodes could result in the loss of other parasites (and/or their DNA), though, so it would be necessary to duplicate sampling efforts. eDNA from nematodes can be recovered from water samples although perhaps with less success than for platyhelminths, which has been attributed to differences in life histories between these parasites (Thomas et al. Reference Thomas, Milotic, Vaux and Poulin2022). Therefore, it, again, may be a matter of optimizing the volume of water filtered in order to have enough nematode eDNA for detection.
Limited resources precluded our testing of other assays, but of the groups we targeted, Nematoda had the most primers available in the literature, with several markers being used to characterize the ‘nemabiome’ or the community of nematodes inhabiting a single host or environmental niche (Avramenko et al. Reference Avramenko, Redman, Lewis, Yazwinski, Wasmuth and Gilleard2015; nemabiome.ca). Such markers include those also targeting the 18S rRNA gene (e.g., Kawanobe et al. Reference Kawanobe, Toyota and Ritz2021) and the ITS2 rRNA gene region (e.g., Redman et al. Reference Redman, Queiroz, Bartley, Levy, Avramenko and Gilleard2019) not tested herein, and others targeting the 28S rRNA gene (e.g., Porazinska et al. Reference Porazinska, Giblin-Davis, Faller, Farmerie, Kanzaki, Morris, Powers, Tucker, Sung and Thomas2009), COI (e.g., Macheriotou et al. Reference Macheriotou, Guilini, Bezerra, Tytgat, Nguyen, Nguyen and and Derycke2019), 12S and 16S rRNA genes (Chan et al. Reference Chan, Chaisiri, Morand, Saralamba and Thaenkham2020). The concern with deviating from frequently used markers is that even if the alternative markers allow better differentiation of species or broader detection of phylogenetic groups, we are currently limited by the number of reference sequences available for taxonomic assignment. Therefore, until databases are populated with parasite sequences, the tradeoff is to utilize a widely used marker for which there are many reference sequences.
The general eukaryotic primers we used (Bråte et al. Reference Bråte, Logares, Berney, Ree, Klaveness, Jakobsen and Shalchian-Tabrizi2010) allowed for detection of Myzozoa, some of which were also detected by the nematode primers. However, besides myzozoans, these primers were inefficient at recovering parasite sequences overall as we were only confident in categorizing 3.3% of the ASVs and 1.8% of the sequences as ‘parasite’. The other generic primer assay that we tested, which targets the ninth variable region (V9) of the 18S rRNA gene (Amaral-Zettler et al. Reference Amaral-Zettler, McCliment, Ducklow and Huse2009), was successful at amplifying our samples; however, we opted to sequence the V4 region as it yields a longer fragment (470 vs. 200 bp). On the plus side, non-target taxa detected included free-living animals such as amphibians, annelids, fishes, and mollusks. Although not an aim of this study, this provides evidence for the potential of using eDNA metabarcoding to also inform the presence of putative hosts during a ParasiteBlitz; this is of particular interest for animals that are cryptic or invasive, helpful in the elucidation of parasite life cycles, and essential to the understanding of parasite population dynamics.
Overall, active water sampling was the most productive in terms of capturing the most parasite groups during this ParasiteBlitz, especially compared to passive water sampling. On the other hand, sediment sampling yielded the most parasite ASVs and sequences, with parasite richness estimates 3–5 times higher than those from active water samples, depending on habitat. The parasite communities from each habitat were distinct based on hierarchical clustering analysis regardless of sample type, except in the case of passive water. The most efficient method for capturing diversity within a particular parasite group, however, did appear to vary: heatmap plots suggested that while sediment may be more efficient for capturing microsporidian and myzozoan diversity, active water sampling may be more efficient for platyhelminths and myxozoans. Additional studies are needed to confirm if this would be the case across multiple localities, especially for Platyhelminthes and Myxozoa, for which we had limited datasets, but these results align well with our current knowledge of the life cycles of digeneans and myxozoans. Passive water filtering has been shown to be effective for detecting free-living biodiversity (e.g., fishes), which is, as mentioned above, particularly informative when in the presence of rare species (Bessey et al. Reference Bessey, Jarman, Simpson, Miller, Stewart, Keesing and Berry2021; Bessey et al. Reference Bessey, Gao, Truong, Miller, Jarman and Berry2022; Chen et al. Reference Chen, Li, Zhao and Yao2024). In our case, passive filtration resulted in 19 ASVs that were not detected using the other methods despite having low DNA yields, and thus, this method has a definite value and should not be disregarded but further optimized for future applications. To our knowledge, the effect of submersion time on DNA recovery from passive filters has not been investigated for parasites. Such time might be critical, contrary to what is known of fish detection, for which submersion time did not affect DNA yield (Bessey et al. Reference Bessey, Jarman, Simpson, Miller, Stewart, Keesing and Berry2021; Bessey et al. Reference Bessey, Gao, Truong, Miller, Jarman and Berry2022) and for which 5 min submersion was enough time to capture fish mitochondrial eDNA; furthermore, species richness estimates were comparable to those of active filtration methods (Bessey et al. Reference Bessey, Gao, Truong, Miller, Jarman and Berry2022). Further optimization of active water and sediment sampling for eDNA metabarcoding of parasites is also necessary. As briefly mentioned above with respect to Microsporidia, the volume of water and membrane size are important considerations for optimal parasite capture. Collection depth of sediments should also be investigated to help guide future eDNA sampling from substrate. Lastly, we were careful to take the eDNA samples a day prior to avoid disturbance of the environment generated during the traditional part of our survey (e.g., netting of fish, deployment of traps). In particular, we wanted to avoid re-suspension of humic substances, especially humic acids, which are known to affect the purity of DNA (Tebbe & Vahjen Reference Tebbe and Vahjen1993) and inhibit PCR (Sidstedt et al. Reference Sidstedt, Rådström and Hedman2020). In hindsight, however, we question if controlled disturbance (e.g., via agitation of a certain area of water/sediment for some length of time) could have resulted in an increased number of taxa recovered from the passive filters, ultimately reducing the number of methods to be used in surveys for parasites.
In optimizing methods for eDNA metabarcoding workflows, it is also important to consider the context dependencies, as different localities vary in their ecological characteristics (e.g., freshwater vs. marine, eutrophic vs. oligotrophic). Most ASVs were unique to a single habitat, with expected overlap between the 2 brackish and the 2 freshwater habitats, respectively. The few ASVs common to all 4 habitats were all microsporidians, possibly because 1) the spores of these parasites are particularly resistant to environmental conditions (Murareanu et al. Reference Murareanu, Sukhdeo, Qu, Jiang and Reinke2021), 2) their host(s) may be tolerant to harsh abiotic conditions such as hypoxia or a wide salinity gradient (Williams et al. Reference Williams, Hamilton, Jones and Bass2018), and/or 3) these parasites have a simple life cycle and may lack host specificity (Willis & Reinke 2022). In our case, the pond and wetland had very low oxygen concentrations (0.16–0.22 mg L–1 and 2.22 mg L–1, respectively), and given the connectivity of the aquatic system studied, ASVs found in the tidal creek could be the result of downstream transport of these parasite spores and/or DNA (Williams et al. Reference Williams, Hamilton, Jones and Bass2018).
Although powerful, incorporating eDNA metabarcoding into a ParasiteBlitz was challenging. Library preparation took more time than expected as previously published metabarcoding primers and protocols did not all work on our samples, and as a consequence, we had a prolonged phase of troubleshooting primers and protocols. This experience taught us that maximum efficacy requires the need for an eDNA ‘team’, expediting this process, and ensuring that the speed component of the ParasiteBlitz is not lost. In the same vein, such a team could involve students in most aspects of the library preparation (e.g., DNA extraction, PCR) in the spirit of involving non-experts in BioBlitzes. It also emphasized the need for the development and subsequent use of more specific and novel parasite metabarcoding primers, especially for parasite taxa not detected using available markers and that have been tested across more localities and sample matrices, thus accelerating the library preparation process. Furthermore, the development of a curated, comprehensive parasite reference sequence database would be extremely useful, as even after taxonomic assignment, manual filtering to keep only parasitic taxa was laborious. The process of determining ‘target’ vs. ‘non-target’ vs. ‘maybe target’ is currently common to most metabarcoding studies, as primers are not specific enough, making it necessary to overshoot sequencing depth (e.g., number of reads per sample) to account for the lack of specificity of the primers used and manually curate and filter the data.
To gain a more comprehensive understanding of parasite diversity within a given ecosystem, sampling techniques should aim to maximize species detection. At the current stage, eDNA metabarcoding is particularly helpful for informing ParasiteBlitzers on which hosts and parasite taxa they might encounter during traditional surveys, thus enabling more effective organization and canalization of efforts. By having this information a priori, traditional surveys could target sampling of certain hosts and thereby expedite the process and potentially result in fewer individuals being sacrificed. Finding potential introduced species further demonstrates the usefulness of this method, as we can now return to this locality with targeted questions about these parasites, as well as inform the appropriate officials in charge of managing and mitigating introduced species.
In conclusion, we have demonstrated the productivity and feasibility of eDNA metabarcoding for parasites from a small set of environmental samples using a multi-marker approach, and its usefulness during ParasiteBlitzes by identifying >1,000 parasite ASVs corresponding to ~600 parasite OTUs across 4 aquatic habitats. We detected sequences from 9 named parasite species: only 2 of which had been reported previously from SC, with some reflecting gaps in knowledge of parasite diversity in SC and/or possible range expansions from the Gulf of Mexico, and others reflecting potential introductions into SC and/or new species lineages. We expect that eDNA metabarcoding of parasites will become more efficient and accurate as more primers are developed, protocols are optimized, and reference databases are populated by parasite sequences. There is a significant need to optimize the workflow from the point of collection to maximize recovery of all parasite groups. The fact that we did not detect some parasite groups in eDNA, despite observing them in the field, highlights inadequacies in primer coverage and reference databases, but also the complementarity of survey methods—none is perfect, and each has a unique and essential role. At this stage of our knowledge and the current status of reference sequences, using traditional methods that include description and molecular characterization of novel species in concert with eDNA metabarcoding during ParasiteBlitzes is necessary.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0022149X25100801.
Acknowledgements
Thanks to the College of Charleston (CofC) Foundation and Dr. Matthew Rutter, Academic Director of Stono Preserve, for access to Stono Preserve. William Roumillat and Pete Meier (CofC) assisted with sampling. Drs. Dominik Büchner, Till-Hendrik Macher, Daniel Grabner, and Kamil Hupało provided technical assistance and advice in the molecular laboratory. Dr. Florian Leese provided sampling equipment and supplies and reviewed a draft of this manuscript. This is contribution # 2506 of the CofC Department of Biology.
Financial support
This work was supported by a Department of Commerce NOAA Federal Award (I.dB. and K.M.H.-S., grant number NA22OAR4170114), a DIOS Incentive Fund Project, Hasselt University (M.P.M.V. and N.K., grant number DIOS/OEYLRODE/2022/001, contract number R-12947); the Special Research Fund (BOF) of Hasselt University (M.P.M.V., grant number BOF20TT06, and N.K., grant number BOF21PD01); a Brain Pool programme for outstanding overseas researchers of the National Research Foundation of Korea (S.G., grant number 2021H1D3A2A02081767); and the Tartar Research Fund, Department of Microbiology, Oregon State University (S.D.A.).
Competing interests
The authors declare none.
Open access
