Carnivore Contact: A Species Fracture Zone Delineated Amongst Genetically Structured North American Marten Populations (Martes americana and Martes caurina)

North American martens are forest dependent, influenced by human activity, and climate vulnerable. They have long been managed and harvested throughout their range as the American marten (Martes americana). Recent work has expanded evidence for the original description of two species in North America — M. americana and the Pacific Coast marten, M. caurina — but the geographic boundary between these groups has not been described in detail. From 2010 to 2016 we deployed 734 multi-taxa winter bait stations across a 53,474 km² study area spanning seven mountain ranges within the anticipated contact zone along the border of Canada and the United States. We collected marten hair samples and developed genotypes for 15 polymorphic microsatellite loci for 235 individuals, and 493 base-pair sequences of the mtDNA gene COI for 175 of those individuals. Both nuclear and mitochondrial genetic structure identified a sharp break across the Clark Fork Valley, United States with M. americana and M. caurina occurring north and south of the break, respectively. We estimated global effective population size (N_e) for each mountain range, clinal genetic neighborhood sizes (NS), calculated observed (H_o) and expected (H_e) heterozygosity, fixation index (F_ST), and clinal measures of allelic richness (Ar), H_o, and inbreeding coefficient (F_IS). Despite substantial genetic structure, we detected hybridization along the fracture zone with both contemporary (nuclear DNA) and historic (mtDNA) gene flow. Marten populations in our study area are highly structured and the break across the fracture zone being the largest documented in North America (F_ST range 0.21–0.34, mean = 0.27). With the exception of the Coeur d’Alene Mountains, marten were well distributed across higher elevation portions of our sampling area. Clinal NS values were variable suggesting substantial heterogeneity in marten density and movement. For both M. americana and M. caurina, elevationaly dependent gene flow and high genetic population structure suggest that connectivity corridors will be important to ensuring long-term population persistence. Our study is an example of how a combination of global and clinal molecular data analyses can provide important information for natural resource management.

Introduction

North American marten (Martes spp.) are a group for which genetic data are changing the understanding of how populations and lineages are structured across large landscapes (e.g., Colella et al., 2019) and provide an important example of how genetic data can inform wildlife management. Marten are actively harvested and managed throughout North America and may be sensitive to forest management (Kirk and Zielinski, 2009; Cushman et al., 2011; Moriarty et al., 2016), human recreation (Slauson et al., 2017), and climate change (Cushman et al., 2011; Wasserman et al., 2012a).

Originally classified as two species — the American (Martes americana) and Pacific Coast (Martes caurina) marten (Merriam, 1890) — further morphometrics resulted in the assignment of M. caurina as a subspecies of M. americana (Wright, 1946). Although some authors continued to refer to two separate groups of marten in North America (Dawson and Cook, 2012), marten were subsequently considered, and managed, under the single species umbrella of M. americana (Colella et al., 2019). Recent molecular work has supported the original morphologically based two species classification of M. caurina and M. americana (Stone et al., 2002; Small et al., 2003; Dawson et al., 2017). Analysis of widely spaced samples from western North America supports the longstanding notion (Wright, 1946) that there is a broad area of contact in the Inland Pacific Northwest including portions of British Columbia, Idaho, Montana, and Washington (Colella et al., 2019).

Lacking from the picture, but crucial to wildlife management, is more detailed sampling to delineate where each lineage occurs and the level of genetic structure and gene flow between and amongst those lineages. Although obtaining the data necessary to calculate these metrics can be challenging (Lindenmayer and Likens, 2010), the inclusion of multiple species in field surveys is becoming a cost effective way to increase the efficiency of field efforts because well-designed multi-species inventory techniques allow for simultaneous data collection for many species (Tucker and Simmons, 2009; Robinson et al., 2017; Lucid et al., 2018a).

From 2010 to 2016 we deployed 734 multi-taxa winter bait stations across a 53,474 km² study area within the known marten contact zone (Colella et al., 2019) from which we collected mammal hair samples. We used marten hair samples from those collections to develop this study for which our objectives were to (1) identify the location of the fracture zone between lineages, (2) measure genetic structure within and between lineages, (3) calculate clinal and global effective population sizes for mountain ranges, and (4) consolidate that information into conservation and management recommendations.

Materials and Methods

Study Area

The study area consisted of a 53,474 km² area containing portions of British Columbia, Idaho, Montana, and Washington (Figure 1). It comprised portions of the Coeur d’Alene (CDA), Purcell, Saint Joe (Joe), Central Selkirk, South Selkirk, West Cabinet, and Valhalla Mountain ranges. Mountain ranges are generally intersected by wide, agriculturally dominated, low elevation valleys. Elevations range from 300 m in valley bottoms in the southern portion of the study area to >3000 m on mountain tops toward the north.

FIGURE 1

Figure 1. Study area map and species range map (inset). Marten detections include DNA, remote camera (US), and track (BC) detections.

Field Methods

We placed a survey grid over the seven mountain ranges in our study area, consisting of 453 5 × 5 km survey cells in the United States (US) and 237 10 × 10 km survey cells in Canada. Grid size in Canada was larger due to remoteness and lack of accessibility of this portion of the study area. We deployed ≥1 winter multi-species bait station per cell with a total of 734 bait stations deployed in 680 cells. Bait stations consisted of a bait (ungulate portion or frozen beaver [Castor canadensis]) attached to a tree, remote camera (all US stations and some Canada stations), and 12.30 caliber bronze rifle bore gun brushes (US) or barbed wire (Canada) positioned concentrically around the tree that acted as hair snags (see Robinson et al., 2017 for details). At Canadian stations, track surveys were completed at each bait station visit (see Kortello et al., 2019 for details).

DNA Extraction, Microsatellite Amplifications, and mtDNA Sequencing

