A multidisciplinary approach to estimating wolf population size for long‐term conservation

The wolf (Canis lupus) is among the most controversial of wildlife species. Abundance estimates are required to inform public debate and policy decisions, but obtaining them at biologically relevant scales is challenging. We developed a system for comprehensive population estimation across the Italian alpine region (100,000 km2), involving 1513 trained operators representing 160 institutions. This extensive network allowed for coordinated genetic sample collection and landscape‐level spatial capture–recapture analyses that transcended administrative boundaries to produce the first estimates of key parameters for wolf population status assessment. Wolf abundance was estimated at 952 individuals (95% credible interval 816–1120) and 135 reproductive units (i.e., packs) (95% credible interval 112–165). We also estimated that mature individuals accounted for 33–45% of the entire population. The monitoring effort was spatially estimated thereby overcoming an important limitation of citizen science data. This is an important approach for promoting wolf–human coexistence based on wolf abundance monitoring and an endorsement of large‐scale harmonized conservation practices.


INTRODUCTION
Quantifying abundance is fundamental for wildlife conservation and management.Reliable information on wildlife population size is key to management and harvest regulation decisions, defining conservation status, and furthering understanding of ecosystems functioning (Williams et al., 2002).However, population abundance continues to be challenging to estimate in nature (Waples & Feutry, 2022), especially when species are elusive, vagile, and distributed over large areas (Blanc et al., 2012).This is because population sampling is never complete, making it necessary to use analytical methods that account for imperfect detection when estimating population size (Kéry et al., 2009).Important advances have been made in understanding imperfect detection in part due to the wealth of capture-mark-recapture data generated by noninvasive genetic sampling (NGS), whereby individual genotypes are identified from DNA extracted from samples (e.g., scat, hair) left behind in the environment (Schwartz et al., 2007).Large carnivore monitoring programs have benefited from this approach by combining noninvasive genotyping with estimations of population size from capture-recapture models (Caniglia et al., 2012;Cubaynes et al., 2010;Marucco et al., 2009).More recently, spatial capture-recapture (SCR) approaches have been used to estimate densities of elusive animals based on DNA sampling (Kéry et al., 2011;Proffitt et al., 2015) and by accounting for multiple sources of available information (Gopalaswamy et al., 2012).An SCR approach was applied for the first time to estimate wolf population densities in Spain and Scandinavia, based on NGS data (Bischof, Milleret, et al., 2020;López-Bao et al., 2018) and evaluation of the effects of group association (Bischof, Dupont, et al., 2020).
The presence of large carnivores, especially wolves (Canis lupus), is highly controversial, and assessment of their protection status or decisions about their management are often accompanied by intense political and public debate (Chapron & López-Bao, 2014).At the beginning of the 21st century, wolf populations were mainly small and highly fragmented due to centuries of persecution.Subsequently, multiple concerted efforts to improve wolf conservation were initiated worldwide (Boitani, 2003), and populations have increased across North America and Europe (Chapron et al., 2014), so successfully that some wolf populations have lost their protected status (Carroll et al., 2020).Regardless of their status, wolf populations continue to be at the center of political debates (Darimont et al., 2018), and this intense public scrutiny makes reliable and up-to-date population size estimates highly sought after by various government offices.Despite the recent methodological advances mentioned above, estimating the size of wolf populations is currently even more challenging due to the recent spatial and numeric expansion of the species.Monitoring wolf populations has become more costly and logistically demanding because it requires investigations over large areas with multipolitical jurisdictions and processing of numerous samples with harmonized protocols.In addition, the politically sensitive nature and high public profile of wolf monitoring put pressure on researchers and institutions to obtain accurate estimates regularly and over short periods, despite often limited available resources.
We developed and implemented an integrated approach to monitoring wolves in the Italian alpine region to develop a longterm and evenly distributed coordinated network of specialized operators for the systematic collection of data over an extensive area highly fragmented by administrative boundaries and to produce estimates of key population metrics with direct rele-vance for conservation and management (e.g., wolf abundance and density number of packs, proportion of mature individuals) and uncertainty associated with these metrics.While citizen science is useful in large-scale monitoring (Dance, 2022;Kosmala et al., 2016), its reliability is often questioned (Conrad & Hilchey, 2011;Kalén et al., 2022).To overcome this problem, we trained and coordinated a network of specialized operators to perform NGS and simultaneously record search effort to control for spatial and temporal variation in sampling intensity.Individual genotypes obtained from the NGS were analyzed with novel SCR models that controlled for spatial and individual variation in detection probability to produce the first comprehensive estimates of key parameters for wolf population status assessment in the Italian alpine region.

