Modelling species distribution from camera trap by‐catch using a scale‐optimized occupancy approach

Abstract

Habitat selection is strongly scale‐dependent, and inferring the characteristic scale at which an organism responds to environmental variation is necessary to obtain reliable predictions. The occupancy framework is frequently used to model species distribution with the advantage of accounting for imperfect observation, but occupancy studies typically do not define the characteristic scale of the modelled variables. We used camera trap data from winter wildlife surveys in the Swiss part of the Jura Mountains to model occupancy of wild boar (Sus scrofa) and roe deer (Capreolus capreolus). We used a three‐step approach: (1) first, we identified factors influencing detectability; (2) second, we optimized the characteristic scale of each candidate explanatory variable; and (3) third, we fit multivariable, multiscale occupancy models in relation to land cover, human presence and topography. Wild boar occupancy was mainly influenced by the interaction between elevation within 2500 m and the proportion of forested areas within a 2500 m, with a nonsignificant additional effect of the interaction between ruggedness within 1900 m and the proportion of forested areas within 2500 m as well as the distance to urban areas. Roe deer occupancy was mainly associated with the interaction between ruggedness within 900 m and the proportion of open landscape within 900 m, with an additional nonsignificant effect of the interaction between elevation within 1500 m and the proportion of open landscape within 900 m as well as the distance to urban areas. Incorporating scale optimization in occupancy modelling of camera trap data can greatly improve the understanding of species‐environment relationships by combining the possibility of occupancy models to correct for detection bias and simultaneously allowing to infer the characteristic scale at which certain factors influence the distribution of the organisms studied.

Introduction

Predicting species-environment relationships is at the heart of ecological research. Species distribution is driven by a complex interplay of abiotic and biotic environmental factors such as climate, topography and resource availability, while being influenced by the presence of conspecifics, competitors and predators (Morrison et al., 2006). Different habitat components can drive habitat selection (i.e. the process by which an animal chooses its habitat; Hall et al., 1997) at different spatiotemporal scales (Mayor et al., 2009). The process is often described using hierarchical orders of selection (Johnson, 1980). Species distribution is considered to be related to environmental conditions at broad scales (first-order selection), while an individual's home range placement within this geographic extent (second-order of selection) and habitat use within the home range (third- and fourth-order of selection) are considered to be driven by conditions at medium to small scales. Implicit to Johnson's hierarchical framework, each ecological process that drives habitat selection at one order can be explicitly linked to a specific ecological neighbourhood at which an organism responds to environmental variation (‘characteristic scale’; McGarigal et al., 2016). Habitat selection and habitat use are thus inherently scale-dependent, and inferring the characteristic scale for each environmental variable is crucial to make robust inferences about the relationship between a species and its environment (Mayor et al., 2009; McGarigal et al., 2016).

