DNA metabarcoding reveals wolf dietary patterns in the northern Alps and Jura mountains

Understanding predator-prey interactions is crucial for wildlife management and human-wildlife coexistence, particularly in multi-use landscapes such as Western Europe. As wolves (Canis lupus) recolonize their former habitats, knowledge of their diet is essential for conservation planning and public acceptance. However, data from such regions is so far scarce, in particular for the Alpine region and the Jura Mountains.  As opportunistic hunters, wolves adapt their diet to the local prey structure, suggesting that abundant wild ungulates are the main prey source. We also expect diet composition to be influenced by region, season, and social status. DNA metabarcoding has emerged as a powerful tool in ecological research, offering high-resolution insights into dietary composition, yet its application in carnivores remains limited. Using a DNA metabarcoding approach, we analyzed 698 wolf scat samples collected in Switzerland from 2017 to 2024. We found red deer (Cervus elaphus) was the dominant prey in most areas, and together with roe deer (Capreolus capreolus) and chamois (Rupicapra rupicapra), accounted for 80.8% of the retrieved sequences. We found similarities in prey selectivity between the Alps and the Jura Mountains, but found significant differences across seasons and between pack and non-pack wolves. This study provides the first detailed dietary analysis of wolves in the Northern Alps and Jura Mountains, offering critical insights for wildlife management. Our findings highlight the importance of wild ungulates in wolf diet and underscore the value of non-invasive DNA approaches as a reliable conservation and biomonitoring tool.

Methods
Scat collection and sampling context

Scat samples were collected opportunistically across Switzerland as part of the national wolf monitoring by trained authorities (Marucco et al. 2023), and had already been genetically assigned to specific wolf individuals by STR (microsatellite) genotyping (Dufresnes et al. 2019). Our analysis is based on 698 samples collected between 02.01.2017 and 07.04.2024 from different regions and seasons. These samples belong to 250 genetically identified wolves (males = 155, females = 93, unknown = 2, pack = 181, non-pack = 69). As the scat samples were not collected systematically or randomly, we have taken into account the context in which they were found (i.e., at a livestock kill, a carcass of a wild prey, or independently of a carcass). We compared the proportions of livestock and natural prey species between the three contexts with Pearson’s Chi-squared test (Pearson 1900). We also tested whether the inclusion of samples found next to carcasses changed the overall result significantly.

DNA extraction

Scat samples were stored in alcohol and kept at 4 °C, until they could be extracted in a purpose-built low-DNA content laboratory at the University of Lausanne. For the DNA extraction, we used the QIAamp DNA Stool kit (Qiagen, Hilden, Germany), which is specifically optimized for scat samples. We followed the manufacturer protocol from step 6, i.e., the protocol was started directly with the addition of InhibitEX to 200mg of wet scat sample, skipping the addition of buffer ASL. Extracted samples were diluted 5-fold before PCR amplification.

Generalist primer and blocking primer

DNA extracts were amplified using a generalist primer pair highly specific for vertebrates (Vert01,(Taberlet et al. 2018)). The barcode targets the 12S mitochondrial rDNA region and amplifies within a range of 56 to 132bp. To limit the amplification of wolf sequences, we added a blocking primer for Canis lupus (5’-CTATGCTTAGCCCTAAACATAGATAATTTTACAACAAAATAATT-C3-3’) to the PCR reaction mix (Jarman et al. 2002). The first 6 bases of the 5’-end of the blocking primer overlap the last 6 bases of the 3’-end of the Vert01_forward primer, inhibiting the amplification of C. lupus target DNA. Given that the most abundant vertebrate DNA on a wolf scat is wolf DNA, adding a blocking primer provides a solution to increase the proportion of prey sequences obtained, and reduce the ratio of redundant wolf sequences in the final sequencing output. We tested the performance of the blocking primer through quantitative-PCR (qPCR), by comparing the amplification success of the scat samples with increasing blocking primer final concentrations (2.5μM, 5μM, 10μM et and 0μM).

