2 Materials and Methods

2.1 Fecal pellet collection and genetic analysis

For each population, we flew 3 surveys to collect fecal pellets during winter (December to March), with sampling occasions spaced approximately one-month apart. Following the aerial survey protocol outlined in Hettinga et al. (2012), aerial transects were systematically flown at 3-km intervals across each entire caribou population range using rotary- or fixed-wing aircraft, or a combination of both aircraft, to locate caribou feeding locations, for a total of 69,070 km flown across the seven ranges (Table 1). Once located, personnel landed at each feeding site and collected fecal samples; this included collecting samples from backtracking on caribou trails. All pellet samples were kept frozen at -20°C until DNA extraction was performed.
In the laboratory, fecal samples were thawed and the mucosal coat surrounding the pellets was removed for DNA analysis. The extraction protocol used to amplify the DNA is outlined in Ball et al. (2007). Following quantification of target caribou DNA, samples were diluted down to a working stock concentration of 2.5ng/ul. We amplified the DNA at 9 variable fluorescently labelled microsatellite loci (FCB193, RT7, RT1, NVHRT16, BM888, RT5, RT24, RT6, OHEQ; Bishop et al., 1994; Cronin, MacNeil, & Patton, 2005; Wilson, Strobeck, Wu, & Coffin, 1997) to generate individual-specific genetic profiles, along with caribou-specific Zfx/Zfy primers for sex identification. The amplification protocol is outlined in Ball et al. (2007). Following amplification, each sample was genotyped on the ABI 3730 DNA Analyzer (Applied Biosystems) Microsatellite alleles were scored with the program GeneMarker v1.91® (SoftGenetics, State College, PA) and followed a protocol documented in Flasko et al. (2017) and McFarlane et al. (2018). Unique individuals were identified using the program ALLELEMATCH (Galpern, Manseau, Hettinga, Smith, & Wilson, 2012). We retained samples that amplified at ≥5 loci and re-amplified apparent unique genetic profiles represented by a single sample using two independent scorers to confirm unique individual identities (Hettinga et al., 2012). An error rate per locus was calculated using these re-amplification results.

2.2 Framework

2.2.1 Empirical SCR modeling

We used a maximum likelihood approach implemented in the R packagesecr (Efford, 2018; R Core Team, 2019) to estimate boreal caribou densities. SCR models are comprised of a submodel for the distribution of animals in the area of study (population density, D ), and a submodel for the detection process, given the detection probability (the intercept of the detection function, g0 ) and given a parameter for\(\ \)scaling the detection function (the spatial extent of an individual’s use of the landscape - \(\sigma\); Borchers & Efford, 2008; Efford et al., 2009). For our empirical data, we treated each survey as an occasion within a single session. We discretized the study area into a 1500 m grid of proximity detectors, (which record the presence of individuals at each detector without restricting movement; Efford et al., 2009). The area of integration for SCR models needs to be large enough such that animals residing beyond the study area have a negligible chance of being detected (Borchers & Efford, 2008; Efford, 2004; Royle & Young, 2008). We therefore defined our state-space with a 15 km2 buffer around all study areas. We ran models for females, males, and both males and females together.
We estimated the parameters of the SCR detection function (g0 and\(\sigma\)) by maximizing the conditional likelihood, and derived density (D ) from the top AICC-ranked models (Anderson, Burnham, & White, 1994; Borchers & Efford, 2008). We used the hazard exponential form of the detection function, because area search data models the cumulative hazard of detection (Efford, 2011). Models assumed that individuals were identified correctly, populations were demographically closed during sampling, and detections were independent and conditional on activity center (Borchers & Efford, 2008; Efford, 2004). We assessed sources of variation on the detection parameters with time and behaviour effects on both g0 and\(\sigma\).

2.2.2 Testing assumptions of homogeneous distribution

Boreal caribou is a non-migratory ecotype of caribou and have relatively small home ranges compared to wide-ranging carnivores such as brown bears (Graham & Stenhouse, 2014; Lamb et al., 2018) and black bears (Whittington & Sawaya, 2015). Boreal caribou group size fluctuates throughout the year; group size is lowest during spring and summer when cows become solitary for calving, increases before the rut, and peaks in early to late winter (Thomas & Gray, 2002). To assess how the distribution of the animals (i.e. clustering) affected the precision and relative bias of our estimates, we simulated different population distributions using three of our empirical datasets (Little Smoky, Cold Lake, and Slave Lake). Different distributions can be used for the simulations including a homogeneous Poisson distribution, inhomogeneous or clustered Poisson distributions (Efford, 2019a). The chosen population distribution should reflect the distribution of the study species. Our empirical data approximated a Neyman-Scott clustered Poisson distribution which was then used for the simulations (Efford, 2019a). To simulate multiple detections in very close proximity, we set the spatial scale (\(\sigma\)) of the 2-D kernel for locations within each cluster to be 1. To simulate varying levels of clustering, we varied the fixed number of individuals per cluster (see Figs S1.1-S1.3). We selected starting values for D, g0 and \(\sigma\) from the empirical model runs (Table 2). We carried out all simulations in the secr R package (Efford, 2018; R Core Team, 2019).

2.2.3 Assessing precision and relative bias of different sampling designs using empirical data

We repeated the empirical population analyses with subsamples of data to explore how reduced sampling intensity affected the relative bias and precision of the density estimates from our empirical study. We rarified the data by reducing the number of sampling occasions and reducing the number of aerial transects flown. For the reduced number of sampling occasions, all possible 2-occasion combinations were run (occasions 1 and 2; occasions 2 and 3; and occasions 1 and 3). Aerial transects were removed from the original spatial field data, keeping either every second or third transect line to emulate sampling strategies of 6 km or 9 km transects. Only the samples collected along the remaining transect lines were retained, and only those detectors along the remaining transect lines were used in the analysis. We used the coefficient of variation (CV) as the metric for precision, and calculated the relative bias (RB = ( \(\hat{D}\) - D)/D ) as the metric for bias (as in Tobler & Powell, 2013; Efford & Boulanger, 2019; Efford & Fewster, 2013; Kristensen & Kovach, 2018). We compared estimates from the reduced datasets (\(\hat{D}\)) to those based on the empirical dataset (D ). We considered models with CV <20% (following Pollock et al. 1990) and relative bias <15% (Otis, Burnham, White, & Anderson, 1978) as favourable outcomes. Models with CV <30% and RB <20% can also be considered favourable (Kristensen & Kovach, 2018), because high precision may be difficult to achieve for rare and low-density species.
We calculated the precision and relative bias of each subsampling scenario. Scenarios were considered precise when CV < 30% and RB < 20%. To determine how the number of captures, number of recaptures, and number of spatial recaptures (recaptures at different locations) influence the precision and relative bias of the estimates, we correlated the precision and relative bias of the estimates with these parameters for each scenario, and then globally.