Camera traps have become a well-established tool to gather occurrence data from many animal taxa and to model their distribution and habitat use (Burton et al., 2015; Caravaggi et al., 2020; O'Connell et al., 2010; Rovero & Zimmermann, 2016). While most camera trap surveys target a focal species, camera traps inherently detect multiple species, and observations of nontarget species provide valuable data for multispecies analyses (Edwards et al., 2018; Mazzamuto et al., 2019). Detectability (i.e. the probability an animal triggers the camera trap and is recorded if it is present) of different species is influenced by specific animal characteristics (e.g. body mass and group size), environmental variables, camera trap model and camera trap setup (Hofmeester et al., 2019). Depending on the research question, variation in detectability resulting from environmental variation or differences in species traits needs to be considered when modelling in order to ensure reliable inferences from camera trap data (Edwards et al., 2018; Hofmeester et al., 2019; Mazzamuto et al., 2019).

Here, we seek to determine the size of the ecological neighbourhood in which environmental variation influences habitat use sensu McGarigal et al. (2016) using occupancy models on camera trap data. Occupancy models are particularly suited to model species distribution and habitat use from camera trap observations (e.g. Niedballa et al., 2015; Oberosler et al., 2017; Rovero et al., 2013; Wevers et al., 2020) as the occupancy modelling framework advantageously uses repeated sampling of detection/nondetection data typical of camera trap surveys to account for imperfect detection (MacKenzie et al., 2017). While resource selection function models (RSFs) applied to presence-only data in a use-available design often explicitly determine the characteristic scale at which an organism responds to each environmental variable (Fattebert et al., 2018; Laforge et al., 2016; McGarigal et al., 2016; Zeller et al., 2017), this is less common in occupancy modelling. Historically, some occupancy studies have incorporated Johnson's hierarchical selection framework into a multiscale approach by contrasting fine- and coarse-scale occupancy models, thereby comparing different orders of selection (e.g. Biggs & Olden, 2011; Hagen et al., 2016; Long et al., 2011; Mordecai et al., 2011; Sunarto et al., 2012). Methods for integrating the characteristic scale explicitly into occupancy modelling (scale-selecting multispecies occupancy models; Frishkoff et al., 2019; Stuber & Gruber, 2020) have recently been developed, but case studies in occupancy modelling that include the characteristic scale are limited (but see, e.g., Frank et al., 2017; Niedballa et al., 2015; Stevens & Conway, 2019). Determining the characteristic scale could be especially important when modelling habitat use of nontarget species for which the survey design was not optimized. As camera trap spacing does not always match a species’ characteristic scale, modelling habitat use at a predefined scale can lead to inaccurate conclusions on which environmental variables are most influential in determining a species’ distribution (Mayor et al., 2009).

To address this research gap, we used by-catch data from a camera trap survey designed to monitor Eurasian lynx (Lynx lynx) in the Swiss Jura Mountains (Foresti et al., 2014; Kunz et al., 2016; Zimmermann et al., 2015) to model wild boar (Sus scrofa) and roe deer (Capreolus capreolus) winter habitat use. In the original study setup, camera trap placement was optimized to enhance lynx detectability for photographic capture–recapture by placing two camera traps opposite each other on forest roads and trails frequently used by lynx. While lynx, like many carnivores, actively use roads and trails, ungulates mostly avoid roads as part of their antipredator response (D’Amico et al., 2016; Harmsen et al., 2010). Hence, we predict that the placement of camera traps on trails or roads to enhance lynx detection negatively affects detectability of both nontarget species (Prediction 1).

Wild boar and roe deer are considered woodland species (Barrios-Garcia & Ballari, 2012; Morellet et al., 2011; Rutten et al., 2019) that have colonized agricultural landscapes (Martin et al., 2018; Morelle et al., 2016; Morellet et al., 2011; Rutten et al., 2019). While both forests and crops provide food resources for wild boar and roe deer, there is a contrast in safety, as the species are exposed in the open landscape but, especially for roe deer, susceptible to lynx predation in the forest. Both wild boar and roe deer are known to adapt their habitat use depending on human disturbance (Bonnot et al., 2013; Fattebert et al., 2017; Fischer et al., 2016; Martin et al., 2018; Morelle & Lejeune, 2015). Roe deer are found to shift between sites with high and low elevation (Coulon et al., 2008; Mysterud, 1999) but tend to avoid rugged terrain (Lone et al., 2014). Global wild boar distribution is limited predominantly by temperature, both high and low (Cordeiro et al., 2018; Markov et al., 2019). In the Jura Mountains, we expect that wild boar occupancy is positively affected by forest availability but negatively affected by elevation and human presence, while roe deer occupancy is positively affected by forest availability but negatively affected by elevation, ruggedness and human presence (Prediction 2).

Instead of relying on the scale determined for the lynx study, we pseudo-optimized the characteristic scale of each variable (McGarigal et al., 2016). Large land cover features influence occupancy on a coarser grain than local features with a less defined edge (Niedballa et al., 2015). In the Jura Mountains, we expected elevation, slope and ruggedness to be influential topographic features and forest, open landscape and urban areas to be influential habitat features driving occupancy of both species. Hence, with regard to scale, we predict that these features drive occupancy on a coarse grain, with the characteristic scale being similar for all features (Prediction 3).

Materials and Methods Study site

The Jura is a subalpine mountainous region that stretches along the border of Switzerland and France. Most of the Jura is covered by forest, both deciduous and coniferous, interspersed with pastures and alpine meadows. Elevation ranges between 484 and 1718 m, with the highest peaks in the south and elevation decreasing towards the north. Annual precipitation ranges between 1300 and 2000 mm (Blant, 2001) and is more pronounced in the West. With increasing altitude, climate becomes harsher, and in winter snow is often found at higher elevation.

Camera trap survey

Xenon white flash camera traps (Cuddeback Ambush, Cuddeback Capture or Cuddeback C1) were deployed in a lynx monitoring programme. At each site, two camera traps were mounted opposite each other at a height of 70 cm for at least 60 days (see Zimmermann & Foresti, 2016, for details). A total of 200 sites were sampled in the period February–April over three consecutive years. In the central Jura region, 53 sites were sampled in 2014, 51 sites in the southern region in 2015 and 96 sites in the northern region in 2016. Camera trap sites were selected using systematic nonrandom sampling: a grid of 2.5 × 2.5 km was superimposed on the area, after which deployment sites were chosen in every second grid cell by identifying forest roads and trails frequently used by lynx to maximize the probability of lynx detections for modelling their abundance and density using photographic capture-recapture (Foresti et al., 2014; Kunz et al., 2016; Zimmermann et al., 2015).

Scale-optimized occupancy modelling

Occupancy models use a detection matrix in which detection and nondetection are recorded for a number of consecutive occasions (i.e. surveys). We pooled detection data from the two camera traps at a site, which was shown to increase detection probability and improve precision of occupancy and detection estimates (Evans et al., 2019; Wong et al., 2019). To reduce nondetections and facilitate modelling, we pooled the detection/nondetection data of three consecutive 24-h periods to define one detection occasion (Shannon et al., 2014). We constructed species-specific detection matrices for wild boar and for roe deer using the CamtrapR package (version 1.2.3, Niedballa et al., 2016). We modelled occupancy using the Unmarked package (version 0.13-2, Fiske & Chandler, 2017) after standardizing explanatory variables to a mean of zero and standard deviation of one using the decostand function in the Vegan (version 2.5-6) package in the R environment (version 3.6.3, R Core Team, 2019).

We followed a three-step approach, in which we combined the two-step process of identification of detection and occupancy variables (Ciarniello et al., 2007; MacKenzie et al., 2017); Niedballa et al., 2015) with a scale-optimization approach borrowed from recent multi-scale RSF modelling (Fattebert et al., 2018; Laforge et al;, 2016; McGarigal et al., 2016; Zeller et al., 2017): (i) we first identified variables influencing detectability; (ii) we then optimized the size of the ecological neighbourhood, or characteristic scale, at which organisms respond to each environmental variable in the occupancy component of univariate occupancy models; and finally, (iii) we built multivariable occupancy models using the identified detection variables(s) and the scale-optimized occupancy variables.

