Vertebrate animals are known to consume other species’ faeces, yet the role of such coprophagy in species dynamics remains unknown, not least due to the methodological challenges of documenting it. In a large-scale metabarcoding study of red fox and pine marten scats, we document a high occurrence of domestic dog DNA in red fox scats and investigate if it can be attributed to interspecific coprophagia. We tested whether experimental artifacts or other sources of DNA could account for dog DNA, regressed dog occurrence in the diet of fox against that of the fox’ main prey, short-tailed field voles, and consider whether predation or scavenging could explain the presence of dog DNA. Additionally, we determined the calorific value of dog faeces through calorimetric explosion. The high occurrence of dog DNA in the diet of fox, the timing of its increase, and the negative relationship between dog and the fox’s main prey, point to dog faeces as the source of DNA in fox scats. Dog faeces being highly calorific, we found that foxes, but not pine martens, regularly exploit them, seemingly as an alternative resource to fluctuating prey. Scattered accounts from the literature may suggest that interspecific coprophagia is a potentially frequent and widespread form of interaction among vertebrates. However, further work should address its prevalence in other systems as well as the implications for ecological communities. Tools such as metabarcoding offer a way forward.

Study area

The study was based in the northwest of the Cairngorms National Park, Scotland (57°09′34″N, 03°51′40″W). The study area, with a temperate oceanic climate, is covered by seminatural (Caledonian) forest and a mix of Scot’s pine (Pinus sylvestris L.) plantations and clear-felled areas (for details see Zalewska et al., 2021). Field voles are the predominant prey species for many carnivores in the area, including red foxes (this study), and undergo cyclic changes in abundance (Lambin et al., 2000). Small mammal surveys indicated low field vole densities for the duration of the diet study (2018 and 2019).

Sampling

To study the diet of mammalian carnivores, scats were collected along transects on unpaved forest roads and trails in three to six sites (see below). Sites were 6 to 11 km² and selected to maximise availability of transects and variability of habitats as determined by the proportion of clear-felled and other herbaceous vs forest habitats, and deer-culling pressure, which introduces carrion in the system. The minimum distance among neighbouring sites’ centroids ranged between 5 and 8 km. Sampling was done from February to April and from May to July in 2018, and again the same periods in 2019. These reflect ecologically relevant periods of food scarcity and the breeding season of prey of interest, to which we refer as winter and spring, respectively. Three sites were visited during the first season, six during the second season, and five during the following two seasons. The first season consisted of two visits, with samples collected in both, while the following seasons had three visits each where scats were cleared during the first one and collected in the following two to maximise collection of fresh samples. On average visits were 21 days apart, ranging from 12 to 48 days. Tracks and trails were scanned by one or two surveyors walking abreast. Scats were identified to probable species in the field by size and shape (Summers et al., 2015). Approximately 1 cm³ of faecal matter was collected from 2,887 samples using disposable wooden sticks and stored into 95% ethanol and then transferred to self-indicating silica or directly stored into silica. Additionally, 298 samples were collected by volunteers. These were frozen whole and 1 cm³ transferred to silica up to 18 months later. A further 232 samples were collected by a pre-existing and ongoing study in one of our study sites. These samples were also frozen whole and collected in a single visit coinciding with the last visit of our sampling seasons. From 3,417 available samples, 2,084 were selected for metabarcoding analysis to maximise spatial and temporal coverage of independent meals. The selected samples included all samples putatively identified in the field as fox (763) or Eurasian badger (Meles meles L.; 85), and a subset of the samples tentatively identified in the field as pine marten (Martes martes L.; 973), least weasel or stoat (Mustela erminea and Mustela nivalis L.; 213), and unidentified samples (50). For additional details on the selection of samples see Text S1 in Appendix 1.

DNA Metabarcoding

We analysed the vertebrate component of the diet of mammalian carnivores through DNA metabarcoding (Shehzad et al., 2012). Total DNA was extracted to a final volume of 400 μl using the NucleoSpin Soil Kit (Macherey-Nagel, Germany) after mixing samples with a sodium phosphate lysis buffer for 15 minutes (Taberlet, Prud’Homme, et al., 2012). The 12S mitochondrial rRNA gene was then amplified by triplicate using the universal vertebrate primer 12SV5 (Riaz et al., 2011). Amplification of fox, marten and badger DNA was only partly blocked through blocking oligonucleotides to allow host identification (Appendix 1 Text S2; Vestheim & Jarman 2008). Extraction and amplification were carried out in dedicated and separate rooms. Sequencing was outsourced to Fasteris (Geneva) and done using a NextSeq 500 (Illumina, USA).

