https://orcid.org/0000-0002-5381-8514
Abstract
In the westernmost Canadian province of British Columbia, fishers (Pekania pennanti) occur in low-elevation forested ecosystems in the central and northern portions of the province, with several substantial mountain ranges and rivers bisecting this distribution. The effect of these geographic features on gene flow and population structuring is unknown and may contribute to fishers differentiating into 2 or more populations within the province. To better understand population structuring within the range of fishers in British Columbia, we collected tissue samples from 491 fishers from throughout the province and evaluated variation of neutral genetic markers using a 385-base pair section of the mitochondrial DNA control region and 13 microsatellite loci. Fishers appeared to be genetically structured as hierarchical stepping-stone populations where an upper hierarchical level of genetic partitioning resulted in 2 populations at the provincial scale, with 3 subpopulations occurring in the central interior region. We observed gene flow of 2 to 5 migrants per generation between the 2 upper populations, which was substantially less than the 90 migrants per generation estimated as needed to maintain genetic homogeneity. This resulted in differentiated upper populations (haplotype FST = 0.272, microsatellite FST = 0.049) characterized by relatively low Ne (Columbian population = 252, 95% CI = 185 to 332; Boreal population [British Columbia only] = 136, 95% CI = 92 to 234). The consequence of low gene flow and resultant population structuring is that the Columbian population likely receives relatively little genetic and demographic support from other populations and, combined with precipitous recent declines in its numbers, puts this population at heightened conservation risk. As a result, swift and effective actions are needed to conserve this vulnerable population of fishers.
Understanding the genetic structure of wildlife populations is a key aspect of revealing past and present demographic and evolutionary processes and has important implications for management and conservation of wildlife (Allendorf et al. 2022). This is especially true for populations that are genetically isolated from each other, as they are prone to increased genetic drift that can lead to population differentiation. Small, isolated populations are also susceptible to inbreeding, which can compromise their ability to adapt to local environments or respond to changing conditions. Furthermore, populations with little connectivity to others may be susceptible to declines because they cannot rely upon migrants to offset reductions in population numbers caused by stochastic events or human-caused threats like habitat change. Fishers (Pekania pennanti) are a case in point.
Fishers are medium-sized carnivores in the mustelid family that occur throughout the Boreal, Carolinian, and Northwestern Montane forests of North America. Although relatively common and locally plentiful in their eastern range (e.g. Koen et al. 2007), fishers are generally much less abundant in the western portion of their distribution (e.g. Lofroth et al. 2010). In the westernmost Canadian province of British Columbia, fishers are at the northwestern limit of their range and an uncommon carnivore of low-elevation forests of central and northeastern portions of the province (British Columbia Conservation Data Centre 2015). Fisher populations in British Columbia, particularly in the central interior region of the province, are of considerable conservation concern due to declines in numbers of 30% to 50% over the past 20 years (BC Conservation Data Centre 2020), increasing habitat threats from logging and wildfire, and low survival and reproductive output (Lofroth et al. 2023). These threats and population limitations contribute to the prediction that fishers may be extirpated in the central interior region within as few as 11 years (Fogarty et al. 2022).
Habitat and distribution research conducted over the past 30 years suggests that fishers are not distributed contiguously throughout British Columbia, as previously thought (e.g. Hagmeier 1956; Proulx et al. 2004). The species occurs over approximately 500,000 km2 within the province, but this distribution is bisected by a >30-km-wide band of high-elevation forests and alpine ecosystems formed by the Rocky, Omineca, and Skeena Mountain ranges (Fig. 1). These mountainous areas are typified by deep, soft snow that hampers the locomotion of fishers (Raine 1983; Krohn et al. 1997), making movement during winter energetically unsustainable. Cavity-bearing trees needed by fishers for reproduction (Weir et al. 2012) are virtually nonexistent in these ecosystems as well. As a result of these factors, mountainous zones are generally unsuitable for fishers and the species is rarely found in these areas in British Columbia (Weir and Lara Almuedo 2010; Davis and Weir 2021). This band of unsuitable mountainous habitat may affect movement of individuals and reduce the exchange of genes within the province.
Distribution of fishers and locations of genetic samples from 491 fishers collected between 1997 and 2017 in central interior (dark gray) and boreal (light gray) regions of British Columbia, Canada. The majority of samples (96.3%) were collected between 2015 and 2017.
Reduced gene flow may contribute to an increased likelihood of genetic drift and population differentiation among fishers in British Columbia. Unlike other species in the Martes clade that show moderate levels of population structuring (e.g. wolverines Gulo gulo: Rico et al. 2016; American martens Martes americana: Kyle and Strobeck 2003, but see Lucid et al. 2020), Fisher populations show substantial substructure at even small distances (Kyle et al. 2001; Wisely et al. 2004). Because of this, minor geographic barriers that impede genetic exchange may exacerbate the tendency for structuring and facilitate population differentiation, particularly in areas where habitat is constrained (Aubry et al. 2004), such as is found in British Columbia.
If gene flow within the range of fishers in British Columbia was restricted for enough time, demes may have formed that experienced sufficient genetic drift to differentiate into distinct and isolated populations. Although fishers that occur in the boreal region of the province are likely part of a larger contiguous population that occurs across the boreal forests of Canada, fishers in the central interior region are separated by >800 km from populations in California and Oregon and >400 km from populations in Idaho and Montana (Lofroth et al. 2010). As a result, fishers in the central interior region are limited to receiving migrants and genes from the boreal region of the province and demographic support and rescue may be unlikely, thereby challenging the persistence of the species in the province.
We wanted to increase our understanding of the magnitude and effects of gene flow and population differentiation among fishers in British Columbia, given the potential restricted genetic connectivity, possible population differentiation, and conservation concern facing fishers in the province. Specifically, we wanted to examine whether patterns of gene flow and structuring within British Columbia were consistent with geographic features that may have affected movement of genes and better understand the level of connectivity and differentiation between fishers from the central interior region of the province and those in the boreal region. Our research used georeferenced tissue samples collected from throughout British Columbia to assess gene flow, genetic drift, and population structuring and to evaluate whether fishers occur as at least 2 discrete populations within the province. Information on genetic structure among these putative populations would be valuable not only to better understand Fisher populations in British Columbia, but also to provide better information for decision-makers upon which to evaluate the status and identify conservation or management actions needed for each population unit.
Materials and methods
Study area
The study area comprised all areas of Fisher occurrence within the province of British Columbia (Fig. 1), with fishers being distributed primarily in the central and north–northeastern portions of the province. The ‘central interior’ region has mixed plateau topography and is dominated by Sub-Boreal Spruce, Sub-Boreal Pine–Spruce, Montane Spruce, and Interior Douglas-Fir biogeoclimatic zones (MacKenzie and Meidinger 2018). The northern and northeastern ‘boreal’ region has similar topography and is dominated by the Boreal White and Black Spruce biogeoclimatic zone (MacKenzie and Meidinger 2018). These regions are separated by approximately 50 km of high-elevation mountain ranges (Weir and Lara Almuedo 2010) comprised of the Engelmann Spruce–Subalpine Fir, Spruce–Willow–Birch, and Boreal Altai Fescue Alpine biogeoclimatic zones (MacKenzie and Meidinger 2018), which are dominated by a deep (e.g. >1.5 m) and soft snowpack much of the year (DeLong 2004). Four narrow (<1.35-km-wide) relatively low-elevation forested passes connect Fisher distribution between the 2 regions.
Sample collection
We collected tissue samples from the hides of legally harvested fishers between 2014 and 2017. We also collected hair samples between 1997 and 2010 from fishers that were handled in the course of radiotelemetry studies or translocation projects (Weir and Corbould 2008; Weir et al. 2013; Lewis 2014), but most samples (96.7%) were collected during the 2015 and 2016 trapping seasons. All capture and handling methods met or exceeded capture and handling guidelines outlined in the protocols for Wildlife Capture and Handling (Resource Inventory Committee 1998), followed American Society of Mammalogists guidelines (Sikes et al. 2016), were approved by a member of the Canadian Council on Animal Care, and were carried out under Wildlife Act permits C076979 and FSJ05-9483.
DNA extraction, amplification, and sequencing
Laboratory work was performed at the US Forest Service Rocky Mountain Research Station, National Genomics Center for Wildlife and Fish Conservation. Genomic DNA was extracted from tissue and hair samples using QIAGEN DNeasy Blood and Tissue kits (Qiagen, Valencia, California) according to manufacturer’s instructions but using modifications for hair samples following Mills et al. (2000).
We amplified the left domain of the control region of mitochondrial DNA (mtDNA; 385 base pairs) using primers L15926 and H16498 (Kocher and Wilson 1991) for a sample of 491 fishers. Reaction volumes of 30 μL contained 50 to 100 ng DNA, 1× reaction buffer (Life Technologies, New York), 2.5 mM MgCl2, 200 μM of each dNTP, 1 μM of each primer, 1 U AmpliTaq Gold polymerase (Applied Biosystems, Foster City, California). The PCR program was 94 °C/5 min, [94 °C/1 min, 55 °C/1 min, 72 °C/1 min 30 s] × 34 cycles, 72 °C/5 min. The quality and quantity of template DNA was determined by 1.6% agarose gel electrophoresis. PCR products were purified using ExoSap-IT (Applied Biosystems, Foster City, California) according to manufacturer’s instructions, and amplicons were sequenced at Eurofins Genomics (Louisville, Kentucky) using standard Sanger sequencing protocols. Resulting DNA sequence data were viewed and aligned with the software Sequencher (Gene Codes Corp., Ann Arbor, Michigan). We deposited the derived sequences to GenBank (accession numbers ON100606 to ON100611).
We also genotyped samples from 467 individual fishers at 13 microsatellite loci that have been described in previous mustelid studies: Mp0059, Mp0144, Mp0175, Mp0197, Mp0200, Mp0247 (Jordan et al. 2007), Ma1 (Davis and Strobeck 1998), Mer022, Mvis020, Mvis072 (Flemming et al. 1999), Ggu101, Ggu216 (Duffy et al. 1998), and Lut733 (Dallas and Piertney 1998). The reaction volume (10 μL) contained 1.0 μL DNA, 1× reaction buffer (Applied Biosystems), 2.0 mM MgCl2, 200 mM of each dNTP, 1 mM reverse primer, 1 mM dye-labeled forward primer, 1.5 mg/mL BSA, and 1 U Taq polymerase (Applied Biosystems). The PCR profile was 94 °C/5 min for an initial denaturing step followed by 45 cycles of 94 °C/1 min, 55 °C/1 min, and 72 °C/30 s. DNA from hair samples was amplified using the multi-tube approach (Eggert et al. 2003). Products were visualized on a LI-COR DNA analyzer (LI-COR Biotechnology). Data were error-checked using the program Dropout (McKelvey and Schwartz 2005), GenAlEx v6.5 (Smouse and Peakall 2012), and MICRO-CHECKER (Van Oosterhout et al. 2004).
Replicate genotyping was performed with all hair DNA samples, with each sample PCR-amplified initially twice across the 13 variable microsatellite loci to screen for allele dropout, stutter artifacts, and false alleles (DeWoody et al. 2006). To minimize genotyping error, if any locus failed to amplify in either replicate or if there was a discrepancy PCR amplification, the sample was repeated twice more. If a genotype was confirmed by this repeat analysis, then it was retained. If a genotype failed again, the sample was assigned a missing score at the failed locus. We reamplified 5% of the tissue samples across loci. We removed from analysis any sample for which amplification failed at one-third of the loci (i.e. 5 loci) or more after attempts to reamplify a locus.
mtDNA analyses
We used DnaSP v4 (Rozas et al. 2003) to identify identical mtDNA sequences, collapse them to unique haplotypes, and calculate the number of variable sites and haplotype diversity. We included 28 homologous P. pennanti sequences from GenBank (accession numbers AB601566, AY143663 to AY143674, GU121228 to GU121237, HQ705176 to HQ705180) in our analysis. We used the R package ‘pegas’ (Paradis 2010) to calculate nucleotide and haplotype diversity for the entire data set and for each geographic region, population, and subpopulation. We used the haplotype data and Arlequin v3.5 (Excoffier and Lischer 2010) to calculate pairwise FST values and conduct analysis of molecular variance (AMOVA), which estimated the pattern of partitioning of variance between populations and within populations, for groupings derived from geographic regions, STRUCTURE clusters, and TESS clusters (see Population clustering, below). We also used the “pegas” package to produce a minimum spanning tree network using the ‘haploNet’ function.
Microsatellite loci analyses
We used a variety of approaches with georeferenced microsatellite data to evaluate the genetic structure of fishers within British Columbia. We used the R package ‘adegenet’ (Jombart 2008) to calculate locus-by-locus observed (Ho) and expected (He) heterozygosity, and we tested if loci met Hardy–Weinberg expectations using an exact permutation test. The per-locus number of alleles and allelic richness were calculated with FSTAT v2.9.3 (Goudet 1995). We also evaluated pairwise linkage disequilibrium (LD) between loci using a Markov chain method in GENEPOP 3.4 (Raymond and Rousset 1995), with 1,000 dememorization steps, 1,000 batches, and 1,000 subsequent iterations. Next, to assess patterns of isolation by distance, we performed a Mantel test that examined the correlation between pairwise matrices of genetic distance and Euclidean geographic distance. We used the “adegenet” function ‘mantel.randtest’ to assess the correlation between Nei’s genetic distance and Euclidean geographical distance for the entire data set and for central interior samples only, with 999 replications for the tests. We then used Arlequin v3.5 (Excoffier and Lischer 2010) to perform AMOVA, where we examined 2 putative geographic populations: those fishers from the central interior (n = 356) and boreal (n = 111) portions of the species range and assessed putative population differentiation using the fixation index metric (FST). Lastly, to further compare genetic similarity between samples and to explore putative population structure, we performed a principal component analysis (PCA) with the R package ‘adegenet’ and calculated the number of private alleles in each population.
Population clustering
We used the program STRUCTURE v2.3.4 (Pritchard et al 2000), which uses aspatial Bayesian assignment to infer population clusters via K-means clustering, to identify genetic clusters that reflected population divisions. We conducted analyses using values of K = 1 to 10, the admixed ancestry model, correlated allele frequencies, information on sampling location, and a burn-in of 100,000 followed by 1,000,000 Markov chain Monte Carlo (MCMC) iterations, with 10 runs for each value of K to estimate LnP(K). We grouped samples based on their location for the calculation of MedMedK, MaxMedK, MedMeanK, and MaxMeanK (sensu Puechmaille 2016, below). We used the Wildlife Management Unit (Ministry of Forests 2023) from which each sample was collected as prior location information to assist the clustering algorithm. The ‘prior location’ option in STRUCTURE is merely a method to group samples that came from the same general location (in our case, Wildlife Management Unit), but these groupings remained aspatial in the Bayesian clustering process. Because sampling among clusters was likely uneven and population structuring may have been hierarchical, we used the methods of Puechmaille (2016) in StructureSelector (Li and Liu 2018) to calculate ΔK and LnP(K), then set a threshold of 0.5 in the estimation of MedMedK, MaxMedK, MedMeanK, and MaxMeanK, which we used to infer the most likely value of K. We used CLUMPAK (Kopelman et al. 2015) to combine the multiple iterations for each K. We considered samples to be confidently assigned to a cluster at an average q-value among replicates of 0.7. We also conducted separate STRUCTURE analyses following the methods above for samples collected from the central interior and boreal region, respectively, to evaluate whether hierarchical structuring occurred in the microsatellite data set.
We used an additional Bayesian clustering algorithm implemented in the software TESS v.2.3.1 (Chen et al 2007), which incorporates spatial information to inform individual ancestry estimates. This method is useful in determining potential genetic barriers within continuous populations. We varied K between 2 and 10 and performed simulation runs with the BYM admixture model and default values, with a total of 100,000 sweeps, 10,000 as burn-in, and 10 runs for each K. To estimate the best K, the deviation information criterion (DIC) was averaged across runs for each K. We selected the point of greatest inflection of the curve of average DIC score across the range of K as the most likely value of K.
Effective population size
We estimated effective population size (Ne) using the LD (i.e. nonrandom allele associations) method as implemented in NeEstimator v2 (Do et al. 2014) with the minimum frequency for included alleles (i.e. Pcrit) equal to 0.05. We calculated Ne for fishers from the central interior and boreal regions separately to investigate how demographic processes varied across space.
Gene flow
We used several approaches to assess gene flow within the distribution of fishers in British Columbia. We first used a Bayesian approach implemented in BayesAss v3.0.5 (Wilson and Rannala 2003) to estimate gene flow, measured as the proportion of individuals per generation that are migrants derived from another population, between the putative central interior and boreal populations. BayesAss is based on a Bayesian MCMC method and attempts to estimate gene flow rates in the most recent generations. We confirmed that the MCMC mixing parameter values of migration rates, allele frequencies, and inbreeding coefficients obtained acceptance rates of between 20% and 60% for proposed changes. We performed 5 replicate runs with 1 × 107 MCMC iterations, a burn-in phase of 1 × 106 iterations, sampling every 1,000th iteration, and initializing each run with a different random number generator seed to compare the posterior mean parameter estimates and verify concordance. Finally, we used Migrate-n v3.6.4 (Beerli and Felsenstein 1999; Beerli 2009) to estimate long-term average rates of gene flow per generation under a coalescent theory model, jointly estimating the mutation-scaled migration rate M (m/μ) and the mutation-scaled effective population size θ (4Neμ). We estimated starting values of θ and M with an FST calculation as detailed in the software and specified the Brownian motion model under the maximum likelihood criteria. We performed 2 identical runs with different starting seeds to assess whether the chains sampled parameter space properly. Additional parameters included 10 short chains for 10,000 generations with a sampling increment of 20, and 3 long chains for 100,000 generations with a sampling increment of 20.
The minimum number of effective migrants needed to prevent genetic drift and maintain long-term genetic homogeneity between populations are affected by factors that reduce the effective population size relative to census population size (Vucetich and Waite 2000). Many of these factors are prevalent in Fisher populations in British Columbia, including an unbalanced sex ratio in adults (1.86 adult females:adult males; Weir 2003), a long generation time relative to the average age of adults (median generation time = 3 years, mean lifespan of adult females = 3.33 years [SD = 1.82, n = 80], mean lifespan of adult males = 3.05 years [SD = 1.29, n = 44]; Weir 2003), fluctuating population sizes among years (Weir and Corbould 2006; Weir et al. 2011), variable reproduction rate among adults (frequency of reproduction varied between 0% and 100% among 50 individuals; data from Lofroth et al. 2023), and a skewed age structure (75% of population are nonadults; Weir 2003). We used the empirical data listed above and followed the methods of Vucetich and Waite (2000) to estimate the number of migrants needed to maintain genetic homogeneity and inhibit genetic drift between regions.
Results
We collected 491 samples of fishers from throughout their range in British Columbia, including 372 samples from the central interior and 119 samples from the boreal regions of the province (Fig. 1). We successfully genotyped 467 individual fishers at 13 microsatellite loci and identified mtDNA haplotypes of 491 individuals. Direct-count repeat genotyping error rates per allele ranged from 0.00 to 0.022 across samples and loci.
mtDNA haplotypes
We observed 6 variable sites and 6 haplotypes (Fig. 2) within our mtDNA samples. All sequences were identical to haplotypes identified by Drew et al. (2003) and were labeled as such in Fig. 2 to provide context within the range of the species. Across all samples, nucleotide diversity was 0.0058 and haplotype diversity was 0.78, with the Drew-7 haplotype being most common (34%) and Drew-1 least common (2.5%).
Haplotype network for 491 fishers collected throughout British Columbia, Canada, between 1997 and 2017. Each line crossing the connection between haplotypes represents 1 substitution. Frequency of occurrence of each haplotype is indicated by radius of circle, with proportion of each haplotype from the central interior region shown in dark gray and boreal region in light gray.
We detected considerable differences in the distribution of haplotypes among regions (Fig. 3). Two haplotypes were widespread across the sampling area (Drew-6, Drew-7), 2 occurred predominately in the central interior region (Drew-4, Drew-9), 1 occurred in, or within 100 km of, the boreal region (Drew-11), with 1 haplotype found exclusively in the central interior region (Drew-1). We observed significant differentiation in haplotypes among geographic regions (AMOVA, P < 0.001), with an FST score of 0.272 between samples collected in the boreal and central interior regions of the province. Differentiation in haplotypes among the STRUCTURE K = 4 clusters was also significant (AMOVA, P < 0.001: FST scores; red–green = 0.178, red–orange = 0.355, red–blue = 0.207, green–orange = 0.175, green–blue = −0.006, orange–blue = 0.156; cluster colors from Fig. 5).
Spatial distribution of 6 haplotypes from 491 fishers collected throughout British Columbia, Canada, between 1997 and 2017. Dark gray parts of the distribution maps show the central interior portion of fisher range; light gray shows the boreal portion. Inset charts show frequency of each respective haplotype among central interior (dark gray) and boreal (light gray) regions.
PCA of genotypes derived from 13 microsatellites for 467 fishers collected from the central interior (CI) and boreal regions of British Columbia, Canada, between 1997 and 2017. The centroid of the point distribution for each region lied outside the 75th percentile inertia ellipse of the point distribution for the other region. Axis 1 explained 8.1% and axis 2 explained 6.2% of the genetic variation, respectively.
Microsatellite summary statistics
We calculated a suite of metrics using microsatellite data for samples collected from throughout the range of fishers in British Columbia. The average number of alleles per locus was 6 and ranged from 2 to 11, 6 loci deviated from Hardy–Weinberg equilibrium (HWE), and 50 of the 78 pairwise comparisons between loci exhibited significant LD. When examined by region, these discrepancies decreased considerably, as we detected 4 loci deviating from HWE and significant LD in only 9 of 78 pairwise comparisons of loci from samples collected in the boreal region. Likewise, we documented 4 loci deviating from HWE and 16 instances of significant LD from central interior samples collected north and west of the Fraser River and 4 loci deviating from HWE and 9 instances of significant LD occurred among samples collected east of the Fraser River. None of the 13 loci that we examined showed consistent deviation from HWE or LD among these 3 geographic areas. Across all samples, the average observed heterozygosity was 0.59 and average expected heterozygosity was 0.64, with the overall inbreeding coefficient (FIS) equaling 0.079 (Table 1).
Genetic diversity measures of 497 haplotypes and 467 genotypes derived from 13 microsatellite loci from samples collected from fishers throughout British Columbia, Canada, between 1997 and 2017, including number of individuals genotyped (n), haplotype diversity (h), nucleotide diversity (π), mean rarefied allelic richness (AR), mean observed heterozygosity (HO), mean expected heterozygosity (HE), and mean inbreeding coefficient (FIS; 95% CI). STRUCTURE and TESS groupings were clusters to which each individual was assigned (color-shape refers to symbols in Figs. 5 and 6), whereas geographic groupings were based on the region from which samples were collected.
| Grouping | n | h | π | AR | HO | HE | FIS |
|---|---|---|---|---|---|---|---|
| All samples | 467 | 0.78 | 0.0058 | 6.36 | 0.59 | 0.64 | 0.079 (0.051 to 0.113) |
| Geographic region | |||||||
| Central interior | 356 | 0.76 | 0.0048 | 5.98 | 0.59 | 0.62 | 0.058 (0.028 to 0.093) |
| Boreal | 111 | 0.66 | 0.0061 | 5.66 | 0.60 | 0.64 | 0.076 (0.028 to 0.15) |
| STRUCTURE K = 4 | |||||||
| Red-circle | 127 | 0.69 | 0.0064 | 5.67 | 0.56 | 0.64 | 0.129 (0.04 to 0.291) |
| Green-square | 120 | 0.71 | 0.0045 | 4.98 | 0.54 | 0.59 | 0.084 (−0.007 to 0.215) |
| Orange-star | 156 | 0.72 | 0.0038 | 5.52 | 0.55 | 0.60 | 0.096 (0.006 to 0.255) |
| Blue-triangle | 64 | 0.71 | 0.0043 | 4.80 | 0.53 | 0.59 | 0.107 (0.017 to 0.293) |
| TESS K = 4 | |||||||
| Red | 126 | 0.71 | 0.0065 | 5.68 | 0.60 | 0.65 | 0.074 (0.03 to 0.128) |
| Green | 134 | 0.73 | 0.0047 | 5.06 | 0.58 | 0.60 | 0.039 (−0.01 to 0.093) |
| Orange | 131 | 0.72 | 0.0038 | 5.59 | 0.60 | 0.61 | 0.017 (−0.015 to 0.051) |
| Blue | 76 | 0.72 | 0.0043 | 4.94 | 0.58 | 0.59 | 0.037 (−0.008 to 0.091) |
| Grouping | n | h | π | AR | HO | HE | FIS |
|---|---|---|---|---|---|---|---|
| All samples | 467 | 0.78 | 0.0058 | 6.36 | 0.59 | 0.64 | 0.079 (0.051 to 0.113) |
| Geographic region | |||||||
| Central interior | 356 | 0.76 | 0.0048 | 5.98 | 0.59 | 0.62 | 0.058 (0.028 to 0.093) |
| Boreal | 111 | 0.66 | 0.0061 | 5.66 | 0.60 | 0.64 | 0.076 (0.028 to 0.15) |
| STRUCTURE K = 4 | |||||||
| Red-circle | 127 | 0.69 | 0.0064 | 5.67 | 0.56 | 0.64 | 0.129 (0.04 to 0.291) |
| Green-square | 120 | 0.71 | 0.0045 | 4.98 | 0.54 | 0.59 | 0.084 (−0.007 to 0.215) |
| Orange-star | 156 | 0.72 | 0.0038 | 5.52 | 0.55 | 0.60 | 0.096 (0.006 to 0.255) |
| Blue-triangle | 64 | 0.71 | 0.0043 | 4.80 | 0.53 | 0.59 | 0.107 (0.017 to 0.293) |
| TESS K = 4 | |||||||
| Red | 126 | 0.71 | 0.0065 | 5.68 | 0.60 | 0.65 | 0.074 (0.03 to 0.128) |
| Green | 134 | 0.73 | 0.0047 | 5.06 | 0.58 | 0.60 | 0.039 (−0.01 to 0.093) |
| Orange | 131 | 0.72 | 0.0038 | 5.59 | 0.60 | 0.61 | 0.017 (−0.015 to 0.051) |
| Blue | 76 | 0.72 | 0.0043 | 4.94 | 0.58 | 0.59 | 0.037 (−0.008 to 0.091) |
Genetic diversity measures of 497 haplotypes and 467 genotypes derived from 13 microsatellite loci from samples collected from fishers throughout British Columbia, Canada, between 1997 and 2017, including number of individuals genotyped (n), haplotype diversity (h), nucleotide diversity (π), mean rarefied allelic richness (AR), mean observed heterozygosity (HO), mean expected heterozygosity (HE), and mean inbreeding coefficient (FIS; 95% CI). STRUCTURE and TESS groupings were clusters to which each individual was assigned (color-shape refers to symbols in Figs. 5 and 6), whereas geographic groupings were based on the region from which samples were collected.
| Grouping | n | h | π | AR | HO | HE | FIS |
|---|---|---|---|---|---|---|---|
| All samples | 467 | 0.78 | 0.0058 | 6.36 | 0.59 | 0.64 | 0.079 (0.051 to 0.113) |
| Geographic region | |||||||
| Central interior | 356 | 0.76 | 0.0048 | 5.98 | 0.59 | 0.62 | 0.058 (0.028 to 0.093) |
| Boreal | 111 | 0.66 | 0.0061 | 5.66 | 0.60 | 0.64 | 0.076 (0.028 to 0.15) |
| STRUCTURE K = 4 | |||||||
| Red-circle | 127 | 0.69 | 0.0064 | 5.67 | 0.56 | 0.64 | 0.129 (0.04 to 0.291) |
| Green-square | 120 | 0.71 | 0.0045 | 4.98 | 0.54 | 0.59 | 0.084 (−0.007 to 0.215) |
| Orange-star | 156 | 0.72 | 0.0038 | 5.52 | 0.55 | 0.60 | 0.096 (0.006 to 0.255) |
| Blue-triangle | 64 | 0.71 | 0.0043 | 4.80 | 0.53 | 0.59 | 0.107 (0.017 to 0.293) |
| TESS K = 4 | |||||||
| Red | 126 | 0.71 | 0.0065 | 5.68 | 0.60 | 0.65 | 0.074 (0.03 to 0.128) |
| Green | 134 | 0.73 | 0.0047 | 5.06 | 0.58 | 0.60 | 0.039 (−0.01 to 0.093) |
| Orange | 131 | 0.72 | 0.0038 | 5.59 | 0.60 | 0.61 | 0.017 (−0.015 to 0.051) |
| Blue | 76 | 0.72 | 0.0043 | 4.94 | 0.58 | 0.59 | 0.037 (−0.008 to 0.091) |
| Grouping | n | h | π | AR | HO | HE | FIS |
|---|---|---|---|---|---|---|---|
| All samples | 467 | 0.78 | 0.0058 | 6.36 | 0.59 | 0.64 | 0.079 (0.051 to 0.113) |
| Geographic region | |||||||
| Central interior | 356 | 0.76 | 0.0048 | 5.98 | 0.59 | 0.62 | 0.058 (0.028 to 0.093) |
| Boreal | 111 | 0.66 | 0.0061 | 5.66 | 0.60 | 0.64 | 0.076 (0.028 to 0.15) |
| STRUCTURE K = 4 | |||||||
| Red-circle | 127 | 0.69 | 0.0064 | 5.67 | 0.56 | 0.64 | 0.129 (0.04 to 0.291) |
| Green-square | 120 | 0.71 | 0.0045 | 4.98 | 0.54 | 0.59 | 0.084 (−0.007 to 0.215) |
| Orange-star | 156 | 0.72 | 0.0038 | 5.52 | 0.55 | 0.60 | 0.096 (0.006 to 0.255) |
| Blue-triangle | 64 | 0.71 | 0.0043 | 4.80 | 0.53 | 0.59 | 0.107 (0.017 to 0.293) |
| TESS K = 4 | |||||||
| Red | 126 | 0.71 | 0.0065 | 5.68 | 0.60 | 0.65 | 0.074 (0.03 to 0.128) |
| Green | 134 | 0.73 | 0.0047 | 5.06 | 0.58 | 0.60 | 0.039 (−0.01 to 0.093) |
| Orange | 131 | 0.72 | 0.0038 | 5.59 | 0.60 | 0.61 | 0.017 (−0.015 to 0.051) |
| Blue | 76 | 0.72 | 0.0043 | 4.94 | 0.58 | 0.59 | 0.037 (−0.008 to 0.091) |
In our assessment of isolation-by-distance patterns, we found that correlations between paired matrices of genetic distance and Euclidean distance differed among areas. When we considered all samples collectively, we rejected the null hypothesis that genetic distance and Euclidean distance were unrelated (Mantel’s r = 0.17; P < 0.001) with a weak positive relationship between genetic distance and Euclidean distance. However, isolation by distance was substantially less prevalent when we considered samples from different geographic regions: boreal region, Mantel’s r = 0.07, P < 0.048; central interior region north of Nechako and Upper Fraser Rivers, Mantel’s r = 0.08, P < 0.031; central interior region south of Nechako and west of Fraser River, Mantel’s r = 0.01, P < 0.335; central interior region east of Fraser River, Mantel’s r = −0.05, P < 0.837 (see Fig. 1 for river locations).
Assessment of putative populations
We used microsatellite data to evaluate whether broad patterns of structuring were consistent with 2 putative populations based upon geographic region. Samples from the central interior region had FIS = 0.058 (95% CI = 0.029 to 0.080), whereas fishers from the boreal region had FIS = 0.076 (95% CI = 0.032 to 0.111). We observed 14 private alleles among fisher samples in the central interior region and 5 private alleles in the boreal region. Using microsatellite loci, we observed significant structuring (AMOVA, P < 0.001) and calculated an FST value of 0.046 between fishers in the central interior and boreal regions of the province. When we used PCA to compare allele frequencies between the 2 geographic regions, we observed that axis 1 of the PCA plot explained 8.1% of the genetic variation within the data set while axis 2 explained 6.2%, which suggested differentiation between the central interior and boreal regions but with considerable overlap. The center of the distribution of genotypes in PCA-adjusted space for fishers from the central interior region did not overlap the 75th percentile inertia ellipse from the boreal region, and vice versa (Fig. 4).
Population cluster analysis
Both the MeanLnP(K) and Puechmaille methods showed the strongest support for K = 4 clusters in the STRUCTURE analysis of the full data set (Fig. 5; Supplementary Data SD1). With this cluster arrangement, 69% of samples in the red-circle cluster were confidently assigned (i.e. q-value ≥ 0.7; color-shapes from Fig. 5), 65% were confidently assigned in the blue-triangle cluster, and 40% and 49% were confidently assigned in the green-square and orange-star clusters, respectively. The correct self-assignment rate was 98% for samples collected in the boreal region/red-circle cluster, whereas within the central interior region, the correct self-assignment rate was 73% for samples collected north of the Nechako and Upper Fraser Rivers/green-square cluster, 90% for samples collected south of the Nechako River and west of the Fraser River/orange-star cluster, and 81% for samples collected east of the Fraser River/blue-triangle cluster (cluster color-shapes from Fig. 5). Only 2 of the 340 samples that were assigned to clusters from the central interior region (i.e. green-square, orange-star, and blue-triangle clusters) were collected in the boreal region. A total of 19 of the 127 samples that were assigned to the red-circle cluster were collected in the central interior region, 95% of which were within 100 km of the transition between boreal and central interior regions (mean = 67.2 km, SD = 34.9 km; Fig. 5). Fixation indices (FST) ranged from 0.026 between the green-square and orange-star clusters to 0.069 between the red-circle and blue-triangle clusters (Fig. 5). Using the Evanno method, the 2-cluster assessment (i.e. K = 2) received the strongest support, with the ΔK for K = 2 (110.8) scoring 6 times higher than the next-best cluster sizes (ΔK = 16.9 for K = 3; ΔK = 6.5 for K = 4).
When we examined samples collected from the central interior region, we detected the strongest support for a cluster size of K = 3 among samples using the Puechmaille method, although MeanLnP(K) reached its maximum at K = 4 with a value of −10,827.9 (SD = 13.8; Supplementary Data SD2). The clusters that were produced with this STRUCTURE analysis were very similar to the 3 clusters from the central interior region estimated from the full-data set STRUCTURE analysis, with 315 of 356 samples (88%) assigned to the same clusters between analyses. Within the boreal region, the maximum MaxMedK and MaxMeanK scores were for a cluster size of 2, although the MedMedK, MedMeanK, and maximal MeanLnP(K) score indicated that a single cluster occurred in this region (Supplementary Data SD3).
We also evaluated population clustering using analyses that incorporated spatial location of the samples. The TESS analysis converged and the DIC curve for the TESS analysis of K = 2 to K = 10 clusters suggested a slight change in curvature at K = 4 or K = 5 (Fig. 6), a result that may occur when fine population structure exists (Caye et al. 2016). Because the aspatial STRUCTURE analysis suggested that fine-scale structure was indeed occurring within Fisher populations in British Columbia, we used TESS spatial analyses at K = 4 clusters to compare aspatial and spatial cluster results. Values of FST between spatial clusters ranged between 0.031 and 0.071 at K = 4, with the higher values between the boreal region cluster and other genetic clusters that occurred in the central interior region (Fig. 6).
Aspatial cluster analysis of 467 genotypes derived from 13 microsatellites collected from throughout the range of fishers between 1997 and 2017 in British Columbia, Canada. Top graph shows the apportioning of ancestry in STRUCTURE for K = 4, sorted by geographic region and sub-region from which samples were drawn. Map shows the location of sample collection for individuals assigned by STRUCTURE to each cluster at K = 4. Putative geographic barrier between the boreal and central interior regions is shown, as are fixation index scores (FST) for comparisons between clusters.
DIC scores of TESS analysis of K = 2 to 10 and spatial pattern of population membership inferred by TESS at K = 4 for fishers collected between 1997 and 2017 in British Columbia, Canada, based on 467 genotypes derived from 13 microsatellites. Fixation index scores (FST) are shown for comparisons between K = 4 clusters.
We observed high congruence in assignment of samples between STRUCTURE and TESS K = 4 clusters. That is, 95% of samples that were assigned to the red-circle STRUCTURE cluster were assigned to the red TESS cluster, 83% of samples that were assigned to orange-star STRUCTURE cluster were also assigned to the orange TESS cluster, 94% of the samples that were assigned to blue-triangle STRUCTURE cluster were assigned to the blue TESS cluster, and 96% of samples that were assigned to the green-square STRUCTURE cluster were also assigned to the green TESS cluster. Average genetic diversity measures for each region and cluster varied (Table 1).
Effective population size
The effective population size (Ne) of the central interior region was 252 (95% CI = 185 to 332) and 136 (95% CI = 92 to 234) for the boreal region of British Columbia (i.e. not including the larger pan-boreal population that likely extended across Canada).
Gene flow
We evaluated gene flow within fishers in British Columbia using several different approaches. Bayesian estimates of gene flow between the 2 geographic regions converged (see Supplementary Data SD4 for tracer output) and were nearly identical at 2.3% of the population per generation (95% CI: 0% to 4.9%) for the central interior and 2.6% (95% CI: 1.2% to 4%) for the boreal regions. We also calculated the number of migrants per generation (Nm) using the θ and M parameter values estimated by Migrate-n (i.e. the number of migrants per generation moving to population x from population y would be represented by the following equation: Nm = [(θx * My → x)/4]) that, similar to our Bayesian results, indicated nearly identical gene flow between the 2 geographic populations at 3 migrants per generation (95% CI = 2.8 to 3.5) for the central interior region and 3 migrants per generation (95% CI = 2.6 to 3.3) for the boreal region. Finally, using the FST indirect measures of gene flow criterion (FST = 1/(4Nem + 1)), an FST = 0.046 between the central interior and boreal regions suggested 5 migrants into each population per generation. Following the methods of Vucetich and Waite (2000), with the mean standard deviation of log counts equal to 0.12 (from Weir and Corbould 2006; Weir et al. 2011) and Ne/Nt estimated at 0.24, we calculated that 90 migrants per generation were required to maintain genetic homogeneity between fishers in the central interior and boreal regions of British Columbia.
Discussion
We observed complex patterns of variation in neutral genetic markers showing that fishers in British Columbia were not a single panmictic population with genes moving widely throughout the province, but rather, we detected significant population structuring, relatively little gene flow between regions, and substantial population differentiation. Our geographic, aspatial, and spatial population cluster analyses consistently showed a complex pattern of hierarchical population structuring that appeared to be linked to geographical features. Relatively restricted gene flow among these populations likely facilitated population differentiation that, when coupled with small population sizes, puts several populations of fishers in British Columbia at elevated conservation risk.
Our clustering assessments showed that fishers in British Columbia were hierarchically structured into several distinct populations. In our STRUCTURE analysis, we identified K = 4 as the most likely cluster size, although we also detected a strong signature of K = 2 using the ΔK method, which has been shown to be indicative of upper-level hierarchical structuring (Evanno et al. 2005). Puechmaille (2016) noted that, in populations with hierarchical stepping-stone genetic structure, ΔK was proficient at identifying the highest level of population structure even though other metrics are needed to detect substructuring adequately (e.g. MaxMedK, MaxMeanK). Hierarchical stepping-stone structure is typified by a contact zone that separates demes (e.g. subpopulations) into 2 upper-level groups (e.g. populations), with demes in each upper-level group primarily exchanging genes with neighboring demes within that group. For example, subpopulations on 1 side of a contact zone would exchange genes more frequently with neighboring subpopulations than with those on the other side of the contact zone. Indeed, we found that cluster assignments of our samples at K = 4 were nested within the upper hierarchical K = 2 cluster assignments. These results, along with our observation of K = 4 being a likely cluster size in the TESS analysis, leads us to believe that fishers in British Columbia were structured as hierarchical stepping-stone populations where an upper hierarchical level of genetic partitioning resulted in 2 populations at the provincial scale, with 3 subpopulations occurring in the central interior region.
Not only was structure evident within the range of fishers in British Columbia, these populations and subpopulations exhibited significant differences in patterns of genetic variation. At the highest level of genetic partitioning, differentiation was pronounced among haplotype data (i.e. FST = 0.272), with significant differentiation between the genetic characteristics of fishers on either side of the Rocky, Omineca, and Skeena Mountain ranges. Differentiation was also evident at this scale among microsatellite allelic data (mean FST = 0.049 among methods), which was similar to values documented in a previous small study of fishers that included British Columbia (FST = 0.06; Kyle et al. 2001) and to “sky island” populations of American martens that were structured by mountain ranges in southeastern British Columbia (FST = 0.048; Lucid et al. 2020). Significant stepping-stone differentiation also occurred among subpopulations in the central interior population, with neighboring subpopulations showing less differentiation than more distant subpopulations (e.g. Fig. 5; red-circle–blue-triangle clusters FST = 0.069 compared to FST of 0.026 between neighboring green-square–orange-star clusters), but differentiation among subpopulations was not as extreme as observed among discrete fisher populations in California (FST = 0.054 to 0.127; Tucker et al. 2013). Differentiation was also evident among mtDNA haplotypes at the subpopulation scale. Several features of the current-day geography of British Columbia likely contributed to the development of this hierarchical structure and pattern of population differentiation.
At a broad provincial scale, high-elevation mountain ranges likely reduced the exchange of genes between higher-level clusters such that significant differentiation occurred between fishers in the central interior and boreal regions. Our estimated bidirectional rate of genetic exchange across the Rocky, Omineca, and Skeena Mountain ranges was ~3 migrants per generation, which although greater than that seen between two highly disjunct Fisher populations in Idaho (0.672 – 0.95 migrants per generation; Lucid et al. 2019), was substantially lower than the estimated 90 migrants per generation needed to preserve genetic homogeneity. Mountain ranges act as barriers to gene flow for other forest-dependent carnivores, like Canadian Lynx (Lynx canadensis), for which the Rocky Mountains separates populations in western Canada (Rueness et al. 2003; Watt et al. 2021) and our results are consistent with past work that has shown that geographic or anthropogenic barriers as narrow as 50 km may impede genetic exchange in fishers (Aubry et al. 2004). Although the low levels of gene flow and the rarity of mis-assignment of samples at this upper level of genetic partitioning suggested that >30 km of unsuitable habitat contributed to differentiation of fisher populations in British Columbia, the Rocky, Omineca, and Skeena Mountain ranges did not appear to be an absolute barrier to gene flow. Specifically, we observed that several samples collected in the Rocky Mountain Trench area north of the upper Fraser River were assigned to the ‘boreal’ population in both the STRUCTURE and TESS clustering, despite occurring on the western side of the mountain ranges. It is unclear as to whether this is the result of recent immigration from the boreal population or whether the boreal population currently extends to this area. More sampling from this area is needed to provide us with a better understanding of the source of this incongruity.
We noted further substructuring and differentiation within the plateau landscapes of central British Columbia, suggesting that other factors contributed to the complex characteristics of the population genetics of fishers in the province. Specifically, large rivers may play a substantial role in substructuring the Fisher population in the central interior region, as both the STRUCTURE and TESS analysis showed consistent geographically discrete genetic clusters whose peripheries aligned with large rivers (Figs. 5 and 6). These rivers were among the largest that occurred within the central interior population and ranged in flow from 239 m3/s (Nechako River) and 814 m3/s (upper Fraser River) to 1,541 m3/s for the main Fraser River (natural long-term mean annual discharge [LTMAD]; Environment and Climate Change Canada 2022). Gene flow of fishers within landscapes has been shown to be affected by the discontinuity of suitable habitat (Tucker et al. 2017) and rivers (Wisely et al. 2004), and the substantial clearing of land, human development, and trapping pressure that occurred along these rivers since European colonization may have recently contributed to structuring within this region. Fishers are capable of crossing medium-sized rivers during winter (e.g. Quesnel River, 99 m3/s; Weir and Harestad 1997), so it is likely that larger rivers were not an absolute barrier but simply reduced gene flow to a sufficient level that structuring and differentiation occurred. The lack of strong within-population structure observed in the boreal region may be indicative of relatively continuous habitat, lower genetic diversity influenced by effective population size and inbreeding, or perhaps because we sampled only a small portion of the much larger pan-Canadian boreal population.
Current-day patterns of genetic variation were also likely affected by British Columbia’s rich and elaborate glacial history that, over the past 20,000 years, has influenced the structure of many cryptic refugia producing complicated genetic histories for many wildlife species. Evidence from the genetic and postglacial biogeography of other habitat specialist mustelids in British Columbia, like American martens, Pacific martens (M. caurina), and North American badgers (Taxidea taxus), may shed light on factors that produced the patterns that we observed in fishers. American martens are believed to have recently expanded from refugia in eastern North America following glacial retreat, whereas Pacific martens show signs of earlier expansion, potentially from multiple refugia along the west coast of British Columbia (Dawson and Cook 2012) and subclade divergence occurred during separation into distinctive refugial populations in western North America (Stone et al. 2002; Schwartz et al. 2020). Likewise, genetic characteristics of North American badgers show a similar pattern where endemic haplotypes developed in cryptic refugia in western British Columbia (Ford et al. 2020). Shafer et al. (2010) found considerable evidence in the literature for multiple refugia off the coast of British Columbia and Alaska, some of which could have housed a relic Fisher population during the last glaciation. Given the distribution of haplotypes that we observed along with the microsatellite results, it is possible that fishers from the central interior region could be remnants of fishers from these glacial refugia, but additional work is needed to evaluate this hypothesis more thoroughly.
Although our population structure and differentiation findings were consistent with the effects that would be anticipated if geographic features affected gene flow, we did not explicitly test these effects. A more thorough assessment of the impact of landscape factors on gene flow within and among populations could be achieved by employing individual-based landscape genetics analyses (e.g. Balkenhol et al. 2020) using more spatially precise genetics samples than ours. These analyses would allow us to better identify the features that affect gene flow, predict gene flow among populations given current and future predicted landscape conditions, and provide natural resource managers with a better understanding of how their land-use decisions affect connectivity of vulnerable populations.
Given that genetic divergence and discreteness are a function of gene flow, population size, and time (Allendorf et al. 2010), the geographic disjunction caused by current-day high-elevation, deep-snow mountain ranges and rivers, along with past glacial history, likely reduced the amount of gene flow for an extended period such that genetic drift and mutation occurred. The biogeographical history of ecosystems that support fishers in British Columbia is not well understood but it is possible that the geographic barrier between the central interior and boreal regions was not always as prominent as it is currently. During the warm dry early Holocene Epoch (10 to 7k before present), dry pine-dominated forests reached well into, if not fully occupied, the current extent of the high-elevation Engelmann Spruce–Subalpine Fir biogeoclimatic zone (Heinrichs et al. 2001) that fishers avoid. During this period, more expansive connectivity corridors between fishers in the boreal and central interior regions may have existed and, as a result, greater genetic exchange may have occurred. This may explain the moderate level of differentiation that we observed in the microsatellite compared to the haplotype data, whereby greater levels of admixture occurred in the past but have since diminished with increased barrier effect as the climate changed to current-day conditions. Alternatively, it is possible that differentiation occurred more recently and that the current-day conditions simply contributed to observed patterns of differentiation. Additional research into this question is needed to better understand the timing of differentiation among the populations and subpopulations.
The high number of migrants needed to preserve genetic homogeneity among Fisher populations in British Columbia, relative to the “one migrant per generation” rule-of-thumb of Wright (1931), arose because several attributes of Fisher populations in the province stray substantially from Wright’s “ideal” population, upon which this rule is based. In British Columbia, Fisher populations are characterized by fluctuations in numbers from year to year (Weir and Corbould 2006; Weir et al. 2011), variance in fecundity (Lofroth et al. 2023), skew in sex ratios (Weir 2003), and overlapping generations. Most of these factors reduce the effective population size Ne relative to the census population size (e.g. Nunney 1991, 1993; overlapping generations moderate this effect), which increases the number of migrants needed to maintain genetic homogeneity and avoid genetic drift (Mills and Allendorf 1996; Vucetich and Waite 2000). These effects have been shown empirically in Brook Trout, where populations are more likely to experience rapid changes in allelic frequency through genetic drift when fewer females contribute to each generation (Nathan et al. 2017).
Although isolation by distance and hierarchical population structuring can be difficult to differentiate (Meirmans 2012), we believe that the patterns of genetic variation within our data were not caused by strong isolation by distance, but rather by true genetic structuring. This is because our sampling and analysis approach allowed us to avoid the pitfalls of identifying false clusters in the presence of isolation by distance (Perez et al. 2018). First, we collected a large number of samples across the entirety of Fisher range in the province, which increased the likelihood that genotypes that we collected were representative of the distribution of genes in the population and helped ensure that STRUCTURE could differentiate between clusters and clines. Second, our observation that the scope and scale of Mantel tests diminished considerably at subregional geographic scales suggests that isolation by distance was not consistent across the entire data set. Most importantly, our aspatial STRUCTURE results were highly congruent with the spatial TESS clustering analyses, which is a powerful indicator of the diminished impact of isolation by distance on the inference of population structure (Perez et al. 2018). All of these factors lead us to be confident that isolation by distance did not substantially affect our ability to discern the population structure of fishers in British Columbia.
In addition to genetic drift and differentiation that we observed in neutral markers, the environmental and ecological conditions in the central interior region of British Columbia may have provided a selective regime where, although unmeasured in this study, local adaptation may have occurred as well. Snow conditions in the central interior region are typically deeper and softer (e.g. average annual snowfall = 316 cm, Sub-Boreal Spruce biogeoclimatic zone; DeLong 2004) compared to the boreal region (e.g. average = 169 cm, Boreal White and Black Spruce biogeoclimatic zone; MacKenzie et al. 2020) and this difference likely affects the ability of fishers to move about their environment (Raine 1983). Furthermore, boreal forests typically support higher diversity and abundances of key prey for fishers than forests in the central interior region (Lofroth et al. 2023). As a result, fishers exploit habitats differently for foraging and moving about their territories in addition to using different habitat features for reproduction and resting (Weir and Corbould 2008; Weir and Lara Almuedo 2010). This likely resulted in local adaptations by fishers in the central interior region that do not occur elsewhere. Evidence of the effect in these differential biotic and abiotic selective pressures can be found in fishers from the central interior region, which have smaller body masses and dimensions, different diets, and half the reproductive output of their counterparts in the boreal region (Lofroth et al. 2023; RDW, unpublished data).
The level of population differentiation among fishers in British Columbia was sufficient that heritable markers could be used to distinguish individuals among regions. The distribution of mtDNA haplotypes was different between the central interior and boreal regions, with 1 haplotype found exclusively in the central interior region (Drew-1; Fig. 3) and 3 occurring predominately in 1 region or the other (Drew-4, Drew-9, Drew-11). Among microsatellite loci, 14 private alleles occurred in fishers from the central interior range portion and 5 private alleles appeared in fishers from the boreal range portion, representing 18% and 8% of the alleles that occurred in each region, respectively. Both spatial and aspatial clustering approaches showed consistent and convergent clustering of genotypes that dependably differentiated fishers from the central interior and boreal regions.
Our genetic results were also consistent with earlier morphometric analyses that identified fishers in the central interior region of British Columbia as belonging to the P. (M.) p. columbiana subspecies (Goldman 1935) which, when considered along with mtDNA analyses that supported subspecies designations (Drew et al. 2003), gives additional support for the validity of the columbiana subspecies. When considering the spatial distribution of the 6 haplotypes within British Columbia (Fig. 2), 4 were restricted or largely restricted by the geographic barrier, which corroborates a putative genetic boundary between the 2 regions. Furthermore, although the Drew-1 haplotype shows a near North America-wide distribution, full mitogenome sequencing showed that this haplotype was comprised of 4 distinct nonsister lineages that sort by subspecies and geographic region, including 1 endemic to the central interior region of British Columbia (Knaus et al. 2011). Additional sampling and a more thorough analysis of the full mitogenome and other genetic measures is needed to evaluate the relationship more fully between the fishers in the central interior region of the province and the P. p. columbiana subspecies.
Given their genetic differentiation (this study) and likely local adaptations, the 2 demes of fishers that occur in British Columbia likely function as demographically and genetically discrete populations. As such, we propose that fishers in the central interior region of the province be described as the Columbian population in reference to their coincidence with the columbiana subspecies, whereas fishers in the boreal forest region be referred to as the Boreal population (British Columbia Conservation Data Centre 2020). Furthermore, the Columbian population may constitute a sufficiently discrete and evolutionarily significant population that it could be considered a Designatable Unit by the Committee on the Status of Endangered Wildlife In Canada.
Conservation implications
Whereas the population of fishers in the boreal region of British Columbia is likely contiguous with the larger Fisher population across the boreal forest and into eastern North America (Proulx et al. 2004), the Columbian population appears isolated and receives little genetic and demographic support from elsewhere. Fishers in this population have low population growth rates (Lofroth et al. 2023) and are extremely susceptible to loss of forest cover through logging (Weir and Corbould 2010); population effects that would be exacerbated by additive trapping mortality (Fogarty et al. 2022). Recent declines in the population to fewer than 511 individuals (British Columbia Conservation Data Centre 2020), increasing habitat and population threats, low effective population size, and evidence of possible inbreeding (i.e. FIS = 0.058), all combine to put this population at heightened conservation risk. These factors, when coupled with low potential for genetic and demographic rescue, are expected to increase the rate of genetic drift and inbreeding such that the ability of the Columbian population to persist into the future is tenuous, particularly given the substructuring that was apparent within this population. As a result, swift and effective conservation actions are needed to preserve this vulnerable population of fishers and avoid its extirpation.
Supplementary data
Supplementary data are available at Journal of Mammalogy online.
Supplementary Data SD1.—MeanLnP(K), ΔK, and MedMedK, MedMeanK, MaxMedK, and MaxMeanK scores for clusters K = 1 to 10 in aspatial STRUCTURE population analysis using genotypes from 13 microsatellites of 467 fishers collected between 1997 and 2017 in British Columbia, Canada.
Supplementary Data SD2.—MeanLnP(K), ΔK, and MedMedK, MedMeanK, MaxMedK, and MaxMeanK scores for clusters K = 1 to 10 in aspatial STRUCTURE population analysis using genotypes from 13 microsatellites of 356 fishers collected between 1997 and 2017 from the central interior region of British Columbia, Canada.
Supplementary Data SD3.—MeanLnP(K), ΔK, and MedMedK, MedMeanK, MaxMedK, and MaxMeanK scores for clusters K = 1 to 10 in aspatial STRUCTURE population analysis using genotypes from 13 microsatellites of 111 fishers collected between 1997 and 2017 from the boreal region of British Columbia, Canada.
Supplementary Data SD4.—Tracer output from BayesAss estimation of gene flow of fishers collected between 1997 and 2017 in British Columbia, Canada.
Acknowledgments
Many thanks to T. Hooft for collecting tissue samples from harvested pelts, C. Engkjer for laboratory assistance, B. Harrower for preliminary data analysis and discussions, and 3 anonymous reviewers for thoughtful reviews and comments on earlier versions of this manuscript.
Author contributions
RDW: Conceptualization, Funding acquisition, Data curation, Formal analysis, Methodology, Writing—original draft, Writing—review & editing, Supervision, Project administration. AMR: Methodology, Formal analysis, Writing—original draft, Writing—review & editing. LR: Writing—original draft, Writing—review & editing, Visualization. KLP: Investigation, Resources. MKS: Investigation, Resources. MKL: Methodology, Formal analysis, Writing—original draft, Writing—review & editing, Supervision.
Funding
Funding for this work was provided by Habitat Conservation Trust Foundation (grant 0-460), Environment and Climate Change Canada, and Ministry of Environment and Climate Change Strategy.
Conflict of interest
None declared.
Data availability
We deposited derived mtDNA sequences to GenBank as accession numbers ON100606 to ON100611.
