Pleistocene–Holocene vicariance, not Anthropocene landscape change, explains the genetic structure of American black bear (Ursus americanus) populations in the American Southwest and northern Mexico

Abstract The phylogeography of the American black bear (Ursus americanus) is characterized by isolation into glacial refugia, followed by population expansion and genetic admixture. Anthropogenic activities, including overharvest, habitat loss, and transportation infrastructure, have also influenced their landscape genetic structure. We describe the genetic structure of the American black bear in the American Southwest and northern Mexico and investigate how prehistoric and contemporary forces shaped genetic structure and influenced gene flow. Using a suite of microsatellites and a sample of 550 bears, we identified 14 subpopulations organized hierarchically following the distribution of ecoregions and mountain ranges containing black bear habitat. The pattern of subdivision we observed is more likely a product of postglacial habitat fragmentation during the Pleistocene and Holocene, rather than a consequence of contemporary anthropogenic barriers to movement during the Anthropocene. We used linear mixed‐effects models to quantify the relationship between landscape resistance and genetic distance among individuals, which indicated that both isolation by resistance and geographic distance govern gene flow. Gene flow was highest among subpopulations occupying large tracts of contiguous habitat, was reduced among subpopulations in the Madrean Sky Island Archipelago, where montane habitat exists within a lowland matrix of arid lands, and was essentially nonexistent between two isolated subpopulations. We found significant asymmetric gene flow supporting the hypothesis that bears expanded northward from a Pleistocene refugium located in the American Southwest and northern Mexico and that major highways were not yet affecting gene flow. The potential vulnerability of the species to climate change, transportation infrastructure, and the US–Mexico border wall highlights conservation challenges and opportunities for binational collaboration.


Abstract
The phylogeography of the American black bear (Ursus americanus) is characterized by isolation into glacial refugia, followed by population expansion and genetic admixture.
Anthropogenic activities, including overharvest, habitat loss, and transportation infrastructure, have also influenced their landscape genetic structure. We describe the genetic structure of the American black bear in the American Southwest and northern Mexico and investigate how prehistoric and contemporary forces shaped genetic structure and influenced gene flow. Using a suite of microsatellites and a sample of 550 bears, we identified 14 subpopulations organized hierarchically following the distribution of ecoregions and mountain ranges containing black bear habitat. The pattern of subdivision we observed is more likely a product of postglacial habitat fragmentation during the Pleistocene and Holocene, rather than a consequence of contemporary anthropogenic barriers to movement during the Anthropocene. We used linear mixed-effects models to quantify the relationship between landscape resistance and genetic distance among individuals, which indicated that both isolation by resistance and geographic distance govern gene flow. Gene flow was highest among subpopulations occupying large tracts of contiguous habitat, was reduced among subpopulations in the Madrean Sky Island Archipelago, where montane habitat exists within a lowland matrix of arid lands, and was essentially nonexistent between two isolated subpopulations. We found significant asymmetric gene flow supporting the hypothesis that bears expanded northward from a Pleistocene refugium located in the American Southwest and northern Mexico and that major highways were not yet affecting gene flow. The potential vulnerability of the species to climate change, transportation infrastructure, and the US-Mexico border wall highlights conservation challenges and opportunities for binational collaboration.