After examining the qPCR results, a blocking primer final concentration of 2μM was deemed sufficient for the first sequencing run. However, many of the samples from this first sequencing run only produced wolf sequences, so these were resequenced with a higher final concentration of blocking primer (4μM) in the subsequent second and third sequencing runs, which mainly contained other additional samples.

DNA metabarcoding and sequencing

Forward and reverse Vert01 primer pairs were combined with an identifying tag, consisting of eight variable nucleotides with at least five bases difference between tags, to identify the sample source of each sequence. Spacers were also added to the 5’-end, i.e. a few random nucleotides of varying length (1 to 3), to increase complexity for enhanced cluster detection during the sequencing run.

PCR reaction mix was as follows: AmpliTaq Gold 360 Master Mix 1x, tagged forward and reverse primers at 0.2μM, BSA at 0.16 mg/ml and Canis lupus Blocking primer at 2-4μM. The final volume per well was 20μL, including 2μL of DNA template. The PCR thermal profile started with denaturation at 95°C for 10 minutes, followed by 40 cycles of amplification. Each cycle was composed of 30 seconds at 95°C, 30 seconds at 49°C and one minute at 72°C, before a final elongation at 72°C for seven minutes. Each scat was amplified in triplicate. For each PCR plate, we added positive controls (an equimolar mix of five different species: Pelophylax ibericus, Rana italica, Timema sp., Cyanoliseus patagonus and Myotis capaccini); extraction controls, to assess the contamination that might have occurred during the DNA extraction part; negative controls, i.e. water, to assess potential contamination during the PCR process; and blank controls, i.e. empty wells, which were used to assess tag-jumps and sequencing errors during the bioinformatics part.

Amplification success was verified on a subset of samples through electrophoresis on a 2% agarose gel. PCR products were then pooled together by plate into a single tube, but amplified positive controls were removed beforehand. After PCR amplification and pooling, amplicons were purified using the MiniElute PCR purification kit (Qiagen, Hilden, Germany). Purified pools were quantified using a Qubit® 2.0 Fluorometer (Life Technology Corporation, USA) and sent to the Fragment analyzer (Advanced Analytical Technologies, USA) to quantify the length and abundance of the amplicons.

Library preparation for sequencing was done following the Tagsteady protocol (Carøe and Bohmann 2020). After adaptors were added, libraries were validated with a fragment analyzer (Advanced Analytical Technologies, USA), to verify the adaptors had correctly attached to the amplicons. Final libraries were quantified through qPCR, normalised and pooled before 150 paired-end sequencing on an Illumina MiniSeq sequencing platform (Illumina, San Diego, CA, USA). Because of the sequential arrival of sequences, sequencing was done in three separate runs, all using a Mid output Kit (Illumina, San Diego, CA, USA).

Bioinformatics

The bioinformatic processing of the raw sequence output and first filtering was done using the OBITools *package (Boyer 2016). In brief, forward and reverse sequences were aligned with a quality score of at least 40. Joined sequences were assigned to samples based on the unique tag combinations attached to the primer pairs. Assigned sequences were then de-replicated, retaining only unique sequences, to greatly reduce the size of the files and the computing time. Afterwards, sequences with less than 10 reads per library were discarded as well as those not fitting the range of the Vert01 barcode lengths. We then calculated pairwise dissimilarities between reads using the *obiclean function. Lesser abundant sequences with single nucleotide dissimilarity were clustered into the most abundant ones. Then, we used the Sumaclust algorithm (Mercier 2013) to further refine the resulting clusters based on a sequence similarity of 97 %. It relies on the same clustering algorithm as UCLUST (Prasad, Madhusudanan, and Jaganathan 2015) and it is used to identify erroneous sequences produced during amplification and sequencing, which are derived from its main (centroid) sequence.

The retained sequences were taxonomically assigned to taxa with a database for Vert01 (Supplementary Information), which was generated using the EMBL database (European Molecular Biology Laboratory). To ensure the accuracy of our automatic taxonomic assignment, these operational taxonomic units (OTUs) were then double-checked using BLASTn on the NCBI database. 

