Basic principles and techniques for wild boar monitoring

Population monitoring is a necessary step in wildlife and habitat management (Engeman et al. 2013, Stephens et al. 2015). In particular, accurate and precise monitoring is fundamental to factual decision-making (Nichols and Williams 2006). Monitoring is equally important in evaluating the performance of control programs and other management actions (Lyons et al. 2008). Population monitoring consists generally of: 1) estimating population density or 2) using relative abundance indexes or quantities related to population abundance (proxies). We expect a tradeoff between the precision and accuracy of these monitoring approaches and their costs, with proxies being less expensive. A sufficent precision of the estimator of population density or index is necessary to design an effective monitoring program (Gerrodette 1987). So, the precision of an estimator must be assessed also for management purposes, otherwise, we cannot infer the variation of population size that a survey is able to detect (Engeman 2005). Common errors in trend analysis include detecting a trend when it does not exist (type I), failing to detect population changes when they exist (type II) and detecting a trend in the opposite direction of the true one (type III) (Nuno et al. 2015). A well-designed survey should be able to reduce these errors to a-priori defined levels, prescribed by the management policy to be implemented, and by the minimum detectable difference to be detected (Skalski et al. 2005, Reynolds et al. 2011).

To reduce costs and complexities of monitoring, wildlife managers often rely on proxies. Anderson (2001) heavily criticized the use of population indexes in wildlife management on the ground of the presence of convenience sampling and uncontrolled impact of covariates. However, later Engeman (2003) contended that under certain circumstances, indexes can be appropriately computed. The use of indexes is common in studies on elusive, wide-ranging and/or low densities species. Nimmo et al. (2015) review evidence that indexes (tracks) in predators and in several other species reflect variation in population abundance; Kojola et al. (2014) successfully used snow track triangles to evaluate wolf Canis lupus abundance in Finland. Index calibration methods are also used routinely to assess lion Panthera leo population trends in Namibia (Dröge et al. 2020). However, recent research (Engeman et al. 2013, Stephens et al. 2015) cast a shadow on the use of proxies, which can be misleading and may introduce severe errors in management actions; for example, when there is not a proportional relationship between the abundance and the proxy (Iannuzzo et al. 2010, Koda and Fujita 2011). Under this circumstance, a proxy cannot demonstrate enough sensitivity at low population densities, and it can exhibit saturation when density is high. Consequently, it becomes useless in these two conditions (scarcity or overabundance), because reliable information is fundamental for management. Further, detectability of the signs used to build the index (e.g. number of observed animals, tracks, pellets) must be considered to avoid the risk of confounding changes in detectability with changes in population abundance (Pollock et al. 2002). Proxies based on the impact of ungulates on vegetation have been also proposed (Morellet et al. 2007), but often the response of vegetation to changes in animal density can be too slow to initiate a quick management response (Tanentzap et al. 2012).

Indexes are often measured only at specific locations and more properly represent site occupancy, or utilization, rather than actual animal density (Putman et al. 2011). Thus, proxies can often be site-specific and then they can hardly be generalized despite the use of rigid and standardised protocols (Williams et al. 2002). Unless calibrated, those indexes/proxies do not provide an estimate of abundance or density that may be required for population management and control, in particular. To be an effective monitoring tool, a population index must be validated at the study area of interest (for instance as done by Corlatti et al. 2016 for the red deer Cervus elaphus, Evans and Rittenhouse 2018, for the black bear Ursus amercanus, and Marchandeau et al. 2006 for the wild rabbit Oryctolagus cuniculus). This literature suggests that often computing a sensible population index may be as complex and expensive as estimating animal density.

Wäber et al. (2013) made a compelling case for the use of density (with other demographic information) instead of indexes by reconstructing the sink–source dynamics of roe deer Capreolus capreolus and Reeve's muntjac Muntiacus reevesi at landscape level. It was concluded that ‘failure to quantify deer numbers and productivity has allowed high density populations to persist as regional sources contributing to range expansion, despite deliberative management programs and without recognition by managers who considered numbers and impacts to be stable'.