| INTRODUC TI ON
The Pleistocene epoch (2.6-0.012 mya) represents a geologic period characterized by massive climatic fluctuations that drove dynamic glacial-interglacial cycles with profound effects on the global distribution and genetic structure of flora and fauna (Hofreiter & Stewart, 2009). Glacial advance contracted species' ranges into refugial pockets of habitat where isolation, selection, and genetic drift resulted in genetic differentiation among populations. Upon glacial recession, species expanded out of their respective refugia into their current distribution resulting in latitudinal patterns of species assemblages, genetic structure, and areas of admixture between formerly isolated populations (Lomolino et al., 1989;Puckett et al., 2015;Shafer et al., 2010). The phylogeographic influence of these glacial-interglacial dynamics has been observed for vagile species like the gray wolf (Canis lupus; Weckworth et al., 2010) and red fox (Vulpes vulpes; Aubry et al., 2009) and for more habitatrestricted species like the woodland caribou (Rangifer tarandus caribou; Klütsch et al., 2012) and American marten (Martes americana; Stone et al., 2002). The location of refugia is highly dependent on the life history of the organism.
The American black bear (Ursus americanus; hereafter, black bear) is a large omnivore endemic to the forests of North America.
Its distribution and genetic structure have been an ebb and flow of isolation and admixture events dictated by glacial tides of the Pleistocene (Puckett et al., 2015). Mitochondrial and nuclear data indicate that black bears were last isolated during the Last Glacial Maximum (LGM) ~26.5 kya and had contracted into three glacial refugia located in Beringia, the Pacific Northwest, and the American Southeast, and a fourth hypothesized refugium in the southwestern United States and northern Mexico (hereafter, the Southwest; Puckett et al., 2015;Varas-Nelson, 2010). As glaciers receded (~20 kya), black bears expanded out of their respective refugia resulting in admixture among populations in west-central and east-central North America and the formation of region-specific subpopulations (Pelletier et al., 2011;Puckett et al., 2015). Black bears across their northern range are genetically diverse and inhabit a large, contiguous landscape with a genetic structure that is consistent with isolation by distance due to female-biased philopatry (Pelletier et al., 2011;Pelletier et al., 2017). Despite support for a fourth glacial refugium in the Southwest (Puckett et al., 2015;Varas-Nelson, 2010), there remains some uncertainty as these inferences are based on limited geographic sampling, particularly in the southeast portion of New Mexico (Onorato et al., 2007;Onorato, Hellgren, Van Den Bussche, & Doan-Crider, 2004). Furthermore, conflicting inferences regarding the genetic structure of black bears have been made in Arizona and New Mexico, where bear populations have been reported as both genetically structured and admixed (Atwood et al., 2011;Varas-Nelson, 2010;Winslow, 2012).
The potential existence of a southwestern refugium for black bears is independently supported by paleoecological reconstruction of woodrat (Neotoma spp.) paleomiddens .
Woodrats collect plant material and pollen blows in or adheres to plants and the macrofossils and pollen become encased in crystallized woodrat urine, or amberat, that form an indurated paleomidden (Spaulding et al., 1990). Investigations of these paleomiddens reveal information about the relative abundance, distribution, and species composition of prehistoric plant communities, which in turn enable assessments of climatic and plant community change   Betancourt et al., 1990, Chp. 21;Holmgren et al., 2006;McAuliffe & Van Devender, 1998;Onorato et al., 2003). Conversely, areas farther north, in northern Arizona and New Mexico and southern Colorado and Utah, contained less hospitable habitat including montane glaciers, tundra, and taiga as well as forests dominated by yellow pine (Pinus spp.), limber pine (P. flexilis), and lodgepole pine (P. contorta). These regions currently harbor more contiguous black bear habitat. As the Holocene aridified, habitats preferred by bears either shifted up in elevation, such as in have been rendered into small, isolated populations more susceptible to genetic drift, eroding genetic diversity (Hooker et al., 2015;Murphy et al., 2017;Murphy et al., 2018). In Florida, major roads have heightened mortality (McCown et al., 2009) and acted as semipermeable barriers, that when coupled with urbanization, fragmented bear habitat, decreased connectivity, and caused appreciable genetic structure among subpopulations (Dixon et al., 2007;Karelus et al., 2017). In the Lower Mississippi Alluvial Valley of Louisiana, human-caused mortality combined with extensive habitat loss and fragmentation forced black bears into a patchwork of small populations isolated by anthropogenic activities; active translocations are underway to help restore bear populations there (Murphy et al., 2018). In several states, roads have been shown to influence movements and patterns of habitat selection, as bears either avoid roads or select areas farther from roads, and in some regions, roads created genetic substructure by acting as filters to bear movement (Cushman & Lewis, 2010;Dixon et al., 2007;Gould et al., 2019;Hiller et al., 2015;Short Bull et al., 2011).
In the Southwest, limited geographic and genetic sampling has obscured the influence of prehistoric and contemporary ecological processes that shape genetic structure and govern gene flow of black bear populations. Our aim was to fill a crucial gap regarding the large-scale population genetic structure of the American black bear by using a suite of microsatellite loci to characterize the genetic profile of 550 individual bears sampled across the Southwest.
We focused on two hypotheses. First, we hypothesized that the current genetic structure could be a consequence of Pleistocene-Holocene vicariance whereby bears occupied forest refugia during the LGM, but then followed changes in the distribution of forests as the Holocene dried and warmed (Pleistocene-Holocene Vicariance Hypothesis). If true, we predicted that bears occupying contiguous forests should be relatively closely related and exhibit little genetic substructure. Bears in the Madrean Sky Islands should be more genetically structured (Atwood et al., 2011;Varas-Nelson, 2010) and should show evidence of gene flow characterized by isolation by resistance. Finally, there should be a pronounced asymmetric pattern of gene flow from south to north. Second, we hypothesized that the genetic substructure of bears could be dominated by anthropogenic activities typifying the Anthropocene (i.e., the concept that we are living in a time when human activities have significant effecs on the global environment). If true, we predicted that the influence of major highways would be manifested by bears being more closely related on the same side of a highway and more distantly related on opposite sides. This should occur irrespective of the intervening habitat matrix. Although highways may not be barriers, they should act as semipermeable filters that influence gene flow (Anthropocene Filter Hypothesis). Expectations from these hypotheses may not be mutually exclusive, but we envision that the strength of evidence gathered through our analysis will expose their relative importance and that the insight gained will illuminate the processes that helped shape the phylogeography and present-day genetic structure of black bears in the Southwest. Our findings will also aid in the conservation and management of black bears by identifying genetically isolated populations and the landscape features promoting or inhibiting genetic connectivity among black bear populations.

