SPECIES IN THE FECES: DNA METABARCODING TO DETECT POTENTIAL GASTROPOD HOSTS OF PARELAPHOSTRONGYLUS TENUIS Consumed By Moose (Alces Alces)

The Grand Portage Band of Lake Superior Chippewa is a federally recognized Indian tribe in northeastern Minnesota, US, and proudly exercises its rights to food sovereignty through subsistence hunting and fishing. Moose are a primary subsistence food used by the Anishinaabeg (people) of the Grand Portage Band historically and presently. Management for and research on maintaining this moose population as a vital subsistence species, thus, sets the context for this article examining pathogen transmission in this culturally important resource.

Multihost wildlife disease transmission systems are theoretically the most likely of all transmission systems to threaten host population viability (McCallum and Dobson 1995; de Castro and Bolker 2005; Weckworth et al. 2020) and among the most difficult systems to understand in enough detail to design interventions (Haydon 2008). Major challenges associated with studying multihost pathogen transmission include difficulty observing contact between hosts and lack of ecologic knowledge regarding one or more hosts (Buhnerkempe et al. 2015). New molecular techniques may overcome these challenges, but remain underapplied in wildlife disease research, particularly in multihost disease systems (Benton et al. 2015).

The technique of DNA metabarcoding opens new avenues for studying multihost pathogen transmission and disease ecology. Metabarcoding matches targeted high-throughput DNA sequencing reads to a curated reference sequence database, in order to identify species in an environmental or fecal sample (Deiner et al. 2017, Ruppert et al. 2019). Metabarcoding is often used in disease ecology to characterize the microbiomes or parasite communities of humans and animals in the diseased and nondiseased state (Gogarten et al. 2020; Mann et al. 2020; Willis and Gabaldón 2020). It also has been used to differentiate between similar-looking vector species, to detect pathogens in hosts, and to determine disease transmission risk (Bohmann et al. 2018; Huggins et al. 2019; Mechai et al. 2021). Despite the demonstrated utility of metabarcoding in understanding disease transmission, it is rarely used to determine multihost transmission pathways.

The persistence of moose (Alces alces), a species of subsistence importance to local Ojibwe tribes, is threatened in the north-central US by the multihost nematode Parelaphostrongylus tenuis (Lankester 2010; Feldman et al. 2017; Debow et al. 2021). In northeastern Minnesota, US, a 58% moose population decline occurred between 2006 and 2017 (Wolf et al. 2021; Giudice 2023), and the population has not recovered. Collaring studies revealed that P. tenuis infection accounted for more adult moose deaths than any other cause, making it a primary driver of the decline (Carstensen et al. 2018). This decline occurred almost entirely within the 1854 ceded territory, an area of over 20,000 km² in northeastern Minnesota, where local Ojibwe tribes retain treaty-reserved rights to preserve their cultures and subsist through moose hunting (Thompson 2017). Therefore, the loss of moose in this region would cause both ecologic and cultural damage.

The multihost life cycle of P. tenuis makes transmission to moose mechanistically complex, and our understanding of slug or snail (hereafter gastropod) intermediate hosts is particularly lacking due to limitations of traditional sampling methods. White-tailed deer (Odocoileus virginianus, hereafter deer) are the natural reproductive host of P. tenuis and shed larvae through feces. Parelaphostrongylus tenuis then infects the intermediary gastropod host, where it develops into an infectious stage before the infected gastropod is consumed by a deer, where the nematode develops to the adult stage and reproduces (Lankester 2001). Higher deer infection rates have been observed in land cover types that harbor more gastropods (VanderWaal et al. 2015); it is hypothesized that cervid ingestion of gastropods is accidental and that consumption rates are proportional to gastropod availability. Infection is generally harmless to deer, but moose, which are susceptible but not a natural host, can suffer severe neurologic damage (Lankester 2010). Previous research identified several gastropod species (Derocerus laeve, Discus cronkhitei, and Succinea ovalis, among others) that are common in moose range and consistently host P. tenuis (Lankester and Anderson 1968; Lankester and Peterson 1996). However, some gastropod host species may be yet to be investigated, because methods used in those studies underestimate both gastropod density and diversity (Hawkins et al. 1998). Furthermore, knowing whether species detected as P. tenuis hosts in previous studies are actually consumed by moose would be a key step toward designing interventions that reduce contact between moose and gastropod hosts (Severud et al. 2023).

We hypothesized that if moose were exposed to P. tenuis through consumption of the intermediate gastropod host, we should be able to detect the DNA of those hosts in moose feces. To determine whether metabarcoding is a viable method for detecting potential gastropod vectors of P. tenuis in the diets of moose and deer, we addressed the following questions: 1) Can metabarcoding detect gastropod DNA in deer and moose feces? 2) What species of gastropods do moose consume?