For wild boar, there are few attempts to apply pellets group counts (Plhal et al. 2014, Fattorini and Ferretti 2020) or to use hunting statistics (Imperio et al. 2010, Davis et al. 2016). Harvest data can be used if, and only if, harvest effort is known (St Clair et al. 2013). Unfortunately, in the case of sport hunting, effort is often unavailable and harvests are unreported. Thus, the suggestion by ENETWILD consortium et al (2018) to use harvest as a broad index of population abundance can be misleading. Spot-light counts are inefficient with wild boar (Focardi et al. 2001). Camera traps can be used to obtain a population index: Bengsen et al. (2011) monitored changes in relative abundance of feral pigs in a tropical rainforest of northern Australia, and Massei et al. (2017) evaluated the performance of a passive activity index from camera traps to monitor wild boar in five English woodlands. However, a meta-analysis by Broadley et al. (2019) strongly suggests that detection rates are likely to confound the variation in density values with that in movement patterns.

We believe that the situation is more promising in cases where methods are aimed at estimating the population density. The application of the random encounter models (REM) (Rowcliffe et al. 2008) can be effective, but it requires an accurate estimate of the daily range in each study area. New approaches based on camera traps have recently been proposed which do not require an estimate of the daily range and are easier than REM to implement by wildlife managers (Howe et al. 2017, Nakashima et al. 2017). Despite this, however, a huge amount of work remains necessary.

Further, non-invasive genetics, capture–mark–recapture (CMR) and nocturnal distance sampling (NDS) have been used with the wild boar. The first two methods use the same family of statistical estimators: a sample of animals is marked, by molecular or physical tags, and recaptured (or resighted) later. Non-invasive genetics was used by Ebert et al. (2012) to estimate wild boar density in the Palatinate Forest (Germany). Their study produced a precise estimate, with a CV estimate around 10%. However, costs can be elevated. CMR estimates, based on ear tagging, were successfully employed in the Preserve of Castelporziano (Italy) (Focardi et al. 2008, Franzetti et al. 2012), at Fort Benning (USA) by Hanson et al. (2008) and in Switzerland by Hebeisen at al. (2008) using radio tags. The method works fairly well but the cost of trapping and resight is high. Further, since wildlife managers usually wish to remove as much specimen as possible, the idea to trap and release wild boar is not always acceptable.

The use of NDS, via FLIR (forward-looking infrared) technologies was first proposed by Gill et al. (1997). NDS was used for Reeves' muntjac and roe deer (Smart et al. 2004, Hemami et al. 2007, Wäber and Dolman 2015) for white-tailed deer (Montague et al. 2017, Haus et al. 2019), fallow deer Dama dama (Focardi et al. 2013), wild boar and roe deer (Morelle et al. 2012). Franzetti et al. (2012) validated the use of NDS as a method to assess the wild boar population of the Castelporziano Preserve (Italy).

A simulation study aiming to compare REM and NDS was recently proposed by Chauvenet et al. (2017). They concluded that density estimators obtained by distance sampling showed a better accuracy than those obtained by REM. In the context of African Swine Fever prevention in Europe, on behalf of the European Food Safety Authority, ENETWILD consortium et al. (2018) reviewed the different methods to assess wild boar populations and concluded that ‘distance sampling with thermography was recommended to estimate wild boar density on a local scale'. Accordingly, we aim to understand whether the promising results obtained by Franzetti et al. (2012) depended on especially favourable conditions or whether it is possible to obtain meaningful results even in areas characterized by lower visibility and animal density than those observed in Castelporziano. In conventional distance sampling, the variance of the density estimator is mainly determined by the encounter rate (Buckland et al. 2001) and thus the larger the sample size, the more precise the estimator. However, the encounter rate depends largely on environmental visibility. Detection probability is higher in open than in close environment and visibility is a metric that accounts for the environmental component of the detection probability, the others being animal size, behaviour and camouflage. Here, we propose a simple method to estimate visibility, which may help wildlife managers to design an effective NDS survey. Our goal is to define guidelines based on experiments and factual statements in order to design efficient NDS surveys of wild boar.