1. Detectability modelling: site-specific and observation-specific covariates

   Imperfect detection—that is, when a species present is not observed—is known to introduce bias in species distribution models (Denes et al., 2015; Gu & Swihart, 2004). Through repeated sampling, occupancy models allow incorporating the underlying detection process at a specific site and relating the detection probability to site-specific or observation-specific parameters (MacKenzie et al., 2017). We included camera trap effort as a possible observational covariate on the detection probability. We defined camera trap effort as the sum of the number of days camera traps were active during a 3-day occasion, with a minimum of 1 day when only one of the camera traps was functioning on one of the 3 days and a maximum of 6 days when both camera traps were functioning for 3 days. We coded camera trap failures (i.e. when the camera trap was not active due to technical problems) as NA in both the effort and detection matrix. We included the year of deployment as a site covariate potentially influencing detectability. However, the three regions in our study site were sampled in three consecutive years; therefore, year of deployment is confounded with survey region, and an effect of year of deployment could indicate either an effect of sampling year or an effect of region. As camera traps were deployed on game trails, hiking trails or forest roads for optimal lynx detection, we added trail type as a site-specific covariate in our detection models. We corrected for possible differences in detection probability due to the varying sensor sensitivity between the three camera trap models used by adding the combination of camera trap models that were deployed at each site as a site covariate. We fitted single season occupancy models with a fixed null model for occupancy (Ѱ) and different possible combinations of the detection (p) covariates. We used model selection to determine the most influential covariates through comparison of AIC values. Even when several models were within ∆AIC <2, we chose to continue with the detection variable with the lowest AIC value.

2. Scale optimization: site-specific covariates

   We deemed mean elevation, mean slope and mean ruggedness to be influential topographic parameters limiting distribution and extracted them from a digital elevation model (EU-DEM). We used CORINE satellite imagery (Copernicus Land Monitoring Service, 2018) to reclassify the available land cover types into urbanized area, forested and open vegetation (lumping crops and meadows). To incorporate the effect of human presence, we included the minimal distance from the camera trap to the nearest urbanized area.

   Spatial scale consists of two components: grain, the spatial resolution measured, and extent, the largest area of investigation (Mayor et al., 2009). Here, we sought to determine the characteristic scale that influences how an organism responds to environmental variation, that is, the optimal grain (Laforge et al., 2016; McGarigal et al., 2016). Extent, in this case, is defined as the size of the study area. Hence, we extracted candidate variables at different grains (Laforge et al., 2016) starting at the smallest grain possible to calculate the proportion of land cover from CORINE and increasing with one pixel in each direction until both wild boar home ranges (mean: 4 km², range: 1–7 km²; Fattebert et al., 2017) and roe deer home ranges (mean: 0.8 km², range: 0.5–1.4 km²; Morellet et al., 2013) were included, with the largest grain corresponding to some of the larger wild boar home ranges measured near the Jura Mountains. We first resampled all candidate variables to a pixel resolution of 100 × 100 m based on CORINE land cover. Subsequently, we calculated for each focal 100 × 100 m pixel the proportion of forest and open landscape and the mean elevation, slope and ruggedness at grains having a diameter of 300, 500, 900, 1500, 1900 and 2500 m with the camera site as the central point. To account for the fact that risk effects (e.g. from predation or human disturbance) persist at different scales for different species, we optimized the scale at which the species respond to the distance to urban areas by using linear distance to urban areas as well as with a varying exponential decay function (Whittington et al., 2011) to model occupancy. We used a base exponential decay (i.e. 1−exp^{(−2×distance)}) as defined by Whittington et al. (2011), which accommodates for the fact that risk is highest at the source and dissipates nonlinearly when moving away from the source. Additionally to the linear distance, we varied the exponential term to mimic a slow decay of risk effects (i.e. risk effects persist over greater distances), an intermediate decay and a fast decay of risk effects (i.e. risk effects persist over smaller distances only).

   To determine the species-specific characteristic scale for each variable, we contrasted univariate occupancy models for each variable at the six predefined grains or four exponential decays. We fitted univariate occupancy models (Ѱ) using the detection covariates that were identified as most important in the previous step (p). We used an information-theoretic framework to determine the best grain for each variable through comparison of AIC values. When several models were within ∆AIC < 2, we chose to continue with the grain with the lowest AIC value as in conventional RSF approaches (Fattebert et al., 2018; Laforge et al., 2016; Zeller et al., 2017).