Sequencing files were analysed using OBITools (Boyer et al., 2016). Sequences found only once in the dataset, containing degenerated bases or that were too short or long (< 60 bp, > 130 bp) were removed first. Molecular taxonomic units (MOTUs) that did not reach 10 reads in at least one PCR were removed too. The remaining were then taxonomically annotated by comparing against a local and a global reference database prepared from the EMBL’s European Nucleotide Archive (http://www.ebi.ac.uk/embl/). Only alignments with at least 95% sequence identity were kept in the final data. Contamination (mostly primate DNA) and tag jumps (i.e., sequences assigned to the wrong PCR during sequencing) were identified at this stage and removed. Samples were assigned to carnivore hosts through a “voting system” wherein a carnivore had to be the most common of all potential hosts in at least two of the three PCR replicates while representing at least 1% of the PCR’s reads. The resulting data was manually curated to address imperfect alignments, assignations to non-native or redundant taxa (e.g., Microtus and Microtus agrestis), or not at species-level. Where multiple MOTUs were assigned to the same taxon and present in the same PCR their reads were combined. Only detections that represented at least 1% average relative frequency of reads across amplifications replicates were used in later analysis (Deagle et al., 2019). For additional details on the laboratory, bioinformatic and manual curation process, see Appendix 1 (Text S2).

Analysis

Both percentage frequency of occurrence (proportion of positive scat samples; % FO) and a modified relative read abundance (RRA), average percentage of reads of a given taxon across positive scat samples, were used to summarise dog and field vole occurrence in the diet of fox and marten.

To test whether foxes exploited dog faeces as an alternative food-source to voles, we modelled the probability of occurrence of dog and field vole in the diet of fox as binomial random variables over time (the four sequential seasons fitted as a categorical variable), and regressed dog probability of occurrence in the diet of fox over field vole % FO, in three generalised linear models. Within each season, variables were aggregated at site level (three to six site levels per season). Sites were considered independent replicates after testing for spatial autocorrelation at sample and site level using the ‘testSpatialAutocorrelation’ function of the ‘DHARMa’ R package (Florian, 2021; R Core Team, 2021). Two models were used for this, one with a binary binomial variable (presence-absence of dog DNA) and sampling season as a categorical covariate and another one with the same variables but a random intercept of site (fitted with ‘glmer’ from the package ‘lme4’; Bates et al., 2022). Sampling season and field vole % FO covariates were fitted in separate models due to strong collinearity. Uniformity and overdispersion of the residuals and presence of outliers were tested using ‘DHARMa’. Quantile deviations were observed for the field vole model but were overall not significant and model assumptions were met. Model estimates are presented alongside their 95% confidence interval.

To address whether dog DNA sequences were artefacts of fox sequences, differences between the sequences of MOTUs assigned to fox (n = 6) and dog (n = 21) were calculated using the ‘adist()’ of the ‘utils’ package. The number of dog and field vole DNA sequence reads per PCR (n = 1,941; 647 fox scats x 3 amplification replicates) were plotted against the number of fox reads per PCR and fitted with 90th quantile regressions over the non-zero component using the ‘rq()’ function of the ‘quantreg’ R package (Koenker, 2021). The number of sequence reads (plus the smallest read count found in the data) were log transformed (base 10) before plotting and fitting the quantile regression. Data management, analysis and visualisation was done with R 4.0.5.

Calorimetry

To obtain a measure of the energetic content of fresh dog faeces, samples from six dogs and households that consumed a range of dry and wet, commercially available, dog food were analysed through calorimetric explosion (Hambly & Speakman, 2005). Samples were weighted before and after drying at 60° C for 14 days to calculate the percentage water content in the samples. They were then homogenised in a blender and compressed into pellets of 0.15-0.25 g and exploded in a Parr 6100 calorimeter using a 1108 oxygen bomb (Parr Instrument Company, USA) after calibration using benzoic acid. Each sample was replicated 2-4 times until the relative standard deviation was less than 1.5% of the mean value. Results are presented per MJ kg^-1 of dry and wet weight and as kilocalories per 100 g of dry and wet weight. Non-combustible mineral residual material was weighed for four of the six dogs and presented as the % weight of the wet and dry pellets. Wet weight of the pellets was back transformed from the dry weight using average water content of this study’s dog faeces (60.5%).

R and/or RStudio.