This study was conducted on the Grand Portage Indian Reservation (GPIR) and public land west of the reservation in Minnesota, US (Fig. 1). The GPIR covers 192 km² and is in Minnesota’s northeastern corner, bordered by Ontario, Canada, to the north and Lake Superior to the east and south. The reservation was established as part of a treaty between several Chippewa tribes and the US in 1854, which also included the ceding of 20,234 km² of land west and south of GPIR to the US government. As part of this treaty, the Grand Portage Band reserved usufructuary rights to hunt, fish, and gather on the lands ceded to the federal government (Thompson 2017). The reservation is located on the transition zone between the boreal and mixed conifer-hardwood forest biomes. Common tree species include sugar maple (Acer saccharum), white pine (Pinus strobus), white cedar (Thuja occidentalis), aspen (Populus spp.), balsam fir (Abies balsamea), and black spruce (Picea mariana). Moose density on the reservation averaged 0.26 individuals/km² (90% confidence interval=0.21–0.34) during the study (Oliveira-Santos et al. 2021).

To maximize the efficiency of moose fecal sampling, we used GPS collars (Vectronic Aerospace, Berlin, Germany) that had been deployed via helicopter capture (Heliwild LLC, Austin, Texas, USA) on 48 moose for other ongoing studies led by the Grand Portage Band. All capture and handling protocols were conducted in accordance with requirements of the University of Minnesota Institutional Animal Care and Use Committee (protocols 1410-31945A, 1601-33318A, 1812-36635A, and 1803-35736A). Details on capture, handling, and collaring are in Oliveira-Santos et al. (2021). Specifically, we mapped coordinates sent by the collars in the previous 24 h using Google Earth (Google Inc., Mountain View, California, USA) and identified locations where a single moose had three or more continuous location fixes that occurred within an approximately 30 m² area. Because this pattern of locations indicated that the moose had spent at least several hours in that area, and therefore was likely to have defecated there (Miquelle 1983), we then searched those locations for samples. In this fashion, we collected a temporally stratified set of fecal samples during June–October in 2017–20. No more than a single sample was collected from a collared animal in a month. Each sample consisted of at least 5–10 individual fecal pellets. We also collected opportunistic samples from uncollared moose and deer and included some opportunistic samples from May. Using collar data or visually assessing moisture content, we ensured samples were less than 2 d old at the time of collection. We recorded GPS coordinates of the location of the sample, the identification number of the animal (if known), the age cohort (calf or fawn, or adult), and date of collection. For samples collected from uncollared animals, we collected from pellet groupings of a different size, shape, or color, to avoid redundantly sampling the same individuals. All samples were collected with a fresh nitrile-gloved hand and placed in a sterile bag (Whirl-Pak, Madison, Wisconsin, USA). After collection, we immediately stored samples on ice until they could be moved to a −20 C storage freezer at the end of the day. In the laboratory, we removed the outer layer of feces using a scalpel to minimize the risk of detecting DNA from gastropods that might have moved across the surface of the feces. We assumed that digesta content would vary between pellets; therefore, we combined three pellets and homogenized them manually in a new sterile bag (WhirlPak). We refroze the samples and stored them at −20 C.

