Natural history bycatch: a pipeline for identifying metagenomic sequences in RADseq data

Sample collection and preservation

We sequenced tissues from multiple sources. We collected tail tissue samples (skin, muscle, and cartilage) in the field from two species of horned lizards, Phrynosoma modestum (six individuals) and Phrynosoma cornutum (nine individuals). We collected the samples using heat-sterilized scissors to avoid contamination, and stored them in RNALater in the field. These samples were collected in southwestern New Mexico in the summer of 2015 (permit number 3,606 to M. Grundler, University of Michigan IACUC protocol number PRO00006234). Samples were stored at ambient temperature in the field (for up to a month) and at −20 °C after being returned to the lab. We also sequenced DNA from night lizard (genus Xantusia) liver tissues (Natural History Museum of Los Angeles County, accession numbers TC1002, TC1003, TC1006, and RLB5221). Collection protocol is not available for the museum samples. The first three samples are Xantusia vigilis liver samples collected in 2012, and the last is a Xantusia riversiana liver collected in 1972. We also collected cloacal swabs (sterile rayon swabs, MW113) from two ribbon snakes, Thamnophis sauritus, and one water snake, Nerodia sipedon, in southeastern Michigan in the fall of 2015 (collected under a Michigan Scientific Collecting Permit 9-16-2015 to I. Holmes). We prevented the swab from coming into contact with the environment, the sampler, or the skin of the animals. We avoided contact with the skin around the cloaca by gently applying pressure to the ventral surface of the animal just anterior to the vent. This pressure slightly everted the cloaca, exposing the mucous membrane and allowing us to insert the swab cleanly. We removed the swab and placed it in a sterile 2 mL vial, then broke the shank so that the cap could be put on. We handled the shank of the swab only above the portion that will be preserved in the vial. The samples were transferred to −20 °C storage within hours of collection. We also sequenced samples from whole digestive tracts of two Sceloporus jarrovi preserved in 95% ethanol and stored at −80 °C (permit number SP673841 to Robert M. Cox). To acquire the samples, we dissected out the total lower intestines. We filled a sterile pipette tip with 100 μL of distilled water. We inserted the pipette tip into the intestine section, and depressed the plunger to force the water through the intestine. We collected the wash in a sterile 1.5 mL vial. The samples were sequenced in the fall of 2015.

Laboratory protocol

We extracted total genomic and metagenomic DNA using Qiagen blood and tissue kits with a 12-hour incubation with proteinase-K prior to the spin column extraction. We used a double-digest RADSeq approach (Peterson et al., 2012), with the enzymes EcoR1 and Msp1 from New England Biolabs. We ligated barcoded Illumina adapters to the sticky ends left by the enzyme cuts, and used a PCR to attach barcoded Illumina primers to double barcode the sequence. We size-selected fragments with genomic inserts between 200 and 300 bp using a Pippin Prep cassette. We performed 125 bp paired-end sequencing the fragments on an Illumina HiSeq platform with V2500 reagents at the University of Michigan Sequencing Core.

Publicly available sequence analysis

We downloaded three doubledigest RADseq datasets from NCBI’s Short Read Archive (SRR1947260 to SRR1947262). They are from the coral snake Micrurus fulvius (Streicher et al., 2016). Details of preservation are reported in the original paper. The authors report that samples were liver, heart, shed skin, or scales, and were preserved in ethanol or stored at −80 F. Details of storage on a per-sample basis were not available. Samples were restricted using the enzymes Sbf1 and Sau3A1, and paired-end sequenced on an Illumina HiSeq 2500 platform (Streicher et al., 2016).

Sequence preparation

We demultiplexed the sequences using pyRAD, removed the low-quality sequences, and clustered reads within samples to 97% identity (pyRAD steps 1–3) (Eaton, 2014). We chose this clustering threshold because many microbial ecologists use a 3% difference in sequences to identify operational taxonomic units. We used the resulting fasta file of clustered sequences for each individual (the pyRAD *.edit file) for all further analyses. Any combination of sequence quality control and clustering programs can be used for this step, for example FastQC or Trimmomatic for filtering (Andrews, 2010; Bolger, Lohse & Usadel, 2014), or vsearch for clustering (Rognes et al., 2016). For the Phrynosoma and Xantusia samples, we continued the pyRAD pipeline to cluster reads across individuals, and used the resulting “*.loci” file for further analyses. In the pyRAD *.loci files, the sequences for each locus are listed in a group, with the individual that provided the sequence identified in the name of that sequence. A standard line break string separates the sequences for each locus. We used a custom R script to take the first sequence for each locus and combine them into a fasta file to be passed to our analysis pipeline (Appendix 1).

Investigating metagenomic sequences

We use NCBI’s discontinuous megablast algorithm to compare all sequences from the *.edit and *.loci files to reference sequences in the online NCBI nucleotide database (Camacho et al., 2009). We use the R package taxize to find the genus and species of each sequence (Chamberlain & Szocs, 2013; Chamberlain et al., 2016). We discard results that aligned to more than one kingdom or phylum with greater than 80% identity. To assess how the threshold for similarity affected the number of sequences that can be identified to phylum, we imposed percent similarity thresholds to the closest matching sequence of 70%, 80%, 85%, 90%, 95%, and 97%. To assess the distribution of parasite sequences across hosts, we screened sequences that clustered across individuals of two horned lizard species, P. cornutum and P. modestum. We performed a similar analysis on samples from two Xantusia species. Finally, we built rarefaction curves using the R package vegan to identify the depth of sampling necessary to identify all genera of metagenomic DNA present in the sample (Oksanen et al., 2018). For each genus for which we have tissue metagenomes (Phrynosoma, Xantusia, and Micrurus) we created a community matrix in which samples are rows and columns are the the number of Chordata sequences, and the number of sequences in each genus of blood parasite. We set a 90% identity match for this analysis.