Data cleaning process continued in R (version 4.0.2). We used the *metabaR *package (Zinger et al. 2021) to refine the sequencing output. The package uses the sample type from each plate, i.e., sample, extraction, negative, blank or positive control, to flag sequences based on their potential to be contaminants. For example, sequences which are more abundant in extraction and negative controls than in samples are labelled as contaminant sequences and removed. The package also compares replicates between them, and discards the ones that are too different from the others, using all the differences between replicates from each sample as reference. We further removed replicates with less than 100 reads in total, as we considered them as low coverage. This process was done separately for each library and allowed us to produce clean datasets, which were directly used for the statistical analysis.

Remaining PCR replicates for each individual scat were then merged. The reads within each PCR replicate were pooled to calculate the presence/absence (PA) of each OTU, the proportional count of OTUs, i.e., Frequency of occurrence (FOO) and proportion of sequences, i.e., Relative read abundance (RRA). Overall, after bioinformatic data cleaning, we obtained vertebrate sequences for all the samples, i.e., 698. Specifically, 665 samples produced at least one non-wolf sequence, while 33 produced only wolf sequences, despite the addition of the blocking primer.

Prey occurrence in dietary data & statistical analysis

We linked the results from the metabarcoding to the seasonal period (spring – early summer = April - July, late summer - autumn = August - November, winter = December - March), geographical region and social status predictors. The attribution of scats to one of the two social statuses was carried out by overlapping the spatial location of the scat with polygons of pack territories present in a given year. Scats collected outside polygons were assigned to non-pack wolves. Non-pack wolves therefore include dispersers, resident single wolves or pairs. Pairs were regarded as packs if they reproduced in the same year.

We ran our statistical analysis regarding the influence of the predictors on diet composition on PA, FOO and RRA of the five most frequent OTUs, later called prey or consumed species. However, we based our results on RRA, as it incorporates detected species in a more quantitative way compared to PA or FOO (P. D. Lamb et al. 2019). Converting read proportion to PA or FOO can have drawbacks, i.e., overestimating the importance of small and rare species because they are counted equally alongside items consumed in greater proportions. RRA may therefore give a more accurate image of the actual consumed proportions (Deagle et al. 2018), although its quantitative accuracy is still an ongoing debate. Each of these metrics was averaged across all samples to provide a mean value per prey item for the general visualization figures (Figure 3, Figure 4, Supporting Information).

We first created a distance matrix with Bray-Curtis dissimilarity with the metaMDS *function of the *vegan *package (Oksanen et al. 2025) and then visualized the dissimilarities with non-metric multidimensional scaling NMDS across the predictors (Buzan et al. 2024; Ricotta and Podani 2017; Wang et al. 2022). We then analyzed dispersion between levels of each predictor using *vegdist *and *betadisper *function with *Bray-Curtis dissimilarity. We statistically tested for differences in dispersion with ANOVA, whereby significance does not presuppose ideal conditions for the further statistical comparison of the groups. \

We assessed differences in the diet composition with a permutational multivariate analysis of variance (PERMANOVA) (Anderson 2001) with *adonis2 *function of *vegan *package, 999 permutations and weights on number of samples per region. Finally, we tested for the influence of specific prey  species on the environmental variables (region, seasonal period, social status) with the similarity percentage method (SIMPER). *Simper *of vegan package compares levels of a categorical variable with a distance matrix as response matrix, in our case the RRA of each detected species per sample. We wrote the results of the diet composition in the language of evidence as a more nuanced approach(Muff et al. 2022).

All data manipulation, visualization, and statistical analysis were performed using R software (R Core Team 2021). Key packages included *tidyverse *for data manipulation and visualization (Wickham et al. 2019), *sf *(Pebesma 2018) and *raster *(Hijmans et al. 2025) for spatial data handling, *tmap *for mapping (Tennekes 2018), and *vegan *was used for ecological and statistical analysis.

For calculation of the consumed biomass we multiplied the RRA of consumed prey species with the minimum biomass required based on the field metabolic rate (Glowacinski and Profus 1997) and an average body weight of 32kg, typical of an Italian wolf (Gazzola, Avanzinelli, et al. 2007). We used RRA, as it offers a more quantitative approach, suitable to the relative consumption of biomass.