3. Scale optimized occupancy modelling

   We checked for collinearity between explanatory variables at their optimal grain by calculating the variance inflation factor (VIF) and correlations with a threshold Spearman rho |r_s| = 0.7 (Dormann et al., 2013). Occupancy models assume that occupancy states remain constant within the study period and that detection histories between sampling sites are independent, that is, temporal and geographical closure (MacKenzie et al., 2017). As the minimal distance between camera traps was smaller than the home range size of the focal species, sampling sites are possibly spatially dependent. Moreover, as both species are highly mobile, geographical closure at a sampling site is unlikely. Hence, we interpreted occupancy estimates as a reflection of both species’ space use (MacKenzie et al., 2017; Wevers et al., 2020). Although pooling data across years can be considered a violation of the temporal closure assumption, it was applied successfully in other occupancy studies (Fuller et al., 2016; Linden et al., 2017). Moreover, sampling occurred in a different part of the Jura each year, and none of the sites were resampled in two different years, so we assumed negligible impact of the violation of temporal closure.

   Based on Prediction 2, we built a priori additive and two-way interaction models (Tables S1 and S2) by combining occupancy variables at their optimal grain (Ѱ) using the detection covariates that were identified as most important in the previous step (p). We used model selection to determine the top model through comparison of AIC values. When several models were within ∆AIC < 2, we used model averaging from the AICcmodavg package (Mazerolle & Mazerolle, 2019) on the top models, showed averaged coefficients (β) to determine the importance and direction of the respective covariates and plotted the response curve of each variable. We predicted occupancy based on the averaged top models and projected the prediction across the study area at the original 100 × 100 m resolution using standardized covariates at their optimal grain for each grid cell. Lastly, we assessed if scale-optimization improved model inferences by contrasting model selection and predictive power (pseudo-R²) of individual models between occupancy models using covariates at their optimal grain versus occupancy models using covariates at one of the six predefined grains.

Results Detectability modelling: site-specific & observation-specific covariates

For wild boar, detection probability was best explained by the Effort + Trail Type model (Table 1). Compared to camera traps placed on forest roads (Reference category), detectability was lower when camera traps were placed on hiking trails (β_hike = −0.71 ± 0.31) and higher when camera traps were placed on game trails (β_game = 0.61 ± 0.32). Detectability was positively influenced by effort (β_effort = 0.30 ± 0.17).

Table 1 Model selection for four detection covariates [Year of deployment, Camera trap Effort, Combination of Camera trap types (Cuddeback Ambush, Cuddeback Capture or Cuddeback C1) and type of trail (forest road, hiking trail and game trail) on which the camera trap was deployed] included in occupancy models with constant occupancy (Ѱ) fit to camera trap data of wild boar and roe deer collected in winters 2014–2016 in the Jura Mountains, Switzerland.

