Comparing minimum number of individuals and abundance from non-invasive DNA sampling and camera trapping in the red fox (Vulpes vulpes)

Applying the most appropriate sampling method is essential for estimating population size. Sampling methods and techniques to estimate abundance may be limited by environmental characteristics, species traits, specific requirements of the techniques, or the economic resources to carry out the sampling. Thus, evaluating multiple sampling methods in monitoring populations is essential for establishing effective conservation strategies. In this study, we compare two of the most commonly used sampling methods with the red fox (Vulpes vulpes) as the type species. On the one hand, we compared the minimum number of individuals (NI) detected by camera trapping, identifying individuals by morphological characteristics with the minimum number of individuals detected by DNA faeces and a set of 16 microsatellites. On the other hand, we estimated abundance by performing an N-mixture model using information from camera-traps to study the relationship between abundance and the minimum number of individuals detected. Results showed that the minimum NI provided by camera trapping was slightly higher than that of DNA faecal genotyping, with 23.66 and 19 individuals, respectively. In addition, abundance and NI detected by camera trapping showed a positive relationship. In contrast, there was a non-significant negative relationship between NI detected by faecal DNA and abundance estimates. Our results suggest using the minimum number of photo-identified individuals as a reliable index to study variation in red fox abundance when other advanced methods cannot be implemented in the study of population size. However, it is necessary to improve the methods of faecal sampling to study the relationship with camera-trap data.