Study areas

The study took place in four different study areas in Italy (Fig. 1).

Castelporziano Preserve

The completely fenced 60-km² Preserve of Castelporziano consists mainly of natural oak woods, with both evergreen (Quercus ilex and Q. suber) and more open deciduous (Q. cerris and Q. frainetto) species and Mediterranean maquis. Mixed or pure forests of domestic pine Pinus pinea and pastures are also present. The climate is Mediterranean. We distinguish three main habitat types: deciduous wood, evergreen wood and Mediterranean maquis. The distributions of habitats and transects are reported in the Supplementary material Appendix 1 Fig. A1.

Casentino National Park

The study was carried out in the Foreste Casentinesi, Monte Falterona e Campigna National Park (360 km²), along the Tuscan-Romagna Apennine ridge. The landscape is characterized by sedimentary rocks, predominantly sandstone intercalated with marl, and by forests, which cover a large part of the area. We distinguish the submontane habitat (300–900 m a.s.l.), where woods are characterized by a high variety of deciduous species, from the montane habitat (900–1600 m a.s.l.). Here, the dominant species is the European beech Fagus sylvatica, usually associated with the sycamore maple Acer pseudoplatanus or with the European silver fir Abies alba, the most abundant conifer within the park. The distribution of habitats and transects are reported in Supplementary material Appendix 1 Fig. A2.

Monte Arcosu

The study was carried out in the Monte Arcosu – Piscinamanna WWF Reserve (SIC ITB041105) (303 km²) in the southern part of the Sardinian isle. The area can be considered montane although highest peaks are slightly over 1000 m a.s.l. The landscape presents an uneven morphology characterized by long and steep valleys that shaped metamorphic and granite rocks. Forests cover a large part of the area, and the vegetation is dominated by holm and cork oak (Quercus ilex and Q. suber). Mediterranean maquis is widespread and consists mainly of Phillirea spp., various species of Erica, strawberry trees Arbutus unedo and common myrtle Mirtus communis. Pastures originated by fire and past livestock presence are dominated by several species of rockroses Cistus spp. The distributions of habitats and transects are reported in Supplementary material Appendix 1 Fig. A3.

Colli Euganei

The study was carried out in the Colli Euganei Regional Park (187 km²). The area is characterized by the highest hills of the Po River plain, which reach a height of 601 m a.s.l.; the landscape has been shaped by sedimentary and volcanic rocks. The main vegetation types are: pseudo-Mediterranean scrub, scattered throughout the area, growing on dry, south-facing volcanic terrain and consist of evergreen species as holm oak, strawberry trees, tree heath Erica arborea; the chestnut forest Castanea sativa grows preferentially on volcanic slopes facing north in cool, deep, siliceous soils; the oak (mainly downy oaks Q. pubescens) forest is mainly found on south-facing slopes on chalky or siliceous soils; dry grasslands are prevalent in the southern part of the area associated to black locust tree Robinia pseudoacacia woods and agricultural crops, especially vineyard and olive groves. The distributions of habitats and transects are reported in Supplementary material Appendix 1 Fig. A4.

Figure 1.

Location of the fhe four study sites in Italy.

Data collection

Thermal imaging equipment can detect the long-wave energy radiated by warm-bodied animals, allowing the detection of animals at night (Havens and Sharp 2015). As the wild boar is more active at night and prefers habitats with dense vegetation cover, thermal imaging can be an effective tool to significantly increase the probability of detecting groups during a DS survey (Supplementary material Appendix 1 Fig. A5A–D). We used a direct view thermal imager (ThermaCAM P640, FLIR Systems Italia, 640 × 380 pixel resolution) sensitive to infrared wavelengths (7.5–13 µm) achieving image magnification with a ×4 lens or a Thermacam PM545, FLIR, Wilsonville, OR, USA with a ×2 lens. Both cameras were able to take pictures of detected animals and to associate them to a voice memo. Cameras were equipped with a coaxial digital laser range finder (Swarovski Optik, LRS 3-1250, Absam, Tyrol, Austria) and an electronic compass (Outback-ES, Riverton, WY, USA). A GPS receiver (GarminTM eTrex Vista, Sciaffusa, Switzerland) recorded track lines and sighting locations. Data collection involved two observers who took turns at sampling. Each transect was walked by one observer alone who stopped every 5–10 m to scan the area around. At the end, the second observer picked the first one up and drove to the next transect, where they exchanged tasks. Observers used headlamps with a red filter when sampling in dense habitat to minimise human-induced flight response and for safe walking.

