The apportionment of dietary diversity in wildlife

Evaluating species’ roles in food webs is critical for advancing ecological theories on competition, coexistence, and biodiversity, but complicated by pronounced dietary variability within species and overlap across species. We combined dietary DNA metabarcoding, GPS tracking, and a machine-learning algorithm to cluster and compare dietary profiles within and among five migratory large-herbivore species from Yellowstone National Park. Interspecific niche partitioning was weak, but statistically significant (PERMANOVA: pseudo-F4,498 = 14.7, R2 = 0.11, P ≤ 0.001), such that some diet profiles from different species were as similar as those from within one species. Instead of affirming species’ identity as a primary determinant of diet composition, we found three statistically different clusters of diet profiles—one concentrated on graminoids and forbs, another on forbs and deciduous shrubs, and a third on gymnosperms—each including samples from all herbivore species. Clusters did not reflect traditional diet classification schemes such as the grazer-browser continuum that is often used to distinguish species by percent grass consumption or use of grassland habitat in African savannas. Instead, clusters in Yellowstone reflected seasonal dietary variation within species that often equaled or exceeded niche differences between species, contributing to our growing understanding of why environmental variability may favor generalist foraging strategies at temperate latitudes whereas specialized grazer and browser guilds appear to predominate in tropical savannas. Data-driven strategies that untangle complex trophic networks without relying on a priori groupings offer promising new insights into wildlife diets, with potential applications in resource management and environmental monitoring.

Methods
We obtained high-resolution diet profiles from fecal samples associated with GPS-collared individuals from five large-herbivore species that migrated over a ~100-km2 area. We processed a total of 697 samples from Yellowstone (inclusive of the 503 samples in this dataset and 194 separate samples from a different sampling scheme) with 38 total extraction blanks to monitor for laboratory cross-contamination. To identify dietary DNA sequences with the plant trnL-P6 marker using an Illumina MiSeq, we developed two reference libraries. The ‘local’ library comprised 334 unique trnL-P6 sequences from 882 vouchered specimens representing 75 plant families. We also built a global library using all trnL-P6 sequences available from ENA in EMBL-EBI using obitools v1.2.13. The global library was built using 26,101 unique trnL-P6 sequences obtained from 223,916 accessions representing at least 641 plant families.

To ensure that both per-base and per-sequence quality scores exceeded Q20, we first ran FastQC on the demultiplexed Illumina sequences obtained for each sample. We used cutadapt to remove primers from forward and reverse reads, allowing a maximum of 2 mismatches and no indels. We assembled sequence reads using the illuminapairedend in obitools using a minimum alignment score of 40. We used the obiuniq command to tally and merge identical sequences across samples, allowing us to quantify relative read abundance (RRA). We discarded sequences that occurred <2 times overall or that were outside the broad range of expected sequence lengths for the marker (<8 bp or >300 bp) using obigrep. We identified putative PCR artifacts as sequences that were highly similar to another sequence (1 bp difference) but with a much lower abundance (0.05%) in the majority of samples in which they occurred and discarded these sequences using the obiclean command. We used our local and global reference libraries to identify each unique plant DNA sequence with the ecoTagcommand in obitools and exported these outputs for further analysis.

Using these outputs, we crosschecked local reference library matches with global reference library matches. At this stage, we eliminated sequences that did not perfectly match a sequence in our local or global reference library (100% identity). We filtered the dataset to include the subset of 503 samples obtained from herds that included GPS collared individuals, as opposed to unaffiliated groups. All samples in this subset achieved >1000 sequence reads that passed quality controls. We thus rarified samples to an even sequencing depth for all downstream analyses and the final dataset comprised 503 samples and included 910 plant taxa; 95.7% of plant taxa were identified to family, 73% to genus, and 44% to species. We used these taxonomic assignments to characterize plant functional groups using growth form classifications of the United States Department of Agriculture (USDA) Plants Database.