Wild boar			Roe deer
Model structure	AIC	∆AIC	Model structure	AIC	∆AIC
˜Effort + TrailType˜1	986.86	0	˜CameraTrapType˜1	2562.87	0
˜TrailType˜1	988.88	2.02	˜Effort + CameraTrapType˜1	2563.49	0.62
˜TrailType + CameraTrapType˜1	989.80	2.93	˜TrailType + CameraTrapType˜1	2563.75	0.88
˜Effort˜1	993.62	6.76	˜TrailType˜1	2576.76	13.90
˜Year + TrailType˜1	995.35	8.49	˜1˜1	2577.20	14.33
˜1˜1	995.74	8.88	˜Effort + TrailType˜1	2577.50	14.64
˜Effort + CameraTrapType ˜1	996.52	9.65	˜Effort˜1	2578.22	15.35
˜CameraTrapType˜1	998.24	11.37	˜Year	2579.80	16.93
˜Year	1012.14	25.27	˜Year + CameraTrapType ˜1	2584.17	21.30
˜Year + CameraTrapType ˜1	1012.14	25.28	˜Year + TrailType˜1	2585.21	22.34
˜˜Year + Effort + TrailType + CameraTrapType˜1	1016.74	29.88	˜Year + Effort + TrailType + CameraTrapType˜1	2588.19	25.32
˜Year + Effort ˜1	1016.86	30.00	˜Year + Effort˜1	2596.26	33.39

The ranking of variables influencing the detection of both species is shown with corresponding AIC and ΔAIC values.

Roe deer detectability was best explained by the Camera trap Type model, although equally parsimonious as the Effort + Camera trap Type model or the Trail Type + Camera trap Type model (Table 1). We kept the Camera trap Type model in all subsequent roe deer occupancy models. Compared to when two Ambush type camera traps (Reference category) were combined, detectability was lower when two C1 type camera traps (β_C1/C1 = −0.76 ± 0.36), an Ambush and C1 type camera trap (β_ambush/C1 = −0.16 ± 0.34), an Ambush and Capture type camera trap (β_{ambush/capture} = −0.46 ± 0.18), an Ambush, Capture and C1 type camera trap (β_{ambush/capture/C1} = −0.83 ± 0.39) were deployed together and higher when two Capture type camera traps (β_{capture/capture} = 0.08 ± 0.16) or a Capture and C1 type camera traps (β_capture/C1 = 0.03 ± 0.33) were deployed together.

Scale optimization: site-specific covariates

For wild boar, the Open, Slope and Ruggedness variables fitted the data best at a grain of 1900 m (Table 2). The Forest and Elevation model fitted the data best at a grain of 2500 m. Urban performed best with a fast exponential decay function.

Table 2 Model selection of univariate occupancy models with six different grain sizes for five explanatory variables of occupancy probability (Ѱ) (open landscape, forest, elevation, ruggedness, slope) and four exponential decays for distance to urban areas. Models were fit to camera trap data of wild boar and roe deer collected in winters 2014–2016 in the Jura Mountains, Switzerland, using previously identified detection covariates (p) and varying occupancy (Ѱ).

Wild boar					Roe deer
Model name	Model structure	Grain	AIC	∆AIC	Model structure	Grain	AIC	∆AIC
Open	˜TrailType + Effort ˜ % Open	1900	982.82	0	˜CtType ˜ % Open	900	2551.82	0
		2500	983.98	1.16		1500	2552.54	0.72
		1500	984.85	2.03		1900	2556.10	4.28
		300	986.17	3.35		2500	2557.27	5.44
		500	987.25	4.43		500	2563.02	11.20
		900	988.55	5.73		300	2564.77	12.95
Forest	˜TrailType + Effort ˜ % Forest	2500	982.36	0	˜CtType ˜ % Forest	1500	2554.00	0
		1900	982.80	0.44		900	2554.21	0.21
		1500	985.14	2.77		1900	2557.02	3.02
		300	986.33	3.97		2500	2557.95	3.96
		500	987.19	4.82		500	2563.29	9.29
		900	988.60	6.24		300	2564.79	10.79
Elevation	˜TrailType + Effort ˜ mean Elevation	2500	980.86	0	˜CtType ˜ mean Elevation	900	2560.91	0
		1900	981.71	0.85		1500	2561.02	0.11
		1500	982.21	1.35		500	2561.28	0.37
		900	982.70	1.84		1900	2561.39	0.48
		500	983.05	2.19		250	2561.42	0.51
		300	983.15	2.29		2500	2561.80	0.89
Ruggedness	˜TrailType + Effort ˜ mean ruggedness	1900	986.52	0	˜CtType ˜ mean Ruggedness	900	2564.55	0
		2500	986.77	0.25		1500	2547.35	0.80
		1500	987.24	0.71		1900	2549.45	2.90
		900	988.44	1.92		2500	2551.10	4.55
		500	988.73	2.21		500	2554.59	8.04
		300	988.75	2.22		300	2559.43	12.88
Slope	˜TrailType + Effort ˜ mean slope	1900	986.58	0	˜CtType ˜ mean Slope	900	2548.56	0
		2500	986.79	0.21		1500	2549.17	0.60
		1500	987.26	0.68		1900	2550.97	2.41
		900	988.47	1.90		2500	2552.35	3.78
		500	988.71	2.13		500	2555.72	7.15
		300	988.79	2.21		300	2559.98	11.42
