The study was conducted in the Coronilla Hill within the “Reserva Natural Salus” (34° 23’37.8168” S; 55° 19’43.05” W) of south-central Uruguay ( Figure 1 ). This is a 13 km² hilly area characterized by a mosaic of woodlands, shrublands, and grasslands. Fecal samples were collected opportunistically from various locations in the study area during 2015-2016, in both cold (June) and warm (January-April) seasons. The field surveys were conducted over three days, and five field assistants were involved in feces collection. Only fresh carnivore feces, identified by their low-humidity signs and morphogenic characteristics, were sampled. These samples were stored in individual plastic tubes in the freezer until DNA extraction to avoid contamination from other environmental sources (Pompanon et al., 2012; De Barba et al., 2014). The three canid species inhabiting the study area are crab-eating fox (C. thous), pampas fox (L. gymnocercus) and domestic dog (Canis lupus familiaris).

We extracted the DNA as soon as the samples arrived at the laboratory (3 to 5 days after collection) using the DNA Stool™ Minikit (Qiagen, Hilden, Germany), following manufacturer’s instructions. The feces were disintegrated using sterile tweezers and scissors, thus taking parts of the interior and exterior of the feces in order to homogenize the material for DNA extractions (De Barba et al., 2014).

Donor species were identified using a three TaqManTM probe real-time PCR assay (Cosse et al., 2017), with species-specific dyes included for crab-eating fox C. thous, pampas fox L. gymnocercus, and domestic dog C. lupus familiaris detection. PCR amplifications were conducted on Rotor-Gene ™6000 (Corbett) with SensiFastprobe™ PCR master mix (Bioline, London, UK) following Cosse et al. (Cosse et al., 2017) procedures.

A fox-blocking oligonucleotide was developed following De Barba et al. (2014) and Vestheim and Jarman (15) to avoid host DNA amplification from feces. The fox-blocking oligonucleotide developed with a phosphate (PHO) 3’ extreme modification (5´-CCACTATGCTTAGCCTTAAACGTAAATGATTYTAT PHO- 3’) was included in the 12S fragment PCR from fecal DNA (Section 2.2), using the primer set 12S-V5, which amplifies a fragment with an average length of 98 bp (Riaz et al., 2011). The PCR consisted of 10µl final volume with 0.04U of Taq polymerase (Invitrogen, USA), 1X Taq buffer, 1.5mM of MgCl, 0.2mg/mL of BSA (Bovine Serum Albumin), 0.05mM of dNTPs, 0.1 µM of forward and reverse primers, 1.6µM of fox blocking primer and a 4-8ng/µl final concentration of DNA. The thermal profile started with a denaturalization step of 94°C for 3 minutes, followed by 35 cycles of 94°C for 30 seconds and 58°C for 90 seconds. We used 40 to 45 cycles in degraded fecal DNA. Positive and negative controls were included in each PCR reaction. A second round of indexing PCR amplifications was conducted to add the Nextera indexes for sequencing on a MiSeq Illumina™ platform. The sequencing library preparation followed the 16S sample preparation guide for the Illumina MiSeq system. Each PCR from the same fox species and taken in the same season shared the same index.

Amplification success was confirmed by visualizing PCR products in 2% agarose electrophoresis using Goodview™ (SBS Genetech, Beijing) nucleic acid stain. PCR products were then purified using 0.8U/µl of Exonuclease I and 0.2U/µl of Thermosensitive Alkaline Phosphatase (FastAP) enzymes (ThermoFisher Scientific™, USA); the thermal profile was 37°C for 1 hour and 30 minutes followed by 5 minutes at 75°C. PCR-purified products were quantified using a Nanodrop TMND-1000 UV-vis Spectrophotometer (Nano-Drop Technologies, Inc., Wilmington, DE, USA) and then pooled in equimolar concentrations. Thus, four pools were generated that included libraries prepared from crab-eating and pampas foxes’ feces collected in warm and cold seasons.

Paired-end (2X 100) sequencing was performed using the Illumina MiSeq platform (Illumina, Inc., San Diego, CA, USA), available at the sequencing facility of the Institut Pasteur, Montevideo. The raw reads obtained were trimmed based on quality (Phred score ≥ 30) using SICKLE software (Joshi and Fass, 2011). PCR primers were removed, applying a locally developed script. Paired-end reads were merged using the PEAR program (Zhang et al., 2014), and the consensus sequences were determined.