DNA extraction from hair samples and microsatellite amplification was performed using standard protocols (Paetkau, 2003; Kendall et al., 2009). We genotyped 241 marten individuals at 16 microsatellite loci (Ma-1, Ma-2, Ma-3, Ma-5, Ma-7, Ma-8, Ma-11, Ma-14, Ma-15, Ma-18, Ma-19 (Davis and Strobeck, 1998), Mer041, Mvis020 (Flemming et al., 1999), MP0055, MP0085, MP0227 (Jordan et al., 2007), plus an un-published ZFX/ZFY gender locus (primers w542: TTC CAG GCA GTA CCA AAC AG and w543: AGG AAA TCA TTC ATG AAT ATC ACT). We used GENEPOP 3.4 option 1 (Raymond and Rousset, 1995) to test each marker for a deficit of heterozygous genotypes. We produced 493 base pair sequences of the mitochondrial cytochrome c oxidase I (COI) for 175 individual marten (GenBank accession numbers MT578095 – MT578269). Macrogen Inc., produced sequence reads using PCR product we provided and primers we designed (forward 5′-CCTATTGATTCGCGCTGA-3′; reverse 5′-GTCTGTAAGTAGTATGGTAATGCC-3′).

Phylogeny Estimation

To place our mtDNA data in geographic context we used published 493 base pair COI sequences (Colella et al., 2019) from 82 specimens (GenBank accession numbers MK320161 – MK320242). Sequence alignment was performed using Geneious 7.1.7 (Biomatters Ltd.) and the subsequent alignment and translation was checked by eye; no indels or stop codons were present. We then estimated the mtDNA gene tree under maximum parsimony (MP) and maximum likelihood (ML) optimization criteria using PAUP^∗ (v. 4.0a, build 166; Swofford, 2003). Maximum Parsimony was performed using the heuristic search algorithm with random addition replicates, tree bisection-reconnection (TBR) branch swapping, and all sites and transformations were weighted equally. For the model-based ML analyses, we used the first (of 4) MP trees to evaluate models of nucleotide sequence evolution via the automodel command in PAUP^∗ under the Bayesian Information Criterion (BIC) and decision theory (DT) (Minin et al., 2003). We then conducted heuristic searches under ML with the Kimura 2 parameter model (Kimura, 1980), TBR branch swapping, ten random-addition replicates, and a random starting tree. Nodal support for the ML tree was evaluated via bootstrap analysis with 500 pseudo-replicates (starting tree built using simple stepwise addition and NNI branch swapping) (Felsenstein, 1985).

Population Genetic Analyses

We assessed genetic structure without a priori constraints using the clustering software STRUCTURE v2.3 (Pritchard et al., 2000). We performed 10 runs at each value of K from 1–5, with 250,000 MCMC iterations after a burn-in of 50,000, using the correlated allele frequencies and the admixture model.

We used GENEPOP 3.4 (Raymond and Rousset, 1995) to calculate observed (H_o) and expected (H_e) heterozygosity in each mountain range, estimate per locus F_ST values (Weir and Cockerham, 1984) between mountain ranges, and non-random associations of alleles within and between loci. We a priori defined populations by mountain range with the exception of the Central and South Selkirks, which we defined as two separate populations because they were separated by a large lake and low elevation valley with considerable human settlement. First, we estimated F_ST using 15 microsatellites (all except Mvis020). Second, to contextualize these observations, we used 9 of 11 markers used in a Canada-wide dataset (Kyle and Strobeck, 2003). We estimated contemporary effective population size (N_e) using the linkage disequilibrium (i.e., non-random allele associations) method (Waples and Do, 2008) as implemented in both Ne Estimator v2 (Do et al., 2014) and LDNe (Waples and Do, 2008) at three levels of lowest allele frequency (0.01, 0.02, 0.05).

Gene Flow Models

The phylogeny estimation and STRUCTURE analyses identified a population break in our study and we used this break to inform our gene flow analyses. We used GeneClass2 (Piry et al., 2004) to test all 48 southern individuals and the 75 northern individuals located closest to the break for evidence of ancestry from each side of the break (Paetkau et al., 1995, 2004).

We partitioned the entire dataset by lineage into northern (M. americana) and southern (M. caurina) populations and evaluated 3 putative models of gene flow and subsequently compared those models to investigate possible restrictions in the directionality of gene flow. We used the coalescent-based program migrate-n (Beerli and Felsenstein, 1999, 2001) to jointly estimate the mutation-scaled effective population size (θ = 4N_eμ) and the mutation-scaled migration rate (M = m/μ; where m and μ are the migration and mutation rates, respectively), which can be used to estimate gene flow by θ ^∗M = 4N_em where N_em is the number of migrants per generation. In model M. americana ⇄ M. caurina, gene flow is permitted in both directions between northern and southern populations. In model M. americana ⟶ M. caurina, gene flow is possible only in the southward direction while in model M. americana ⟵ M. caurina, gene flow is only in the northward direction. For each analysis in migrate-n, we used Bayesian inference for both microsatellite and mtDNA data. Starting values of θ and M were estimated with an F_ST calculation and the following uniform priors were specified: for microsatellites, θ, minimum = 0.0, maximum = 50.0, delta = 1.0 and M, minimum = 0.0, maximum = 500.0, delta = 50.0; for mtDNA, θ, minimum = 0.0, maximum = 0.1, delta = 0.001 and M, minimum = 0.0, maximum = 1000.0, delta = 100.0. Analyses were conducted using four replicates, four static-heated chains (1.0, 1.5, 3.0, 100,000.0) ran for 1,000,000 generations, sampling every 100, after 10,000 were discarded as burn-in. We inspected posterior distribution plots of estimated M and θ values (bin number = 1500) to assess convergence. We converted estimates of M and θ to the fraction of immigrants in population i that came from population j in the last generation via Nm = θ_i ^∗ M_{j → i} (this value was divided by 4 for the diploid, biparentally inherited microsatellites). In order to compare the different gene flow models, we obtained marginal likelihoods for each model (Kass and Raftery, 1995; Beerli and Palczewski, 2010), which can be used to evaluate multiple models when based on the same data using Bayes factors (Bloomquist et al., 2010). In the case of the 15 microsatellite loci, information from all loci was combined into a global estimate by Bezier approximation of the thermodynamic scores. The probability of a certain model was then obtained by dividing the exponent (on the base of e) log likelihoods by the sum of all exponent log likelihoods (Beerli and Palczewski, 2010).

Spatial Patterns of Genetic Diversity

To capture potentially clinal genetic diversity patterns we used Spatial Genetic Diversity (sGD, Shirk and Cushman, 2011, 2014) to calculate local measures and map continuous gradients of genetic diversity across the study area. We used land cover and elevation to create a resistance surface for use in sGD. We expected less-forested habitat to be increasingly resistant to marten movement. Therefore, we assigned forest a resistance of 1 and non-forest areas were assigned incrementally higher resistance with increasing contrast to forest habitat (e.g., Vergara et al., 2017). The highest resistance value of 20 was given to water bodies, urban areas, and ice fields. Elevation resistance was defined using an inverted Gaussian function as has been done previously for marten in the study area (Wasserman et al., 2010). Also, similar to Wasserman et al. (2010) we used a dispersal distance of 15 km through average resistance values on the landscape to create a neighborhood search radius of 40,000 cost units to define genetic neighborhoods. We calculated allelic richness (Ar), observed heterozygosity (H_o), inbreeding coefficient (F_IS), and genetic neighborhood size (NS) in a 15 km diameter surrounding each sampled individual.

Results

Detections and Distribution

Marten were the most commonly detected mammal at bait stations and were well distributed across the study area with the exception of the Coeur d’Alene Mountains (CDA, Figure 1). We detected marten at only 12% of bait stations in the CDAs, but at 63% of other US bait stations.

Phylogeny Estimation

Within the 493 COI nucleotides there are 10 concordant polymorphic sites that differentiate two haplotypes common to individuals from north of the Clark Fork Valley (North1 and North2) from two haplotypes common to marten south of Highway 200 (South1 and South2) (Figures 1, 2). Among the 143 COI sequences from north of Highway 200, 134 had the most common ‘North1’ haplotype while 9 (7 Purcells, 1 Valhalla, 1 South Selkirks) had a ‘North2’ variant that differed at three sites. Among the 32 sequences from south of the Clark Fork Valley, 22 individuals had a standard ‘South1’ haplotype and 6 (all from the Joe) had a variant ‘South2’ that differed at three sites. No southern haplotypes were detected north of the Clark Fork Valley but 4 southern individuals (3 Joe, 1 CDA) had the North1 haplotype (Figure 1). The North2 haplotype had a variant that differed at 2 non-diagnostic sites (94 and 268) but that also aligned with South1 and South2 at one diagnostic site (325). In the phylogeny (Figure 2), both North1 and North2 haplotypes grouped with sequences Colella et al. (2019) identified as M. americana with 88% bootstrap support. South1 and South2 grouped with sequences Colella et al. (2019) identified as M. caurina with 86% bootstrap support.

FIGURE 2

Figure 2. Bootstrap consensus tree of 493bp COI segment generated for this study and those of Colella et al. (2019).

Population Genetic Analyses

We produced genotypes for 16 microsatellite loci for 241 samples that assigned to 235 marten individuals. We removed the six repeat individual samples leaving us with 235 marten individuals (164 male, 61 female, and 10 unknown gender). Five individuals had incomplete genotypes, with a maximum of 2 loci missing per individual. Forty eight individuals were sampled south of the Clark Fork Valley and 187 to the north. Marker Mvis020 showed a large and highly significant (P < 0.0001) deficit of heterozygous genotypes. Of the 39 heterozygotes at Mvis020 we had identified 37 as female. This suggests that Mvis020 is on the X-chromosome (and that we misidentified 2 females as male), so we removed it from our population genetics analyses. None of the other 15 markers showed significant heterozygote deficits with a Bonferroni correction.

Genotypic differentiation (Genepop Option 3) was significant between all 21 pairs of mountain ranges. These results indicate that our population boundaries (valleys between mountain ranges) generally corresponded as genetic barriers to movement (Figure 3). After correcting for the number of tests (120 pairs of markers tested) we found no evidence of non-random associations of alleles between markers with a Bonferroni correction (“linkage disequilibrium”; Genepop Option 2).

FIGURE 3

Figure 3. Geneclass assignment plot with mixed ancestry individual ScotchmanV1A (black square) shown between northern individual and southern individual clusters.

H_o and H_e were similar across the study system with 15-locus averages ranging lowest in the Valhallas (H_o = 0.51, H_e = 0.54) to highest in the Central Selkirks (H_o = 0.60, H_e = 0.60). The Valhallas has the greatest heterozygote deficit (0.03) of any mountain range. The mean H_o-H_e deficit across all mountain ranges was 0.00, and averages of H_o and H_e across mountain ranges were and were similar for M. americana (H_o = 0.57, H_e = 0.58) and M. caurina (H_o = 0.57, H_e = 0.56) (Table 1).

TABLE 1

Table 1. Observed heterozygosity (H_o) and expected heterozygosity (H_e) for each locus by mountain range, species, and full study area.

The 10 STRUCTURE runs at each value of K converged strongly at each K, with likelihoods peaking at K = 4. The largest lnP(D) increase (1,600) occurred at K = 2 identifying a major split at the Clark Fork Valley. Smaller, less discrete, genetic breaks were identified at K = 3 and K = 4, most notably identifying a sharp split between the South Selkirks and Valhallas. The new split introduced between K = 4 and K = 5 divided 1 of the 4 earlier ancestries in such a way that no individual’s ancestry could be traced predominantly to either new group, thus failing to identify evidence of five genetic populations within the study region (Figures 4, 5).

FIGURE 4

Figure 4. Apportioning of ancestry for 235 M. americana and M. caurina in the highest likelihood runs in STRUCTURE at K = 2 (bottom panel) through 5 (top panel).

FIGURE 5

Figure 5. Estimated log likelihood of K for 50 runs of STRUCTURE (10 each at K = 1–5).

The K = 2 model assigned mixed ancestry to several individuals immediately on either side of the Clark Fork. Most notably, the Cabinet sample ScotchmanV1A was estimated to have roughly equal ancestry north and south of the river, suggesting that it may have one parent from either side (Figure 1). Four of the 9 CDA individuals had substantial estimates of northern ancestry, in the range expected of second or third generation hybrids (3/4 or 7/8 southern ancestry). There was no evidence of recent mixed ancestry in the Joe, where a mean of 99% and minimum of 95% of ancestry was assigned to the southern group (gray in Figure 4).

The 15-locus analysis produced F_ST values similar to (range 0.03–0.32, mean = 0.17), but slightly higher than the 9-locus analysis (range 0.03–0.30, mean = 0.15; Table 2). F_ST values confirmed the Clark Fork Valley split identified by Structure at K = 2. Values across the Clark Fork Valley range from 0.21 to 0.32 with smaller values among northern (0.03–0.09) and southern (0.03) mountain ranges.

TABLE 2

Table 2. F_ST values between mountain ranges.

Both estimation methods for contemporary Ne (NeEstimator and LDNe) and the three allele frequencies used produced similar results (Table 3). Although confidence decreased with lower allele frequencies used, N_e values remained generally similar across tests for most mountain ranges. The CDA (Ne = 31.2) and West Cabinets (Ne = 39.8–45.8) had the lowest and most consistent Ne values across allele frequencies. The remaining mountain ranges had higher and less consistent Ne values across allele frequencies.

TABLE 3

Table 3. Estimates of contemporary effective population size (N_e) and 95% confidence intervals obtained from individuals genotyped at 15 microsatellite loci by mountain range.

Gene Flow Models

Of the 123 individuals tested in GeneClass2 for ancestry from either side of the break, only ScotchmanV1A (Figure 1) had a genotype that was inconsistent with pure ancestry on the side of the Clark Fork where it was sampled (p < 0.0001) (Figure 3).

The evaluation of three models tested (M. americana ⇄ M. caurina, M. americana ⟶ M. caurina, and M. americana ⟵ M. caurina) showed evidence of gene flow occurring under each model but gave most support (probability > 0.99) to M. americana ⟶ M. caurina for both mtDNA and microsatellite data. These results indicate gene flow occurring primarily, but not exclusively, from north to south. Estimates of the fraction of immigrants were relatively low for all models. The M. americana ⟶ M. caurina, model estimated the fraction of immigrants at 1.09 (0.031–5.15; mtDNA) and 0.89 (0.0–3.79; microsatellites), M. americana ⇄ M. caurina at 1.18 (0.03–5.7, mtDNA) and 1.99 (0.0–6.89 microsatellites), and M. americana ⟵ M. caurina at 1.05 (0.01–4.8, mtDNA) and 1.82 (0.0–6.04, microsatellites).

Spatial Patterns of Genetic Diversity

Allelic richness was highest in the Joe (mean = 5.42) followed by South (mean = 4.66) and Central (mean = 4.62) Selkirks (Table 4). The Joe samples all showed high levels of allelic richness while the South and Central Selkirk ranges showed clustering of high allelic richness in the central portions of each range with a notable drop in allelic richness toward the northwestern portion of the Central Selkirks (Figure 6). The lowest allelic richness was in the CDAs (mean = 3.81) and Valhallas (mean = 3.84).

TABLE 4

Table 4. Mean values by mountain range of clinal allelic richness (A), observed heterozygosity (H_o), inbreeding coefficient (F_IS), and genetic neighborhood size (NS) of the sGD analysis.

FIGURE 6

Figure 6. Spatial genetic distance (sGD) analysis results. Clockwise from top left: Allelic richness (Ar), observed heterozygosity (Ho), inbreeding coefficient (F_IS), and genetic neighborhood size (NS).

Observed heterozygosity had little variation and values were similar (0.52–0.60) to those calculated in GenePop based on a priori population definitions. F_IS values ranged from −0.02 to 0.04 with mountain ranges averaging lowest in the CDA (mean F_IS = −0.02) and Joe (mean F_IS = −0.01) and highest in the West Cabinets (mean F_IS = 0.03) and Valhallas (mean F_IS = 0.04). Genetic neighborhood sizes averaged lowest for the CDAs (mean NS = 31.20) and West Cabinets (mean NS = 48.81), was highest in the Joe (mean NS = 106.62), and there was a cluster of high NS values in the US portion of the South Selkirks (Figure 6 and Table 4).

Discussion

Species Fracture Zone and Taxonomic Alignment

We observed a larger genetic break across the approximately 5 km wide Clark Fork Valley (F_ST = 0.21–0.30) than was observed by Kyle and Strobeck (2003) between Golden, British Columbia (northern portion of our study area) and the island of Newfoundland (F_ST = 0.21, Kyle and Strobeck, 2003). In addition to lying 4,000 km to the east, Newfoundland has not had a terrestrial connection with the rest of North America since the end of the Pleistocene glaciations (Pielou, 2008). During the last glacial maximum, approximately 20,000 years before present (Clark et al., 2009), the Clark Fork Valley along with the CDA Mountains were part of the large, 300 km long, Lake Missoula (Pielou, 2008). North of Lake Missoula, our study area was largely engulfed by the Cordilleran Ice Sheet (Pielou, 2008) with the remaining southern portion of our study area being ice free and part of the Pacific Northwest Refugium (Shafer et al., 2010). M. caurina are among the 150 known species or species complexes that have disjunct mesic forest inland and coastal populations in Pacific Northwest (Carstens et al., 2004) due to desertification of the Columbia Plateau (Shafer et al., 2010). During last glacial maximum many species persisted in ice free Pacific Northwest refugia (e.g., Lucid and Cook, 2004) while maintaining genetic connectivity to populations persisting in coastal refugia (Sawyer et al., 2019). The fossil record indicates marten persisted disjunctly during last glacial maximum in western and eastern North America with the eastern group expanding west from the late Pleistocene to late Holocene (Graham and Graham, 1994). Our findings are consistent with the hypothesis (Stone et al., 2002) that M caurina persisted in ice free Pacific Northwest refugia while M. americana persisted in eastern refugia and expanded west as glaciers retreated.

The only individual we identified as likely to have a parent from either side of the break (ScotchmanV1A, Figure 1) was sampled north of the Clark Fork Valley, 8 km from the river itself. However, we identified northern mitochondrial haplotypes in 4 of 32 individuals sampled south of the break, including 3 individuals from the Joe, relatively distant from the break, whereas we observed no southern haplotypes to the north. Colella et al. (2019) found a similar asymmetry in the distribution of mitochondrial haplotypes, suggesting that gene flow from north to south has been taking place for some time. The preponderance of evidence suggests M. caurina likely occupied the ice-free Joe during last glacial maximum and M. americana colonized the northern portion of our study area from the east as glaciers retreated. Marten would not have been able to establish populations in the CDA during last glacial maximum due to the periodic Lake Missoula floods but, after a terrestrial connection was available, gene flow likely initiated and was predominate from M. americana (north) to M. caurina (south). This glacial history has lead us to today where, for practical management purposes, the Clark Fork Valley serves as a geographic boundary between the two lineages with the M. caurina lineage occurring to the south and M. americana lineage occurring to the north of the valley.

Effective Population Size

Effective population sizes (N_e), or number of breeding individuals, are relatively modest considering the size of the mountain ranges. Estimates of clinal genetic neighborhood sizes (NS), also number of breeding individuals, range from 13.9 to 149.1 within a 15 km radius of individual marten. Estimation of clinal NS eliminates the assumption of panmixia from breeding individual estimates and has been shown to estimate larger numbers of breeding individuals than N_e (Shirk and Cushman, 2014). However, developing global estimates from NS values has limited feasibility and since being introduced as a concept (Wright, 1946), NS has been used sparingly by researchers leaving few examples from which to extrapolate relevance. In contrast, N_e is a widely accepted as a variable to quantify a population’s capacity to resist loss of genetic diversity through drift (Shirk and Cushman, 2014). Therefore, it was important to use a variety of techniques to provide opportunities for context and comparison of results.

We used seven (NS plus 2 software program and 3 allele frequencies for N_e) techniques to estimate number of breeding animals for each mountain range. Although our N_e confidence intervals were generally quite large, comparing them to the NS estimates show similar general patterns of abundance across the study area with moderate numbers (hundreds) estimated for most mountain ranges and depressed numbers in the CDAs (N_e = 31.2, mean NS = 31.2) and West Cabinets (N_e = 39.8–45.8, mean NS = 48.81).

Inter-Specific Gene Flow

The individual of recent mixed ancestry is direct evidence that contemporary gene is occurring. At first glance, this appears contrary to the Migrate results showing dominate gene flow from north (M. americana) to south (M. caurina). However, all models predict less than one migrant per generation across the break. We suspect our sampling window was too short and sample sizes too small to detect the few first generation hybrids that may be on the landscape and that the longer term process is predominantly north to south, as the Migrate models suggest. This asymmetry may reflect a higher likelihood for M. americana individuals to disperse due to proximate local population conditions or stochastic events.

M. americana Genetic Structure North of the Fracture Zone

The Valhallas stand out as the most genetically distinct M. americana mountain range with the highest northern F_ST values. Compared to Valhalla, West Cabinet (F_ST = 0.09) values are higher than South Selkirk values (F_ST = 0.03) showing an expected lowering in value with reduced distance. However, the adjacent South Selkirks (F_ST = 0.07) show a relatively high level of genetic disjunct with the Valhallas. N_e estimates varied widely for this range (82.6–872.7) but the relatively low mean allelic diversity (3.81) along the highest mean F_IS value (0.04) in the study area suggest the lower N_e estimates may be more reflective of reality.

The West Cabinets are surprisingly disjunct (F_ST = 0.06) from the South Selkirks given the MacArthur Wildlife Corridor purportedly connects the two ranges (Davidson, 2003). The Purcell-West Cabinet F_ST value (0.05) is slightly lower despite no “official” wildlife corridor connecting them to the West Cabinets. Although modeling work has suggested this corridor is important for grizzly bear and wolverine (Gulo gulo) movement (Cushman et al., 2006; Schwartz et al., 2009), field studies have provided suggestive evidence that fisher (Pekania pennanti) may not use the corridor (Lucid et al., 2019). The high level of marten genetic disjunct casts further doubt on the permeability MacArthur Wildlife Corridor by mammals.

In our study area, M. americana avoid crossing un-forested areas (Cushman et al., 2011, but see Mowat, 2006) and gene flow is primarily a function of elevation (Wasserman et al., 2010). That the mean detection elevation for M. americana in the US portion of our study area was 1,467 m (Lucid et al., 2016) is consistent with the Wasserman et al. (2010) finding of a minimum resistance to gene flow at 1,500 m. This result does not apply to marten continentally as marten are broadly distributed across low elevation boreal forests in Canada (Kyle and Strobeck, 2003), rather this result is limited to our study area and, perhaps, other mountainous portions of marten range of similar latitude. Recognizing climate change may shift this ideal elevation upward suggests marten corridor conservation efforts in our study area may be most effective at elevations at or higher than 1,500 m.

Genetic Structure Along the Fracture Zone

The lowest population estimates in the study area occur along the fracture zone in the CDAs (Ne and NS = 31.2) and West Cabinets (Ne = 39.8–45.8, mean NS = 48.81). F_IS values show disparity, however, with the lowest study area values in the CDAs (mean F_IS = −0.02) and the second highest values in the West Cabinets (mean F_IS = 0.03). Although suggestive of a depressed West Cabinet and genetically robust CDA population, the interpretation of positive F_IS values can be less than straight forward. First, the signal created by null alleles, unrecognized population structure within mountain ranges, or a Wahlund effect (Wahlund, 1928) created by recent immigration could all contribute to increased F_IS estimates. Second, the genetic break across the Clark Fork Valley could lead to differential detectability of Wahlund effects. In a scenario where individuals move between populations with a history of high connectivity (i.e., low F_ST), immigrants have genotypes that are consistent with the population into which they move, and thus their addition to the new population does not create a departure from Hardy-Weinberg expectation. By contrast, when immigrants move between populations with a history of strong isolation (i.e., high Fst, such as the Clark Fork Valley), they have genotypes that are inconsistent with origins in their destination population, and thus create a detectable Wahlund effect. The difference between these scenarios is not the current rate of movement, but the historical rate. Given the possibility of such an artifact, we cannot rule out that the identification of relatively high (but statistically non-significant) F_IS in the West Cabinets reflects greater historical separation, rather than contemporary differences in population processes relative to other mountain ranges that we sampled.

Low numbers of marten in the CDAs are supported both by low field survey detection along with low N_e and NS estimates. Therefore, the negative F_IS values are most likely a Wahuld effect as described in scenario 2 above. Although adjacent areas have been modelled as high quality marten habitat (i.e., Wasserman et al., 2012b), the relatively low elevation CDAs may provide lower quality habitat. Assuming relatively low marten density in the West Cabinets is the correct interpretation of the data, it could plausibly be driven by competition or predation by fisher, which are relatively abundant in the West Cabinets but largely absent from the rest of the study area (Lucid et al., 2019). Marten are currently managed as a single biological unit, M. americana, across the state of Idaho with unlimited state-wide harvest (Idaho Department of Fish and Game (IDFG), 2019). Recognizing the statistical imperfections of our results, the prudent management approach would be to re-evaluate harvest regulations in the CDAs and West Cabinets.

M. caurina Genetic Structure South of the Fracture Zone

The Saint Joe Mountain Range had low F_IS values (mean = −0.01), the highest allelic richness (mean = 5.42), and some of the highest population estimates in the study area (mean NS = 106.62, N_e 138.6–763.9). These metrics suggest the Joe hosts a relatively large and genetically robust marten population. Our M. cqurina sample size (n = 48) is lower than M americana (n = 187) primarily due to the low population size in the CDAs and that the geographic area surveyed is smaller in the southern portion of our study area. Future work to determine how the Joe population is structured with other M. caurina populations to the east and south would be a valuable extension of this study.

Conclusion

Conservation and Management Implications

We identified genetic exchange across a strong genetic break that likely has roots in Pleistocene events, leading to the possibility of substantial local adaptation on either side of the break. Hybridization has been proposed as a conservation problem, particularly in areas where marten augmentation has occurred, with a detrimental effect to M. caurina (Colella et al., 2019). However, the natural erosion of genetic barriers between lineages lead us to recommend natural hybridization not be considered a conservation or management problem within our study area.

The strength of the genetic break that we identified shows that historical connectivity across the Clark Fork Valley has been low enough to produce more or less complete demographic independence: the trajectory of these populations will depend to a much greater extent on local birth and death than on immigration or emigration across the Clark Fork Valley. On this basis we suggest that these populations, which we identify as M. americana and M. caurina, should be treated as separate entities for monitoring programs and management decisions such as developing harvest recommendations.

While there is little doubt that historical connectivity across the Clark Fork Valley has been low, we identified evidence of interbreeding across this break in several individuals sampled near the boundary. This raises the possibility of a recent increase in exchange between highly distinct genetic groups, which may have impacts on fitness that are relevant to management. Genetic assessment of individual ancestry is uniquely suited to non-equilibrium scenarios of suddenly increased connectivity (Paetkau et al., 2009). We therefore recommend ongoing sampling within a band spanning perhaps 20 km on either side of the contact zone, combined with genetic assessment of ancestry to identify individuals likely to have dispersed across the Clark Fork Valley, as well as their offspring (F1s). Ideally, if a number of F1 individuals are detected, it would be possible to conduct more detailed assessments of survival and reproductive output of such individuals relative to individuals whose recent ancestry is from a single side of the genetic break.

In the meantime, we recommend that augmentation or translocation efforts should move marten only from areas with robust local populations from the same genetic group. Prior to augmentations in areas where marten density may be low, such as the CDAs and West Cabinets, in depth study should be undertaken to determine reasons behind depressed populations. Furthermore, harvest in these ranges should be more closely monitored and harvest limits considered.

Long term persistence of marten populations in our study area will be dependent on the maintenance of at least some gene flow amongst mountain ranges. Although we would not dissuade corridor conservation at any elevation, marten in the region that we studied would benefit most from corridor conservation efforts focused on areas of elevation 1,500 m or higher. Monitoring marten population genetic structure over time will be a crucial aspect in evaluating effectiveness of these measures and, as we have demonstrated, inclusion of marten into regional multi-taxa survey programs is an effective way to develop hair and genetic data collections.

Data Availability Statement

The datasets presented in this study are publicly available. The raw microsatellite data are available in the Supplementary Material section. The COI sequence data are available at https://www.ncbi.nlm.nig.gov/genbank (accession numbers MT578095-MT578269).

Ethics Statement

Ethical review and approval was not required for the animal study because all samples were collected with non-invasive methods. Author institutions do not require approval. However, we do follow the ethics policy of the American Society of Mammalogists.

Author Contributions

ML conceived the overall study design, coordinated the data analysis, and wrote the manuscript. ML, LR, and SE conceived of the field survey design and oversaw the sample collection and database management in the United States. GM, DH, and AK conceived of field survey design and oversaw sample collection and database management in the Canada. DP oversaw and SG performed the laboratory analysis. SC, SG, JS, AR, and DP analyzed the data. LS contributed to the study design conception, preparation of spatial datasets, interpretation of the data, and the figure preparation. All authors read and approved the final manuscript.

Funding

This research was supported by the following grants: Columbia Basin Trust, Fish and Wildlife Compensation Program, U.S. Fish and Wildlife Service Wildlife and Sport Fish and Restoration Program Competitive State Wildlife Grant (F12AP01101), Idaho Panhandle National Forest (10-CS-11010400-023), U.S. Forest Service Rocky Mountain Research Station (10-JV-11221633-100 & 08-CS-11221633-194), the Wolverine Foundation, Zoo Boise Conservation Fund, and the National Science Foundation (DEB-1457726).

Conflict of Interest

DH was employed by the Seepanee Ecological Consulting Company. AK was employeed by the Grylloblatta Ecological Consulting Company. DP and SG were employed by the company Wildlife Genetics International. LR was employed but they Idaho Department of Fish and Game during the course of this research. However, LR now is employeed by the Rainforest Ecological Consulting Company.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This US portion of this project was part of the Multi-species Baseline Initiative (https://idfg.idaho.gov/baseline). We thank employees of the following organizations for field assistance: Bureau of Land Management, Coeur d’Alene Tribe of Indians, Idaho Department of Lands, Idaho Panhandle National Forest, Kalispel Tribe of Indians, Seepanee Ecological Consulting, and Washington Department of Fish and Wildlife. We thank Friends of Scotchman Peaks Wilderness, Idaho Conservation League, and Selkirk Outdoor Leadership and Education for facilitating the field work of hundreds of volunteer citizen scientists, Potlatch Corporation and Hancock Forest Management for allowing access to privately held land. We also thank multiple wildlife technicians and volunteers for perseverance in the face of adversity. Kootenay Valley Helicopters Ltd., provided safe and efficient field transportation. B.C. Provincial Parks and B.C. Ministry of Forests, Lands, and Natural Resource Operations provided generous in kind support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00735/full#supplementary-material

References

Allendorf, F. W., Luikart, G., and Aitken, S. N. (2013). Conservation and the Genetics of Populations. Oxford: Blackwell.

Google Scholar

Beerli, P., and Felsenstein, J. (1999). Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics 152, 763–773.

Google Scholar

Beerli, P., and Felsenstein, J. (2001). Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proc. Natl. Acad. Sci. U.S.A. 98, 4563–4568. doi: 10.1073/pnas.081068098

PubMed Abstract | CrossRef Full Text | Google Scholar

Beerli, P., and Palczewski, M. (2010). Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics 185, 313–326. doi: 10.1534/genetics.109.112532

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloomquist, E. W., Lemey, P., and Suchard, M. A. (2010). Three roads diverged? Routes to phylogeographic inference. Trends Ecol. Evol. 25, 626–632. doi: 10.1016/j.tree.2010.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunsfeld, S. J., Miller, T. R., and Carstens, B. C. (2007). Insights into the biogeography of the Pacific Northwest of North America: evidence from the phylogeography of Salix melanopsis. Syst. Bot. 32, 129–139. doi: 10.1600/036364407780360094

CrossRef Full Text | Google Scholar

Brunsfeld, S. J., and Sullivan, J. (2005). A multi-compartmented glacial refugium in the northern Rocky Mountains: evidence from the phylogeography of Cardamine constancei (Brassicaceae). Conserv. Genet. 6, 895–904. doi: 10.1007/s10592-005-9076-7

CrossRef Full Text | Google Scholar

Carstens, B. C., Stevenson, A. L., Degenhardt, J. D., and Sullivan, J. (2004). Testing nested phylogenetic and phylogeographic hypotheses in the Plethodon vandykei species group. Syst. Biol. 53, 781–792. doi: 10.1080/10635150490522296

PubMed Abstract | CrossRef Full Text | Google Scholar

Cartens, B. C., Brunsfeld, S. J., Demboski, J. R., Good, J. M., and Sullivan, J. (2005). Investigating the evolutionary history of the Pacific Northwest mesic forest ecosystem: hypothesis testing within a comparative phylogeographic framework. Evolution 59, 1639–1652. doi: 10.1111/j.0014-3820.2005.tb01815.x

CrossRef Full Text | Google Scholar

Clark, C. D., Hughes, A. L., Greenwood, S. L., Spagnolo, M., and Ng, F. S. (2009). Size and shape characteristics of drumlins, derived from a large sample, and associated scaling laws. Q. Sci. Rev. 28, 677–692. doi: 10.1016/j.quascirev.2008.08.035

CrossRef Full Text | Google Scholar

Colella, J. P., Wilson, R. E., Talbot, S. L., and Cook, J. A. (2019). Implications of introgression for wildlife translocations: the case of North American martens. Conserv. Genet. 20, 153–166. doi: 10.1007/s10592-018-1120-5

CrossRef Full Text | Google Scholar

Cushman, S. A., McKelvey, K. S., Hayden, J., and Schwartz, M. K. (2006). Gene flow in complex landscapes: testing multiple hypotheses with causal modeling. Am. Nat. 168, 486–499. doi: 10.1086/506976

PubMed Abstract | CrossRef Full Text | Google Scholar

Cushman, S. A., Raphael, M. G., Ruggiero, L. F., Shirk, A. S., Wasserman, T. N., and O’Doherty, E. C. (2011). Limiting factors and landscape connectivity: the American marten in the Rocky Mountains. Landscape Ecol. 26, 1137. doi: 10.1007/s10980-011-9645-8

CrossRef Full Text | Google Scholar

Davidson, D. K. (2003). “Innovative partnerships that address highway impacts to wildlife habitat connectivity in the Northern Rockies,” in Proceedings of the International Conference on Ecology and Transportation (ICOET 2003), New York, NY.

Google Scholar

Davis, C. S., and Strobeck, C. (1998). Isolation, variability, and cross-species amplification of polymorphic microsatellite loci in the family Mustelidae. Mol. Ecol. 7, 1776–1778. doi: 10.1046/j.1365-294x.1998.00515.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dawson, N. G., Colella, J. P., Small, M. P., Stone, K. D., Talbot, S. L., and Cook, J. A. (2017). Historical biogeography sets the foundation for contemporary conservation of martens (genus Martes) in northwestern North America. J. Mammal. 98, 715–730. doi: 10.1093/jmammal/gyx047

CrossRef Full Text | Google Scholar

Dawson, N. G., and Cook, J. A. (2012). Behind the Genes: Diversification of North American martens (Martes americana and M. caurina). Biology and Conservation of Martens, Sables, and Fishers: A New Synthesis. Ithaca, NY: Cornell University Press, 23–38.

Google Scholar

Demboski, J. R., and Sullivan, J. (2003). Extensive mtDNA variation within the yellow-pine chipmunk, Tamias amoenus (Rodentia: Sciuridae), and phylogeographic inferences for northwest North America. Mol. Phylog. Evol. 26, 389–408. doi: 10.1016/s1055-7903(02)00363-9

CrossRef Full Text | Google Scholar

Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., and Ovenden, J. R. (2014). NEESTIMATOR v2: re-implementation of software for the estimation of contemporary effective population size (N_e) from genetic data. Mol. Ecol. Resour. 14, 209–214. doi: 10.1111/1755-0998.12157

PubMed Abstract | CrossRef Full Text | Google Scholar

Felsenstein, J. (1985). Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39, 783–791. doi: 10.1111/j.1558-5646.1985.tb00420.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Flemming, M. A., Ostrander, E. A., and Cook, J. A. (1999). Microsatellite markers for American mink (Mustela vison) and ermine (Mustela erminea). Mol. Ecol. 8, 1352–1354.

Google Scholar

Frankham, R., Bradshaw, C. J., and Brook, B. W. (2014). Genetics in conservation management: revised recommendations for the 50/500 rules, Red List criteria and population viability analyses. Biol. Conserv. 170, 56–63. doi: 10.1016/j.biocon.2013.12.036

CrossRef Full Text | Google Scholar

Frankham, R., Brook, B. W., Bradshaw, C. J., Traill, L. W., and Spielman, D. (2013). 50/500 rule and minimum viable populations: response to Jamieson and Allendorf. Trends Ecol. Evol. 28, 4–187.

Google Scholar

Franklin, I. R. (1980). “Evolutionary change in small populations,” in Conservation Biology: An Evolutionary-Ecological Perspective, eds M. E. Soulé and B. A., Wilcox (Sutherland, MA: Princeton University Press), 135–149.

Google Scholar

Graham, R. W., and Graham, M. A. (1994). “Late quaternary distribution of Martes in North America,” in Martens, Sables, and Fishers: Biology and Conservation, eds S. W. Buskirk, A. S. Harestad, M. G. Raphael, and R. A., Powell (Ithaca, NY: Cornell University Press), 26–58.

Google Scholar

Idaho Department of Fish and Game (IDFG) (2019). Upland Game, Furbearer, and Turkey Seasons and Rules. Boise, Id: Idaho Department of Fish and Game.

Google Scholar

Jamieson, I. G., and Allendorf, F. W. (2012). How does the 50/500 rule apply to MVPs? Trends Ecol. Evol. 27, 578–584. doi: 10.1016/j.tree.2012.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Jordan, M. J., Higley, J., Matthews, S. M., Rhodes, O. E., Schwartz, M. K., Barrett, R. H., et al. (2007). Development of 22 new microsatellite loci for fishers (Martes pennanti) with variability results from across their range. Mol. Ecol. Res. 7, 797–801. doi: 10.1111/j.1471-8286.2007.01708.x

CrossRef Full Text | Google Scholar

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. J. Am. Statist. Assoc. 90, 773–795.

Google Scholar

Kendall, K. C., Stetz, J. B., Boulanger, J., Macleod, A. C., Paetkau, D., and White, G. C. (2009). Demography and genetic structure of a recovering grizzly bear population. J. Wildl Manag. 73, 3–16.

Google Scholar

Kimura, M. (1980). A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16, 111–120. doi: 10.1007/bf01731581

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirk, T. A., and Zielinski, W. J. (2009). Developing and testing a landscape habitat suitability model for the American marten (Martes americana) in the Cascades mountains of California. Landscape Ecol. 24, 759–773. doi: 10.1007/s10980-009-9349-5

CrossRef Full Text | Google Scholar

Kortello, A., Hausleitner, D., and Mowat, G. (2019). Mechanisms influencing the winter distribution of wolverine Gulo gulo luscus in the southern Columbia Mountains, Canada. Wildlife Biol. 2019, 1–13.

Google Scholar

Kyle, C. J., and Strobeck, C. (2003). Genetic homogeneity of Canadian mainland marten populations underscores the distinctiveness of Newfoundland pine martens (Martes americana atrata). Canad. J. Zool. 81, 57–66. doi: 10.1139/z02-223

CrossRef Full Text | Google Scholar

Lande, R., Barrowclough, G. F., and Soulé, M. E. (1987). Viable Populations for Conservation. Effective Population Size, Genetic Variation and Their Use in Population Management. Cambridge: Cambridge University Press, 87–123.

Google Scholar

Lindenmayer, D. B., and Likens, G. E. (2010). The science and application of ecological monitoring. Biol. Conserv. 143, 1317–1328. doi: 10.1016/j.biocon.2010.02.013

CrossRef Full Text | Google Scholar

Lucid, M., Robinson, L., and Ehlers, S. (2016). Multi-Species Baseline Initiative Project Report: 2010–2014. Boise, ID: Idaho Department of Fish and Game.

Google Scholar

Lucid, M. K., and Cook, J. A. (2004). Phylogeography of Keen’s mouse (Peromyscus keeni) in a naturally fragmented landscape. J. Mammal. 85, 1149–1159. doi: 10.1644/brb-218.1

CrossRef Full Text | Google Scholar

Lucid, M. K., Ehlers, S., Robinson, L., and Cushman, S. A. (2018a). Beer, brains, and brawn as tools to describe terrestrial gastropod species richness on a montane landscape. Ecosphere 9:e02535. doi: 10.1002/ecs2.2535

CrossRef Full Text | Google Scholar

Lucid, M. K., Rankin, A., Sullivan, J., Robinson, L., Ehlers, S., and Cushman, S. (2019). A carnivores’ oasis? An isolated fisher (Pekania pennanti) population provides insight on persistence of a metapopulation. Conserv. Genet. 20, 585–596. doi: 10.1007/s10592-019-01160-w

CrossRef Full Text | Google Scholar

Merriam, C. H. (1890). Descriptions of twenty-six new species of North American mammals. N. Am. Fauna 1890, 1–55.

Google Scholar

Mills, L. S., and Allendorf, F. W. (1996). The one-migrant-per-generation rule in conservation and management. Conserv. Biol. 10, 1509–1518. doi: 10.1046/j.1523-1739.1996.10061509.x

CrossRef Full Text | Google Scholar

Minin, V., Abdo, Z., Joyce, P., and Sullivan, J. (2003). Performance-based selection of likelihood models for phylogeny estimation. Syst. Biol. 52, 674–683. doi: 10.1080/10635150390235494

PubMed Abstract | CrossRef Full Text | Google Scholar

Moriarty, K. M., Epps, C. W., and Zielinski, W. J. (2016). Forest thinning changes movement patterns and habitat use by Pacific marten. J. Wildlife Manag. 80, 621–633. doi: 10.1002/jwmg.1060

CrossRef Full Text | Google Scholar

Mowat, G. (2006). Winter habitat associations of American martens Martes americana in interior wet-belt forests. Wildlife Biol. 12, 51–62.

Google Scholar

Mowat, G., and Paetkau, D. (2002). Estimating marten Martes americana population size using hair capture and genetic tagging. Wildlife Biol. 8, 201–210.

Google Scholar

Nielson, M., Lohman, K., Daugherty, C. H., Allendorf, F. W., Knudsen, K. L., Sullivan, J., et al. (2006). Allozyme and mitochondrial DNA variation in the tailed frog (Anura: Ascaphus): the influence of geography and gene flow. Herpetologica 62, 235–258. doi: 10.1655/0018-0831(2006)62[235:aamdvi]2.0.co;2

CrossRef Full Text | Google Scholar

Paetkau, D. (2003). An empirical exploration of data quality in DNA based population inventories. Mol. Ecol. 12, 1375–1387. doi: 10.1046/j.1365-294x.2003.01820.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Paetkau, D., Calvert, W., Stirling, I., and Strobeck, C. (1995). Microsatellite analysis of population structure in Canadian polar bears. Mol. Ecol. 4, 347–354. doi: 10.1111/j.1365-294x.1995.tb00227.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Paetkau, D., Slade, R., Burden, M., and Estoup, A. (2004). Genetic assignment methods for the direct, real−time estimation of migration rate: a simulation−based exploration of accuracy and power. Mol. Ecol. 13, 55–65. doi: 10.1046/j.1365-294x.2004.02008.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Paetkau, D., Vázquez-Domínguez, E., Tucker, N. I., and Moritz, C. (2009). Monitoring movement into and through a newly planted rainforest corridor using genetic analysis of natal origin. Ecol. Manag. Restor. 10, 210–216. doi: 10.1111/j.1442-8903.2009.00490.x

CrossRef Full Text | Google Scholar

Pielou, E. C. (2008). After The Ice Age: The Return of Life to Glaciated North America. Chicago, IL: University of Chicago Press.

Google Scholar

Piry, S., Alapetite, A., Cornuet, J. M., Paetkau, D., Baudouin, L., and Estoup, A. (2004). GENECLASS2: a software for genetic assignment and first-generation migrant detection. J. Hered. 95, 536–539. doi: 10.1093/jhered/esh074

PubMed Abstract | CrossRef Full Text | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959.

Google Scholar

Rankin, A. M., Wilke, T., Lucid, M., Leonard, W., Espíndola, A., Smith, M. L., et al. (2019). Complex interplay of ancient vicariance and recent patterns of geographical speciation in north-western North American temperate rainforests explains the phylogeny of jumping slugs (Hemphillia spp.). Biol. J. Linnean Soc. 127, 876–889. doi: 10.1093/biolinnean/blz040

CrossRef Full Text | Google Scholar

Raymond, M., and Rousset, F. (1995). GENEPOP on the Web (Version 3.4). Available online at: http://wbiomed.curtin.edu.au/genepop/ (accessed March 20, 2020).

Google Scholar

Robinson, L., Cushman, S. A., and Lucid, M. K. (2017). Winter bait stations as a multispecies survey tool. Ecol. Evol. 7, 6826–6838. doi: 10.1002/ece3.3158

PubMed Abstract | CrossRef Full Text | Google Scholar

Sawyer, Y. E., MacDonald, S. O., Lessa, E. P., and Cook, J. A. (2019). Living on the edge: exploring the role of coastal refugia in the Alexander Archipelago of Alaska. Ecol. Evol. 9, 1777–1797. doi: 10.1002/ece3.4861

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, M. K., Copeland, J. P., Anderson, N. J., Squires, J. R., Inman, R. M., McKelvey, K. S., et al. (2009). Wolverine gene flow across a narrow climatic niche. Ecology 90, 3222–3232. doi: 10.1890/08-1287.1

CrossRef Full Text | Google Scholar

Shafer, A. B., Cullingham, C. I., Cote, S. D., and Coltman, D. W. (2010). Of glaciers and refugia: a decade of study sheds new light on the phylogeography of northwestern North America. Mol. Ecol. 19, 4589–4621. doi: 10.1111/j.1365-294x.2010.04828.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirk, A. J., and Cushman, S. A. (2011). sGD: software for estimating spatially explicit indices of genetic diversity. Mol. Ecol. Resour. 11, 922–934. doi: 10.1111/j.1755-0998.2011.03035.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirk, A. J., and Cushman, S. A. (2014). Spatially-explicit estimation of Wright’s neighborhood size in continuous populations. Front. Ecol. Evol. 2:62.

Google Scholar

Shirk, A. J., Wasserman, T. N., Cushman, S. A., and Raphael, M. G. (2012). “Scale dependency of American marten (Martes americana) habitat relations [Chapter 12],” in Biology and Conservation of Martens, Sables, and Fishers: A New Synthesis, eds K. B. Aubry, W. J. Zielinski, M. G. Raphael, G. Proulx, and S. W., Buskirk (Ithaca, NY: Cornell University Press), 269–283. doi: 10.7591/9780801466076-015

CrossRef Full Text | Google Scholar

Slauson, K. M., Zielinski, W. J., and Schwartz, M. K. (2017). Ski areas affect Pacific marten movement, habitat use, and density. J. Wildlife Manag. 81, 892–904. doi: 10.1002/jwmg.21243

CrossRef Full Text | Google Scholar

Small, M. P., Stone, K. D., and Cook, J. A. (2003). American marten (Martes americana) in the Pacific Northwest: population differentiation across a landscape fragmented in time and space. Mol. Ecol. 12, 89–103. doi: 10.1046/j.1365-294x.2003.01720.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Soulé, M. E. (1980). “Thresholds for survival: maintaining fitness and evolutionary potential,” in Conservation Biology: An Evolutionary-Ecological Perspective, eds M. E. Soulé and B. A., Wilcox (Sunderland, MA: Sinauer Associates), 151–169.

Google Scholar

Stone, K. D., Flynn, R. W., and Cook, J. A. (2002). Post-glacial colonization of northwestern North America by the forest-associated American marten (Martes americana. Mammalia: Carnivora: Mustelidae). Mol. Ecol. 11, 2049–2063. doi: 10.1046/j.1365-294x.2002.01596.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Swofford, D. L. (2003). PAUP^∗: Phylogenetic Analysis Using Parsimony (^∗and Other Methods). Version 4. Sunderland, MA: Sinauer Associates, Inc.

Google Scholar

Thompson, I. D. (1991). Could marten become the spotted owl of eastern Canada? Forestry Chronicle 67, 136–140. doi: 10.5558/tfc67136-2

CrossRef Full Text | Google Scholar

Tucker, N. I., and Simmons, T. (2009). Restoring a rainforest habitat linkage in north Queensland: Donaghy’s Corridor. Ecol. Manag. Restor. 10, 98–112. doi: 10.1111/j.1442-8903.2009.00471.x

CrossRef Full Text | Google Scholar

Vergara, M., Cushman, S. A., and Ruiz-González, A. (2017). Ecological differences and limiting factors in different regional contexts: landscape genetics of the stone marten in the Iberian Peninsula. Landscape Ecol. 32, 1269–1283. doi: 10.1007/s10980-017-0512-0

CrossRef Full Text | Google Scholar

Wahlund, S. (1928). Zusammensetzung von populationen und korrelationserscheinungen vom standpunkt der vererbungslehre aus betrachtet. Hereditas 11, 65–106. doi: 10.1111/j.1601-5223.1928.tb02483.x

CrossRef Full Text | Google Scholar

Waples, R. S., and Do, C. H. I. (2008). LDNE: a program for estimating effective population size from data on linkage disequilibrium. Mol. Ecol. Resour. 8, 753–756. doi: 10.1111/j.1755-0998.2007.02061.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wasserman, T. N., Cushman, S. A., Schwartz, M. K., and Wallin, D. O. (2010). Spatial scaling and multi-model inference in landscape genetics: Martes americana in northern Idaho. Landscape Ecol. 25, 1601–1612. doi: 10.1007/s10980-010-9525-7

CrossRef Full Text | Google Scholar

Wasserman, T. N., Cushman, S. A., Shirk, A. S., Landguth, E. L., and Littell, J. S. (2012a). Simulating the effects of climate change on population connectivity of American marten (Martes americana) in the northern Rocky Mountains, USA. Landscape Ecol. 27, 211–225. doi: 10.1007/s10980-011-9653-8

CrossRef Full Text | Google Scholar

Wasserman, T. N., Cushman, S. A., and Wallin, D. O. (2012b). Multi Scale Habitat Relationships of Martes Americana in Northern Idaho, USA. Fort Collins, CO: U.S. Department of Agriculture.

Google Scholar

Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. doi: 10.1111/j.1558-5646.1984.tb05657.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wright, S. (1946). Isolation by distance under diverse systems of mating. Genetics 31, 39–59.

Google Scholar

Young, M. K., McKelvey, K. S., Pilgrim, K. L., and Schwartz, M. K. (2013). DNA barcoding at riverscape scales: assessing biodiversity among fishes of the genus Cottus (Teleostei) in northern Rocky Mountain streams. Mol. Ecol. Resour. 13, 583–595. doi: 10.1111/1755-0998.12091

PubMed Abstract | CrossRef Full Text | Google Scholar