Sampling design of survey to control effort
Wolves have reestablished in the Western Alps by naturally expanding from the Apennines in the 1990s (Fabbri et al., 2007) and have been reconnected with the Dinaric wolf population in the Eastern Alps since 2014 (Marucco et al., 2022).The wolf historical presence index is described in Figure 1a for 2007-2016.The spatial data in a 10 × 10-km grid are derived from 2 periods: 2007-2011 by Chapron et al. (2014) and 2012-2016 by Kaczensky et al. (2021).Wolf sporadic presence was assigned a value of 1 and wolf permanent presence was assigned a value of 3 in each period and in each cell.We summed the 2 grids from the 2 periods to create an index of historical wolf presence that ranged from 1 to 6 (Figure 1a).These data were used to inform the modeling process (Appendix S1).Wolves are now expanding across the entire Alps and toward foothills.
From 2020 to 2021, we conducted the first comprehensive and coordinated wolf survey in the Italian alpine region to lay the foundation for long-term wolf monitoring.Sampling schemes and protocols were created with the goal of estimating abundance (number of individuals and packs) and species distribution by searching for and collecting environmental DNA left by wolves (scat, urine, hair).The presence of the species was estimated across the entire region of approximately 100,000 km 2 by systematically searching 10 × 10 km cells where signs of presence had been found the previous year and by opportunistically searching all other cells in the Alpine region (Figure 1b).Systematic sampling was based on transects searched once a month from October to April in areas where wolf packs have been detected previously (intensive survey area) and transects searched once every 2 months where only sporadic wolf occurrence had been reported (extensive survey areas) (Figure 1b).No sampling was conducted during the breeding season (May-June), and survey lengths were short, relative to the life history of the species under study, to satisfy the SCR closed population assumption (Dupont et al., 2019;Royle et al., 2014).
This sampling strategy required training and effective coordination of a large network of operators available to travel transects in the entire Alpine region.Twenty-six online and practical training events were held during fall 2020.Overall, 1513 specifically trained operators from 160 institutions or associations constituted the Wolf Network, which operated continuously and extensively throughout the region collecting data with a dedicated mobile app and contributing to opportunistic sampling and systematic surveys along predefined transects on specific dates.The majority of operators (72%) were employees of institutions concerned with wildlife management, and the remaining were trained volunteers.We quantified search intensity and transect operator's experience and used these quantifications in the modeling framework to estimate detection probability (details in Appendix S1).From October 2020 to April 2021, 1776 transects were searched in the area (Figure 1b), for a total of 40,725 km, and there were 1-10 repetitions per transect.

NGS and relatedness analyses
A total of 1918 noninvasive DNA samples collected from 2020 to 2021 were identified to species and individually genotyped.Mitochondrial DNA was used to distinguish wolf samples from other carnivores, including domestic dogs, following Marucco et al. (2009).Wolf samples were individually genotyped and sexed using laboratory-specific sets of 8-16 microsatellite loci.Microsatellite data were checked for quality and potential allelic dropout and false alleles.The detailed laboratory genetic methods are in Appendix S2.We evaluated pack pedigrees and structure for paternal and maternal relationships with spatial information and the programs CERVUS 3.0 (Kalinowski et al., 2007) and ML-RELATE (Kalinowski et al., 2006) to test individual relatedness to other individuals in the putative packs.The relatedness analyses were conducted among individuals in groups sampled within a 150-km radius (Caniglia et al., 2014) or between adjacent areas.Hence, based on additional field data (marking sites, camera traps, reproductive sites, information on dead wolves, etc.), we tested the hypothesis that genotypes belong to the same pack.In some areas, a relatedness analysis was not conducted because of low numbers of detected individuals or because pack reconstruction was not possible-in these cases, the genotype was set as undetermined (NA).The relatedness analyses allowed us to assign genotypes to packs and to the 3 social classes: reproductive individuals (i.e., the breeding pair), offspring, and other (i.e., an unrelated immigrant in a pack, a disperser, or a solitary individual) (details in Appendix S2).