Urban	˜TrailType + Effort ˜ Distance to Urban	Exponential Fast	988.66	0	˜CtType ˜ Distance to Urban	Linear	2560.78	0
		Exponential	988.73	0.08		Exponential Fast	2561.96	1.17
		Exponential Slow	988.84	0.18		Exponential	2562.70	1.92
		Linear	988.85	0.19		Exponential Slow	2564.36	3.58

In the model structure column, detection covariates are shown on the left-hand side of the formula [i.e. Camera trap Effort, Combination of Camera trap types (Cuddeback Ambush, Cuddeback Capture or Cuddeback C1) or Type of trail (forest road, hiking trail and game trail) on which the camera trap was deployed] while occupancy covariates are shown on the right-hand side of the formula. The ranking of variables influencing the occupancy of both species is shown with corresponding AIC and ∆AIC values.

In roe deer, we found that the Open, Elevation, Slope and Ruggedness models fitted best at a grain of 900 m, while Forest fitted best at a grain of 1500 m (Table 2). Urban performed best with a linear decay. In both species however, several variables were equally parsimonious with one or more scales (Table 2).

Occupancy modelling: site-specific covariates

Correlations exceeded threshold values for forest and open landscape (VIF_{WildBoarForest} = 7.76, VIF_WildBoarOpen = 6.74; VIF_{RoeDeerForest} = 3.88, VIF_RoeDeerOpen = 3.95; |r_s| = 0.88) and for ruggedness and slope (VIF_{WildBoarRuggedness} = 160.45, VIF_{WildBoarSlope} = 162.09; VIF_{RoeDeerRuggedness} = 169.83, VIF_RoeDeerSlope = 170.48; |r_s| = 0.99). As we expected ruggedness and elevation to be more important features than slope in driving habitat use of both species (Coulon et al., 2008; Mysterud, 1999), we excluded slope from further analyses, and hence, our Topography model only included elevation and ruggedness in the analysis. As both forest and open landscape were considered important habitats in driving habitat use of wild boar and roe deer, we kept both but ran separate models for each covariate.

Occupancy of wild boar was best explained by the Forest × Elevation interaction model (Table S1), although equally parsimonious as the Topo × Forest + Urban model. Averaged β estimates indicated a significant effect of the interaction between forest and elevation but not of the interaction between forest and ruggedness nor the distance to urban areas (Table 3). Bivariate and response plots (Fig. 1) illustrated modelled effects on wild boar occupancy. Spatial projection (Fig. 3) showed that wild boar occupancy was highest in the northern part of the study area.

Table 3 Averaged parameter estimates, standard error and confidence intervals of variables (i.e. proportion of forest, proportion of open landscape, mean ruggedness, mean elevation and distance to urbanized areas) of top-ranked (i.e. within 2 ∆AIC) wild boar occupancy models (Forest × Elevation and Topo¹ × Forest + Urban) and top-ranked roe deer occupancy models (Open × Ruggedness, Topo¹ × Open and Topo¹ × Open + Urban). Models were fit on camera trap data of wild boar and roe deer in winters 2014–2016 in the Jura Mountains, Switzerland.

Wild boar				Roe deer
Variable	Averaged β estimates	se	CI	Variable	Averaged β estimates	se	CI
Forest	0.43	0.22	−0.01, 0.87	Open	0.54	0.24	0.07, 1.02
Elevation	−0.70	0.24	−1.17, −0.22	Elevation	−0.45	0.27	−0.98, 0.08
Ruggedness	−0.03	0.14	−0.58, 0.39	Ruggedness	−1.33	0.37	−2.04, −0.62
Elevation × Forest	0.54	0.24	0.07, 1.02	Elevation × Open	−0.01	0.29	−0.59, 0.57
Ruggedness × Forest	0.09	0.19	−0.13, 0.77	Ruggedness × Open	−1.10	0.39	−1.86, −0.33
Urban	0.10	0.21	−0.08, 0.83	Urban	−0.11	0.19	−0.49, 0.27

¹The ‘Topo’ variable includes elevation and ruggedness.

Generate image description

Figure 1. Response surfaces showing the effect of the interaction between forest and elevation, the interaction between forest and ruggedness and the distance to urban areas on predicted occupancy (Ѱ) based on the averaged top models of wild boar in the Swiss Jura Mountains, winter 2014–2016. Optimal scale is indicated for each variable as subscript. The grey shaded area in the urban response plot represents confidence intervals.