Staff had been trained at identifying and recognising the different species as well as at recording precise measures for about one month, first working with an expert trainer and then alone. Observers were trained in large enclosures to collect many measurements and the quality of the collected data was evaluated retrospectively, using pictures and voice notes.

We recorded radial distance and compass bearing to the centre of the group detected, group size, animal activity (feeding, moving, lying or standing still) and fleeing response (yes/no) to the presence of the observer. Observations of flushing animals were ignored unless their original location was recorded or was evident from the heat radiating from the ground where they had been lying. Data were collected between 19:00 h and 06:00 h, paying attention to allocate the same effort in each part of the night (split in three sections: 18:00–22:00, 22:00–02:00, 02:00–06:00 h). To minimise the risk of non-independent detections, transects in close proximity were surveyed in different nights. Perpendicular distances were computed using ArcMap (ESRI, Redlands, CA, USA).

The sampling design is summarized in Table 1. For each study area we denoted two or three habitat classes as appropriate. Each survey was replicated more than once, depending on the available budget and on the length of the research contract. Transects were selected among available paths in order to have enough transects to estimate the variance of the encounter rate and evenly cover the study area, while also ensuring that each transect was long enough to reduce the number of transects with no observations.

Further, for each habitat listed in Table 1 we estimated the level of environmental visibility. We established a set of random points in the different habitats along the used transects. Visibility was estimated using a laser range finder positioned on a tripod at eye level. The range, i.e. the distance of reflection of the laser beam, was then estimated at the four cardinal directions, which are at random with respect to the direction of the transect line. In the 1.5% of sampled points, because of the peculiar topography, some of the measured distances exceeded the capability of the range finder (up to 1000 m) and were removed from the analyses.

Distance sampling analysis

Data were analysed using Distance 6.0 (Thomas et al. 2010). We used only the conventional distance sampling (CDS) analysis with per-habitat post-stratification. We did not use more complex modelling approach (cf. Buckland et al. 2004) preferring to have a quick and simple routine to assess the size or density of populations.

First, we performed some preliminary analyses to select the truncation distance, w. Because the truncation distance may influence the CV of density estimators through variations of the encounter rate (number of observed groups per km of transect), contrary to the approach usually employed in CDS, in this study it was necessary to set fixed truncation distance in order to have detection probability and effective strip width (ESW) estimates comparable among habitats and study areas. Upon inspection of the detection functions, g(x), we selected three truncation distances, w, of 50, 75 and 100 m. Comparison of different w values can be informative on the precision attainable by our NDS surveys. Because of this, the density estimation procedures could not follow Buckland et al. (2001), but results should be considered in a comparative framework, focusing on the attainable precision.

Second, we selected the best model between half-normal, uniform and hazard-rate with cosine adjustment, using AIC for model selection. The mean cluster size was estimated using a size-biased regression method in which log(cluster size) is regressed on the estimated detection probability, g(x).

In CDS, the variance of the density estimator is estimated taking into account three sources of variability, which allow highlighting the possible occurrence of critical survey issues: 1) the variance of the encounter rate, empirically estimated from differences among transect lines, 2) the variance of cluster size and 3) the variance of the detection probability (Buckland et al. 2001).

Simulations

In this study the actual density of the surveyed populations was not known, but we investigated bias and precision originating from our survey design through simulations, using a code derived from La Morgia et al. (2015) and La Morgia and Focardi (2016).