SCR modeling and estimate of population size
We built a Bayesian SCR model (Efford, 2004(Efford, , 2004;;Royle & Young, 2008) to estimate wolf population size and map density from NGS data.An SCR model estimates the distribution of individual activity centers (ACs) from repeated individual spatial detections.The distribution of detected individual ACs over the available habitat (S) can be modeled as an inhomogeneous binomial point process (Illian et al., 2008;Zhang et al., 2022) with intensity calculated as where I h is the point process intensity in habitat grid cell h, X h is the vector of covariate values for habitat grid cell h, and β is the vector of associated coefficients.In our analysis, S was defined as the 5-km resolution grid covering the entire Italian alpine region (Figure 1c) surrounded by a 30-km area that allowed for the possibility that individuals with ACs outside the searched area could be detected within it.
To model the spatial variation in wolf density across S, we considered the additive and linear effect of the following 5 spatial covariates: 1, historical wolf presence (considering that current wolf density is likely the result of the recolonization history of the species [Figure 1a]); 2, percent forest cover; 3, percent low natural vegetation; 4, percent bare rock; and 5, human population density (considering that wolves often select for forested habitat and bushlands and avoid bare rocks and human disturbance) (Falcucci et al., 2013;Marucco & McIntire, 2010).A detailed description of each covariate is in Appendix S1.
We implemented a reversible-jump Monte Carlo Markov chain (RJMCMC) approach to select among the covariates (Green, 1995;O'Hara & Sillanpää, 2009).Following this procedure, each regression coefficient was associated with a binary indicator parameter (w) that determines the probability of inclusion of each covariate and thus its importance: (2) To select the most important covariates, we estimated the posterior inclusion probabilities (PIPs) for each coefficient in isolation and the posterior model probabilities, which were calculated for each possible combination of the covariates considered.
We used data augmentation (Royle et al., 2013) to estimate population size, whereby available but undetected individuals were added to the data set of encounter histories and a Bernoulli state variable z i was used to describe the state of individual i: z i ∼ Bernoulli ψ , where ψ is the probability of an individual from the augmented pool of individuals belonging to the population.Population size (N) was then obtained by summing over the vector z 1∶M : where M is the total number of individuals considered in the augmented pool (M >> N).
We also modeled individual sex and social status: where ρ is the proportion of males in the population and θ sex is a sex-specific vector representing the proportion of individuals in each social status category for each sex ( ∑  sex = 1).Bischof, Dupont, et al. (2020) found that wolves form packs; therefore, spatial association between individuals (i.e., aggregation and cohesion), if ignored, can lead to unreliable inferences.However, results of the aforementioned study show that group association only causes noticeable bias in density estimates in extreme cases (i.e., when the population is configured into very few groups of very large size and high levels of coordination in the space use of group members).This was not the case in our study population; groups (packs) consisted of 2-10 individuals.The computational cost of accounting for grouping in the wolf model (e.g., Emmet et al., 2021) would far outweigh its benefits in our study, given the primary focus on density estimation.
To account for spatial variation in individual detection probability, we used the common half-normal function (Borchers & Efford, 2008), where the probability for individual i to be detected at detector j (p ij ) decreases as distance between that individual's AC s i and detector j (d ij ) increases: and where p 0 sex i status i is the sex-and status-specific baseline detection probability,  sex i status i is the sex-and status-specific scale parameter (Royle et al., 2013), X j is a vector of detector-specific covariates, and  is the vector of associated coefficients.We considered the effects of the following 4 spatial covariates explained spatial variation in p 0 : 1, transect length (total length in kilometers of all transects searched during the sampling season); 2, transect operator's experience, because we expected p 0 to increase as the number of kilometers searched in increased and as the training experience of the operator increased; 3, snowfall (mean accumulated surface snow in the sampling season), because we expected the presence of snow to increase the chances of detecting and collecting genetic material; and 4, human population density, because we expected high detection probabilities in high human population density areas, even if no systematic transect was present.A detailed description of each covariate is in Appendix S1.
For the detection model, we used a 5-km grid covering the entire Italian Alps region and used grid cell centers to define detector locations.Following Milleret et al. (2018), we divided each detector grid cell into 25 subcells of 1 × 1 km and used a partially aggregated binomial detection model.The encounter frequency of individual i at each detector j (y ij ) was then modeled as a binomial distribution with a maximum sample size (K) of 25: We used a single detection model for the collection of systematic and opportunistic samples.Instead of considering separate observation models, we focused on modeling the variation in baseline detection probability among detectors with multiple covariates, and specifically covariates of the systematic search effort.Under this model, we made the assumption that a wolf sample can be collected opportunistically anywhere in the detector grid with a probability varying according to the presence of snow and the local human population density and that the systematic search effort increases this probability to collect a wolf sample according to the length of search track recorded and experience of trackers.
We fitted the SCR model with package nimbleSCR 0.1.3(Bischof, Turek, et al., 2020), NIMBLE 0.10.2 (de Valpine et al., 2017), and R 4.1.3(R Development Core Team, 2022).We ran 4 chains of 50,000 iterations each and discarded the first 10,000 samples as burn-in, resulting in a total of 160,000 Monte Carlo Markov chain (MCMC) samples per model to draw inferences from.We assessed convergence based on the potential scale reduction value for all parameters and mixing of the chains with trace plots (Brooks & Gelman, 1998).For mapping density and extracting abundance estimates, we thinned the posterior samples by 10 and based the maps on 16,000 MCMC samples.All the details on the SCR models are in Appendix S1.

Wolf presence data and NGS
In total, 5636 samples were collected (Figure 1c).After removing samples that failed initial quality checks, 1918 of these were analyzed with microsatellites.Following removal of genotypes suspected of having allelic dropouts or false alleles, we obtained 745 genotypes, of which 95% were from scat samples, 1% were from hair and urine samples, and 4% were from saliva samples (Figure 1c; Appendix S2).These 745 genotypes were associated with 449 different wolf individuals (222 females, 213 males, 14 NA), which represents the minimum number of individuals in the study area.The average recapture rate was 1.66 (SD 1.13) detections per individual and 1.73 (1.23) and 1.62 (1.05) detections for females and males, respectively.We verified the hypothesis that multilocus nuclear genotypes belonged to the same pack based on the relatedness analysis and determined the social status of individuals (Appendix S2).

Abundance and density of wolves in the Italian alpine region
Mean wolf abundance for the entire study area during the winter of 2020-2021 was estimated as 952 individuals (95% credible interval [CrI] 816-1120) (Figure 2).Based on the location of individual ACs predicted by the model, we estimated that 684 (95% CrI 600-782) individuals belonged to the west-

Alpine region
Reproductive *Small deviations between the total estimate and the sum of abundance estimates from the constituent subregions may arise due to rounding.
ern and 268 (95% CrI 201-354) to the eastern part of the population (Table 1).We estimated that between 24% and 34% of the population consisted of reproducing individuals, 55-67% were classified as offspring, and 8-13% were classified as others (Table 1).Because individual social status was modeled as a categorical state variable, it was possible to estimate the presence of 135 breeding units or packs (95% CrI 112-164).In addition to abundance estimates, the SCR model was able to estimate determinants of spatial variation in wolf density (Figure 2).The RJMCMC showed strong support for the influence of the historical presence covariate with a PIP of 1.00; for the influence of low natural vegetation with a PIP of 0.97; and for the influence of human population density with a PIP of 0.79.Percentage of bare rock (PIP bare rocks = 0.41)  ).The RJM-CMC analysis also revealed that percent bare rock and percent forest had correlated and opposite effects.The best model estimated a negative effect of percent bare rock (β bare rocks = −0.91 [95% CrI −1.59 to −0.24]), and the second-best model indicated that wolf density was positively associated with percent forest instead (β fores = 0.73 [95% CrI 0.20-1.27])From 45% to 60% (mean = 53%) of individuals estimated by modeling were undetected by NGS.Baseline detection probabilities (p 0 ) of both sexes were higher for reproductive individuals than for pups and others (Figure 3).Baseline detection probability also varied substantially across space (Figure 4; Appendix S1).We detected a strong positive relationship with search effort (β transect length = 0.20 [95% CrI 0.16-0.26])followed by operator experience (β experience = 0.13 [95% CrI 0.10-0.16]),average amount of snowfall (β snow = 0.23 [95% CrI 0.07-0.40]),and, to a lesser extent, the logarithm of human population density (β logpop = 0.06 [95% CrI −0.07 to 0.18]).Additionally, baseline detection probability was slightly higher in the western part (β west = 0.20 95% CrI [−0.08 to 0.49]), but this effect was not significant.Sex and social status also had a pronounced effect on the scale parameter of the detection function; the largest value was estimated for male others (σ male other = 29.6 km [95% CrI 22.2-40.4])and the smallest for reproductive males (σ reproductive male = 2.4 km [95% CrI 2.1-2.8]).Female others also had large-scale parameters (σ female other = 13.5 km [95% CrI 10.6-17.5])compared with reproductive females (σ reproductive female = 3.7 km [95% CrI 3.2-4.1]),but to a lesser extent than males.Finally, both female and male pups had comparable scale parameters (σ female pup = 3.0 km [95% CrI 2.5-3.5] and σ male pup = 2.8 km [95% CrI 2.3-3.4]).

Estimate of key parameters for wolf population status assessment
Ours is the firstnsive estimation of wolf population size in the entire Italian alpine region that took a systematic and intensively coordinated NGS approach.This survey overcame the challenges of relying on the cooperation of highly fragmented administrative jurisdictions and sampling at a large scale.This first comprehensive estimate of 816-1120 individuals is crucial for wolf management and will serve as a robust baseline for future changes in population trajectory.The SCR models accounted for 2 important challenges associated with largescale population-level estimates: not all individuals present in a study area are detected (Kéry & Schaub, 2012) and individuals that reside primarily outside the surveyed area may be detected within it (Bischof et al., 2016).On the western side of the Italian Alps, the latter is especially important because the wolf population is connected to the French Alps (Louvrier et al., 2018).The SCR models were developed with the intent to relax the geographic closure assumption typical of traditional CR models by explicitly allowing for movements of individuals around their latent and estimated ACs (Royle et al., 2014).Thanks to this explicit link between the estimated number of individuals and geographic area (Efford, 2004), we were able to precisely estimate the density for our region and define the population size.We estimated that 45-60% (mean = 53%) of individuals remained undetected during NGS.Hence, we confirmed that minimum counts based on the number of genotypes are misleading when used to reveal variation in abundance over space and time, whereas SCR estimates that account for varying detection probability should be considered the baseline for future population size monitoring.Finally, we obtained information at a large scale (landscape and population level) and with a high grain (local and individual level).This is rarely the case in ecological studies, which, for logistical reasons, usually have to trade one off for the other (grain vs. extent), 1989).
Our model also provided the first estimate of the number of reproductive units with credibility intervals, thanks to the modeling of the individual social status as a categorical state variable (Appendix S1).We estimated the presence of 135 packs (95% CrI 112-164), which represents the first populationlevel estimate of wolf reproductive units with an associated level of credibility worldwide.This metric is particularly relevant for population management because the number of packs is considered biologically meaningful for the persistence of a wolf population (Reinhardt et al., 2019) because the ecology of the species is based on the social and territorial dynamics of packs (Mech & Boitani, 2003).We also estimated that mature individuals accounted for 33-45% of the population (i.e., number of reproductive individuals and other, excluding offspring).The possibility of estimating the proportion of reproductively mature individuals in the population is also highly relevant because it is the basis for the International Union for the Conservation of Nature Red List assessments based on criterion D (IUCN Standards and Petitions Committee, 2022) and is generally used to evaluate wolf conservation status.Unaccounted individual heterogeneity in detection or space use that is linked with age class or sex could potentially lead to significant bias (Cubaynes et al., 2010;Dupont et al., 2023).However, thanks to the relatedness analysis, we were able to account for differences between sexes and among individual status in space use and detectability, therefore accounting for individual heterogeneity (Cubaynes et al., 2010).
We report estimates for the overall area and distinguished between portions of the population documented in the western and eastern part of the alpine regions (Table 1) because connection among these portions is still limited and naturally differentiated (Figure 1).The 2 portions of the population face different recolonization momentum.In the west, recolonization began 30 years ago.Wolves now completely occupy mountainous areas with densities up to 15 wolves/100 km 2 and are now expanding farther in lowland areas.In the east, the recolonization process is recent, beginning in 2014, and although expansion is evident (Figure 1a), large mountain areas are still being recolonized.Overall, our population abundance estimates showed the magnitude of recovery following the wolf alpine extirpation by the early 1900s.