Roe deer occupancy was best explained by the Topo × Open, Open × Ruggedness and Topo × Open + Urban models (Table S2). Averaged β estimates indicated a significant effect of the interaction between ruggedness and open but a nonsignificant effect of the interaction between elevation and open or the distance to urban areas. Bivariate and response plots (Fig. 2) illustrated modelled effects on roe deer occupancy. Spatial projection (Fig. 3) showed that roe deer occupancy was high throughout the Jura, except on the mountain ridges.

Generate image description

Figure 2. Response surfaces showing the effect of the interaction between forest and elevation, the interaction between forest and ruggedness and the distance to urban areas on predicted occupancy (Ѱ) based on the averaged top models of roe deer in the Swiss Jura Mountains, winter 2014–2016. Optimal scale is indicated for each variable as subscript. The grey shaded area in the urban response plot represents confidence intervals.

Generate image description

Figure 3. Spatial projection of wild boar (top) and roe deer (bottom) probability of presence, based on the averaged top occupancy models, across the Jura Mountains, Switzerland, winters 2014–2016. Camera trap sites are indicated with open circles for nondetections and opaque circles for detections of the respective species.

Modelling occupancy of either species at each of the six a priori defined grains without performing scale optimization resulted in different models being retained based on AIC for several a priori chosen grains (Wild boar—Table S3, Roe deer—Table S4). Additionally, the most important factors driving occupancy at the optimal grain appeared nonsignificant when modelled at some of the six a priori grains (Wild boar—Table S5, Roe deer—Table S6). Moreover, we found that our final scale-optimized models performed better than models at some a priori chosen grains based on R² (Wild boar—Table S7, Roe deer—Table S8).

Discussion

We used a three-step approach to model scale-optimized occupancy of roe deer and wild boar. Scale optimization showed that the scale at which variation in the studied environmental covariates influences wild boar and roe deer habitat use corresponded to the species’ respective home range sizes. Occupancy of both species was best explained by a combination of land cover types, topography and distance to urbanized areas. Our results illustrate that modelling occupancy of both species at a priori chosen scales without scale optimization would influence the ecological interpretation on which environmental variables are most influential in driving the species’ habitat use. We argue that incorporating scale optimization as performed routinely with RSFs can improve inferences from occupancy models by accounting for the characteristic scale at which environmental factors drive habitat selection, thereby revealing the appropriate scale at which to model habitat use.

Detectability

As expected (Prediction 1), camera trap-placement targeted at maximizing lynx observations (i.e. placing camera traps on forest trails or roads) had an impact on the detectability of both nontarget species, but using an occupancy framework allowed us to account for these survey-related impacts on detectability. Additionally, consistent with the general findings of Hofmeester et al. (2019), detectability of wild boar and roe deer in our study was influenced by general survey protocol features (i.e. sampling effort and camera trap type). Our results suggest a possible avoidance of roads for wild boar (D’Amico et al., 2016) and emphasize the effect that combining camera trap types with varying sensor sensitivity can have on detectability. While we considered roe deer detectability to be primarily influenced by the type of camera traps deployed on a site, detection models retained within AIC <2 indicated an additional effect of effort and trail type.

Optimized occupancy

We expected that wild boar habitat use would be mostly influenced by elevation, forest availability and human disturbance (Prediction 2). We found that wild boar scale-optimized occupancy was mostly influenced by the interaction between forest availability and elevation, with an additional nonsignificant effect of the interaction between forest availability and ruggedness and the distance to urban areas (Fig. 1). Predicted wild boar occupancy was lowest at high elevation sites with low forest availability but increased with increasing forest availability or decreasing elevation and peaked at high elevation sites with dense forests as well as low elevation sites with sparse forests. As our surveys were carried out in winter, avoidance of high elevation sites without forest could possibly be explained by lower temperatures (Markov et al., 2019), more snow (Thurfjell et al., 2014) and less food availability, while the use of forests at higher altitudes illustrates wild boar as a typical woodland species, exploiting the canopy cover while foraging for resources such as beech-nuts (Ballari & Barrios-Garcia, 2014). At low elevation, greater use of sites with a higher availability of open habitat could be explained by wild boar exploiting left-overs from agricultural crops in the open landscape while hiding in forest fragments (Ballari & Barrios-Garcia, 2014; Rutten et al., 2019).

For roe deer, we expected their habitat use to be mainly driven by elevation, ruggedness, forest availability and human disturbance (Prediction 2). We found that roe deer scale-optimized occupancy was influenced mostly by the interaction between ruggedness and the availability of open landscape, with an additional nonsignificant effect of the interaction between elevation and the availability of open landscape and the distance to urban areas (Fig. 2). Predicted roe deer occupancy was highest at flat open terrain and decreases with increasing ruggedness. While wild boar forages on left-over roots of crops, roe deer mainly use the aboveground part of crops (Morellet et al., 2011). As roe deer is one of the main prey species of lynx in the Jura Mountains, its preference for open areas could be explained by a predator avoidance strategy, although lynx presence would more likely influence fine scale (i.e. Johnson’s [1980] third order) within habitat use rather than home range placement or species distribution (Gehr et al., 2017; Lone et al., 2014; Molinari-Jobin et al., 2002).