Limitations and caveats

The major limitation of our approach is that it relies on public databases to determine the taxonomic identity of sequences. However, public databases do not accurately reflect the diversity of metagenomic taxa. For example, we recover relatively few Archaea sequences from our hind-gut samples. We hypothesize that this reflects the relative lack of Archaea genetic sequences in GenBank to compare against, rather than an absence of Archaea from our samples. However, the number of publicly available, taxonomically identified reference sequences is quickly increasing, so this source of bias should be reduced in the future. Second, our identifications are based on randomly restricted DNA samples, rather than widely-accepted barcode sequences. These sequences can’t be corrected for copy number variation, and we have little to no ability to determine whether two different sequences represent different individuals or whether they are two separate samples of the genome of a single individual. Due to the inherent stochasticity of the ddRAD approach, this method should not be used to quantify the relative or absolute abundance of metagenomic communities. Finally, the sequences in this paper are from an Illumina HiSeq platform. Negative controls are not recommended on this platform, as running one leads to uneven ratios in barcode sequences, which can damage the sequencing quality for the entire run. The lack of negative controls is one reason that this approach should be considered exploratory, rather than as a method for quantifying microbial load in specimens. Well-designed primer sets can account for all of these problems, and should be used to answer questions about relative abundances of metagenomic lineages or community structure.

We note some unexpected taxonomic identifications in our sequences. Specifically, we note that some of the sequences align with highest percent identity with arthropod sequences, or with Streptophyta. These sequences may represent sequencing error that alters a highly conserved host or metagenomic sequence to erroneously align more closely with a non-target sequence. Alternatively, sequences could be labelled incorrectly in the NCBI database, or they could be contamination in genome assemblies (Orosz, 2015). Regardless of the source, this taxonomic error indicates that our method should be used for exploratory purposes only.

Ecology of metagenomic sequences

All of the putative parasite taxa that we can identify with 97% certainty are known parasites of vertebrate hosts. Protopolystoma xenopodis is known from African clawed frogs, in which it attaches to the kidney and feeds on blood, thereby potentially releasing its own DNA into the host’s bloodstream (Theunissen, Tiedt & Du Preez, 2014). Diphyllobothrium latum is known from the digestive tracts of a range of vertebrates, including mammals and fish (Wicht et al., 2009; Schurer et al., 2016). Nippostrongylus brasilensis and Strongyloides stercoralis are nematodes known from mammals. Their lifecycle begins with free-living juveniles that find a host and bore into the bloodstream through the skin. The juveniles migrate to the lungs, where they develop into adults before entering the digestive tract to breed (Koutz & Groves, 1953; Haley, 1961). Soboliphyme baturini is known from mammals, and infects the stomach (Zarnke et al., 2004). Elaeophora elaphi occurs in red deer, where it lives in the portal vein near the heart (Carrasco et al., 1995). While the sequences we detected are probably not the same species as their closest match, they should be closely related, and are likely to have similar life histories. All of the life histories here indicate that parasite DNA could plausibly be shed into the bloodstream.

We found many families of bacteria that are known from vertebrate guts, and some that have not previously been recorded. Actinobacteria, Acidobacteria, Cyanobacteria, Fusobacteria, Spirochaetes, Synergistetes, and Tenericutes have been reported from wild snake hindguts (Colston, Noonan & Jackson, 2015). Plancomycetes and Verrucomicrobia have been found in the guts of wild apes (Yildirim et al., 2010). Euryarchaeota, Deinococcus-Therums, Thermatogae, and Fibrobacteres have been found in dog gut microbiomes (Swanson et al., 2011). Chloroflexi has been recorded from human guts (Campbell et al., 2014), and Synergistetes has been recorded from the gut of young calves (Li et al., 2012). To the best of our knowledge, the bacterial phyla Ignavibacteriae and Armatimonadetes, and the fungal phylum Entomophthoromycota have not previously been reported from vertebrate hind microbiomes. A species in Entomophthoromycota has been found in a cyst in the esophagus of a rat snake, Elaphe obsoleta (Dwyer et al., 2006). Other Entomophthoromycota are pathogens of invertebrates (Gryganskyi et al., 2012), indicating that they are capable of invading multicellular hosts. Ignavibacteriae is a sister phylum to Bacteroidetes and Chlorobi, both known from gut microbiomes (Podosokorskaya et al., 2013). The phylum has been sequenced from wastewater, indicating that it can survive in organic waste (Meng et al., 2015). Armatimonadetes is primarily a soil phylum, but also participates in plant rhizobial communities (Tanaka et al., 2012), and has been found in mosquito salivary glands (Sharma et al., 2014) and decomposing swine manure (Tuan et al., 2014). The three new phyla are all reasonable candidates for the gut microbiome, as they are known to occur within multicellular host tissues. However, caution should be exercised because both the known and novel taxa we identified from the gut samples can also be found in environmental samples. Further study, using carefully selected barcode loci, should be undertaken before these taxa are considered an established part of the gut microbiome.