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 km2 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 (Ne) for each mountain range, clinal genetic neighborhood sizes (NS), calculated observed (Ho) and expected (He) heterozygosity, fixation index (FST), and clinal measures of allelic richness (Ar), Ho, and inbreeding coefficient (FIS). 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 (FST 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.

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 2 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

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 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 2 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.

Study Area
The study area consisted of a 53,474 km 2 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.

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).

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 pseudoreplicates (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 nonrandom 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., nonrandom 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. 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 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 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, 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.

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)

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).  Table 1).
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).
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).

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
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).
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).

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

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 Estimates above the diagonal are based on 9 of the 11 markers used by Kyle and Strobeck (2003). Estimated below the diagonal are based on all 15 microsatellites analyzed for this study.
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 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.  Estimates were done using both *LDNe and **NeEstimator and three sets of lowest allele frequency used (0.05, 0.02, 0.01).

M. americana
Frontiers in Genetics | www.frontiersin.org 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 . 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 , 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 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 . 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.

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 . 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.