The reference database was built using the available sequences of the 12S rRNA mitochondrial gene in GenBank (www.ncbi.nlm.nih.gov) in December 2020, including complete vertebrate mitochondrial genomes and partial sequences. The final dataset comprised 306,362 sequences from 28,270 different vertebrate species. Complete taxonomy information for each sequence included in the reference database was retrieved from Genbank. The reference database used is available at Zenodo (zenodo.org/records/10731091), and the raw reads are accessible from the SRA database under BioProject accession PRJNA1082626.

The pipeline described was completed for each experiment and was implemented through locally developed scripts. A detailed description of the bioinformatics workflow was included as supplementary data ( Supplementary Figure S1 , Supplementary Data 1 ). First, non-redundant consensus sequences (consensus from paired-end reads) were used for taxonomic identification based on Blast against the local reference database. Thus, each non-redundant consensus sequence was used as query and, the most similar reference sequences (best hits) were defined using the highest blast score. Second, taxonomic information for each best hit was retrieved and further considered. Finally, the list of Molecular Operational Taxonomic Units (MOTUs) comprises species derived from the best hits described in the studied region with a minimum of ten supporting reads. Additionally, higher taxonomic levels such as Genera, Family, or Order were assigned when matched species in the database do not occur in the study area, but there is a species of the same taxon in it. Following Turon et al. (2020), we used networks to analyze observed patterns of variations in generated reads. The reads obtained from the same MOTU in each pool were aligned using MUSCLE (Edgar, 2004). Then, networks were built using the minimum spanning method (Bandelt et al., 1999) in PopArt (Leigh and Bryant, 2015), as described in Supplementary Data 2 . Variants with low abundance and only one-point mutations from highly abundant central variants are classified as non-biological (Moreno et al., 2016; Turon et al., 2020). These variants, which arise due to PCR bias, sequencing errors, or degradation (Turon et al., 2020; Zizka et al., 2020), occur punctually in different reads and consequently are highly dispersed in the network and have low relative frequencies. Conversely, highly abundant sequences in central network positions (centroids) are considered genuine biological sequences.

To further support our distinction between erroneous sequences and real biological intra-MOTU variants, we analyzed the pattern of observed mutations in conserved and non-conserved positions. Highly conserved positions within the 12S rRNA sequence were identified by referring to available rRNA secondary structures on the Rfam website (rfam.xfam.org) (Kalvari et al., 2020) ( Supplementary Data 3 ). We reasoned that these specific sequence positions remain stable due to the constraints imposed by purifying selection on the structural integrity. Therefore, variations observed within conserved regions were deemed more likely to be errors, and were utilized to estimate the frequency of noisy variants per site. We estimated the overall mean p-distance for conserved and non-conserved positions using MEGA X (Tamura et al., 2011). Finally, we compared the network reconstructions (see above) and mean p-distance for conserved and non-conserved positions for the most frequent food species identified in each pool.

The Blast search results show a distribution of identity values of best hits ranging from 93 to 100% ( Supplementary Figure S2 ). Mammals showed higher read and variant counts across all cases. During the cold season, crab-eating fox feces revealed the presence of cavy (Cavia aperea) and cow (Bos taurus). A wide variety of toads taxa were detected in the crab-eating fox feces. Specifically, during the warm sampling period, the genus Rhinella was detected, while Pseudis, Leptodactylus, and Physalaemus were identified only in the cold season samples ( Table 1 ). Several bird food items were also detected in crab-eating fox feces. The food items with the highest number of reads during cold and warm seasons were Cavia aperea and Lepus europaeus, respectively ( Table 1 ). On the other hand, in pampas fox feces, the food items with the highest number of reads were Cavia aperea during the cold season and Mus sp. during the warm seasons sampling ( Table 1 ). Other items identified in pampas fox feces include three armadillo species, all of which are present in the study area. Additionally, only one amphibian species and one Passeriform bird MOTU were identified in pampas fox feces during the cold season ( Table 1 ).

Most food items consumed by crab-eating foxes were mammals, both in cold and warm seasons (78.5 and 60.3% of reads, respectively). Also, in this species, amphibians represented 16.6% of the reads in the cold season, and birds 30.7% in the warm season ( Figure 2 ). In pampas fox, mammals were the substantial majority of food items consumed, 99.1% of reads in the cold season and 96.4% in the warm season ( Figure 2 ). For both fox species, rodents constituted the most prevalent order among mammals. During the cold season, 67.0% of reads belonged to the Caviidae family, while in the warm season, 51.8% were attributed to the Muridae family ( Figure 2 ).

