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 (H_o), Allelic Richness (AR), and expected heterozygosity (H_e). 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 km²), but with high H_o (0.81 ± 0.05), AR (7 ± 0.78) and H_e (0.71 ± 0.03), despite a signal of recent bottlenecks. In both species, peripherality was not correlated with H_o 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 H_o. 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).

Study area

The 30,000km² 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 (Rhododendron spp.), dwarf birch (Betula nana), and willow (Salix herbacea). Mammals inhabiting the area include barren-ground caribou, moose (Alces alces) but at very low densities, grizzly bears, wolves (Canis spp), red and arctic foxes (Vulpes spp.), 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² 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_o, and Nei’s unbiased expected heterozygosity (H_e; (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 secr 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 (Efford, 2019, Efford, 2004). We used the mask.check command in secr 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². Models were fitted using the secr.fit 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 secr 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 < 2) with the best pooled models (dAICc < 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.