| Study area
We conducted our study in the southwestern United States (Arizona, Colorado, New Mexico, Texas, and Utah) and northern Coahuila de Zaragoza, Mexico (Figure 1). Orography and climate vary drastically throughout the Southwest with elevation ranging from 21 m at the southwest corner of Arizona to 4155 m-high peaks in southern Colorado. The desert and grassland communities receive the majority of the ~100-300 mm of annual precipitation during the July to October monsoon season with mean monthly maximum temperatures for July from 1961 to 1990 ranging from ~16 to 40°C (Davey et al., 2006(Davey et al., , 2007a(Davey et al., , 2007b. The forest communities receive the majority of the ~300 to 1250 mm of annual precipitation during the winter with mean monthly minimum temperature for January from 1961 to 1990 ranging from ~ −20 to −4°C (Davey et al., 2006(Davey et al., , 2007a(Davey et al., , 2007b (Omernik & Griffith, 2014). Biotic communities at higher elevations and latitudes consist of Petran Subalpine and Petran Montane conifer forests transitioning to mid-elevation Great Basin Conifer and Madrean Evergreen woodlands (Brown, 1994).
Large expanses of low-elevation valleys are composed of biotic communities such as Plains and Great Basin Grassland, Semidesert Grassland, and the Great Basin, Chihuahuan, and Sonoran desert scrub. These areas comprise a low elevation "sea" that is not typically used by black bears, isolating them on montane 'sky islands' (Brown, 1994;Hellgren et al., 2005;Olson et al., 2001).

| Sample and marker selection
We collected genetic samples from individual black bears through hunter harvest, live-capture, noninvasive genetic sampling, and vehicle collisions. We attempted to sample ≥ 25 individuals from each mountain range to obtain an adequate representation of allele frequency and diversity within each assumed subpopulation (Hale et al., 2012). Despite there being evidence that bears from west Texas and northern Coahuila de Zaragoza, Mexico, have genetic signatures more similar to bears in the southeast United States (Onorato et al., 2007; Durnin et al., 2007;Ostrander et al., 1993;Paetkau et al., 1998;Paetkau & Strobeck, 1995;Taberlet et al., 1997). Wildlife Genetics International in Nelson, British Columbia, Canada, generated all genotypes. Detailed laboratory protocols for microsatellite amplification and assignment error can be found in Paetkau (2003) and Gould et al. (2018). We obtained permits under the Convention on International Trade in Endangered Species (Export Permits 12US86418A/9, 13US19950B/9, 13US199551B/9, 14US43944B/9, 15US61420B/9, 15US69493B/9, 15US69502B/9) to export samples to Canada for analysis. Our research was authorized by the New Mexico Department of Game

| Describing genetic structure and estimating gene flow
We conducted all genetic analyses in program r version 3.5.2 unless otherwise specified (R Core Team, 2018). We tested for linkage disequilibrium (LD) using the r package genepop and null alleles and deviations from Hardy-Weinberg equilibrium (HWE) using the r package popgenreport (version 3.0.0; Adamack & Gruber, 2014; version 1.0.5; Rousset et al., 2017). We conducted a significance test for null alleles by assessing if a bootstrapped 95% confidence intervals (CI) for each locus overlapped zero, whereby overlap would indicate that the frequency of null alleles does not differ significantly from zero. We applied a Bonferroni correction (α = 0.05 divided by the number of F I G U R E 1 Distribution of genetic samples and subpopulations of American black bears (Ursus americanus) in the American Southwest and northern Mexico. geneland identified 6 and 14 subpopulations using the uncorrelated (polygons) and correlated (symbols) allele frequency models, respectively. The 6 larger subpopulations are named clusters.
pairwise comparisons for each test) of α ≤ 0.0005 (LD) and α ≤ 0.003 (HWE) to reduce the likelihood of a false-positive significance test, and we generated allele-frequency statistics for each locus in popgenreport. For each identified subpopulation, we quantified genetic diversity using unadjusted private alleles (A P ) and private alleles using rarefaction (A PR ), which accounts for differences in sample size among subpopulations, using hp-rare v1.0 (Kalinowski, 2004(Kalinowski, , 2005. We also quantified genetic diversity by estimating expected (H E ) and observed (H O ) heterozygosity and allelic richness using rarefaction (A R ) with the r package diversity (version 1.9.9; Keenan et al., 2013). We used diversity to calculate deviations from random mating (F IS ) and genetic differentiation among subpopulations (F ST ) along with their 95% CI based on 1000 bootstrap iterations. We classified values of F ST from 0.05-0.14, 0.15-0.24, and ≥0.25 as moderate, high, and very-high differentiation, respectively, and considered differentiation to be biologically meaningful if the lower CI was ≥0.05 (Hartl & Clark, 1997).
We used two Bayesian clustering programs to characterize population structure, geneland and structure (Guillot et al., 2005;Pritchard, 2000). Both programs use multi-locus genotypes to infer the number of genetic subpopulations (K) maintaining both Hardy-Weinberg and linkage equilibrium. geneland, however, uses spatial data to infer the spatial boundaries that separate the K subpopulations (Guillot et al., 2005). Because geneland has been shown to outperform other Bayesian clustering methods in detecting barriers to dispersal in fewer generations for species with higher dispersal abilities (Blair et al., 2012;Safner et al., 2011), we based our inferences on geneland. Results from the structure analysis were similar and are available, along with the methods, in Appendix S1. We performed 10 independent runs using the uncorrelated and correlated allele frequency models. We used both models because the former is less sensitive to departures from model assumptions while the latter is more apt to detect subtle genetic differentiation (Guillot et al., 2008;The Geneland Development Group, 2018). We varied K from 1 to 31 (the maximum number of sampling locations +1) and then used the model with the highest mean posterior probability to select K and assigned individuals to the population in which their estimated proportion of ancestry (i.e., the Q-value) was the greatest.
We optimized all models using 500,000 Markov Chain Monte Carlo iterations, 1000 burn in, a 100-iteration thinning interval, an uncertainty of 2 km for GPS coordinates, and a maximum rate of 1650 nuclei for the Poisson-Voronoi tessellation (three times the number of individuals). We implemented our analysis in the r package geneland using program r (version 3.4.4; R Core Team, 2018; The Geneland Development Group, 2018). After assessing population structure, we again assessed for LD, null alleles, and deviations from HWE, and if these tests failed then we reassessed for population structure until tests were not significant.

| Environmental variables
Available food resources should influence habitat selection, as black bears must accumulate large-fat stores for both hibernation and reproduction (Costello et al., 2003) and food for bears in arid environments is tied to precipitation (precip; Zlotin & Parmenter, 2008).
Black bears are forest obligates and have evolved morphological and behavioral adaptations associated with exploiting forest stands (Herrero, 1972) and they require thermal refugia (Lara-Díaz et al., 2018) because they are susceptible to hyperthermia (Sawaya et al., 2017). We modeled these features using canopy height (canopy), percent canopy (percan), and water bodies (water) as canopy provides thermal cover and water is necessary for thermoregulation, especially if bears crossed more inhospitable land cover such as desert. Male black bears have been shown to use less rugged areas (Costello, 2010;Johnson et al., 2015;Onorato, Hellgren, Van Den Bussche, & Skiles Jr., 2004), so we used a Terrain Ruggedness Index (TRI) to represent potential movement corridors. Linear-water features (streams) contain food, escape, and thermal cover, and are travel corridors (Atwood et al., 2011;Johnson et al., 2015). Roads can elicit negative behavioral and genetic effects and can influence bear distribution (Dixon et al., 2007;Gould et al., 2019), so we assessed their effect by estimating road density (rd. density). Interstates and highways (rd. major), which if not acting as barriers, may inhibit gene flow by heightening mortality rates (Little et al., 2017). We did not explore the influence of other anthropogenic activities such as agriculture or human settlements because the former is uncommon and sparsely distributed across New Mexico while the latter is correlated with road density.
We determined the spatial extent of the environmental variables, and subsequent resistance surfaces, by buffering all sample locations by 61 km based on the maximum-dispersal distance for black bears in the Sangre de Cristo and Mogollon mountains, New Mexico (Costello, 2010 Hansen et al., 2013). We obtained canopy height (continuous covariate) at a 1 km resolution using the National Aeronautics and Space Administration EARTHDATA Spatial Data Access Tool (https://daac.ornl.gov). We obtained location data for streams (binary covariate) and water bodies (binary covariate) from the National Hydrography Dataset (https:// www.usgs.gov/core-scien ce-syste ms/ngp/natio nal-hydro graphy).
We derived TRI (continuous covariate) using a National Elevation Dataset 30 m digital elevation model (www.natio nalmap.gov) and the Benthic Terrain Modeler in arcMap. We calculated road density (km/25 km 2 ; continuous covariate) and major roads (categorical covariate) using data from Open Street Map (www.opens treet map.org). For major roads, we used three classifications: interstate highways, state highways, and county roads. We resampled each resistance layer to a 5 km resolution using bilinear interpolation to reduce the computational intensity of the optimization process We created and manipulated all resistance surfaces using arcMap v10.4.1 (Environmental Systems Research Institute, Redlands, CA, USA). We assessed correlation among covariates using a Pearson's correlation coefficient of r ≥ |0.60|. We found a correlation between canopy height and percent canopy (r = 0.72) and removed the latter from our analysis.

| Generating the resistance surface
We used the r package resistancega to optimize resistance surfaces, assess the effect of pairwise-effective distance on pairwisegenetic distance, and conduct model selection while accounting for non-independence among the pairwise data (version 4.1-0.2.1; Peterman, 2018;Peterman et al., 2014). resistancega optimizes resistance surfaces using a genetic algorithm (a process based on the theory of natural selection) that eliminates the subjective assignment of resistance values by expert opinion and the limited exploration of the optimized parameter space (Peterman, 2018;Peterman et al., 2014). The optimization process begins with the selection of parameter values that control the transformation, shape of the transformation, and resistance value for a continuous surface, or if a categorical surface, the assignment of values to each resistance level. After each iteration, pairwise effective distances among all individuals are calculated and a linear mixed-effects model is then fit to the data where effective distance is used to predict genetic distance among individuals (Clarke et al., 2002;Peterman, 2018;Peterman et al., 2014). The relative support for the combination of parameter values at each iteration is assessed using an objective function from the mixed-effects model, and once the objective function can no longer be improved, surface optimization is completed.
We quantified effective landscape distance using randomwalk commute times in the r package gdistance (version 1.2-2; Van Etten, 2017). We quantified pairwise-genetic distance using the individual-based metric proportion of shared alleles (D ps ) in the r package adegenet (version 2.1.1; Jombart, 2008). We applied a monomolecular and Ricker transformation along with their inverse, reverse, and inverse-reverse forms to each continuous-resistance covariate to explore the functional relationship between each covariate and resistance to movement. We constructed our model using a maximum of three covariates due to computational intensity and assessed the relative support among resistance surface models using Akaike's Information Criterion adjusted for small sample size (AICc) with models >2 AICc units from the top model being discounted (Burnham & Anderson, 2002;Hurvich & Tsai, 1989) and by calculating the model weight (w i ). We conducted the optimization process in program r (version 3.4.4; R Core Team, 2017) using the Bridges high-performance computing system at the Pittsburgh Supercomputing Center (Nystrom et al., 2015;Towns et al., 2014).
To explore how geographic distance vs. landscape distance affects pairwise-genetic distance, we used AIC C to rank the fit of linear mixed-effects models using Euclidean distance, the top-ranked resistance surface, and a model that combined both.

| Relative degree and direction of gene flow
We investigated asymmetric gene flow by estimating relative migration among the estimated subpopulations using the divMigrate function in the r package diversity where maximum relative gene flow is set at 1 and minimum at 0 (Keenan et al., 2013;Sundqvist et al., 2016). We calculated G ST , a measure of population differentiation and an analog of F ST , for network plots, conducted 1000 bootstrap iterations to generate 95% CIs to evaluate if asymmetric gene flow was significant and chose to display connections ≥0.50.

| Describing genetic structure
We genotyped 550 (285M:265F) individuals from 28 localities (Appendix S2: Table S1). We found a moderate percentage of null alleles (4-12%) across all loci in this total sample (Appendix S2: Table S2). We found 82% of the pairwise comparisons among loci (n = 105) for LD to be significant (P < 0.0005 after Bonferroni correction) and all loci were out of HWE (p < .003 after Bonferroni correction; Appendix S2: Tables S3-S4). These metrics suggest that genetic structuring may occur among black bears across the Southwest.  Table S1). We found the CXX20 locus to be non-randomly associated with G10B and G10H in the BM subpopulation (Appendix S3: Tables S2-S4). The G10U and G10L loci were out of HW proportions in the ECPSRM and DMS, respectively (Appendix S3: Table S5). The presence of null alleles and linkage disequilibrium suggests these loci may not accurately represent genetic structure and diversity while deviations from HW proportions suggest there could be additional genetic structure that was not detected under the uncorrelated allele frequency model.
Allelic richness was lowest in the SIS (A R = 3.91) and highest in the DMS (A R = 5.43); the TP had the second highest allelic richness despite small sample size (Table 1). The number of private alleles using rarefaction was lowest in DMS (A PR = 0.09) and highest in the TP (A PR = 1.79; Table 1; Appendix S3: Table S6). Observed heterozygosity ranged from 0.42 to 0.64 was slightly lower than H E (0.44-0.62) for all regional subpopulations except for the TP (Table 1). The F IS estimates suggested deviations from random mating within the ECPSRM and DMS subpopulations, but along with H O being lower than H E for both subpopulations, it is more likely that a Wahlund effect, rather than nonrandom mating, is occurring, which indicates greater substructure within these two regions (Wahlund, 1928). Overall, the two most isolated subpopulations, BM and the TP, displayed the highest levels of genetic differentiation compared to all other subpopulations (Table 2).
the correlated allele frequency model identified 14-genetic clusters that closely tracked the sampled mountain ranges (Appendix S2: Figure S1). The BM and TP subpopulations from the regional results remained while the larger clusters were broken down into 12 subpopulations ( Figure 1). We did not find evidence of null alleles (Appendix S4: Table S1-S8). All loci within each respective subpopulation were in HWE (Appendix S4: Table S9). Because there was no discernable pattern of null alleles, LD, or HW disequilibrium for ≥1 locus at ≥1 subpopulation we retained all loci in our analyses (Morin et al., 2010). TA B L E 1 Number of individuals (N), private alleles (A P ), private alleles using rarefaction (A PR ), allelic richness using rarefaction (A R ), observed (H O ) and expected (H E ) heterozygosity, and a measure of deviations from random mating (F IS ) and its 95% confidence interval (LCI and UCI) based on 1000 bootstrap iterations for American black bear (Ursus americanus) subpopulations in the American Southwest and northern Mexico.

| Landscape features regulating gene flow
The top-ranked resistance-surface model was well supported (w i = 1.00), substantially outperformed the second-ranked model (ΔAIC c = 47.18), and included canopy, precipitation, and TRI (Appendix S5: Tables S1 and S2). The transformations that best represented the relationship between canopy, precipitation, and TRI with resistance to movement were the inverse monomolecular, inverse Ricker, and monomolecular, respectively, indicating that resistance decreased as canopy increased, decreased as precipitation increased until the covariate reached moderate levels at which point resistance started to increase, and increased as TRI increased (Table 4). More simply, resistance was lowest in areas with higher forest canopy, higher levels of precipitation, and less rugged landscapes. Precipitation contributed the most to the top-ranked resistance surface (58%) followed by canopy (40%) with a small contribution from TRI (2%). This top-ranked model received considerable support when compared to Euclidean distance alone, suggesting isolation by resistance better explained the observed-genetic pattern than isolation by distance (  (Table 4). Our analysis did not show support for any resistance-based models that included road density or major roads (Appendix S5: Tables S1 and S2).

| Relative degree and direction of genetic connectivity
The directional relative migration network-clustered populations in the northern part of our study region along the Colorado-New  Figure 2; Appendix S3: Figure S1; Appendix S3: Table   S7). The mountain range subpopulations exhibited a similar pattern with subpopulations from the central portion of the study area clustering together and asymmetric gene flow in a northward direction (Figures 1 and 3; Appendix S4: Figure S1; Appendix S4: Table S11).

| DISCUSS ION
Our study further supports the hypothesis that the Southwest

| The influence of forest refugia on genetic structure
Our analyses supported the hypothesis that the current genetic structure is a consequence of Pleistocene-Holocene vicariance whereby bears occupied forest refugia during the LGM, but then followed changes in the distribution of forests as the Holocene dried and warmed (Pleistocene-Holocene Vicariance Hypothesis). The Southwest likely served as a refugium for black bears during various periods in the Pleistocene and Holocene when habitats in more northerly latitudes were dominated by more cold-adapted plant species that black bears do not typically use .
Thus, for much of the late Pleistocene and into the early Holocene, the dominant paleovegetation community of the region was a TA B L E 3 Estimated pairwise genetic differentiation (F ST ) and their 95% confidence intervals based on 1000 bootstrap iterations for mountain range subpopulations identified by geneland using the correlated allele frequency model for American black bears (Ursus americanus) in the American Southwest and northern Mexico. piñon-juniper-oak woodland that black bears inhabited, similar to the plant community selected by black bears in the Sky Islands today (Onorato et al., 2003). The isolated Sky Island mountain ranges, currently inhabited by black bears, were most likely functionally connected by this piñon-juniper-oak woodland (Van Devender, 1990a).
Black bears are omnivorous, but vegetation, fruits, and nuts comprise 70-90% of the diet, supplemented with insects and vertebrates (Delgadillo Villalobos et al., 2019). In spring, they feed on grasses and other vegetation, in mid-late summer on soft mast, such as berries, and in late summer-fall prior to hibernation they forage on hard mast, such as acorns and piñon pine nuts (Beck, 1991;Costello et al., 2001;Onorato et al., 2003). In lower elevations within the Southwest, they also feed on sotol (Dasylirion spp.), yucca (Yucca spp.), and prickly pear cactus (Opuntia spp.; Delgadillo Villalobos et al., 2019). Although found in semiarid shrublands, black bears are primarily a forestadapted species and forests are important habitats across their range (Evans et al., 2017;Gould et al., 2019;Onorato et al., 2003). Thus, it stands to reason that black bears would track the abundance of their main food over the short term, which would explain contemporary movements and dispersal patterns and populations would track the distribution of their primary habitat over the long term, which would explain species distribution and population genetic structure.

| The influence of transportation infrastructure on genetic structure and gene flow
Our analyses did not support the hypothesis that interstate highways are limiting the movement of black bears across the Southwest Roads, urbanization, and interstate highways can negatively influence carnivore populations and the size of the interstate (e.g., number of traffic lanes) and relative traffic flow may be contributing factors as well (Riley et al., 2014;Serieys et al., 2015). Perhaps one of the most extreme cases has occurred in the Santa Monica Mountains of southern California where the morass of urbanization and grand thoroughfares has restricted population size and caused degradation in genetic variation in mountain lions (Puma concolor; Riley et al., 2014). Although bear resource use is negatively affected by roads (Gould et al., 2019), the effect of roads on bear movements and gene flow varies across their range. Roads had little effect on movements in remote areas such as in Idaho (Cushman et al., 2006) but had major impacts on movement patterns and genetic structure in more heavily urbanized areas such as Florida (Dixon et al., 2007;McCown et al., 2009). The highways in the Southwest receive less TA B L E 4 Model selection results for two optimization runs derived using Akaike's Information Criterion corrected for small sample size (AIC c ) comparing the top-ranked resistance surface optimized using linear mixed-effects models with maximum likelihood population effects parameterization to models composed of Euclidean distance (Distance Only) and Euclidean plus the top-ranked resistance surface (Top resistance surface + Distance).  (Blair et al., 2012;Epps et al., 2005;Safner et al., 2011).

| The scale of population genetic structure in Southwestern black bear populations
Regionally, subpopulation boundaries followed the distribution of three major ecoregions (Omernik & Griffith, 2014) (Lomolino et al., 1989). Atwood et al. (2011) also found genetic substructure among black bear populations along the US-Mexico border within the Sky Islands region.
There was a relatively high degree of genetic differentiation when comparing Boulder Mountain, Utah and the Trans-Pecos region in west Texas and northeast Mexico to the other Southwest subpopulations. We believe this differentiation is due to genetic isolation rather than incomplete sampling. The origination of the Boulder Mountain population is more enigmatic and more information is needed to determine if it is a product of eastward expansion by populations from the Pacific Northwest refugium or through the expansion and isolation of populations from the north (Lackey et al., 2013;Malaney et al., 2018;Puckett et al., 2015).
We had a small sample from the Trans-Pecos region, but that subpopulation had the highest allelic richness, the highest observed heterozygosity, many private alleles, and some private alleles occurred at high frequency, indicating a period of isolation and genetic differentiation followed by little to no connectivity (Slatkin, 1985). The existence of private alleles in the Trans-Pecos region is most likely a product of isolation followed by bears recolonizing west Texas from the Sierra Madre Oriental (Onorato et al., 2007;Onorato, Hellgren, Van Den Bussche, & Doan-Crider, 2004 (Harris, 1987(Harris, , 1989(Harris, , 1993(Harris, , 2003Messing, 1986;Saunders, 1977;Slaughter, 1975). These caves could be sampled for ancient environmental DNA or ancient DNA could be amplified from the fossils themselves to further our understanding of the distribution of refugia and the movements and genetic structure of bears in the American Southwest and northern Mexico.

| The influence of landscape resistance and geographic distance on gene flow
Our Our resistance-based models of gene flow revealed that areas with higher canopy cover and precipitation, essentially forested habitats, were associated with higher rates of gene flow and within these areas, gene flow was facilitated by less rugged areas. Gene flow was also affected by geographic distance. Thus, populations connected by contiguous forest had high gene flow (e.g., Eastern Colorado Plateau and Southern Rocky Mountains, Datil-Mogollon Section, Mexican Highlands, and Sacramento Section) and populations separated by desert and isolated by distance had lower gene flow (e.g., Trans-Pecos and Boulder Mountain). Because female black bears are highly philopatric, the effect of distance on gene flow may be governed by their behavior, but also mediated by long-distance dispersal events by male bears with male bears also having been shown to select against ruggedness (Apps et al., 2006;Johnson et al., 2015;Lara-Díaz et al., 2018;Pelletier et al., 2011).
So, it appears that habitat most likely acts as a conduit (e.g., forest) or filter (e.g., desert) to bear movement and that geographic distance plays an important role owing to intersexual differences in movement behavior.

| Conservation implications
American black bears in the Southwest occupy a naturally frag- Climate change is contributing to a rise in aridity and temperature in the Southwest and has led to increases in insect outbreaks, intense droughts, and catastrophic wildfires resulting in substantial tree mortality reducing the distribution and quality of bear habitat over the long term (Gould et al., 2019;Thorne et al., 2018;Williams et al., 2010). Increasing human development and population growth is likely to increase human population density and traffic rates resulting in higher rates of road mortality, which could lower genetic connectivity and enhance fragmentation, heightening the extinction risk for some subpopulations (Dixon et al., 2007;Ernest et al., 2014;Riley et al., 2014).

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Microsatellite data used in this study are available through USGS ScienceBase, https://doi.org/10.5066/P91COLPR.