Comparison of seven DNA metabarcoding sampling methods to assess diet in a large avian predator

Description

Notes

Methods

We collected diet samples from Buteo lagopus throughout their North American nonbreeding range in the conterminous United States and Canada during the nonbreeding season (November – March) from 2020 – 2023. We sampled diet from captured birds via: 1) talon swabs, 2) beak swabs, 3) cheek (i.e., buccal cavity) swabs, 4) cloacal swabs, 5) cloacal fecal loops, and 6) fecal samples. We also collected 7) fecal samples from uncaptured free-ranging birds. In total, we collected 592 samples from 189 individuals including 113 captured and 76 uncaptured Buteo lagopus.

All samples were extracted in a laboratory dedicated to low quality and quantity DNA samples.  No forms of high-quality DNA were handled or stored in this laboratory. We collected DNA onto three "substrates" as described above: nylon bristle swabs (beak and talon samples), foam swabs (cheek and cloacal samples) and homogenized feces (fecal samples, fecal loops). Each substrate had its own DNA extraction protocol, but the protocol was the same for each substrate regardless of the source of the DNA. DNA was extracted from nylon bristle swabs (beak and talon samples) using the Qiagen DNeasy Blood and Tissue Kit.  Samples were vortexed for 15 s to dislodge prey remains cells from the bristles. The DNA/RNA Shield (1 mL) was then split between two new 2 mL microcentrifuge tubes (500 mL in each tube). After overnight incubation, extraction volumes were scaled up to match the amount of starting material. The entire extraction volume for a sample (~3.0 mL) was then spun through one spin column to combine the extraction back into one tube. The extraction was then completed according to the manufacturer's protocol. DNA was extracted from foam swabs (cheek and cloacal samples) also using the Qiagen DNeasy Blood and Tissue Kit, as has been done in previous diet metabarcoding studies for cloacal swabs in sharks (Clark et al., 2023) and for esophageal and cloacal swabs in Tyto alba (barn owl; Elmore et al., 2023). The swab and 600uL of DNA/RNA shield were transferred to a new 2 mL microcentrifuge tube. Like the bristle swab, extraction volumes were scaled up to match volume of starting material (~1.8 mL) and the entire extraction volume was spun through a spin column. The swab tip was carried through to the first centrifugation step and placed in the spin column to ensure all DNA present made it onto the filter. DNA was extracted from fecal samples using the QIAamp Fast DNA Stool Mini Kit.  Samples were vortexed then 200 mL of homogenized sample was used in the extraction. An extraction negative containing no sample was included in each extraction to monitor for potential contamination.

We targeted the V5 region of the 12S Ribosomal RNA gene. We prepared 14 libraries. All samples were run in duplicate. Additionally, each library contained two negative controls (nanopure water instead of DNA) and two positive controls, a low-concentration equimolar mix of tissue-derived DNA from two vertebrate species unlikely to occur naturally in our samples: Ovis canadensis (bighorn sheep) and Puma concolor (cougar).

We used a two-step metabarcoding library preparation, which consisted of two rounds of PCR. Two-step library preparation methods are popular, have been used for a variety of diet studies (Goldberg et al., 2020; Bourbour et al., 2021), and are generally modified versions of Illumina's amplicon sequencing guidelines (Illumina, 2013). The first round of PCR, the amplicon PCR, amplified the target gene region using primers that also contained overhangs for Illumina sequencing adapters, which were then added in the second round of PCR, the index PCR. We used custom-designed fusion primers for the amplicon PCR. The primers were based on those designed by Riaz et al. (2011) for vertebrates but had small additional heterogeneity spacers (one base pair) to increase complexity, which can improve sequencing yields (Jensen et al., 2019). Additionally, we synthesized two versions of the reverse primer, with different 5' initial bases to better capture prey richness. 

All sequencing was performed at the University of Oregon Genomic and Cell Characterization Core Facility (UO GC3F) in Eugene, Oregon. Libraries were run in three rounds: round one contained our pilot, round two contained libraries 2-7 and round three contained libraries 8-14. Each library was afforded approximately 15% of a PE 150-bp Illumina NovaSeq 6000 SP sequencing lane. We did not use a host-blocking primer because 1) we wanted to maximize the recovery of avian prey taxa that are not well differentiated at 12S, and 2) designing and validating a primer specific to rough-legged hawk was outside the scope of this project. We anticipated high host read counts, particularly within cheek, cloacal, and fecal samples and consequently increased our mean reads per library to maximize prey read recovery.

Bioinformatics

Samples were demultiplexed by the UO GC3F. Next, adapters and primers were trimmed using cutadapt (Martin, 2011). From there, samples followed the pipeline found here: https://github.com/sckieran/Metabarcoding_Pipeline_Waits, including all scripts used. Briefly, forward and reverse reads were merged with PEAR (CC licensed: https://cme.h-its.org/exelixis/web/software/pear/). Next, unique reads (ASVs, no error correction) were pulled from each sample using FASTX-collapser (open source: http://hannonlab.cshl.edu/fastx_toolkit/index.html) and singletons and doubletons were removed to reduce noise. We then filtered sequence data following the general approach outlined in De Barba et al. (2014). First, we removed low frequency noise by relative read abundance and read count. We removed any sequence shorter than 50 bp after merging, and any sequence that was not present at a rate above 0.1% the total read count in any sample. In total, we removed 120,319 ASVs before BLAST. However, to ensure that our filters were appropriately conservative, we separately BLASTed all ASVs removed in this way. A summary of the removed ASVs, including mean and median read count and number of samples that ASV was found in, can be found in the Supplemental Material (section 3; Table S2) along with further details about the filtering methods. Overall, there were no taxa removed from any sample in our final dataset due to the pre-BLAST ASV filters, suggesting that we successfully reduced noise without impacting downstream detections. We built a custom local BLAST-formatted reference database that included all potential prey taxa and several common contaminants and controls. We queried each unique sequence in our filtered samples against our BLAST database and assigned it the lowest unambiguous taxonomic rank present in the top-scoring BLAST hit by bitscore. We then removed ASVs matching at <98% identity. Within each sample, we removed ASVs present in <2 PCR replicates and built consensus ASV profiles between replicates by taking the mean of ASV read counts. We then examined PCR negative read counts and filtered each sample to either: 1) remove ASVs with fewer read counts than the maximum observed in any PCR negative sample (filter1; read count threshold = 4,379); or 2) remove sequences with fewer read counts than the maximum non-host or non-human taxa in any PCR negative sample (filter2; read count threshold = 449). We took additional steps to ensure the geographic and taxonomic accuracy of taxonomic assignments (see Supplemental Material, section 3 for additional details regarding taxonomic assignment and sequence filtering).