The sequence networks obtained in almost all cases showed a typical star shape, with a unique most frequent central sequence variant and less frequent sequence variant located in the periphery separated by few changes from the central one ( Figure 3 , Supplementary Data 2 ). Our analyses reveal no difference in the estimated variability (measure as overall mean p-distance) between the structural conserved and the non-conserved sequence regions. Moreover, no main differences can be observed between the network generated using both conserved and non-conserved sequences ( Figure 3 ). This result suggests that intraspecific variation observed in food items comes from erroneous reads. Taken together, we conclude that there is only one real biological variant identified in C. aperea in the cold season on both fox species. Also, in the warm season, only one natural variant could be considered in L. europaeus on crab-eating fox and only one in B. taurus on pampas fox ( Figure 3 ).

The experimental design used in this study for carnivore species detection and diet assessment proved to be sensitive enough to allow qualitative trophic evaluations (Buglione et al., 2018). We also present a novel and strong blocking oligonucleotide for pampas and crab-eating fox for the 12S rRNA gene fragment. The use of this kind of oligonucleotides is an issue that needs to be addressed since the predator DNA is much more abundant than the consumed taxon DNA. Despite the advantages of the speed, accuracy, and potential of the HTS approach for dietary studies, some considerations should be taken into account (Pompanon et al., 2012). The taxonomic resolution of the DNA fragment chosen should have broad taxonomic coverage, suitable discrimination power, and high amplification efficiency. The latter is usually linked to short amplicon length (100-250bp) (Pompanon et al., 2012). In this sense, the DNA fragment chosen here fulfills all these conditions and has been previously successfully used in several carnivores and omnivorous metabarcoding feeding studies in which they evaluated its usefulness and efficiency (Ficetola et al., 2009; Roemer et al., 2009; Shehzad et al., 2012; De Barba et al., 2014). Taxonomic bias or underrepresentation of local species in the reference database may affect the analysis, as can be observed in the distribution of MOTUs among the identified food items ( Figure 2 ). Thus, it is mandatory to improve and increase the quality and taxonomy coverage of the sequence database, including local species data. Meanwhile, the definition of MOTUs at higher taxonomic levels may help explore diet composition based on available data on public databases (Pompanon et al., 2012). Several DNA sequences for the 12S rRNA gene from vertebrate species inhabiting the study area were not found in the GenBank database. However, in some cases, the identification of taxon at the Genera level was possible (e.g., Pseudis and Colaptes). A similar approach was used when two or more species in the database displayed equal similarity scores due to the insufficient resolution of 12S rRNA for species assignment. In some cases, the identification at the Genera level was enough since only one species of the group is present in the study area (e.g., Cavia aperea).

Both false-positive and false-negative results are frequent in diet metabarcoding approaches (Taberlet et al., 2012; De Barba et al., 2014; Corse et al., 2019). To prevent these methodological shortcomings, we used universal PCR primers, blocking oligonucleotide primers fox species and included negative controls, and as a result, relatively low-frequency taxa were detected (De Barba et al., 2014; Corse et al., 2019).

In line with published diet studies that rely on morphological characteristics of taxon remains from the feces of both fox species, this study also identifies mammals (rodents and hares) as the most common vertebrate food items. Cavies (Cavia aperea) were detected in almost all seasons and consumed by both fox species. The only exception was seen during the warm season in crab-eating fox ( Table 1 ), where minor consumption of meat animal proteins may occur because large amounts of fruits are accessible and there are less energetic requirements (Di Bitetti et al., 2009). Also, in accordance with our results, published trophic analyses concluded that rodents, particularly Cavies, were the main taxon consumed by these foxes (Farias and Kittlein, 2008; Di Bitetti et al., 2009; Bossi et al., 2019). The brown hare, Lepus europaeus, is also one of the most abundant food items identified, and species from this Genera have been previously described as important prey for both fox species (García and Kittlein, 2005; Farias and Kittlein, 2008; Bossi et al., 2019).

Interestingly, some dietary composition differences were detected between fox species ( Table 1 ). As reported by Bossi et al. (2019), several amphibian taxa remains were present, especially in crab-eating fox feces, while different armadillo species are specific prey items for pampas fox. Birds were found as food items in both fox species, but only the Passeriform Order was identified in pampas fox. Moreover, during the warm season, Mus sp. was found in pampas fox feces in the study area near houses; this is not surprising, given that this rodent species frequently inhabits areas near human settlements.