Introduction
The recent acceleration of biodiversity loss urgently requires the development of monitoring programmes to understand trends and spatial patterns of the abundance of wildlife species. Among these, the order Carnivora is among the most threatened groups worldwide (Gittleman et al. 2001). Carnivores are emblematic species with an important ecological role; however, they also cause conservation conflicts with human and their activities (Linnell and Strand 2000;Prugh et al. 2009;Ritchie et al. 2012). Many carnivore species are elusive and solitary mesopredators (hereafter mesocarnivores), with nocturnal habits, large territories and low population densities. These characteristics make non-invasive survey methods highly suitable for studying and monitoring these species. Camera trapping and DNA-genotyping are currently the most relevant non-invasive sampling methods applied to monitor mesocarnivores and analyse their abundance assisted by different analytical techniques (Srbek-Araujo and Chiarello 2005;Trolle et al. 2007;O'Connell et al. 2011;Sollmann et al. 2013;Rodgers et al. 2014). However, these methods are unreliable in all scenarios because environmental characteristics and species traits may limit them, as well as the specific requirements of the techniques to estimate species abundance or the economic resources available to carry out the sampling. Therefore, because different sampling methods and techniques can be applied to estimate mesocarnivore abundance in different regions or seasons, it is necessary to assess the relationship between these methods and how they can be selected in terms of their sampling effort, economic costs, precision and accuracy of abundance estimates (Silveira et al. 2003;Gaidet-Drapier et al. 2006;Gompper et al. 2006;Balme et al. 2009).
Abundance can be assessed using population size estimates or indices of relative abundance. Among the former, capture-recapture methods are widely used by biologists to estimate parameters of wildlife population (Pollock 1976) using the capture histories of identified individuals to draw the detection probability. Camera trapping and non-invasive DNA sampling are among the most widespread techniques used to identify individuals (Karanth and Nichols 1998;Jackson et al. 2006;Mondol et al. 2009;Galaverni et al. 2012). However, other techniques such as N-mixture methods can be used when individuals cannot be individually identified. N-mixture models use data from spatially replicated count surveys to effectively estimate population sizes while accounting for the detection process (Royle 2004;Ficetola et al. 2018;Kidwai et al. 2019;Costa et al. 2020). However, these methods imply a great effort and resources as they require a high number of spatial replicates to account for imperfect detection (instance >20 sites; Kéry and Schaub 2012) and a high probability of detection (Royle 2004;Veech et al. 2016) to obtain reliable abundance estimates. Under these limitations for estimating abundance, relative abundance indices are helpful and convenient for studying population size. Indices of relative abundance are positively correlated with the population size of the species (Caughley and Sinclair 1994), so changes in index values reflect a change in actual abundance (Romesburg 1981;Anderson 2003). Among relative abundance indices faeces counts along transects and capture rates by camera-traps are the most important. However, using relative abundance indices should be taken with caution because the monotonic relationship with true abundance may not be constant and positive over time and across habitats. In addition, different indices may show different sensitivities to the same habitat characteristics; therefore, relative abundance estimates will depend on the index used. In addition, the indices do not control the prob-ability of detection of a species as do the technique mentioned above (CR and N-mixture methods), which may influence abundance estimates. Moreover, we should consider that the efficient use of camera trapping and DNA faeces sampling might be limited under some circumstances that could affect abundance estimation. For example, environments with high humidity and temperature poorly preserve DNA in faeces (Murphy et al. 2007). Temperature differences can also lead to reduced detection of animals and undesirable detection biases in camera-traps triggered by 'heat-in-motion' (Meek et al. 2015).
Identifying relative abundance indices that follow a constant relationship with abundance without depending on the sampling method used, habitat characteristics or probability of detection is of great interest when other advanced methods cannot be implemented in the study of population size. Here, we compare the relationship between camera trapping and DNA faecal sampling to study population size during the same spatio-temporal settings. We compared the minimum number of individuals (NI) identified by each sampling method in several areas covering different habitat characteristics. We propose the minimum NI as a relative index of abundance and study the relationship of this index with abundance estimation. In particular, we used N-mixture models (Royle 2004) to estimate abundance using cameratraps. Previous studies focused on carnivore species with spotted fur or animals marked to facilitate the identification of individuals using camera-traps (Karanth and Nichols 1998;Silveira et al. 2003;Soisalo and Cavalcanti 2006;Mosquera et al. 2016). We used one of the most abundant mesocarnivore species in Europe and particularly in the Iberian Peninsula (Lloyd 1980), the red fox (Vulpes vulpes). Although the red fox lacks a distinct fur pattern, all individuals show distinctive characteristics (i.e. body and tail coat; muzzle, head and ears shape; leg and paw marking; injuries) that allow individual for identification (Sarmento et al. 2009;Dorning and Harris 2019). Furthermore, territory marking by foxes' character facilitates faeces sampling (Goszczyński 1990). These characteristics make the fox a model species to test non-invasive sampling techniques in carnivores. The abundance of red foxes is not well-documented in all areas of the Iberian Peninsula. However, some studies have tested different non-invasive sampling techniques to determine its abundance Jiménez et al. 2019a;Jiménez et al. 2019b). The red fox has an critical ecological role in seed dispersal (Campos and Ojeda 1997;Juan et al. 2006), the transmission of rabies disease (Chautan et al. 2000), and its impact on game species (Beja et al. 2009). In Spain, the red fox is a game species with annual hunting quotas. However, the quota implementation follows non-scientific criteria based on abundance estimations. Also, the lack of selective techniques for fox control may harm other endangered carnivores without an effective control of fox abundance (Virgós and Travaini 2005), which highlights the importance of the species for the carnivore guild and associated conservation conflicts.
Using red fox as model species, we tested the following questions: (1) are both methods correlated in detecting the same minimum NI? (2) is the minimum NI index correlated with red fox abundance estimates by both methods? (3) which method is the most economical for estimating red fox abundance in the same spatial and temporal environment? Our applied aim is to provide methodological information to assist decision-making in research requiring estimates of the relative abundance of mesocarnivores.

Study area
The red fox is continuously distributed in the Autonomous Community of Madrid (Spain). Thus, we selected seven areas in this region for our study: Carabaña (C) and Villarejo de Salvanés (V) in the southeast, La Berzosa (B) and San Mames (SM) in the north and Robledo (R), Quijorna (Q) and Pelayos de la Presa (P) in the southwest (Fig. 1).
The southeast part of the mesomediterranean floor has a temperature range of 5.7-7.5 °C in winter and 20.7-25.2ºC in summer and annual precipitation of around 390 mm per year. Here, rainfed crops, irrigated, and pasture predominate. Oak (Quercus ilex rothundifolia and Quercus coccifera) forests on calcareous soils result in a fragmented landscape characterized by a mixture of open spaces, low vegetation (https://www.comunidad.madrid/ servicios/urbanismo-medio-ambiente/parque-regional-sureste), also with a strong human influence (Pascual et al. 2010).
The northern part of our study area lay between the meso-and supramediterranean floors characterized by milder summers and colder winters than the southeastern area; precipitations are also more abundant in this area. The landscape consists of Sclerophyllous vegetation (vegetation with hard leaves, short internodes and leaf orientation parallel or oblique to direct sunlight), oaks with gum rockrose (Cistus spp) and abundant granite boulders. Like- Fig. 1 Sampling areas, cameras, and trail distribution in the study areas. Carabaña (C), Villarejo de Salvanés (V), La Berzosa (B), San Mames (SM), Robledo (R), Quijorna (Q) and Pelayos de la presa (P). The dots show the camera locations, and the lines show the faecal trails (red is trail 1; yellow is trail 2, and green is trail 3) wise, scrublands of broom (Cistus spp) and thyme (Thymus vulgaris) with pine woods are present from an altitude of 1400 m a.s.l. The southwest region of Madrid is characterized by a great variety of landscapes, reliefs and plant species, with predominant coniferous forests, especially stone pine (Pinus pinea). On the slopes, there are holm oaks (Quercus ilex), cork oaks (Quercus suber), other oaks (Quercus spp), and junipers (Juniperus communis), which are typical species of Mediterranean forests. Rainfall is abundant throughout the year except in summer, and winters are cold with occasional frosts and snowfalls.
The sampling period took place in 2017 (B, C, SM and R) and 2018 (V, Q and P) to cover different seasons: B, C and P in summer, V and Q in spring, SM in autumn, and R in winter. In this way, differences in detectability and marking that depend on seasonality were avoided (Ralls et al. 2010).

Camera trapping sampling and faeces collection sampling
At each site, we placed eight to ten cameras (Scoutguard SG562C LED White model), spaced approximately 450-600 m apart, covering as large an area as possible to maximize the number of individuals photographed and reduce the likelihood of unsampled foxes (similar to Sarmento et al. 2009). We used ArcGIS 10.2 (ESRI Inc., Redlands, California, USA) to create a minimum convex polygon (MCP) using the camera locations. We placed sardines and a commercial lure (HAGOPUR® Premium Attractant Fox) at the camera-trap sites to increase the probability of detection of foxes (Heinlein et al. 2020;Sebastián-González et al. 2020). The cameras operated for 35 days and were checked every seven days to replenish baits, collect the photographs and check the battery (see Martin-Garcia et al. 2022 for additional information on the photographic sampling design).
We collected faeces in 21 trails of 1 km distributed in seven areas (i.e. three trails per area) (Supplementary Table S1). We cleaned each trail on the first day of camera placement to ensure an exclusive collection of faeces deposited within the sampling period. After that, each trail was sampled every two weeks three times during camera placement to increase the probability of detecting faeces from all individuals in the population. Trails were also inspected the day after the cameras were checked. We chose trails based on proximity to the camera placement areas to increase the probability of detecting the same individuals by both methods (Fig. 1). Fresh faeces that could potentially belong to foxes due to morphology and odour were collected by the same operators. All the operators had broad experience in recognizing carnivore faeces. We placed faeces in 96% alcohol for the first 12 h and then stored them on silica gels at 2-4ºC until they were processed (Nsubuga et al. 2004).

DNA extraction and amplification
We extracted DNA from the faeces using a QIAamp DNA Stool Kit (Qiagen) following the manufacturer´s protocol. DNA purity and quantification were determined with a Nano-Drop® 2000 spectrophotometer and Qubit® 3.0 fluorometer Quantitation Kit (Invitro-gen™). For species verification, we amplified a short region (120 bp) of the mitochondrial (mt) DNA "ND1" gene. Samples were BLAST using the BLASTx option at the NCBI platform to identify foxes. PCR reactions were performed in a total volume of 20 µL (the protocol for species verification is explained in detail in Supplementary material S1). All amplifications were carried out using filter tips in separate rooms (pre-and post-PCR), and negative controls/blanks were included in all amplifications to avoid contamination. PCR products were run and visualized on a 1.5% agarose gel using gelgreen (BIOTIUM). PCR products were sequenced at Macrogen (Netherlands).

Photo-identification and microsatellite genotyping
Individuals were identified based on traits such as body size, age range and, the appearance of the tail, spotting at specific points and other diagnostic features. We followed Sarmento et al. (2009) and Dorning and Harri (2019) to select traits to assist identification ( Table 1; see also Fig. 2 for examples of individual identification). We did not account for seasonal changes in coat for each individual. The sampling time was not long enough to appreciate these seasonal differences, and we did not repeat sampling in the same area during different seasons.
Cubs were not individually identifiable due to their juvenile fur lacking sufficient distinctive features. However, juveniles tend to move together, so we decided to consider the maximum number of juveniles appearing together in the same photo as the minimum number of juveniles in the population and include it in the study.
We reviewed photos by a second observer to control for photo-identification bias and reduce overestimation (Foster and Harmsen 2012;Ferreras et al. 2017;Johansson et al. 2020). A third observer re-analyzed fox photographs when the first and second observers disagreed on the number of foxes identified. We considered the minimum number of individuals as the mean of the number of individuals identified by each observer. We also estimated the standard error and confidence intervals of the number of individuals photo-identified.
In microsatellite genotyping, allele sizes were manually scored relative to an internal size standard (GeneScan GS600LIZ) in GeneMapper version 5 (Applied Biosystems ™). Each locus was amplified three times to minimise genotyping errors. Genotypes were accepted as reliable if: (1) a heterozygote was observed at least twice in two independent reactions, and/ Source of variation Body coat Condition and colour of coat on body, belly and chest (black or white patches or stripes). Head/face/ muzzle shape Size of the head; length of the muzzle; the fullness of cheeks; marks on top of nose; colour and shape of fur markings on each side of the muzzle. Ears Length; shape (rounded or pointed); colour and texture of the coat on the inside of the ears; mottling on the back of the ears; wounds on the rims of ears Leg and paw markings Height and shape of black socks on the legs; black or white coat on the front of the thighs; coat colour on the inside of the legs; mottling of the coat and presence of white markings. Tail coat Coat condition; colour and pattern; shape and size of the dark patch around the supracaudal gland.

Tail shape
Length; thickness; straightness; tip shape (pointed, rounded, tapered, tufted, curly or flattened) Injuries Sarcoptic mange infection; bites; scars and deformities. Table 1 Features selected for fox identification or (2) a homozygote was observed at least twice independently and failed amplification of the third replicate (Taberlet et al. 1996;Frantz et al. 2003;Flagstad et al. 2004). We defined ambiguous genotypes as (1) genotypes with only one amplified replicate, (2) genotypes that failed all replicates, (3) genotypes that generated different alleles in each replicate, or (4) genotypes with two identical homozygous replicates and a different third. These genotypes were annotated as missing alleles (000000) for the corresponding marker. We used FreeNA program to test the null allele frequency (Chapuis and Estoup 2009). Gimlet 3.4 (Valière 2003) was used to evaluate heterogeneity observed (Ho) and expected (He); Hardy-Weinberg equilibrium (H-W equilibrium) and probability of identity (P ID ) to evaluate the reliability of the microsatellites used. The P ID among genotypes (Kalinowski et al. 2007) is the most widely used statistical method to quantify the power or ability of molecular markers to distinguish between two individuals. We tested P ID , P ID (sib) , multilocus P ID and Multi-locus P ID (sib) . The P ID is the probability that a single unrelated individual has this genotype (individuals are randomly mated); P ID (sib) is the probability that a single full sibling has this genotype (sister-sibling only population), and the multi-locus P ID is the P ID calculated over several loci by sequentially multiplying the P ID value over the loci (considering that the loci are independent) (Scandura et al. 2001). The multi-locus P ID (sib) is the P ID (sib) calculated over several loci. We considered that a multi-locus P ID (sib) less than 0.01 (Mills et al. 2000) or a multi-locus P ID between 0.001 and 0.0001 was sufficiently sensitive for identification and avoided underestimation (Waits et al. 2001). The presence of null alleles was checked with the FreeNa program. We also tested the effect of removing markers with a null allele frequency above 0.30 (Dakin and Avise 2004;Huang et al. 2016).
We used Cervus 3.0.7 (Kalinowski et al. 2007) and Gimlet 3.4 to create genotype profiles for all the samples and test the consensus in the number of individuals provided by both programs. Cervus 3.0.7 identifies samples with identical genotypes for the specified number of loci. From our set of 16 microsatellites, we scored individuals as the same if they had identical genotypes for at least eight or more common loci. This grouping method with matching samples ensures a conservative number of identified individuals by minimizing individuals created through erroneous, multi-locus genotypes (Mondol et al. 2009). Also, we allowed a maximum of two mismatching loci to control genotyping errors and increase success in genotype assignment (Kalinowski et al. 2007). We used Gimlet 3.4 to reconstruct the consensus genotypes. Finally, we tested the genotype reassortment function with the assumption that missing alleles are considered distinctive alleles.

Model selection: abundance estimation
We ran N-mixture models to estimate fox abundance from camera trapping without individual identification and accounting for the influence of imperfect detection. We used unmarked package (Fiske and Chandler 2011) in R software (R Core Team 2022) using the pcount function. This function estimates abundance in a hierarchical model. The actual abundance is estimated from the local variation of abundance (λ) at i sites using j temporal counts controlling for detection probability. There are two linked processes for estimating abundance (Kéry and Schaub 2012): a) Abundance process (λ): Fitted by a Poisson distribution with mean λ and the variation of local abundance at site i λ ∼ P oisson (λ) b) Observation process ( p ): Fitted by a binomial distribution of the observed counts ( y i,j ) of individuals at each site (i ) in each temporal replicate ( j ) with a probability of detection.
We counted each independent fox trapping event in each camera-trap per occasion (24 h). We considered more than one fox capture per occasion when we detected several foxes together in the same independent fox capture event. We first ranked three models to study abundance and detection probability among sampling areas: (1) constant detection probability among areas with a variation of abundance estimates among areas (i.e. p (.) ~ λ (site)), (2) variation in detection probability and variation of abundance estimates among areas (i.e. p (site) ~ λ (site)), and (3) variation in detection probability among areas with constant abundance estimates among areas (i.e. p (site) ~ λ (.)). We compared the performance of the Poisson, zero-inflated Poisson, and negative binomial distributions for each model. For model selection, we used Akaike's Information Criterion (AIC) (Akaike 1974) corrected for small sample size (AICc) (Burnham and Anderson 2002). We run a chi-square test in Nmix.gof.test function of package AICcmodavg (Mazerolle and Mazerolle 2017) to assess the goodness-of-fit and overdispersion of the selected model. We then estimated the posterior distribution of detection and abundance ( λ ) using empirical Bayes random effects (ranef) methods from the unmarked package. We used a parametric bootstrap approach with a simulation of 5000 bootstrap samples for each fit assessment. We obtained the mean abundance of each area, the standard error and the confidence interval.

