DOI: 10.5061/dryad.p8cz8w9vz
title: Density and genetic diversity of grizzly bears at the northern edge of their distribution

description: Species at the periphery of their range are typically limited in density
 by lower habitat quality. As a result, the Central-Marginal Hypothesis
 (CMH) predicts a decline in genetic diversity of populations towards the
 periphery of a species’ range. Grizzly bears (Ursus arctos) once ranged
 throughout most of North America but have been extirpated from nearly half
 of their former range, mainly in the south. They are considered a species
 at risk even in Canada’s remote North, where they occupy the northernmost
 edge of the species’ continental distribution in a low-productivity tundra
 environment. With climate change, one of their main food items in the
 tundra (caribou), which has always shown yearly fluctuations, is
 declining, but simultaneously, grizzlies appear to be expanding their
 range northward, in tundra environment. Yet, a lack of population density
 estimates across the North is hindering effective conservation action. The
 CMH has implications for the viability of peripheral populations, and the
 links between population fluctuations, potential bottlenecks and genetic
 diversity need to be determined to contribute to species’ conservation.
 Using non-invasive genetic sampling from 2012 to 2014, and autosomal DNA
 genotyping (via-microsatellites), we estimated bear density using a
 spatial capture-recapture framework and analysed genetic diversity using
 observed heterozygosity (Ho), Allelic Richness (AR), and expected
 heterozygosity (He). We compared our findings to other studies that used
 comparable methodologies on this and a related species (Black bears; Ursus
 americanus). We found densities of grizzly bears that were low for the
 species but characteristic for the region (5.9 ± 0.4 bears/1000 km2), but
 with high Ho (0.81 ± 0.05), AR (7 ± 0.78) and He (0.71 ± 0.03), despite a
 signal of recent bottlenecks. In both species, peripherality was not
 correlated with Ho but was negatively correlated with density. We suggest
 that the apparent growth of this expanding population of grizzlies offsets
 the negative impacts of recent bottlenecks on Ho. Indigenous Knowledge
 provides historical context (on the order of centuries – e.g., arctic
 large mammal fluctuations, grizzly bear bottlenecks) for the current bear
 population dynamics (on the order of decades – e.g., climate change,
 northern grizzly bear expansion).