Scale optimization

Instead of choosing a scale a priori to model habitat selection, we pseudo-optimized a set of relevant scales to tease out the characteristic scale of each variable (McGarigal et al., 2016). While we modelled occupancy as a measure of habitat use, the scale optimized approach enabled us to identify the characteristic scale at which these environmental variables drive habitat use. Wild boar occupancy models performed better when landscape variables were extracted at grains between 1900 to 2500 m, while roe deer occupancy models performed better when landscape variables were extracted at grains between 900 to 1500 m (Prediction 3). Our results suggest that the environmental features we considered in our models drive habitat use at a scale corresponding to an average home range size of each species, therefore being factors that influence home range choice, analogous to Johnson (1980) second order of selection (i.e. home range placement at the landscape scale). Our results thus support general guidelines for sampling occupancy at the home range scale of the study species.

We found that our scale-optimized models performed better than models fitted with a priori chosen grains, either larger or smaller. Importantly, occupancy modelled at a priori defined grains would lead to different ecological interpretations of the species’ habitat use. We did find that, for several variables, the best ranking grain did not fit the data significantly better than one or more other grains, but the best ranking grains were similar in scale, for example, roe deer models at grains of 900 and 1500 m. This possibly indicates that the differences between grains were simply not large enough, as shown by comparison of mean, standard error, minimum and maximum (Table S9). We argue that modelling habitat use at an arbitrary predefined grain can lead to inaccurate conclusions on which environmental variables and scales are most influential in determining a species’ distribution, consistent with Mayor et al. (2009).

We chose to cap our range of grains to the grain of the original sampling grid designed for lynx, which corresponded to the upper limit of measured wild boar home ranges in the Jura Mountains. The characteristic scale for several of the variables studied corresponded to the largest grain measured, possibly indicating that they might be affecting species occurrence at an even larger scale (Jackson & Fahrig, 2015). Moreover, by performing scale optimization with univariate models only, we could not account for the fact that interactions between variables could influence the scale at which these variables drive habitat use nor for the fact that a variable could influence habitat use at multiple scales (Frishkoff et al., 2019). Being less complex, our approach is an easy to use extension of the typical two-step occupancy approach using Unmarked. Hence, by introducing scale optimization borrowed from the RSF literature, our approach provides a simple, yet effective improvement to modelling habitat use that is intuitive to implement for single species, single season occupancy modelling.

The main focus of this study was to incorporate the characteristic scale in modelling the habitat use of wild boar and roe deer in an occupancy framework. We see added value of scale optimization in modelling by-catch data from camera trap surveys designed to target one specific species. By demonstrating that the optimal grain corresponds to the respective home range of wild boar and roe deer, we suggest that, in our case, we modelled the home range placement process (Johnson, 1980—second order). Optimizing the characteristic scale of each candidate variable prior to occupancy modelling therefore ensured that, regardless of the spacing and placement of the camera traps in the original sampling design, habitat use at the camera sites was modelled at the relevant scale for the factors driving occurrence of these two nontarget species in the study area.

Conclusion

Being able to predict which environmental factors drive species distribution is essential for the design of evidence-based management and conservation policies. By firstly determining the characteristic scale of candidate environmental variables and subsequently incorporating these in a multiscale occupancy model, we can reach an optimal match between the studied environmental variables and a species’ requirement. We showed that large land cover features influenced wild boar and roe deer occupancy in our study at relatively coarse grains and that habitat type and topography seemed to be the key drivers in habitat use at the home range scale. Specifically, we found that the winter distribution of wild boar was primarily dependent on a trade-off between elevation and forest habitat. Roe deer occupancy was mostly directed away from rugged terrain while being modulated by the availability of open landscape. We argue that by incorporating the characteristic scale into occupancy models when analyzing camera trap data, we improved our habitat selection inferences by combining the strength of occupancy modelling in correcting for detection bias with the ecological relevance of scale optimization borrowed from the RSF literature.

Acknowledgments

This work makes use of data and/or infrastructure provided by the foundation Carnivore Ecology and Wildlife Management (KORA) and funded by the Federal Office for the Environment. We thank especially D. Foresti, F. Kunz as well as civilians and volunteers for their help during field work and data entry and R. Bürki for checking the data integrity. This work would not have been possible without the support of the wildlife managers of the cantons of Baselland, Bern, Jura, Neuchâtel and Vaud and the help of their game wardens. JW is funded by a BOF-mandate from Hasselt University.