Using the QIAamp PowerFecal Pro Kit (Qiagen, Hilden, Germany), we extracted DNA from each fecal sample following the manufacturer’s protocol and using a QIAcube instrument (Qiagen). We quantified the DNA extract concentration with the Qubit dsDNA Broad Range Assay Kit (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Following quantification, we ran samples on a 1% agarose electrophoresis gel and visualized the extracted DNA with a gel imager (Thermo Fisher Scientific) to confirm adequate molecular weight for downstream procedures. After DNA extraction, we submitted samples to the University of Minnesota Genomics Center (UMGC; Minneapolis, Minnesota, USA). At UMGC the extracted DNA was cleaned with the AxyPrep Mag PCR Clean-Up Kit (Corning Inc., Corning, New York, USA) to remove PCR inhibitors. Using primers designed to amplify cytochrome c oxidase subunit I (COI) from a variety of invertebrate phyla, including Mollusca (LCO1490: 5′-ggtcaacaaatcataaagatattgg-3′ HC02198: 5′-taaacttcagggtgaccaaaaaatca-3; Folmer et al. 1994), UMGC performed PCR on a 710-base pair (bp) fragment of the COI gene. Specific reagent concentrations and thermocycling conditions were described in Gohl et al. (2016). Following PCR, UMGC cleaned the product with AMPure XP beads (Beckman Coulter, Brea, California, USA) at a concentration of 1.8× the sample to remove impurities. Then, UMGC assessed the success of the PCR reaction with Agilent TapeStation (Agilent Technologies, Santa Clara, California, USA) and only proceeded with library preparation on samples that had bands from 600–1,200 bp. During library preparation, UMGC added index sequences and Nextera adapter sequences to the PCR product following standard protocols. They used AMPure XP beads at a concentration of 0.65× to size select for fragments >300 bp. Aiming for ∼75,000 reads per sample, UMGC sequenced the samples with the Illumina MiSeq (Illumina, San Diego, California, USA) with a paired-end 300-bp run. All raw Illumina data generated herein are publicly available via the National Center for Biotechnology Information (GenBank) Sequence Read Archive (under project accession no. PRJNA872187).

After data generation, UMGC demultiplexed reads and removed index and adapter sequences with bcl2fastq v2.20 (Illumina). Because paired-end reads could not be merged due to a lack of overlap, we only used forward reads in our metabarcoding pipeline (Leff et al. 2015). We removed primer sequences with USEARCH v8.1 (fastx_truncate; Edgar 2010). Using FastQC software v0.11.5 (Andrews 2010), we then examined the quality scores for each sample. We truncated all reads to 260 bp with USEARCH to remove low-quality tails. Next, we quality-filtered our reads with USEARCH by setting the maximum number of expected errors to one, which corresponds with the most probable number of errors in a filtered read being zero (fastq_maxee 1; Edgar and Flyvbjerg 2015). We used VSEARCH v2.3.4 to dereplicate reads and used the USEARCH cluster-otus program to create 97% similar operational taxonomic units (OTUs; v2.3.4 Rognes et al. 2016). We discarded all OTUs that were based on a single read. Using a threshold of 97% similarity between the reference and query, we globally aligned the remaining OTUs to a reference database of COI sequences of gastropods. The reference database (Ratnasingham and Hebert 2007) contained all Gastropoda COI sequences deposited in the Barcode of Life Database (BOLD) as of 10 November 2021. This database included at least one species representative of 83% of the surveyed gastropod species in our study area (Lankester and Peterson 1996; Cyr et al. 2014; Severud et al. 2023) and at least one genus representative of the remaining 17%. We also tested a 90% similarity threshold between the reference and the query to minimize the chance of excluding species that only had a genus representative in the database.

Upon identifying aligned OTUs, we retrieved OTUs that were representative of the raw forward read and corresponding reverse read (535 bp total). We used the Basic Local Alignment Search Tool (BLAST) to compare the aligned sequence to all sequences in the National Center for Biotechnology Information database (National Center for Biotechnology Information 2022).

After obtaining a putative identification of the aligning reads with BLAST, we performed a phylogenetic analysis that explored the relatedness between our gastropod sequences and other published COI sequences within the genus to confirm our species identities. Specifically, we trimmed the primers from the forward and reverse reads and performed a multiple alignment using MAFFT v7.305 (Katoh et al. 2002), including sequences from GenBank and BOLD to represent the genetic diversity within the genus and outgroups. Using MEGA software (v10.0.5; Kumar et al. 2018), we found the best evolutionary model for each alignment. We used RAxML v8.2.9 with 1,000 bootstraps to generate a maximum likelihood phylogenetic tree (Stamatakis 2014). Finally, we visualized phylogenetic trees and examined bootstrap support for species identifications using FigTree v1.4.0 (Rambaut 2010).

We collected 258 moose fecal samples (n=9 in 2017, n=20 in 2018, n=114 in 2019, and n=115 in 2020; Fig. 1) and 48 deer samples (n=3 in 2017, n=38 in 2018, and n=7 in 2019). All deer samples and 54 moose samples were from individuals probably not radiocollared. Our final dataset included 22 samples from May (one moose, 21 deer), 57 samples from June (50 moose, seven deer), 68 samples from July (63 moose, five deer), 56 samples from August (45 moose, 11 deer), 41 samples from September (all moose), and 56 samples from October (54 moose, two deer). We successfully extracted fecal DNA from 297 samples (249 moose and 48 deer). We amplified detectable PCR product from 265 samples (234 moose, 31 deer), and after sequencing, our depth range was 1,653–138,697 reads per sample (mean depth=26,644 reads per sample; SD=36,183 reads). After quality filtering, the remaining read depths ranged from 5 reads per sample to 121,015 reads per sample, with an average of 12,820 reads per sample (SD=30,227 reads). Using quality-filtered reads, we generated an average of 171 OTUs per sample (SD=154 OTUs), with a range of 2–818 OTUs per sample. Eleven OTUs from 11 samples aligned to our database with >97% similarity. However, after performing a BLAST search to verify the species identities given by the reference database alignment, we found that one of the reference sequences labeled as Gastropoda actually represented human contamination (sample identification QHAK608-21 in BOLD). Excluding the samples aligning to the human-contaminated sequence, we observed sequences from five fecal samples that aligned to our reference database for gastropods. Lowering the alignment stringency to 90% did not result in any additional gastropod alignments.

Using BLAST, we confirmed the presence of Punctum minutissimum in two moose fecal samples collected on 14 September 2018 (535 bp; 100% shared identity; E-value=4e-132; 25 reads) and 5 July 2017 (535 bp; 100% shared identity; E-value=8e-134; two reads); Zonitoides arboreus in one moose fecal sample collected 17 September 2019 (535 bp; 100% shared identity; E-value=1e-131; seven reads); and Helisoma spp. from two moose fecal samples collected on 20 June 2019 (535 bp; 99.64% shared identity; E-value=1e-137; two reads) and 24 June 2018 (535 bp; 100% shared identity; E-value=2e-139; 22 reads; Table 1). The closest species level match for the Helisoma spp. detections was Helisoma trivolvis (89.45% shared identity; E-value=5e-91). We did not detect gastropods in the deer fecal samples.

Phylogenetic analyses confirmed our BLAST species identifications (Figs. 2–4). Bootstrap support for classification of our reads as Punctum minutissimum and Z. arboreus was 100% (Figs. 2, 4). The phylogenetic analysis did not provide further clarity on the species identity of the reads we classified as Helisoma spp. but provided strong support for our genus classification (bootstrap support=100%).

Our findings confirm that moose consume gastropod species known to become infected by P. tenuis. Zonitoides arboreus is a well-documented P. tenuis host (prevalence=0.4–12%; Lankester and Peterson 1996; Lankester 2001) and a common terrestrial snail in our study area. We also detected moose consumption of two additional gastropod species that have not been screened extensively for the ability to harbor P. tenuis.

We also demonstrate that fecal metabarcoding can be used to gain novel insight on interactions between hosts of a multispecies pathogen transmission system. Understanding the complex interactions in multihost pathogen transmission systems is a difficult but critical step in designing effective interventions (Haydon 2008). In particular, documenting cross-species transmission events in a field setting is crucial to transmission modeling efforts but may be impossible to observe with traditional methods (Buhnerkempe et al. 2015). Molecular techniques remain underapplied in wildlife epidemiology (Benton et al. 2015), and our study exemplifies how increasingly-accessible high-throughput sequencing can be leveraged to observe cryptic interactions between moose and potential P. tenuis gastropod intermediate hosts.

We did not determine if the Z. arboreus snail we documented carried P. tenuis. However, it is likely that moose eat enough forage to be exposed to infected Z. arboreus, even at low prevalence levels. Lankester and Peterson (1996) hypothesized that even at prevalence rates of 0.08%, deer fawns may ingest one infected gastropod every 51 d. Given significant differences in body size, moose must consume more than deer to meet daily energy requirements (Moen et al. 1997), thus, moose may occasionally consume infected Z. arboreus.

Understanding where moose consume Z. arboreus and other intermediate hosts is important for designing effective interventions such as altered timber management practices that change gastropod distributions and assemblages (Severud et al. 2023). Based on previous gastropod surveys in our study area (Lankester and Peterson 1996; Cyr et al. 2014; Severud et al. 2023), Z. arboreus is less common than other known P. tenuis hosts (e.g., Discus cronkhitei, Derocerus laeve, and S. ovalis) and rarely climbs vegetation (McCoy and Nudds 1997; Hawkins et al. 1998; Nekola 2003), making it less likely to be accidentally consumed by moose that are browsing live foliage. Moose may encounter Z. arboreus in forest types that contain paper birch (Betula papyrifera), because moose and deer spend significant time there, and those forest types support high Z. arboreus densities (Kearney and Gilbert 1978).

Punctum minutissimum, the other terrestrial gastropod we detected, is underrepresented in studies of P. tenuis prevalence despite its high abundance on the landscape (Nekola 2003). This may be a result of biased sampling methodology, as many previous studies used wet cardboard sampling methods that might not attract this species (Lankester and Peterson 1996; McCoy and Nudds 1997; Hawkins et al. 1998). Although Punctum minutissimum has not been screened extensively for its ability to host P. tenuis (Lankester and Anderson 1968; Platt 1989; Cyr et al. 2014), it shares the ground-dwelling trait that may make Z. arboreus a common host (Pilsbry 1939). Punctum minutissimum and moose both have documented preferences for recently logged or disturbed areas (Kralka 1986; Timmermann and McNicol 1988; Severud et al. 2023), so further research is warranted to better understand the P. tenuis host potential of Punctum minutissimum.

Helisoma, the only aquatic genus of snail we detected, is common in our study area and not known to host P. tenuis, but the potential role of aquatic gastropods in P. tenuis transmission has received little research attention. Exposure of aquatic snails to P. tenuis is possible, as L1 larvae leave deer fecal pellets and settle on the bottom of a water body when submerged (Duffy et al. 1999). Larvae can survive in water for several months (Lankester 2018), and water may even provide a dispersal mechanism (Lankester and Anderson 1968). Laboratory experiments found that L1 larvae could infect and mature into L3 larvae in the Lymnaea genus of aquatic gastropods (Anderson 1963), but the Helisoma genus remains unexamined for its competency as a P. tenuis host. Moose use wetlands and forage on aquatic vegetation most intensely in northeastern Minnesota during June and July (Van Ballenberghe and Peek 1971; Bump et al. 2017), which is when we detected Helisoma spp. consumption. Because moose might consume numerous aquatic gastropods during this time, the potential role of aquatic gastropods in the P. tenuis transmission cycle merits further research.

Based on previous hypotheses about P. tenuis transmission, the limited number of gastropod DNA detections in moose and deer feces is surprising. Lankester and Peterson (1996) hypothesized that deer fawns consume approximately 50 gastropods a day in autumn, yet we only detected five gastropod consumption events in moose and none in deer over 265 samples. In addition, we did not detect the gastropod species most common in our study area based on timed searches (Arion spp., 68% of all species; Severud et al. 2023). Unlike previous studies, our study was specifically designed to observe consumption events. Our limited detections indicate either that consumption events are less common than previously assumed, that species relative abundance on the landscape does not correlate directly with cervid ingestion, or that the sensitivity of our method for detecting consumption events is low.

Measuring and improving the sensitivity of this technique precludes recommending wider implementation and would inform future studies. We would expect DNA to degrade through the digestive system to some extent, which potentially would limit detection, particularly when consumption is a rare event or of low quantity. Thus, more work is needed to estimate the sensitivity of detection of the metabarcoding approach, potentially by feeding captive deer or moose known quantities and species of gastropods. If the test has low sensitivity, developing gastropod-specific primers that amplify a smaller fragment of DNA might greatly improve sensitivity (King et al. 2008). Importantly, this metabarcoding method potentially offers a way to directly observe the previously unobservable phenomenon of gastropod consumption, so pursuing a better understanding of the sensitivity and optimizing the protocol is worthwhile.

Insight about which gastropod species moose are eating, the host competency of those species, and their ecology could be used by forest and wildlife managers to design habitat modifications or other interventions that reduce the density of those gastropods in areas commonly used by moose. Although such actions alone may not prevent local extinction of moose, it would give managers another tool in a currently limited toolbox to understand P. tenuis transmission. The techniques presented in this article offer the potential to better understand a variety of gastropod-borne diseases, which are a major threat to human and animal health worldwide but often poorly understood (Giannelli et al. 2016). These techniques may also be used to better describe the life cycles of trophically transmitted parasites that can have important consequences for wildlife populations (LoGiudice 2003).

We thank the Grand Portage Band of Lake Superior Chippewa, Department of Biology and Environment, particularly Tyler Walters, Yvette Ibrahim-Chenaux, Roger Deschampe, Edmund Isaac, Eric Robley, Kimberly Teager, Madeline Grunklee, Krishna Woerheide, and Heather Fox, for assistance in collecting this extensive sample set. We further thank Maya Sarkar, David Myhre, Charlena Bergman, Suzanne Stone, and the University of Minnesota Genomics Center for help with preparing and sequencing these samples. We are grateful to Laramie Lindsey and Evan Kipp for providing guidance during the bioinformatic analysis and to Katherine Marchetto and William Severud for providing ecologic context for our findings. Funding for the research herein was provided by the Minnesota Environment and Natural Resources Trust Fund, as recommended by the Legislative-Citizen Commission on Minnesota Resources, the University of Minnesota College of Veterinary Medicine Minnesota Agricultural Experiment Station funds (MIN-62-099), the Van Sloun Foundation, the University of Minnesota Genomics Center, the Great Lakes Restoration Initiative, the Grand Portage Band of Lake Superior Chippewa, start-up funds awarded to P.A.L., the University of Minnesota College of Veterinary Medicine, and generous private donations.