description: Study area The
 30,000km<sup>2</sup> study area, centred at N 64.2°, W 110.0°,
 was located in the Southern Arctic (Coppermine River Upland Ecoregion,
 CRU) and Taiga Shield (Takijuq Lake Upland Ecoregion, TLU) Ecozones
 (Figure 1) (Ecological Stratification Working Group Canada 1995). Its
 northern limit extended to the border of NT and Nunavut, and its southern
 limit to the tree line. The mean annual temperature of the CRU ecoregion
 was -7.5 ˚C, and -10.5 ˚C for the TLU ecoregion (Ecological Stratification
 Working Group Canada 1995). Average temperatures at the time of the study
 were likely higher than those cited above, due to a changing climate (Post
 et al. 2009). The tundra area is considered semi-arid, with mean annual
 precipitation ranging from 200-300 mm, with mostly continuous permafrost
 characterized by a rolling landscape of uplands, lowlands, and plateaus
 (Ecological Stratification Working Group Canada 1995). Eskers, created
 through glaciation and composed of stratified sand and gravel, are found
 throughout the landscape and form most of the relief in this area.
 Lowlands are generally a mosaic of sedge dominated wetlands composed of
 fens and bogs as well as lakes and rivers. Vegetative cover is typically
 dominated by heath and shrub species such as Labrador tea
 (<em>Rhododendron</em> spp.), dwarf birch (<em>Betula
 nana</em>), and willow (<em>Salix herbacea</em>).
 Mammals inhabiting the area include barren-ground caribou, moose
 (<em>Alces alces</em>) but at very low densities, grizzly
 bears, wolves (<em>Canis </em>spp), red and arctic foxes
 (<em>Vulpes </em>spp<em>.</em>), and small
 rodents. Grizzly bear sampling approach
 Non-invasive genetic sampling of
 bears can be conducted using natural rub sites (e.g., trees), or
 human-made objects usually placed within predetermined grid cells
 delineated prior to detector deployment (Woods et al. 1999, Karamanlidis et al. 2010,
 Dumond et al. 2015, Boulanger et al. 2018). The spatial organization of a
 detector array for use in capture-recapture studies depends on the
 movement ecology of species (Sollmann et al. 2012, Royle et al. 2014). We
 divided the study area into 221 square grid cells of 144
 km<sup>2</sup> each. This area was based on a rough
 approximation of 14-day home ranges of barren ground female grizzly bears,
 which was based on previous studies (McLoughlin et al. 1999) (Fig. 1).
 Owing to the lack of natural rub sites in the tundra, we constructed hair
 snare posts from wooden boards wrapped in barbed wire and fastened
 together into a tripod shape so that posts could be transported by
 helicopter and deployed in the field (Dumond et al. 2015) (Figure 2). Hair
 snares were deployed near the centre of each grid cell, while avoiding
 locations on lakes or waterbodies, which also influenced logistics of
 access. They were baited with a non-reward scent lure corresponding to the
 seasonal availability of food sources to attract grizzly bears, as
 identified during a Traditional Ecological Knowledge (TEK) workshop held
 in the community of Lutsel K’e, NT in 2015 (Jessen 2017), and based on
 previous studies (Gau et al. 2002b). In early and late summer, rotten
 cow’s blood and fish oil were used while bergamot, raspberry, and
 cranberry oil were used in the mid-summer (Supplementary Table
 S1). We conducted six
 sampling occasions (rounds of site visits) per session (year), each
 lasting 10-14 days to minimize hair sample degradation from weather
 exposure, which increases with time (Dumond et al. 2015)(Supplementary
 Methods; Supplementary Table S1). Hair samples, i.e., clumps of hair
 captured by a single barb, were collected during each visit, and placed in
 labelled paper envelopes which were stored in a cool, dry place. Sampling
 of the northern half of the study area was carried out from mid-June to
 mid-September of 2012 and 2013, and sampling of the southern half of the
 study area was identically conducted from mid-June to mid-September of
 2013 and 2014 (Fig. 1). Laboratory and statistical
 analyses of genetic data A
 high-quality set of hair (either ≥ 30 underfur or ≥ 2 guard hair roots)
 was chosen from each hair sample and analyzed for species confirmation,
 sex, and genotype using established techniques, including established
 genotyping error-checking protocols (Paetkau 2003). We used a ZFX/ZFY
 gender marker, plus 8 microsatellite markers (G10B, CXX110, G1D, G10H,
 G10J, G10M, G10P, MU59) (Paetkau et al. 1999). Using 5 loci is considered
 sufficient for accurately detecting individuals in brown bears (Waits et
 al. 2001), and typically 7-8 loci are used in studies of brown bear
 populations (Boulanger et al. 2001,
 Dumond et al. 2015).
 All genetic analyses were conducted by Wildlife Genetics
 International (WGI) in Nelson, British Columbia. Scoring errors and the presence of null
 alleles were detected using MICROCHECKER v2.2.3 (Van Oosterhout et al.
 2004). We also tested markers for deviations from Hardy-Weinberg
 Equilibrium (HWE) and for Linkage Disequilibrium (LD), using the exact
 probability test in Genepop v4.2 (Rousset 2008) (Supplemental Methods S2).
 Allelic richness (AR), H<sub>o</sub>, and Nei’s unbiased
 expected heterozygosity (H<sub>e</sub>; (Nei and Roychoudhury
 1974)) were used as measures of genetic diversity. Both were used because
 allelic richness measures the number of alleles in a population
 standardized by sample size and is a measure of the raw amount of
 variation at loci, while expected heterozygosity accounts for both the
 number of alleles and the evenness of allele frequencies. AR was
 calculated using the rarefaction method implemented in FSTAT v2.9.4
 (Goudet 2003), and He was obtained using Genetix v4.05.2 (Belkhir et al.
 2004). We used the program NeEstimator to estimate
 effective population size (Ne) with its
 linkage disequilibrium single‐sample estimator without lowest allele frequency
 restriction (Do et al.
 2014). We
 also tested for heterozygosity excesses and signs of a genetic bottleneck
 using the program Bottleneck (Cornuet and Luikart 1996, Cristescu et al.
 2010). Its infinite allele model (IAM), two-phase model (TPM), and
 stepwise mutation model (SMM) were all applied. It should be noted that
 all are legitimate, although imperfect models of mutations, and
 simulations results are inconclusive on which is the most appropriate for
 autosomal microsatellite data (i.e., our data) in particular (Shriver et
 al. 1993). Ultimately these models, for which no consensus currently
 exists on which is best, allow assessing genetic signatures of population
 bottlenecks or of population expansions; however, these cannot be
 distinguished from genetic signatures resulting from natural fluctuations
 in density. Estimating density using spatial
 capture-recapture We
 developed spatial capture-recapture encounter histories of individual
 bears based on the hair samples (Woods et al. 1999, Mowat and Strobeck
 2000, Boulanger et al. 2018). We used spatial capture-recapture (SCR)
 models to estimate density. SCR models join an observation model that
 describes the decreasing probability of detecting an individual as a
 function of the distance between a detector location and an animal’s home
 range centre, to a spatial point process model that describes the
 distribution of animal home range centres on a landscape (Efford 2004,
 Royle and Young 2008, Efford and Fewster 2013). We estimated bear density in the study area
 for the years 2012, 2013, and 2014 separately, using the
 <em>secr</em> package version 4.5.5 (Efford 2022) in the program R,
 Version 4.1.3 (R Core Team
 2022). To account for peripheral individuals with home range
 centres outside the detector array, we applied a 48 km buffer zone,
 approximately double the greatest Root Pooled Spatial Variance for males
 or females for 2012, 2013, and 2014 (Eﬀord, 2019, Efford, 2004). We used
 the mask.check command in <em>secr</em> to ensure that the
 buffer size was appropriate. The study area and the buffer zone summed
 together made up an area of approximately 85,000
 km<sup>2</sup>. Models were fitted using the
 <em>secr.fit</em> command that maximizes model likelihood
 through integration over the unknown locations of individuals’ home range
 centres. We used a halfnormal Gaussian encounter probability function in
 all density estimates. We classified the hair snares as “proximity”
 detectors that allow for multiple, independent animal detections at a
 single detector in each sampling occasion, but only one detection per
 individual per sampling occasion. SCR models explicitly account for the
 distance and spatial layout of the sampling stations, even if not fully
 evenly spaced (Royle and Young 2008, Efford and Fewster 2013). We excluded
 water bodies from usable habitat using the polygon “inland lakes and
 rivers” from the Canadian 2011 Census boundary files (Open Canada 2022)
 and the habitat mask feature from the <em>secr</em>
 package. For each year, we
 conducted model selection using Akaike’s Information Criterion corrected
 for small samples (AICc), to find the most parsimonious fitting model
 (Anderson and Burnham 2004). Densities of grizzly bears might differ by
 sex and SCR models can account for sex either by conducting analyses
 separately for females and males or by conducting analyses for all animals
 pooled but including sex-specific covariates. To determine the best model
 structure for density estimation, we proceeded in three steps.
 (1) To determine which parameters
 should be sex-specific, we first compared 8 “sex-specific” models with all
 combinations of the following covariates: no sex covariates, sex-specific
 densities, sex-specific baseline detection probabilities, and sex-specific
 sigmas (Supplementary Table S2). Sigma is the spatial scale parameter
 (also called the movement parameter) in SCR models, and in terrestrial
 wildlife species is related to the radius (r) of 95% home range estimates
 (e.g., from telemetry studies) through the equation sigma = r/2.45 (Sun et
 al. 2014). (2) We then
 compared 12 “pooled” models without sex specific covariates with
 combinations of the following detection probability covariates: null model
 (no covariates), individual local behavioural response (bk), individual
 global behavioural response (Bk), linear time (T), factorial time (t), and
 lure type (lure). We tested for behavioral responses because grizzly bears
 tend to revisit rub trees and other hair snag sites that they previously
 visited (Morehouse and Boyce 2016, Lamb et al. 2018). We did not fit
 models that included both lure type and time covariates because lure type
 was partly correlated with time (Supplementary Table S3). (3) To determine the model to be used for
 density estimation, we combined the best sex-specific models (dAICc
 &lt; 2) with the best pooled models (dAICc &lt; 2), also including
 interactions of detection covariates with sex, to arrive at a suite of
 competing “combined” models. We then used the best combined model for
 parameter estimation.