Coordinated network of operators for long-term wolf monitoring and conservation practice
By conducting the first large-scale wolf survey in the Italian alpine region, we created a large-scale network of specialized operators that are now available for long-term biodiversity conservation in practice.By training operators to collect and record search effort data, we were able to overcome the main challenge of citizen science data, that is, the lack of formalized quantification of spatial-temporal heterogeneity in search effort (Kalén et al., 2022).Moreover, we mainly engaged institutional personnel to overcome the high administrative fragmentation in conservation practices by various government offices.Although the creation of this network was time consuming, results demonstrate its usefulness, and future wolf abundance monitoring plans and analyses should benefit from the infrastructures put in place, decreasing the effort needed for future population assessment.The creation of this network has also been key to the involvement of practitioners, raising awareness and spreading knowledge at the local scale.This effective participatory approach can serve as the foundation for long-term and reliable wolf abundance monitoring, favoring wolf-human coexistence through knowledge dissemination, and provide an important asset for long-term biodiversity conservation.The multidisciplinary approach we adopted, which involved not only the concerted efforts of several jurisdictions through field sampling organized by wolf biologists, but also the coordination and harmonization between genetic laboratories, allowed an unusually high level of interaction and discussion among several disciplines (i.e., wildlife biology, conservation genetics, ecological modeling).This effort resulted in the adoption of an effective sampling design and modeling analysis, which took into account the biology and distribution of the species, local needs by authorities, and optimized survey efforts and relatedness analyses, while minimizing genotyping errors, all of which contributed critically to this analysis of the wolf population in the Italian Alps.

