Quantifying spatiotemporal occupancy dynamics and multi‐year core‐use areas at a species range boundary

Many species face large‐scale range contractions and predicted distributional shifts in response to climate change, shifting forest characteristics and anthropogenic disturbances. Canada lynx (Lynx canadensis) are listed as threatened under the U.S. Endangered Species Act and were recently recommended for delisting. Predicted climate‐driven losses in habitat quality and quantity may negatively affect the north‐eastern Minnesota lynx population, one of the six remaining resident populations in the contiguous United States. We develop a large‐scale monitoring protocol and dynamic occupancy modelling framework to identify multi‐year core‐use areas and track spatiotemporal occurrence at the southern periphery of the species range.


| INTRODUC TI ON
Many species face large-scale range contractions and predicted range shifts in response to climate change, shifting forest characteristics and anthropogenic disturbances (Gooliaff & Hodges, 2018;Karanth, Nichols, Hines, Karanth, & Christensen, 2009;Lenoir & Svenning, 2015). Populations near a species range periphery may be particularly susceptible to extinction, while also providing valuable opportunities for species conservation (Channell & Lomolino, 2000;Lenoir & Svenning, 2015). For example, range peripheries may represent suboptimal habitats that are regularly colonized but suffer from low persistence probabilities and possible loss of connectivity with larger source populations (Kéry, Guillera-Arroita, & Lahoz-Monfort, 2013). Conversely, populations near range peripheries can provide important sources for long-term conservation due to increased speciation potential and greater survivorship when core populations undergo strong declines (Channell & Lomolino, 2000;Squires et al., 2013). Characterizing metapopulation dynamics (e.g. occurrence, colonization, persistence and turnover) is therefore crucial to understanding the importance of populations at range peripheries and the spatiotemporal factors affecting range expansion (i.e. colonization) and contraction (i.e. failure to persist).
Occupancy models, and more specifically, dynamic occupancy models provide a flexible framework to investigate spatiotemporal dynamics in species occurrence, colonization and persistence using detection/non-detection data that are often less intensive to collect relative to count or mark-recapture data (Kéry et al., 2013;MacKenzie, Nichols, Hines, Knutson, & Franklin, 2003;Royle & Kéry, 2007). Dynamic occupancy methods extend single-season approaches by directly modelling the colonization and persistence processes underlying seasonal variation in occupancy, creating an efficient approach to assess species distributional changes, metapopulation dynamics (colonization, extinction and persistence) and spatiotemporal factors associated with these processes (Frey, Strong, & McFarland, 2012;MacKenzie et al., 2003;Yackulic, Nichols, Reid, & Der, 2015). Metrics derived from dynamic occupancy models, such as turnover probabilities and identification of high use areas, are especially important in understanding species dynamics at the periphery of their range. For example, sites near a range periphery may display high turnover rates due to low persistence but high colonization probabilities due to immigration from a core area (Kéry et al., 2013). Conversely, multiple processes may lead to low turnover rates, including sites that are no longer within a species range (e.g. sites that are never re-colonized) or sites that are consistently used (e.g. "core-use areas" that always persist). A primary benefit of dynamic occupancy models is the ability to disentangle processes leading to low turnover probabilities and subsequently identify sites that are robust to changes across multiple years and important to local persistence (Donaldson, Bennie, Wilson, & Maclean, 2019).
Identification of core-use areas (i.e. areas of intense use) is a recurrent objective in many habitat, resource selection and spaceuse studies as their identification may form the basis for targeted conservation efforts (Burdett, Moen, Niemi, & Mech, 2007;Manly, McDonald, Thomas, McDonald, & Erickson, 2002;Powell, 2000;Vander Wal & Rodgers, 2012). Telemetry data are commonly used to identify core-use areas, specifically 2nd or 3rd order resource selection of marked individuals (e.g. individual home range location or use of habitat components within a home range; Johnson, 1980;Manly et al., 2002). Telemetry data provide a wealth of information on individuals; however, telemetry studies are often expensive, invasive and difficult to implement in studies of low-density species of conservation concern (Burdett et al., 2007). Alternatively, occupancy methods consider the plot as a sampling unit and provide inferences to landscape-level processes and 1st order resource selection (geographic range; Johnson, 1980) without the need for marking individuals (Manly et al., 2002;study design I). Incorporating concepts of core-use areas with multi-year dynamic occupancy models provides a conceptual approach to identify multi-year, landscape-scale core-use areas based on population-level dynamics. Quantifying multi-year core-use areas in this context provides a useful tool to identify sites vital to conservation planning and prioritization due to population-level resiliency and consistent use of specific sites across multiple years.  (Koen, Bowman, & Wilson, 2015;McKelvey, 2000;Murray, Steury, & Roth, 2008;Schwartz, Mills, McKelvey, Ruggiero, & Allendorf, 2002). While lynx populations at their southern population-level mechanisms driving species persistence and key areas for conservation protection.