We validated our distance sampling surveys using four different population sizes for each study area. The chosen population sizes were selected to reasonably encompass the actual population size: two populations were smaller than the one estimated during the study (–30% and –25% in size), and one larger population (+20% in size). The survey effort was set as equal, double e triple to the one applied in each study area (Supplementary material Appendix 1 Fig. A1–A4). The rationale of this simulation exercise was that if our sampling designs did not introduced systematic biases with an increase in effort and population size, we should have observed an increase in the accuracy and the precision (%CV) of the estimates produced by the simulations.

Simulations were implemented in R (< www.r-project.org>). To simulate wild boar populations, we first generated cluster sizes from a look-up table of observed clusters specific to each study area. Second, we generated the spatial distribution of independent and identically distributed clusters via a homogeneous Poisson process (applying the function rpoint of the R package spatstat; Baddeley et al. 2015).

For each study area, each population size and each effort set, we performed 500 iterations. To simulate the detection process, we began by calculating the distances between the samplers (transects) and animal clusters. The simulation of the observer's activity was based on a random variate, u, uniformly distributed in (0,1), and the cluster was considered detected if u < g(x).

Given the simulated data, we finally performed separate DS analyses for each set of simulated observations according to standard CDS methodology and using functions provided by the mrds package (Laake et al. 2018). Truncation distance was set to 100 m (i.e. equal to the truncation distance adopted in the analyses (cf. Fig. 3). We used a hazard rate function to model animal detections as a function of distance. The goodness-of-fit of each model was tested through a χ² test, excluding the results of runs if the test yielded p ≤ 0.1. The standard deviations of both bias and CV were empirically calculated.

Table 1.

NDS in the different study areas and habitats in Italy. For the four study areas and habitat types, we report the area surveyed, the number of transects, their mean length and the range (min–max), the survey effort and the starting and ending dates of surveys.

Figure 2.

Box plot of the distribution of the visibility (abscissae, log scale) recorded in each habitat of four study areas in Italy (ordinates).

Survey experiments

In Fig. 2, we report the visibility values recorded in the different study areas and habitats computed as the mean of the distances recorded by the laser range finder. There were important differences in environmental visibility among areas (Wilcoxon test, χ²₃=473.4, p < 0.0001): Casentino (submountain) exhibited the lowest (median score 144) and Castelporziano the highest visibility (median score 955). Even more striking were the differences among habitats within areas (Table 2). Open areas were characterized by visibility values >50 m while dense forests by values <10 m. Note that these distributions are extremely leptokurtic with maximal values observed in maquis (Castelporziano, kurtosis = 121.3, Monte Arcosu, kurtosis = 58.8).

Figure 3 shows the estimated detection functions for the considered habitats. In Casentino N.P. we observed a quite irregular distribution of detections in both habitats, which suggest poor estimates of the detectability. In general, it was noted that habitats with irregular histograms are the ones less sampled (Table 1). On the other hand, Castelporziano is characterized by suitable detectability histograms, especially for the woods. Also, the woodland habitat in Colli Euganei R.P. presented a proper statistical behaviour, while the distribution of detections is more irregular in the fields. Monte Arcosu presented a well-fitting distribution for woodlands and poorly fitting distributions for maquis and pastures. Figure 3 also reports the mean detection probability, which ranges from 0.56 of the submontane habitat of the Foreste Casentinesi N.P. to a low 0.17 of the maquis in Monte Arcosu. In general, note that open areas were not characterised by high detection probabilities as one would expect because of terrain roughness and the presence of bushes and other obstacles that can limit actual detectability. Maquis was always characterised by a smaller detectability than that recorded in woods.

Table 2.

Estimated visibility in the different habitats of four study areas in Italy. We report sample size, visibility, the degrees of freedom used to estimate between habitat difference, χ² and significativity.

Figure 3.

Detection probability plot (red line) and histogram of detections (blue line) in different habitats of four study areas in Italy, w = 100 m. Always data were not significantly different from the estimated g(x) as evaluated by the Kolmogorov–Smirnoff test. The worst fitting (D = 0.084, p = 0.30) was observed for the woodland of Monte Arcosu where the distribution is more spiked near 0 than expected.

We report animal densities and confidence intervals in Fig. 4. In Castelporziano, densities were the highest in all habitats while in Casentino N.P. densities were the lowest. In the other areas there were intermediate density values, ranging between 10 and 20 boar km^–2. There was no significant difference among habitats within the same study area (Table 3). For comparison, we also reported the same analysis run for the encounter rates, which showed, instead, significant differences among habitats within the same study area except for the Casentino N.P.

Figure 4.

Density estimates (±95% confidence limits) in different habitats of four study areas in Italy, w = 100 m.

The effect of the encounter rate on the precision of the density estimates is reported in Fig. 5. As expected, the larger the encounter rate, the higher the precision obtained (F_4,22 = 11.96, p < 0.0001), although significant differences among study areas resulted from the analyses. This relationship was stronger in Colli Euganei, Monte Arcosu and Casentino N.P. than in Castelporziano (F_3,22 = 5.69, p = 0.005). To note that the coefficients of the linear regression between CV estimates and the log of the encounter rates were significant for Colli Euganei (t = 3.2, p = 0.005) and Monte Arcosu (t = 6.1, p < 0.0001) but not for Casentino (t = 0.12, p = 0.90) and Castelporziano (t = 1.05, p = 0.30).

The estimated CV values of the density estimators adopted in the different habitats surveyed are reported in Fig. 6. Half of the surveys were characterised by a CV estimate <30%. As expected from previous analyses, the smallest CV estimates were obtained in Castelporziano and Monte Arcosu, while the estimator for Casentino was characterised by 30 < %CV < 50. Surveys carried out in Colli Euganei produced heterogeneous results. The mode of %CV estimates is around 25%. The three plots of Fig. 7 evidence clearly the effect of truncation. Usually, the smaller the w the better the %CV estimate. For instance, in Monte Arcosu, shortening w from 100 to 75 m reduced the CV estimate from 95% to 55%. This well illustrates the concept that a wise choice of the truncation distance can significantly improve estimation, especially when the sample size is small.

Table 3.

Between-habitat comparison of encounter rate and density. We report the value of the t statistics, the Satterthwaite degrees of freedom to control for inequalities of variances and test probability. Statistics were performed using pre-computed results from DISTANCE 5.0.

The three factors which might influence the precision of the density estimators (the cluster size, the detection probability and the encounter rate) had a variable weight in the different study areas. This is clearly shown in Fig. 7. The estimator of cluster size played a minor role in all study areas. Indeed, we always found a small estimate of mean cluster size (Monte Arcosu, 1.32 ± 0.04; Castelporziano 1.59 ± 0.09; Casentino N.P. 1.79 ± 0.08; Colli Euganei 3.70 ± 0.54). The variance of the encounter rate estimator dominated the overall variance of wild boar densities in Castelporziano and in Casentino N.P. while in the other two areas the contribution given by the estimator of the encounter rate and of the detection probability was comparable.

Simulations

The simulation results are reported in Fig. 8. As expected, the CVs improved with the effort, with values <20% when the sampling effort was tripled. Interestingly, the variability around the average CV is small, indicating a low possibility that a given survey is characterised by a large CV estimate. The bias values also decreased as a function of effort, although there was a tendency toward negative bias. In Castelporziano and in the Casentinesi N.P., the bias was negligible, while in Colli Uganei R.P. it was larger (respectively 15–20%) and in Monte Arcosu it was beyond 20%, even when the sampling effort was tripled. Inside each study area, the CVs and the bias converged with the highest sampling effort, while with lower efforts their values were much more scattered, which was expected.

Figure 5.

%CV of density as a function of the log(encounter rate). On the left Casentino (red), Colli Euganei (green) and Monte Arcosu (grey), On the right the plot for Castelporziano (black). w = 50, 75 and 100 m.

Figure 6.

The distribution of the %CV in the different habitats surveyed in Casentino N.P. (red), Colli Euganei (green), Monte Arcosu (grey) and Castelporziano (black) for three different w values: 50, 75 and 100 m.