Pearson correlation
We used a Pearson's correlation test to quantify the proportionality and similarity between (a) minimum photo-identified and genotyped NI, (b) minimum photo-identified NI and abundance estimates, (c) minimum genotyped NI and abundance estimates. We performed a logarithmic transformation when data were not normally distributed. We followed Prion's classification to establish correlation ranges (Prion and Haerling 2014).

Individual identification: photo-identification and microsatellite genotyping
We obtained 309 photos of foxes. We identified a total of 23.66 individuals with a mean of 3.38 individuals per area (Table 2). Indentified foxes per each observer were described in Supplementary Table S5.
We collected 77 faeces along the trails (Table 2). According to the molecular species identification, mtDNA sequence ("ND1" gene) identified 69 samples of Vulpes vulpes. Based on GIMLET analyses, six samples (three from C and three from Q) could not be assigned to a specific genotype and were eliminated from the analyses, resulting in 63 faecal genotypes.
Consensus genotypes of sample types with 16 microsatellite loci (N = 63) revealed the presence of 19 individuals ( Table 2). All loci included in the study were polymorphic. We found 3-12 alleles per locus and an average of 7.06 annotated alleles per locus. The mean allelic dropout rate was 24% (loci = 16; SD = 0.149) among loci and 35% (samples = 69; SD = 0.235) among samples. The single locus probabilities were combined to obtain the total probability over the 16 loci, assuming the independence of different loci. Results of Ho, He and H-W equilibrium of 16 microsatellites are shown in Supplementary Table S6. The different identity probabilities for the 16 microsatellites loci were multi-locus P ID = 1.98 × 10 -14 and multi-locus P ID (sib) = 3.43 × 10 − 06 . Based on the identified probability, the 16 loci considered were sufficient to distinguish with 99% certainty between sibling red foxes (Supplementary Table S6 ). We obtained null allele values of less than 0.30 in almost all loci, except for INU030, which had a null allele frequency of 0.33 (Supplementary Table  S6). We repeated the analyses after removing INU030 and obtained the same results regarding the same number of individuals and sample genotypes. However, the multi-locus P ID increased slightly to 1.16 × 10 -13 . Therefore, we decided to retain this microsatellite in subsequent analyses (Huang et al. 2016). The number of individuals identified was consistent between the two software tools used (Cervus 3.0 and Gimlet 3.4), with a single exception in the P location population, where we detected two individuals by Gimlet 3.4 but could not obtain results with Cervus 3.4. The number of assigned faeces to each individual is described in Supplementary Table S5.