K E Y W O R D S
Canada lynx, core area, detection probability, dynamic occupancy, range dynamics, site occupancy, turnover periphery are of high management and conservation importance, successful monitoring in these regions is notoriously difficult due to the challenges of monitoring low-density, cryptic and highly variable populations (Barber-Meyer, Ryan, Grosshuesch, Catton, & Malick-Wahls, 2018;Murray et al., 2008;Squires, McKelvey, & Ruggiero, 2004;Squires, Olson, Turner, DeCesare, & Kolbe, 2012).
We developed a snow track survey protocol and dynamic occupancy modelling framework to investigate lynx occupancy, colonization, persistence and to identify multi-year landscape-level core-use areas. More broadly, our framework to identify multi-year core-use areas using detection/non-detection data provides a flexible and transparent approach for prioritizing key sites for conservation protection. We collected lynx detection/non-detection data from 2014 to 2019 in Superior National Forest and designated critical habitat in north-eastern Minnesota, USA. Our study area is at the southern boundary of the lynx range and covers one of the six geographic regions supporting lynx populations in the contiguous US (USFWS, 2017). Lynx occupancy and abundance are highly variable in this region, which highlights the complexities of developing efficient monitoring protocols to identifying spatial patterns and core-use areas across broad spatial and temporal scales. As lynx were recently recommended for delisting from the ESA, it is important to improve our understanding of lynx habitat use and metapopulation dynamics particularly in north-eastern Minnesota where fewer protections and predicted habitat change may increase the rate of population decline (US Fish and Wildlife Service, 2017).

| Study site and data collection
The study area encompassed 22,100 km 2 covering the North-eastern  (Barber-Meyer et al., 2018;Moen, Burdett, & Niemi, 2008) and other regions where lynx are less likely to occur (e.g. Voyageurs National Park [Moen et al., 2008] Figure 1). Although not stratified by habitat type, observers specifically surveyed a range of habitat conditions and regions where expected lynx occurrence ranged from low to high ( Figure 1; Supplement 1). We note, however, that sites in roadless areas of the Boundary Waters Canoe Area Wilderness (north-central region of the study area) were underrepresented in this study ( Figure 1).
For analysis purposes, we overlaid a 5 × 5 km grid across the study area resulting in 884 25-km 2 cells ( Figure 1). Grid cells were selected to match other ongoing surveys in Superior National Forest, specifically we followed grid cell definitions of the North American Bat Monitoring Program (NABat, Loeb et al. 2015), which roughly corresponded to half the size of female lynx winter home ranges (95% fixed kernel estimate = 44 km 2 ; Burdett et al., 2007). Each grid cell was considered a "site" and encounter histories were generated for each cell based on detection (1) or non-detection (0)