Despite the smaller number of crab-eating fox feces sampled, many more taxa were detected in their diet, even those previously reported in other areas of its geographic distribution (Di Bitetti et al., 2009; Bossi et al., 2019). Although the limitation of the pooling strategy used here and the consequences of the relatively low number of samples included, the metabarcoding approach represents an improvement from the macroscopic trophic analysis. In the latter, the comparison of the results from different studies is challenging, mainly due to differences in taxonomic level prey resolution (Di Bitetti et al., 2009; Pompanon et al., 2012).

Cow DNA was identified in three of the four pools analyzed. The only exception was crab-eating fox during the warm season, the pool with a reduced number of collected feces ( Table 1 ). This observation is not surprising since the study area is a Reserve with a popular picnic zone with barbecues. Thus, foxes could be consuming leftovers because predation attacks on cows are improbable (Farias and Kittlein, 2008). Moreover, livestock consumption was earlier related to their scavenger habits (Farias and Kittlein, 2008; Bossi et al., 2019). The presence of livestock in fox’s diet needs to be further analyzed and clarified because canids are persecuted by predation of domestic animals (mainly sheep) and are recurrent targets of poisoning and hunting by ranch owners (García and Kittlein, 2005; Bossi et al., 2019). Note, however, that no sheep DNA was found among food items in our study. This result is highly significant for conservation biology because these fox species are considered livestock predators in the countries where they occur. Unfortunately, additional studies in other areas dedicated to sheep breeding, including methodological approaches that distinguish carrion consumption from predation, are necessary. For this reason, future feeding metabarcoding DNA studies should include new sampling and molecular methods to elucidate trophic interactions such as secondary predation, scavenging, coprophagy, or parasitism.

The study of the genetic variability of food species was achieved through network reconstruction ( Figure 3 ). This methodology utilizes the frequency and network position of detected variants to discriminate between biological variability and variability generated from erroneous reads (see above). The development of new methods to distinguish between both sources of observed variation is much needed to improve the accuracy of metabarcoding based on HTS. This is particularly important in environmental samples, where DNA is typically of low quality and degraded. Therefore, PCR cycles or sequencing depth need to be increased to achieve amplification success (Adams et al., 2019). To identify genuine genetic variations, it is advisable to incorporate positive controls, utilize multiple replicates, and employ high-quality filtering (Adams et al., 2019). Additionally, the reduction of erroneous sequences can be achieved through data denoising by in-silico methods (Tsuji et al., 2020; Zizka et al., 2020). Unfortunately, low abundant real biological variants could also be lost (Elbrecht et al., 2018). In this work, variations observed in sequence sites that are part of the conserved secondary structure of rRNA helices among vertebrate species, are presumed to be primarily errors. The null or marginal differences found in the estimated p-distances between structurally conserved and non-conserved sequences suggest that most of the variation observed is mainly random and comes from erroneous reads ( Supplementary Data 3 ). Thus, we can state that only one variant was detected in all cases; the most frequent and central sequence in each network reconstruction should be considered the only real biological variant ( Figure 3 ). There are non-mutually exclusive explanations for this lack of variability in the consumed taxon identified. First, our study was focused on a small area, and the number of samples per pool was low. Therefore, the low genetic variability observed may characterize the populations analyzed in the study area. Second, the DNA fragment used here, which allows species and taxon discrimination, may be too conserved for intraspecific genetic analysis (Riaz et al., 2011). This is a primary approach, further analyses, including non-coding intraspecific variable DNA fragments, should be done to confirm the efficiency of this strategy. In addition, expanding the database of sequences of neotropical species for the 12S rRNA enables us to better understand natural variability through traditional population genetic analysis methods.

This work is the first vertebrate diet composition analysis on Neotropical canids conducted using a metabarcoding approach. A sensitive and accurate methodology was applied for crab-eating fox and pampas fox diet characterization to achieve this goal. Our results agree with previous trophic analyses in these fox species based on macroscopic methods, validating the methodological strategy. Moreover, a first approach for detecting consumed species sequence variants was presented based on network reconstruction using non-coding DNA fragments. The extensive use of these techniques will allow going further in the understanding of South American wild canids feeding ecology and in the knowledge about their role in ecosystem regulation. Further analyses are needed to validate the approach for estimating intraspecific diversity through sequence analysis of prey DNA in generalist predators’ feces. In this regard, it is necessary to compare these results with more traditional population genetic methods.