Model selection: abundance estimation
Model selection resulted in two top-ranked models (lowest AICc). The model with lower AICc indicated a constant detection probability and variation in abundance estimation between areas. The second model maintained a constant abundance estimate and a variable detection probability between areas (Table 3).
In both models, the Poisson distribution was more supported than the zero-inflated Poisson and negative binomial. Although the two top-ranked models were close (i.e. ∆AICc < 2), we decided to select the first top-ranked model because it explained our predictions for

Discussion
Evaluating multiple sampling methods in monitoring carnivore populations is essential for establishing effective conservation strategies (Caughley and Sinclair 1994;Sadlier et al. 2004;Barea-Azcón et al. 2007). We used camera trapping and DNA faecal genotyping to estimate the minimum number of foxes. We estimated fox abundance by implementing N-mixture models to assess and compare the relationship between abundance and minimum NI provided by camera trapping and DNA faecal sampling. First, we found that the estimation of minimum NI provided by camera trapping was slightly higher than that of DNA faecal genotyping. Second, there are indications that areas with more NI identified were those with higher fox abundance, following a positive relationship between abundance and NI detected by camera-trapping. However, we also found a non-significant negative relationship between NI detected by faecal DNA and abundance estimates. The comparison of camera trapping and DNA faecal sampling methods showed a slight variation in the number of the identified foxes. We identified three more individuals with camera-traps than DNA faecal sampling (minimum NI was identical in Carabaña and Pelayos). However, despite this slight variation, we found no significant correlation between the minimum NI calculated from photo-identified and faecal genotyped. Comparisons between both non-invasive sampling methods to estimate the population abundance of carnivore species are relatively common in the literature (Mondol et al. 2009 Felis silvestris)). Galaverni et al. (2012) found concordance between the two sampling methods in the number of identified wolves, which supported their complementarity. However, the minimum number of wolves identified was slightly higher with the use of DNA than with the cameras. In particular, Mondol et al. (2009) found non-significant differences between sampling methods in tigers, although they detected three more individuals using camera-traps than DNA methods.

