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.