| Dynamic occupancy modelling
We parameterized the dynamic occupancy model as a Bayesian state-space model with a submodel for the latent occupancy process and second submodel for the detection process conditional on the occupancy state (Royle & Kéry, 2007). We chose this formulation to more easily derive associated metrics such as turnover and core-use area. Following notation of other state-space dynamic occupancy models (e.g. Broms, Hooten, Johnson, Altwegg, & Conquest, 2016;Royle & Kéry, 2007), we let y itj denote the detection (1) or non-detection (0) at site i ∈ 1,2, … ,N sites, in year t ∈ 1,2, … ,T years, during survey j ∈1, 2, …, J it surveys at site i in year t and z it denotes the true occupancy status at site i in year t. The number of surveys at site i in year t, J it , can vary and may be zero for some sites and years (Royle & Kéry, 2007).
Dynamic occupancy models estimate four processes: (a) initial occupancy ( i1 ), (b) colonization ( it ; i.e. probability an unoccupied site becomes occupied), (c) persistence ( it ; i.e. probability an occupied site remains occupied), and (d) detection probability (p itj ; Table 2). We consider each winter a primary period with secondary periods consisting of repeated surveys within a winter. Between primary periods, local occupancy status may change as a function of site-specific colonization and persistence probabilities. The occupancy state is assumed constant within a primary period. Due to the likely violation of the geographic closure assumption within a winter sampling period, occupancy probability should be interpreted as the probability that lynx used a site during the winter sampling period (MacKenzie et al., 2018).
We investigated four spatial habitat covariates that we expected to affect lynx occupancy, persistence and colonization (Aubry, Koehler, & Squires, 2000;Fuller, Harrison, & Vashon, 2007; Hoving  Figure S1. Submodels for persistence and colonization also included fixed year effects to address annual variation in these processes (Table 2). Overall, this formulation allows initial occupancy to spatially vary as a function of habitat characteristics and persistence and colonization probabilities to vary both spatially and temporally. While more complex spatiotemporal interactions likely influence lynx dynamics in this region, we limited our modelling approach to these primary covariates of interest due to the large number of parameters, relatively short duration of our study (5 years), and the limitation of having habitat covariate data that only varied spatially but not temporally ( Figure S1).  Observers also recorded four covariates expected to affect In this way, the proportion of the N sites that are occupied is simply defined as the set of sites that were used by at least one lynx during the winter.
To investigate spatial patterns in turnover probabilities and identify multi-year core-use areas we derived two additional metrics.
First, site-and year-specific turnover probabilities were calculated as the probability a site changed occupancy status between years ( it ; MacKenzie et al., 2018, p. 362).
Site-specific multi-year average turnover probabilities ( i ) were then monitored as the average site-specific turnover probability across all years. Second, we derived the number of years a site was occupied interpreted as the probability a site qualified as a core-use area. We identified sites with w i > 0.50 as probable core-use areas. This approach to identify multi-year core-use areas is easily adapted across study systems, while providing a transparent, replicable definition of core-use areas from explicit metapopulation dynamics of occupancy, colonization and persistence.

| Implementation
We fit the dynamic occupancy model in a Bayesian framework using JAGS (Plummer, 2003) and the jagsUI package (Kellner, 2015) accessed through R version 3.5.1 (R Development Core Team, 2018).
We used Normal(0,10) priors for all logit-scale regression parameters and Beta(1,1) priors for the back-transformed intercept parameters (Table 2). A sum to zero constraint was enforced for all categorical regression parameters, effort was log-transformed, and all other covariates were standardized to have mean zero and variance of one.
We ran six parallel MCMC chains. Each chain contained 5,000 adaptation iterations, followed by 5,000 burn-in iterations, and 50,000 posterior iterations thinned by 5. Chain convergence was visually  (Gelman et al., 2013). We report results as posterior medians and 2.5th and 97.5th quantiles (95% CRI).

| RE SULTS
Surveys occurred each winter during 2014-2019 and covered >17,000 km across 431 of the 884 sites (Table 1; Figure 1). The number of sites surveyed per year ranged from 242 to 304 and the number of sites with detections ranged from 39 to 74 (Table 1) We identified 55 sites classified as probable core-use areas in the eastern and north-eastern region of the study area ( Figure 5). Coreuse areas were relatively contiguous in the eastern region, consisting of sites with both high average snowfall and high per cent conifer forest ( Figure 5). While turnover probabilities were also near zero for sites in the south and west, this region was dominated by consistently low occupancy probabilities (Figures 4 and 5).
Detection probability was positively related to survey effort and date ( Figure 6). Detection probability declined as snow condition deteriorated, however differences between Good and Fair snow conditions were not significant ( Figure 6). Similarly, detection probability was higher in sites with Good road conditions ( Figure 6). In general, detection probability remained <0.30 when snow conditions were Poor but exceeded 0.50 under Good snow conditions and substantial effort (e.g. >15 km surveyed per grid cell; Figure 6).

| D ISCUSS I ON
Lynx populations in the contiguous US are expected to decrease in abundance and become more disjunct due to climate-driven losses in habitat quality and quantity (US Fish and Wildlife Service, 2017).
Understanding the relationships between habitat and demography is essential for species management and conservation, particularly at the periphery of a species range. We found strong support for relationships among local habitat characteristics and lynx occupancy, multi-year core-use areas from detection/non-detection data provides a useful tool for conservation planning and management without the need for capturing or handling individuals. Additionally, the focus of these methods on 1st order landscape-level resource selection is crucial for large-scale conservation planning that aims to identify areas where populations are highly resilient, robust to change and disproportionately contribute to regional persistence (Donaldson et al., 2019).
Landscape characteristics are known to affect lynx abundance, survival and reproduction with lynx predominantly occurring in young coniferous forests that support high snowshoe hare abundances (Aubry et al. 2000;Fuller et al., 2007;Hoving et al., 2005;Vashon et al., 2008). Lynx are also more abundant in areas of deep and persistent snow (Hoving et al., 2005), where they are believed to outcompete other terrestrial predators for access to hare populations (Carroll, 2007). We found that average snowfall was positively related to initial occupancy, persistence and colonization of lynx.
While parameter estimates describing habitat relationships are not directly comparable between occupancy and telemetry studies, it is reasonable to expect similar patterns between the two. Lynx-habitat relationships at the southern periphery based on telemetry data have shown increased lynx presence in areas with more coniferous and regenerating young forests (McCann & Moen, 2011;Moen et al., 2008), which is consistent with our results that lynx were more likely to occupy and persist in areas with increasing per cent conifer forests. Occupancy and telemetry data also support similar landscape-level trends associating lynx with mid-successional forests, while mature forests are used but may not be preferred (Bayne, Boutin, & Moses, 2008;Gooliaff & Hodges, 2018;Holbrook et al., 2019;Hoving et al., 2005;Murray et al., 2008).
Snow conditions are commonly suggested as a factor affecting lynx range limits especially at the southern periphery (Gooliaff & Hodges, 2018;Hoving et al., 2005). In this study, areas with higher average snowfall were associated with higher occupancy, colonization and persistence probabilities. Additionally, our dynamic occupancy approach separated the effects of broad-scale snowfall patterns on lynx  , 2018). Our ability to tease apart the interacting effects of snow conditions and forest characteristics was limited by the relatively short duration of our study (5 years) and limited data on regional snowfall. Understanding if forest management practices may buffer lynx from the possible negative effects of climate change (e.g. reduced snow cover) would be particularly useful for lynx conservation efforts and will hopefully be explored by future studies with access to longer time series.
Our study demonstrates an efficient multi-year, landscape-scale monitoring programme for a species of conservation concern, a perpetual challenge for a variety of programmes, especially those on federal lands (Noon, Bailey, Sisk, & Mckelvey, 2012). Specific to our lynx study system, future monitoring could benefit from increased spatial coverage of surveys (including roadless areas in the Boundary Waters Canoe Area Wilderness), stratified sampling designs and exploration of spatial autocorrelation (Bled, Nichols, & Altwegg, 2013;MacKenzie et al., 2018).
While dynamic occupancy methods are applicable to any number of species, spatially and temporally replicated track surveys are a particularly promising approach to monitor species across broad spatial and temporal scales since they do not require physical capture or marking, utilize indirect observations (e.g. tracks, scat) and are cost-effective across large geographic areas and multiple years (Karanth et al., 2011;Louvrier et al., 2019;Whittington, Heuer, Hunt, Hebblewhite, & Lukacs, 2014). Additionally, our incorporation of landscape-scale resource selection and core-use area concepts with dynamic occupancy models provides a framework to identify population-level mechanisms driving species persistence across heterogenous landscapes. These methods improve our capacity to assess species distributional changes, quantify population-level spatial dynamics and enhance the contribution of monitoring programmes to inform conservation planning and prioritization.

ACK N OWLED G EM ENTS
We are grateful to the numerous surveyors through the years espe-

DATA AVA I L A B I L I T Y S TAT E M E N T
Canada lynx detection/non-detection data, covariates and R script to run the dynamic occupancy model are available at https://doi. org/10.5061/dryad.66t1g 1jzq.

B I OS K E TCH
The shared interests of our research team include the application of field studies and quantitative methods to better understand spatial patterns in wildlife demography that inform management decisions and conservation efforts.