Table 4
Comparison between the minimum number of individual photo-identified (NI photo-identified), the minimum number of individual faecal genotyping (NI faecal genotyping) and the abundance. Estimation of fox abundance in each area was conducted using Poisson models. The following statistics are shown: mean abundance (abundance mean); standard error of abundance (abundance S.E); lower and upper interval coefficient of abundance ( We identified fewer individuals using DNA faecal sampling than camera trapping (especially in San Mamés and La Berzosa). This difference could be explained by limitations in the faeces sampling design due to the use of faeces for communication between foxes. Carnivores mark their home ranges and territories with faeces deposited in a non-random distribution (Kruuk 1978;Macdonald 1980;Gorman 1990;Soler et al. 2009). Accounting for this behaviour, we chose random but well-delimited paths (trails and roads) that covered the entire area sampled by cameras. Fox faeces often mark the boundary of territories or sites with critical resources (Macdonald 1985;Barja et al. 2001). Because resource and territory marking occurs unevenly between individuals, marks of specific individuals in a given area may be over-represented (Gorman and Trowbridge 1989). Moreover, defecation rates differ between individuals (Cavallini 1994), particularly between males and females, adults, and juveniles (Goszczyński 1990;Peterson et al. 2002;Ralls et al. 2010;Fawcett et al. 2012), which can influence the number of individuals detected. Defecation rate also vary according to season and diet (Andelt and Andelt 1984;Goszczyński 1990). Other effects, such as changes in marking behaviour caused by the spatial distribution of roads (Vilà et al. 1994;Barja and List 2014;Zaman et al. 2019) may also lead to increased non-uniform distribution of faeces. In our study area of Villarejo, in contrast to other locations, the number of individuals identified by DNA in faeces was higher than using cameras. These differences suggest some limitations of the camera trapping method. Although 32 photos were obtained for identifying individuals, we obtained 15 unidentified photos of either new or previously detected individuals. Camera traps could alter red fox behaviour, specially in females due to the sounds and flashes they make (Meek et al. 2014). Predation risk or anthropogenic disturbance (Lucherini et al. 1995) might also increase trap avoidance, thus influencing individual detection. These potential behavioural changes might impact the more elusive individuals more effusively, with the consequent bias in their detection.
Although we found a non-significant correlation between the genotyped and photo-identified minimum NI, the minimum photo-identified NI and the abundance estimation showed a positive relationship. This result suggests the possibility of using NI as a straightforward index to explain variations in species abundance. We validated individual red fox photoidentification according to previous studies (Sarmento et al. 2009;Dorning and Harris 2019 but see, Güthlin et al. 2014) to obtain a minimum NI. However, we recommend at least three observers to reduce the bias in over-or under-estimating the number of individuals (Foster and Harmsen 2012;Ferreras et al. 2017;Johansson et al. 2020). Photo identification also has limitations for juveniles or cubs. Identifying juveniles becomes difficult when individuals are alone. Therefore, we consider the minimum number of juveniles as the total number of juveniles together in the same photo. However, we acknowledge this approach might be susceptible to underestimating juvenile populations because we considered only the minimum number of juveniles. Nevertheless, we assume that adding juveniles to the study generated less bias when comparing the minimum NI identified and the abundance estimated by both sampling methods. The identification of individuals has become a standard method for the abundance estimation of animal populations, such as using the Spatial Capture-Recapture (SCR) method (Efford 2014;Royle et al. 2014;Rodgers et al. 2014;Wegge et al. 2019). However, CR methods for abundance estimates require a sufficient number of recaptures to obtain accurate estimates (Otis et al. 1978). SCR approaches need at least 20-25 recaptures, including spatial recaptures to correctly describe the movement , or alternatively using movement data from telemetry tagged individuals (Jimenez et al. 2019) to improve the estimate. Due to the data limitations, we used N-mixture models to estimate fox abundance.
The N-mixture model estimates abundance using count data without needing individual identification or reference to the effective trapping area. Other studies focused on the reliability of the N-mixture models to estimate abundance. Basile et al. (2016) found that N-mixture models and SCR methods yielded a similar estimation of the abundance of the short-toed treecreeper (Certhia brachydactyla). Ficetola et al. (2018) obtained limited differences between N-mixture models and capture-mark-recapture to estimate the abundance of small vertebrates. Also, to avoid violating assumptions of N-mixture caused by double counting of a single sampling occasion (Link et al. 2018), we only considered more than one fox capture per occasion (24 h) when we detected several foxes together in the same capture event. Our results between photo-identified NI and abundance showed that NI results are helpful when the minimum number of identified individuals in the population is not sufficient, or the number of captures of identified individuals is not enough to produce reliable abundance estimates (Otis et al. 1978). NI may also be recommended when sufficient temporal and spatial replicates are unavailable or the assumptions of N-mixture models are not met (e.g. independence in the case of gregarious animals). Thus, Martin-Garcia et al. (2022) found that the minimum photo-identified NI might not be biased by detection probability, thus obtaining the same predictors as the N-Mixture models to explain abundance patterns.
We found no significant correlation between minimum NI genotyped and abundance. Including more faeces could help identify a potential relationship between the minimum number of individuals identified by DNA faeces sampling, NI photo-identified, and abundance (Wegge et al. 2019;Lindsø et al. 2022). Future research should include random transects out of existing trails to increase the number of faeces and detections of individuals. Implementing a random transect design that covers many landscapes with different compositions and configurations (Güthlin et al. 2012) could reduce the bias caused by some individuals marking more intensively along the trails. In this vein, we could better refine and compare the relationship between estimated abundances of faecal DNA and camera trapping sampling methods (Rodgers et al. 2014) and between NI genotypes and abundance. In addition, new state-of-the-art genomic approaches based on SNPs approach (e.g. RAD-seq) may increase the accuracy of DNA amplification decreasing the loss of samples and consequently improving abundance estimates (Andrews et al. 2016;De Barba et al. 2017;Erwin et al. 2021). Another timely method is the SNP genotyping method based on high-throughput real-time PCR technique known as Dynamic Array™ by Fluidigm® that is mainly used for degraded samples such as faeces and ancient DNA studies samples (Kraus et al. 2015).
In our experience, using camera trapping was cheaper than DNA sampling methods to study red foxes. Our budget for the camera trapping method was 2756 euros (including cameras, baits and placement in the study areas) compared to 7500 euros for the faecal DNA method (including fieldwork, sample shipment and protocol optimization). Several aspects should be considered in terms of budget. Firstly, a small pilot optimization study beforehand could help to reduce costs in future analyses. The optimization process with fewer samples is helpful to check the protocol used and obtain preliminary results avoiding using all the valuable samples. Secondly, using other less costly genetic techniques. Single nucleotide polymorphism (SNPs) is frequently used as a new genotyping method for individual and sex identification (Parker et al. 2021;Buchalski et al. 2022;Lopez-Bao et al. 2020). To date, developing a SNP panel for individual identification is considered an efficient and cost-effective method to simultaneously genotype hundreds of individuals (Carroll et al. 2018). Lastly, researchers should consider an extra budget for the replacement of the cameras in case of loss or failure during the study. Consequently, the budget for camera trapping may increase depending on the number of cameras needed for our research.
Our research highlights the importance of correctly selecting sampling methods for abundance studies. Researchers can adjust the broad choices in sampling methods to their available funds and logistics. Different methods can perform differently and provide different results; thus, it is required to identify first the costs and limitations of the potential techniques for our specific research objectives and the species to study. Our results suggest the minimum photo-identified NI is a reliable index for studying abundance variation when other methods are unavailable. In contrast, it is necessary to improve the methods of faeces sampling to estimate population size and to explore its relationship with camera trap data. Sampling designs with transects away from existing trails will increase the probability of finding more faeces. We should mention that both methods were compared to study abundance over a short period. On the other hand, DNA sampling in faeces could be useful to identify individuals over a longer period, over years and seasons when this would be very difficult with photos (Bellemain et al. 2005). In addition, assessments of genetic diversity, population substructure, gene flow, paternity, and heritability are easy to evaluate with DNA stool genotyping (Pilot et al. 2014;Zanin et al. 2016). However, this was not the main objective of this study. Regarding our model species, the red fox, camera trapping methods can be a reference for future red fox management actions. In our study, the validity of using the camera trapping method to estimate fox population abundance is also motivated by its lower cost when compared to the faecal DNA genotyped method. However, further research on the cost-effectiveness of new genetic methods is encouraged.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.