FIGURE 1
FIGURE 1 (a) Wolf historical presence index for each grid cell in the Italian Alps from 2007 to 2016 at a 10 × 10-km resolution (data from 2 documented periods: Chapron et al. [2014] for 2007-2011 and Kaczensky et al. [2021] for 2012-2016), (b) wolf noninvasive genetic sampling in the same region (black lines, search transects; light green, cells sampled only opportunistically [i.e., trained operators were present but no transects were established]; bright green, cells sampled once every 2 months; dark green, cells sampled once every month), and (c) wolf noninvasive genetic samples collected in the winter of 2020-2021 (red) and genotyped samples included in the spatial capture-recapture analysis (blue).

FIGURE 2
FIGURE 2 (a) Wolf density in the Italian alpine region estimated with the spatial capture-recapture model fitted to the noninvasive genetic data collected in winter 2020-2021 and (b) spatial determinants of wolf density in the Italian alpine region (i.e., effects plots for wolf historical presence [WHP]): percentage of low natural vegetation (herbaceous), human population density (human pop), and percentage of bare rock and sparse vegetation according to the best model based on weights from reversible-jump Monte Carlo Markov chain (shading, credible intervals).

FIGURE 3
FIGURE 3 Scale parameter (σ) and baseline detection probability (p 0 ) of the detection function by sex and social status (reproductive [RI], offspring, or other) for wolves in the Italian alpine region in winter 2020-2021 as estimated by the spatial capture-recapture model fitted to the noninvasive genetic data (point, median; shading, density plot; width, frequency; whiskers, quantiles).

FIGURE 4
FIGURE 4 Spatial determinants of wolf baseline detection probability in the Italian alpine region in winter 2020-2021 as estimated by the spatial capture-recapture model fitted to the noninvasive genetic data (shading, credible intervals).

TABLE 2
Model selection for the determinants of wolf density in the Italian alpine region for winter 2020-2021.
Abbreviations: forest, percent forest cover; herbaceous, percent low natural vegetation; pop, human population density; rock, percent bare rock; WHP, historical wolf presence.a Model weights based on posterior model probabilities for the different covariate combinations.Only models with weight >0.01 are shown.