Estimating spatially variable and density‐dependent survival using open‐population spatial capture–recapture models

Abstract Open‐population spatial capture–recapture (OPSCR) models use the spatial information contained in individual detections collected over multiple consecutive occasions to estimate not only occasion‐specific density, but also demographic parameters. OPSCR models can also estimate spatial variation in vital rates, but such models are neither widely used nor thoroughly tested. We developed a Bayesian OPSCR model that not only accounts for spatial variation in survival using spatial covariates but also estimates local density‐dependent effects on survival within a unified framework. Using simulations, we show that OPSCR models provide sound inferences on the effect of spatial covariates on survival, including multiple competing sources of mortality, each with potentially different spatial determinants. Estimation of local density‐dependent survival was possible but required more data due to the greater complexity of the model. Not accounting for spatial heterogeneity in survival led to up to 10% positive bias in abundance estimates. We provide an empirical demonstration of the model by estimating the effect of country and density on cause‐specific mortality of female wolverines (Gulo gulo) in central Sweden and Norway. The ability to make population‐level inferences on spatial variation in survival is an essential step toward a fully spatially explicit OPSCR model capable of disentangling the role of multiple spatial drivers of population dynamics.

E14 crosses the Scandinavian Peninsula from east to west and constitutes a clear barrier for the wolverine population, as little movement of individuals in and out of that area are observed (Gervasi et al. 2016). In total, we used 2443 genotyped samples of female wolverines collected during seven consecutive winters (2014-2020; Appendix S4 Table S1-2). The Scandinavian wolverine population is transboundary, with two strikingly different harvest regimes across the national border between Sweden and Norway due to differences in conflict potential with freeranging semi-domestic sheep and legal regulations (Bischof et al. 2020b). Sweden has a very restrictive harvest of wolverines, while Norway has been much more liberal (Gervasi et al. 2019). Wolverines in Norway are exposed to an annual license hunt as well as an extensive damage control program. Thus, we also obtained recovery locations and genetic identification data from 118 females legally culled during the study period. Another eight individuals were found dead due to other causes (e.g., poaching, car collision; Table S3).

OPSCR model
To estimate survival probabilities of wolverines, we built a Bayesian hierarchical OPSCR model to estimate spatially-explicit causes of mortality as described in paragraph 2.1 of the main text, but with some customization linked to the population and monitoring regime of the wolverine in Scandinavia (Ergon and Gardner 2014, Bischof et al. 2016, 2020b).
The demographic model We used the same demographic model as described in paragraph 2.1.2 of the main text. Individuals considered alive could either die from culling and transition to zi,t=3 with probability ht, or die from other causes and transition to zi,t = 4 with probability wt, so that zi,t ~ dcat(0, Φt, ht, wt ), where Φt = 1−ht −wt. All legal culling mortality events were reported, but most other mortality remains cryptic. Imperfect detection of non-culling mortality prevents further breakdown of estimates by cause-specific mortality, such as natural, traffic, and poaching deaths.
We created a binary covariate (Country) with value 0 for habitat cells in Norway and 1 for habitat cells in Sweden. We modeled individual mortality hazard rates as a function of individual activity center location. We modelled spatial variation in culling and other causes of wolverine mortalities as an additive function of country (Norway = 0, and Sweden=1) and local density at time t (β , β ): For the illustration of the method, we did not account for other sources of variation in mortality. However, it could be possible to estimate temporal variation in mortality, for example by using annually varying slope parameters ( , ) or intercepts ( , ).

The movement model
We used an inhomogeneous point process to model the distribution of individual activity centers (ACs) with a spatial intensity (where s is a vector of x and y spatial coordinates of ACs) (Zhang et al. 2020). We discretized the habitat into a grid of 20 x 20 km habitat cells to allow the placement of individual AC , as a function of a spatial covariate X (main text, Figure 3). The initial individual AC locations , were conditional on X: where , is the value of the spatial covariate at , and β the slope parameter describing the relationship between the habitat covariate and density. Here, X was defined as the average number of known wolverine dens as a proxy for wolverine density (Bischof et al. 2020b. For t>1, the probability density of , , was conditional on the spatial covariate X and the Euclidean distance to , : where τ is the standard deviation of a bivariate normal distribution centered on , . Under this specification, movement is described as an isotropic Gaussian random walk weighted by the spatial covariate X (Ergon and Gardner 2014, Bischof et al. 2020b, Zhang et al. 2020, and τ regulates the distance that individuals are likely to move between years. Such movement feature has been shown to help distinguish between mortality and emigration (Ergon andGardner 2014, Gardner et al. 2018). 2) Average distance from the nearest road (Roads): the distance from each detector to the closest road (1:100,000, the Swedish mapping, cadastral and land registration authority;

The observation model
N50 kartdata, the Norwegian Mapping Authority,). This variable represents accessibility, which we predicted to facilitate detectability.
3) Yearly average percentage of snow cover in each detector grid cell (MODIS at 0.1 degrees resolution, www.neo.sci.gsfc.nasa.gov, accessed 2019-10-11) between December 1-June 31 (Snow). As wolverine NGS during winter relies heavily on the presence of snow, we predicted that greater snow cover increases detectability. 4) Indicator of whether an individual was detected or not during the previous monitoring season (PrevDetection). Previous detection could be expected to positively influence the probability of being detected at subsequent occasions (Gervasi et al. 2014).

5) Country and year. We estimated yearly baseline detection probabilities ( )
for each country to account for differences between the national monitoring regimes. where , 2 is an indicator function used to condition detection on the individual being alive. This design allowed us to reduce the number of detectors j involved in the calculation of , , while retaining as many binary detections as possible (Milleret et al. 2018). In addition, we added a 60 km buffer (>6 , (Sun et al. 2014)) around the detector grid to allow the placement of AC, and therefore the movement of individuals in and out of the trapping grid , Bischof et al. 2020b).
Dead recovery model We followed the procedure described in paragraph 2.1.4 in the main text.
We modelled dead recoveries within the entire habitat (main text, Figure 3). Only the observation of dead recoveries that were caused by culling were modelled using eqn 9 in the main text. This means that we did not use the spatial information contained in the 8 dead recoveries caused by other causes than culling. Dead recovery information from all causes other than culling was only used to set the individual state of individuals zi,t to 4 when the mortality event occurred.
We performed several assessments of SCR/OPSCR model robustness throughout the development process that led to the models used in the analysis. Based on simulations, we assessed robustness to nonindependence between observations (Bischof et al. 2020a), spatial aggregation of detection information (Milleret et al. 2018), population closure (Dupont et al. 2019), ignoring spatial heterogeneity in detection probabilities (Moqanaki et al. 2020) or misspecifying the detection function (Dey et al. 2021). However, while numerous goodness-of-fit (GOF) tests for non-spatial capture-recapture models exist (Pradel et al. 2005), there is a paucity of robust GOF tests for Bayesian SCR and especially OPSCR (Dey et al. 2021).  Figure S1. Yearly abundance of female wolverines within the entire habitat from 2014 to 2020 (main text, Figure 3). Violins show the posterior distribution and points the median estimates.