https://cpw.cvlcollections.org/files/original/a9dfac8e8bdbf1a6004f021f6c8ca711.pdf c4119107ba38090606d6aff6d9b8fb49 PDF Text Text The research in this publication was partially or fully funded by Colorado Parks and Wildlife. Dan Prenzlow, Director, Colorado Parks and Wildlife • Parks and Wildlife Commission: Marvin McDaniel, Chair • Carrie Besnette Hauser, Vice-Chair Marie Haskett, Secretary • Taishya Adams • Betsy Blecha • Charles Garcia • Dallas May • Duke Phillips, IV • Luke B. Schafer • James Jay Tutchton • Eden Vardy �Received: 18 September 2018 | Revised: 11 August 2019 | Accepted: 20 September 2019 DOI: 10.1111/mec.15261 ORIGINAL ARTICLE Urbanization impacts apex predator gene flow but not genetic diversity across an urban‐rural divide Daryl R. Trumbo1 | Patricia E. Salerno1 | Kenneth A. Logan2 | Mathew W. Alldredge3 | Roderick B. Gagne4 | Christopher P. Kozakiewicz5 Nicholas M. Fountain‐Jones6 Kevin R. Crooks 8 | Simona Kraberger4 | | Meggan E. Craft6 | Scott Carver5 | Holly B. Ernest7 | Sue VandeWoude 4 | 1,9 | W. Chris Funk 1 Department of Biology, Colorado State University, Fort Collins, CO, USA Abstract 2 Apex predators are important indicators of intact natural ecosystems. They are also Colorado Parks and Wildlife, Montrose, CO, USA 3 Colorado Parks and Wildlife, Fort Collins, CO, USA 4 Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, CO, USA 5 Department of Biological Sciences, University of Tasmania, Hobart, TAS., Australia 6 Department of Veterinary Population Medicine, University of Minnesota, Saint Paul, MN, USA 7 Department of Veterinary Sciences, University of Wyoming, Laramie, WY, USA 8 Department of Fish, Wildlife, and Conservation Biology, Colorado State University, Fort Collins, CO, USA 9 Graduate Degree Program in Ecology, Colorado State University, Fort Collins, CO, USA Correspondence Daryl R. Trumbo, Department of Biology, Colorado State University, Fort Collins, CO, USA. Email: daryl.trumbo@gmail.com Funding information Directorate for Biological Sciences, Grant/ Award Number: 1413925 and 723676 sensitive to urbanization because they require broad home ranges and extensive contiguous habitat to support their prey base. Pumas (Puma concolor) can persist near human developed areas, but urbanization may be detrimental to their movement ecology, population structure, and genetic diversity. To investigate potential effects of urbanization in population connectivity of pumas, we performed a landscape genomics study of 130 pumas on the rural Western Slope and more urbanized Front Range of Colorado, USA. Over 12,000 single nucleotide polymorphisms (SNPs) were genotyped using double‐digest, restriction site‐associated DNA sequencing (ddRADseq). We investigated patterns of gene flow and genetic diversity, and tested for correlations between key landscape variables and genetic distance to assess the effects of urbanization and other landscape factors on gene flow. Levels of genetic diversity were similar for the Western Slope and Front Range, but effective population sizes were smaller, genetic distances were higher, and there was more admixture in the more urbanized Front Range. Forest cover was strongly positively associated with puma gene flow on the Western Slope, while impervious surfaces restricted gene flow and more open, natural habitats enhanced gene flow on the Front Range. Landscape genomic analyses revealed differences in puma movement and gene flow patterns in rural versus urban settings. Our results highlight the utility of dense, genome‐scale markers to document subtle impacts of urbanization on a wide‐ranging carnivore living near a large urban center. KEYWORDS effective population size, gene flow, genetic diversity, landscape genomics, Puma concolor, urbanization 1 | I NTRO D U C TI O N al., 2017; Theobald, 2005). Habitat fragmentation due to urbanization can have important impacts on predator movement, dis- Urbanization is a major threat to biodiversity, and in particular to ease, and survival (Carver et al., 2016; Fountain‐Jones et al., 2017; apex predators with broad home ranges (Cohen, 2003; Crooks et Markovchick‐Nicholls et al., 2008). This reduced connectivity can 4926 | © 2019 John Wiley & Sons Ltd wileyonlinelibrary.com/journal/mec� Molecular Ecology. 2019;28:4926–4940. �| TRUMBO et al. 4927 lead to smaller, more isolated populations, where less gene flow Zeller, Vickers, Ernest, & Boyce, 2017). Furthermore, the substantial and genetic diversity, as well as smaller effective population sizes area requirements of large carnivores such as pumas can enhance (Ernest, Vickers, Morrison, Buchalski, & Boyce, 2014; Riley et al., their role as umbrella species, whose protection also benefits co‐ 2006; Vandergast, Bohanak, Weissman, & Fisher, 2007) ultimately occurring species through broadscale habitat preservation (Thorne, cause local and regional extirpations through environmental and Cameron, & Quinn, 2006). demographic stochasticity and inbreeding depression (Allendorf, The southern Rocky Mountains in western Colorado, USA sup- Luikart, & Aitken, 2013). Moreover, increased human recreational port natural habitats with high puma densities, as well as many rural activities in wildlife habitats associated with nearby urbanization can and urban human developments (Hornocker & Negri, 2009). The change wildlife movement patterns and habitat usage, exacerbating Western Slope of the Rocky Mountains primarily consists of large the impacts of fragmentation (Lewis et al., 2015; McKinney, 2002). areas of contiguous public wildlands with an abundant prey base for As human populations continue to expand worldwide, urban areas pumas, interspersed with small rural and exurban developments, in- are becoming larger and more extensive on the landscape. However, cluding the Uncompahgre Plateau region near the town of Montrose we do not fully understand how urbanization affects natural eco- (Western Slope Study Area; Figure 1). In contrast, the Front Range is systems near wildland‐urban interfaces (Magle, Hunt, Vernon, & a rapidly urbanizing, major metropolitan area on the Eastern Slope of Crooks, 2012; Radeloff et al., 2005). the Continental Divide, where urbanization is spreading from lower Large carnivores are important indicators of intact natural eco- elevation areas in and around the Denver Metropolitan Area into systems, as they require an abundant and sustainable prey base, as adjacent wildland habitats in the foothills of the Rocky Mountains. well as high habitat connectivity to support their broad home ranges Pumas continue to persist near this wildland‐urban interface, in- (Sergio et al., 2008; Sergio, Newton, Marchesi, & Pedrini, 2006). cluding adjacent to the city of Boulder on the western edge of the However, understanding the effects of urbanization on large carni- Denver Metropolitan Area (Front Range Study Area; Figure 1; Lewis vores is difficult due to their low population densities and secretive et al., 2015; Moss, Alldredge, Logan, & Pauli, 2016). From 2010– nature (Hornocker & Negri, 2009; Logan & Sweanor, 2001; Riley et 2017, Colorado was the eighth fastest growing U.S. state by popu- al., 2006). Camera traps, radiotelemetry, and GPS collars provide lation (577,829 residents added) and the sixth fastest by percentage valuable information on animal home ranges and population sizes (11.5% population growth; US Census Bureau, 2017), with most of (e.g., Blecha, Boone, & Alldredge, 2018; Lewis et al., 2015), but these this growth occurring along the eastern edge of the Front Range. studies are expensive, time consuming, and can only monitor a small Thus, comparative studies of puma movement and gene flow in one fraction of the total population for limited time periods. Population of the most populous states in the midcontinental USA, which also and landscape genetics can provide additional, complementary supports a robust puma population, can provide insight into the ef- techniques for a more detailed understanding of wildlife popula- fects of urbanization on this important apex predator. tions (Balkenhol, Cushman, Storfer, & Waits, 2016; Epps, Wehausen, Here, we tested how different landscape factors, including ur- Bleich, Torres, & Brashares, 2007; Lowe & Allendorf, 2010). Genetic banization, enhance or restrict gene flow and genetic diversity in a studies provide an indicator of functional landscape connectivity large apex predator across an urban‐rural divide in Colorado, USA. through measures of gene flow, effective population sizes of breed- A large sample of pumas were utilized from (a) the rural Western ing individuals, and cost‐efficient monitoring of genetic diversity Slope and (b) the more urbanized Front Range (n = 130; 76 in the across broad geographic areas (McRae, Beier, Dewald, Huynh, & Western Slope, 54 in the Front Range; Figure 1). We used double Keim, 2005; Solberg, Bellemain, Drageset, Taberlet, & Swenson, digest restriction site associated DNA sequencing (ddRADseq) to 2006). Moreover, recent high‐throughput sequencing technologies genotype pumas at 12,444 single nucleotide polymorphism (SNP) enable the genotyping of many more thousands of loci than previ- loci to evaluate the potential differences in gene flow, effective ously possible, providing higher power to detect the often subtle population sizes, genetic diversity, and population structure in population genetic structure of wide‐ranging species such as large these two different landscapes. We tested landscape genomic hy- carnivores (Holderegger, Kamm, & Gugerli, 2006; Luikart, England, potheses by correlating key landscape factors with puma genetic Tallmon, Jordan, & Taberlet, 2003). distance measures. We hypothesized that pumas in the more ur- Pumas (Puma concolor; other common names include mountain banized Front Range would have (a) smaller effective population lions, cougars, panthers, catamounts) are a large, apex predator sizes, (b) lower levels of genetic diversity, and (c) more landscape with one of the broadest latitudinal ranges of any terrestrial carni- factors related to urbanization that restrict gene flow, relative to vore, spanning western North America, Central America, and South the rural Western Slope landscape. America (Hornocker & Negri, 2009). Pumas are sensitive to urbanization, requiring broad‐scale landscape connectivity to persist, and are thus useful indicators for monitoring the effects of urban fragmentation (Beier, 1995; Crooks, 2002; Maletzke et al., 2017). Given sufficient habitat area and landscape connectivity, however, pumas 2 | M ATE R I A L S A N D M E TH O DS 2.1 | Samples and sequences can still persist within and adjacent to urban systems (Blecha et al., Puma blood and tissue samples were collected as part of ongo- 2018; Lewis et al., 2015; Riley et al., 2014; Wilmers et al., 2013; ing monitoring efforts by Colorado Parks and Wildlife in both the �4928 | TRUMBO et al. (a) (b) (c) (d) (e) (f) (g) (h) (i) F I G U R E 1 Study area in the Western Slope and Front Range of the southern Rocky Mountains of Colorado, USA. Landscape genomic analyses included 76 pumas from the Western Slope and 54 pumas from the Front Range (white circles). Resistance surfaces, shown for the Western Slope, represent alternative hypotheses of the effects of landscape variables on puma dispersal and gene flow (red = high gene flow, blue = low gene flow) for: (a) percent impervious surface cover (negative effect on gene flow), (b) land cover (forested, open‐natural, and developed: positive, neutral, and negative effects on gene flow), (c) percent tree canopy cover (positive effect), (d) vegetation density (positive effect), (e) river and stream riparian corridors (positive effect), (f) roads (negative effect), (g) minimum temperature of the coldest month (negative effect), (h) annual precipitation (positive effect), and (i) topographic roughness (positive effect). We also tested isolation by geographic Euclidean distance. Land cover base maps show forests (green), shrub and grasslands (tan), urban areas (red), agriculture and ranchlands (brown and yellow), and alpine tundra (grey). Note that impervious surface shows very little spatial variation on the Western Slope (a) because there is very little impervious surface present in this rural region [Colour figure can be viewed at wileyonlinelibrary.com] Western Slope and Front Range regions of the southern Rocky sampling represents a large proportion of the resident pumas pre- Mountains of Colorado, USA (Figure 1; Carver et al., 2016; Lewis et sent in both regions during the sampling period, as Lewis et al. (2015) al., 2015). Samples were collected from 2005–2014 on the Western estimated 14.4 (SE 1.6) and 14.7 (SE 1.3) resident pumas occupying Slope and 2007–2013 on the Front Range. Western Slope samples the Western Slope and Front Range study areas at a single time consisted of 36 males and 40 females, and Front Range samples con- point, from motion camera and telemetry data collected in 2009 and sisted of 23 males, 30 females, and one puma of unknown sex. Our 2010, using mark‐resight analysis. �| TRUMBO et al. 4929 Genomic DNA was extracted from tissue or blood using QIAGEN was the number of allele mismatches between replicate pairs di- DNeasy Blood &Tissue kits (QIAGEN Inc.). We genotyped a total of vided by the number of loci, and SNP error rate was the proportion 76 individuals from the Western Slope and 54 individuals from the of SNP mismatches between replicate pairs. Front Range using the ddRADseq protocol described in Peterson, Locus, allele, and SNP error rates, as well as the number of SNPs Weber, Kay, Fisher, and Hoekstra (2012) and sequenced on Illumina retained, were largely consistent across Stacks parameter settings HiSeq2500 and 4000 machines (Illumina) using 100 bp single‐end (Table S1), with the exception of the minimum number of identical, sequencing at the University of Oregon Genomics Facility (gc3f. raw reads required to create a stack (–m). Increasing –m to seven uoreg​on.edu). We tested nine different combinations of restric- reduced allele and SNP error rates, but not the locus error rate, tion enzymes on puma samples for digestion efficiency and eval- while also reducing the number of SNPs retained by up to 40%. uated the size ranges of fragment distributions using an Agilent For our final de novo assembly, we balanced retaining SNPs with Tapestation 2200 (Agilent Genomics). We chose the digest enzymes reducing locus, allele, and SNP error by using the parameter set- EcoRI‐HF (6b recognition) and NlaIII (4b recognition) and a target tings –m = 3, –M = 4, –n = 4, and max_locus_stacks = 3 (Table S1). fragment size range of 300–400 bp (excluding adapters). We used a In practice, the Pearson's correlation was very high between ge- Blue Pippin with a 2%, internal standard, 100–600 bp gel cartridge netic distances calculated with the final parameter settings versus (Sage Science) for size selection and a biotinylated P2 adapter with the most conservative (i.e., –m = 7) parameter settings (r = 0.997 DynaBeads (Peterson et al., 2012) to purify the polymerase chain re- for proportion of shared alleles). We then exported the SNP ma- action (PCR) template for the final enrichment. PCR was performed trix with the for 12 cycles and five reactions were tested for each pool of indi- retaining SNPs that were present in at least 20% of individuals viduals. We initially genotyped 16 individuals multiplexed into an by population, and retaining a single random SNP per locus. This Illumina 2500 HiSeq lane to estimate maximum multiplexing based matrix was further filtered for missing data in Plink v.1.07, first by on a target of >12× coverage per locus. After assessment of locus locus, then by individual, and then by minor allele frequency (MAF) coverage, we proceeded to multiplex 48 and 66 individually‐bar- using multiple combinations of thresholds for reducing missing coded samples on Illumina 2500 and 4000 HiSeq lanes, respectively, data in the matrix (see github.com/pesal​e rno/PUMAg​e nomics). using the Peterson et al. (2012) flex adaptors. After evaluating missing data from SNP matrices, we retained the populations program in Stacks (Catchen et al., 2013), matrix with a more stringent locus filter (excluding loci missing 2.2 | Bioinformatics pipeline and filters >25% individuals) and a less stringent filter on minor allele frequency (excluding loci with MAF < 0.01). We additionally filtered We evaluated read quality for each sequencing lane using Fast QC loci that were found at position 95 (the last position of our reads) (bioin​forma​t ics.babra​ham.ac.uk) and assembled our SNP data due to a higher number of SNPs present in this position, suggest- set de novo using Stacks 1.41 (Catchen, Hohenlohe, Bassham, ing increased error rates due to low sequence quality towards the Amores, & Cresko, 2013). Details on Stacks code and parameter end of the sequencing read. To compare landscape resistances settings used are on the GitHub repository; github.com/pesal​ with putatively neutral loci, we used a Principal Components erno/PUMAg​e nomics. We demultiplexed and filtered sequenc- Analysis (PCA) to identify loci showing strong signatures of se- ing reads using the program lection relative to neutral background genomic variation with the v process _ radtags in Stacks . Due to sensitivity of downstream genotyping with different Stacks pa- program PCA dapt (Luu, Bazin, & Blum, 2016). We found 12, puta- rameter settings (Mastretta‐Yanes et al., 2015; Paris, Stevens, & tively adaptive, outlier loci using a false discovery rate of 10%, so Catchen, 2017), we incorporated individual sample replicates in we filtered these outliers out for downstream landscape genomic library preparations. In each library, we included three within and analyses to avoid confounding neutral demographic patterns with three between library replicates, which were used for estimating patterns generated by loci under selection. genotyping error rates for different combinations of parameters used to construct loci with the denovo _ map. pl pipeline in Stacks . We ran 11 different de novo assemblies, varying four different 2.3 | Population genomics and structure Stacks parameters one at a time while leaving the other param- Population genomic statistics were calculated for the two sampling eters at their default values (Mastretta‐Yanes et al., 2015), that regions, the Western Slope and Front Range (Figure 1). Observed affect locus, allele, and SNP error rates and the number of loci and expected heterozygosity (Hobs and Hexp), nucleotide diversity genotyped, consisting of (a) minimum number of identical, raw- (π), inbreeding coefficient (FIS), and population genetic differen- reads required to create a stack (‐m), (b) number of mismatches tiation (FST ) were calculated using the populations program in Stacks allowed between loci when processing a single individual (‐M), (c) with SNP loci that passed previous filters, excluding a single indi- number of mismatches allowed between loci when building the vidual (sample_1382) that did not pass the 75% missing data thresh- catalog (‐n), and (d) maximum number of stacks at a single de novo old. We estimated allelic richness (Ar) using hp‐rare 1.0 (Kalinowski, locus (‐max_locus_stacks) (Table S1). Locus error rate was calcu- 2005), which corrects for variance in sample sizes using rarefac- lated as the number of loci present in only one of the samples of a tion. Two complementary, individual‐based genetic distances were replicate pair divided by the total number of loci, allele error rate calculated: proportion of shared alleles distance (Dps; Bowcock et �4930 | al., 1994) using the TRUMBO et al. adegenet r v. 3.3.3 package and relatedness dis- tance (r; Smouse & Peakall, 1999) using the PopGenReport r pack- Two genetic distance measures were used as response variables in landscape genomic analyses: proportion of shared alleles distance age. We then calculated mean genetic distance among individuals (Dps; Bowcock et al., 1994) and relatedness distance (r; Smouse & for each region, corrected for geographic distance (i.e., genetic Peakall, 1999). Environmental resistances among individuals were distance per km), as individuals that are further apart are expected calculated using Circuitscape (McRae, 2006) for each landscape to have higher genetic distances due to neutral isolation by dis- resistance surface (McRae, 2006; Row, Knick, Oyler‐McCance, tance population processes (Balkenhol et al., 2016; Wright, 1942). Lougheed, & Fedy, 2017). Circuitscape resistances are a useful tool Effective population sizes (Ne) were estimated using the linkage in landscape genetics because they summarize all potential move- disequilibrium method in NeEstimatorv.2.01 (Do et al., 2014), using ment pathways simultaneously, as opposed to least cost paths that the correction for chromosome number (Waples, Larson, & Waples, evaluate only a single idealized pathway, and thus assume the study 2016), which has been shown to be a robust method for inferring Ne organism has complete knowledge of the landscape and always using SNP data sets and large sample sizes (Waples, 2016; Waples chooses the ideal pathway (Balkenhol et al., 2016; McRae, 2006). et al., 2016). We evaluated overall genetic structure as well as ge- Landscape variables were tested for multicollinearity, both prior netic differentiation among the two sampling sites (Western Slope to and after calculating environmental resistances in Circuitscape, and Front Range) using PCA and Discriminant Analysis of Principal to ensure Pearson's r correlations < .7 and variance inflation factor Components (DAPC) in the (Jombart, 2008; (VIF) scores <5 in final landscape genomics models, as collinearity Jombart, Devillard, & Balloux, 2015) and Admixture ancestry analy- can cause instability in parameter estimation in regression models sis (Alexander, Novembre, & Lange, 2009). Individual membership (Tables S3 and S4; Dormann et al., 2012; Row et al., 2017; Warren, probabilities of DAPC were based on the first 30 principal compo- Glor, & Turelli, 2010). r package adegenet nents and the first discriminant function. We used the function as- Two complementary methods were used to estimate the ef- signplot identify individuals that were putative migrants or admixed fects of environmental resistances on genetic distances: multiple based on the individual DAPC assignment probabilities. We used the regression on distance matrices (MRDM; Legendre, Lapointe, & find.clusters command in adegenet and minimized cross validation error in Admixture to estimate the number of populations (i.e., K). Casgrain, 1994) using permute v.3.4 and maximum likelihood of population effects (MLPE; Clarke, Rothery, & Raybould, 2002; Row et al., 2017; van Strien, Keller, & Holderegger, 2012) using the 2.4 | Landscape genomics lme 4 r package. MRDM is a permutational, distance matrix‐based approach that has been traditionally used in landscape genetic Geographic Information Systems (GIS) data were collected for dif- analyses, whereas MLPE is a newer linear mixed effects modelling ferent landscape factors that we hypothesized would affect puma technique that models pairwise comparisons as a random effect dispersal and gene flow in Colorado. Table 1 provides details on GIS and environmental resistances as fixed effects (Balkenhol et al., data sources, spatial resolution, and ecological justification for each 2016). Recent evaluations of landscape genetic approaches found landscape factor. Study area extents were calculated and landscape linear mixed effects modelling using MLPE to be more accurate, variables were compared across regions by buffering individual although both approaches performed well (Shirk, Landguth, & data points by a typical female puma dispersal distance of 34.6 km Cushman, 2017). Therefore, we included the traditional MRDM (Logan & Sweanor, 2001), dissolving overlapping buffers, and calcu- approach as well as MLPE in order to utilize multiple, complemen- lating zonal statistics within each region (Western Slope and Front tary techniques for inferring associations between landscape fea- Range) using ArcGISv. 10.1 (ESRI, Redlands, California). Landscape tures and gene flow. For MRDM and MLPE, genetic distances were data were converted into resistance surfaces using the Reclassify the response variable and environmental resistances were explan- and Raster Calculator tools in ArcGIS. The following hypothesized atory variables. Additionally for MLPE, a random effect matrix relationships of landscape factors with puma gene flow were mod- of individual comparisons was included to control for the nonin- eled: percent impervious surface cover (negative effect on gene dependent, pairwise structure of the data, and landscape resis- flow), land cover (forested, open‐natural, and developed: positive, tances were standardized to units of standard deviation centered neutral, and negative effects on gene flow), percent tree canopy on the mean (Row et al., 2017; van Strien et al., 2012). Models cover (positive effect), vegetation density (positive effect), river and were ranked using the Bayesian information criterion (BIC), and stream riparian corridors (positive effect), roads (negative effect), top models within five BIC units are reported (Richards, 2015). minimum temperature of the coldest month (negative effect), annual precipitation (positive effect), topographic roughness (positive effect), and elevation (positive effect). Additionally, we included an isolation by geographic distance model, which would be supported if none of the landscape variables had an effect on gene flow except 3 | R E S U LT S 3.1 | Genotyping and filtering SNP matrices for straight line, Euclidean distance between individuals (Balkenhol Initial Stacks processing retained a single random SNP per 95 bp read et al., 2016; Wright, 1942). Table S2 describes methods and justifica- and SNPs present in at least 20% of individuals by population, result- tion for converting raw landscape variables to resistance surfaces. ing in a matrix of 98,813 SNPs. These SNPs were further filtered in �Climate Min. temp. Ann. precip. Minimum temperature of the coldest month Mean annual precipitation Veg. density Enhanced vegetation index Riparian River and stream riparian corridors Tree cover Roads Road corridors Percent tree canopy cover Imperv. Percent impervious surface cover Vegetation Land cover Land cover: forested, open‐ natural, and developed Land cover Geo. dist. Code Isolation by geographic distance Landscape variable Mean annual precipitation accumulation (mm) calculated from 1970–2000 weather station data, interpolated between stations Mean annual minimum temperature of the coldest month (°C) calculated from 1970–2000 weather station data, interpolated between stations Density of vegetation calculated from chlorophyll reflectance in visual and near‐infrared spectra Percentage of tree canopy cover River and stream riparian corridors, with 50 m buffers on each side Roads, with 50 m buffers on each side Percentage of impervious surface Multiple land cover categories collapsed into 3 costs of movement: forested (lowest), open natural areas (medium), and developed (highest) Euclidean, straight‐line distance between individuals Description Global Climate Data (world​ clim.org/bioclim; Hijmans et al., 2005), 1 km Global Climate Data (world​ clim.org/bioclim; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), 1 km Moderate Resolution Imaging Spectroradiometer (modis. gsfc.nasa.gov), 250 m National Land Cover Database (mrlc.gov/nlcd2​011.php; Homer et al., 2004), 30 m National Hydrography Data set (nhd.usgs.gov), 30 m Colorado Department of Transportation (dtdap​ ps.color​adodot.info/otis), 30 m National Land Cover Database (mrlc.gov/nlcd2​011.php; Homer et al., 2004), 30 m National Land Cover Database (mrlc.gov/nlcd2​011.php; Homer et al., 2004), 30 m No environmental data; model assumes only distance affects gene flow, 30 m Data source, spatial resolution Model of isolation by straight‐line distance (Wright, 1942) Forested habitats provide the most cover for hunting and dispersal, open natural areas are intermediate, and developed areas are the least suitable habitat for dispersal (Crooks, 2002; Lewis et al., 2015) Human development results in increased noise, lights, and hunter access, limiting dispersal (Ernest et al., 2014; Maletzke et al., 2017; Riley et al., 2006) Roads increase mortality, noise, lights, and hunter access, limiting dispersal (Maletzke et al., 2017; Newby et al., 2013; Riley et al., 2006) River and stream riparian corridors provide vegetative and topographical cover for dispersal, as well as water sources attracting prey species (Dickson et al., 2005; Hilty & Merenlender, 2004; Naiman, Decamps, & Polluck, 1993) Low tree canopy limits cover for ambush predation and concealment, and restricts dispersal (Blecha et al., 2018; Logan & Sweanor, 2001; Sweanor et al., 2000; Warren et al., 2016) Low vegetation density limits cover for ambush predation and concealment, and restricts dispersal (Blecha et al., 2018; Hilty & Merenlender, 2004; Sweanor et al., 2000; Warren et al., 2016) Low minimum temperatures and high snowfall, found at high elevation mountain ridgelines (e.g., alpine tundra habitats) restrict hunting, breeding, and dispersal (Hornocker & Negri, 2009) Dry habitats with low precipitation accumulation limit prey species for hunting and vegetative cover, restricting dispersal (Logan & Sweanor, 2001; McRae et al., 2005) ArcGIS Reclassify tool, Circuitscape ArcGIS Spatial Analyst ArcGIS Spatial Analyst ArcGIS Analysis Tools, Spatial Analyst ArcGIS Analysis Tools, Spatial Analyst ArcGIS Spatial Analyst ArcGIS Spatial Analyst ArcGIS Spatial Analyst ArcGIS Spatial Analyst | (Continues) Ecological justification Calculation Environmental variables used for landscape genomic analyses, data sources, spatial resolution, and ecological justification Distance Category TA B L E 1 TRUMBO et al. 4931 �Elev. Elevation Elevation calculated from digital elevation models Topographic complexity based on variance in elevation within a moving window Description National Elevation Data set (lta.cr.usgs.gov/ned) National Map Tool (viewer.natio​ nalmap.gov), 30 m National Elevation Data set (lta.cr.usgs.gov/ned) National Map Tool (viewer.natio​ nalmap.gov), 30 m Data source, spatial resolution Increased forest cover and reduced human development found at higher elevations enhances hunting, breeding, and dispersal (Hornocker & Negri, 2009) ArcGIS Spatial Analyst 11,889 11,958 Western Slope Front Range 54 indiv. 76 indiv. Ngen 0.243 0.241 Hobs 0.263 0.272 Hexp 0.0028 0.0029 π 1.89 1.94 Ar 0.084 0.118 FIS 0.024 FST 0.51 (0.18) 0.29 (0.06) D PS/km (SE) 0.17 (0.07) 0.10 (0.02) r/km (SE) 48.1 (46.3–49.9) 81.2 (77.6–84.8) Note: Population genomic measures are observed heterozygosity (Hobs), expected heterozygosity (Hexp), nucleotide diversity (π), allelic richness (Ar), inbreeding coefficient (FIS), genetic differentiation among populations (pairwise FST ), mean genetic distance among individuals corrected for geographic distance (D PS and r per km) with standard errors (SE), and effective population size (Ne) with 95% confidence intervals (CI) based on parametric bootstrapping. Area (km2) Ne (95% CI) Steep, topographically‐complex canyons and mountain slopes provide cover for hunting and dispersal (Dickson et al., 2005; Hornocker & Negri, 2009) Ecological justification Geomorphometric and Gradient Metric Toolbox (Cushman, Gutzweiler, Evans, & McGarigal, 2010), ArcGIS Spatial Analyst Calculation Study areas (km2), number of individuals genotyped (Ngen), and population genomic parameter estimates from the Western Slope and Front Range of Colorado Topo. rough. Code Topographic roughness Landscape variable (Continued) Region TA B L E 2 Topography Category TA B L E 1 4932 | TRUMBO et al. �| TRUMBO et al. 4933 Plink by removing loci that were present in less than 75% of individu- estimates low or high (Waples, 2016; Waples et al., 2016). Ne re- als, which resulted in a matrix of 20,355 SNPs. Only a single individual mained consistently higher in the Western Slope, although it dif- was excluded based on our >75% missing loci per individual threshold. fered between the earlier and later sampling periods, indicating the After excluding SNPs present in the 95th sequencing base position and population may be expanding (Table S5). We were able to differenti- with minor allele frequency <0.01, we retained 12,456 SNPs. PCAdapt ate between the Western Slope and Front Range regions based on detected twelve outlier loci, putatively under selection, while account- PCA, DAPC, and Admixture ancestry analysis, and K = 2 was the best ing for population structure (K = 2). After removing these putatively supported value of K by minimizing cross validation error (Figure 2; adaptive loci, the final neutral data set contained 12,444 SNPs (Table Alexander et al., 2009). The proportion of correct individual assign- S1; github.com/pesal​erno/PUMAg​enomics). ment to populations based on DAPC was high for most individuals in both the Western Slope (0.98) and the Front Range (0.96), and the 3.2 | Population genomics and structure The two study areas encompass similar geographic extents: 2 2 assign plot identified putative migrants and admixed individuals in both regions (Figure 2b). However, the Admixture ancestry analysis showed more admixed individuals, particularly in the Front Range 11,889 km for the Western Slope and 11,958 km for the Front (Figure 2c). We also analyzed both regions separately for population Range (Table 2). Measures of genetic diversity (Hobs, Hexp, π, Ar,) and substructure (Figure S1), and there was no signature of population inbreeding (FIS) were similar for the Western Slope and Front Range differentiation within the Western Slope or Front Range. (Table 2). However, the effective population size (Ne) was smaller, mean genetic distances among individuals (D PS/km and r/km) were higher, and there was more admixture in the more urbanized Front 3.3 | Landscape genomics Range (Table 2, Figure 2). We also calculated Ne using subsets of The Front Range has more urban development than the Western individuals (i.e., pre‐ and post‐2010 individuals in the Front Range, Slope, with more impervious surface cover and a higher density of pre‐ and post‐2011 individuals in the Western Slope), because mul- roads (Figure 1, Table 3, Table S2). The Front Range also has more tiple overlapping generations may bias effective population size tree canopy cover, higher vegetation density, and higher annual (a) Principal components analysis (c) Admixture ancestry analysis (b) Discriminant analysis of principle components Western slope Western slope Front range Front range F I G U R E 2 Population structure from (a) Principal components analysis (PCA), (b) Discriminant analysis of principal components (DAPC), and (c) Admixture ancestry analysis. Individuals assigned to the Western Slope and Front Range are green and blue, respectively. K = 2 was most supported in Admixture ancestry analysis using cross validation error [Colour figure can be viewed at wileyonlinelibrary.com] �4934 | TA B L E 3 TRUMBO et al. Habitat differences between the Western Slope and Front Range of Colorado Western Slope Front Range Landscape data Min Max Median Mean Elevation (m) 1,453.5 4,362.9 2,354 2,418.0 552.5 1,474.5 Tree canopy cover (%) 0 100 20 29.9 31.4 0 Impervious surface (%) 0 100 0 0.5 4.1 0 Minimum temp. coldest month (°C) −20.2 −9.5 −13.3 −13.9 2.8 Annual precipitation (mm) 208 1,137 458 483.3 171.5 Enhanced vegetation index −1,806 8,955 4,634 4,434.9 1,858.3 Topographic roughness 0 27,924.6 11 53.1 129.6 Std Dev Min Max Median Mean Std Dev 4,347.1 2,365 2,374.9 629.3 100 32 35.0 33.6 100 0 4.0 13.5 −19.9 −8.3 −12.7 −12.6 2.8 359 1,006 452 496.4 121.9 −1,969 9,132 5,416 4,957.8 1,800.9 0 20,067.0 25 56.2 100.8 Landcover 1 10 1 3.2 3.0 1 10 1 4.2 3.7 Roads 1 10 1 1.7 2.5 1 10 1 2.7 3.6 Riparian areas 1 10 10 9.4 2.3 1 10 10 9.4 2.3 Note: Units are percent cover for impervious surface and tree canopy cover; resistance values for landcover, river and stream riparian corridors, and roads; °C for temperature; millimeters for precipitation; meters for elevation; and unitless measurements based on chlorophyll reflectance and variance in elevation, respectively, for enhanced vegetation index and topographic roughness. precipitation than the Western Slope (Table 3), probably due to the stream and river riparian corridors, roads, and minimum temperature high desert habitats (i.e., the Colorado Plateau ecoregion) in the of the coldest month; and for the Front Range included the same Western Slope being drier than the grassland and shrub habitats landscape variables plus impervious surface cover. found at lower elevations of the Front Range (i.e., the Great Plains ecoregion; McMahon et al., 2001). Landscape genomic patterns of pumas were different in the rural Western Slope compared to the more urbanized Front Range, Prior to running Circuitscape, landscape raster surfaces were with the exception of geographic distance being supported in both largely uncorrelated (i.e., Pearson's r < .7), with the exception of el- regions (Tables 4 and 5). In the Western Slope, tree canopy cover evation, which was positively correlated with annual precipitation was consistently positively correlated with gene flow in MRDM and and negatively correlated with minimum temperature of the cold- MLPE models. In addition, low minimum temperatures of the coldest est month in both regions, and vegetation density, which was nega- month were negatively correlated gene flow, and riparian habitats tively correlated with annual precipitation in the Front Range (Table were positively correlated with gene flow, in three of the top MLPE S3). After Circuitscape analyses, environmental resistance variables models (Tables 4 and 5). In contrast, in the Front Range, tree canopy showed more collinear relationships than raw raster surfaces (Table cover and percent impervious surface cover were negatively asso- S4), probably due to Circuitscape resistances being higher for indi- ciated with gene flow in the top MLPE models (Table 5). Because viduals separated by larger geographic distances (McRae, 2006). the relationship between tree cover and gene flow was the opposite Therefore, we removed landscape variables from both regions that of what we hypothesized in the Front Range, we also inverted the were strongly correlated with many other variables, until all VIF tree cover resistance surface (i.e., higher tree cover = higher resis- scores were less than 10 (Row et al., 2017). Variables retained were tance), reran Circuitscape and MLPE analyses, and higher tree cover geographic distance, river and stream riparian corridors, roads, im- still showed significant negative correlations with gene flow in this pervious surface cover, tree canopy cover, vegetation density, and region. minimum temperature of the coldest month. However, vegetation density was still correlated with geographic distance in both regions, and impervious surface was correlated with geographic distance in 4 | D I S CU S S I O N the Western Slope (Table S4). We removed these variables as well, resulting in Pearson's r correlations less than .7 for all explanatory The apex predator puma (P. concolor) persists in many urbanized variables, and final model VIF scores of 4.07 in the Western Slope regions throughout its range, yet the localized effects of recent and 3.54 in the Front Range. Thus final MRDM and MLPE models for urban sprawl remain unclear. Here, we compared patterns of the Western Slope included geographic distance, tree canopy cover, landscape genomic connectivity and genetic diversity of pumas �| TRUMBO et al. 4935 across two regions that span an urban‐rural divide in Colorado, with urbanization on the Front Range, population‐level genetic USA. Landscape genomic connectivity patterns differed between diversity and inbreeding measures were similar to those on the regions, such that genetic distances were higher and urbaniza- rural Western Slope. This suggests that recent urban sprawl in the tion (i.e., percent impervious surface cover) restricted gene flow Colorado Front Range has not yet had a large impact on the ge- in the more urbanized Front Range, whereas forest and riparian netic diversity of pumas. This is in contrast to more isolated puma cover were most important for enhancing gene flow on the rural populations in other highly urbanized landscapes such as south- Western Slope. Despite finding reductions in gene flow associated ern California and Florida, which exhibit reduced genetic diversity TA B L E 4 Multiple regression on distance matrices (MRDM) landscape genomic results from the Western Slope and Front Range of Colorado Region Genetic distance Landscape factors Direction of effect r2 p Western Slope Dps tree cover + .07 0.001 r geo. dist. − .03 0.001 Dps geo. dist. − .05 0.001 r geo. dist. − .04 0.001 Front Range Note: Response variables were individual‐based genetic distances, i.e., proportion of shared alleles (Dps) and relatedness (r). Explanatory variables, after removing correlated variables, were the geographic (Euclidean) distance model (geo. dist.), percent impervious surface cover, percent tree canopy cover, river and stream riparian corridors, roads, and minimum temperature of the coldest month. Forward selection followed by backward elimination was performed, with 1,000 random permutations of the dependent distance matrix per step, using Bonferroni‐corrected p‐to‐enter and p‐to‐remove alpha values of .05. Standardized beta coefficients were used to assess the direction of effect of each landscape variable on gene flow. Only univariate models were supported. TA B L E 5 Colorado Maximum likelihood of population effects (MLPE) landscape genomic results from the Western Slope and Front Range of r2 Region Genetic distance Landscape factors Direction of effect Western Slope Dps Tree cover + .13 0 Tree cover + .14 2.5 Min. Temperature − Tree cover + .17 0 Geo. Dist. − Tree cover − .16 3.4 Tree cover + .16 3.4 Riparian areas + Geo. Dist. − .16 3.7 Tree cover + .16 3.8 .12 0 .13 4.1 .14 0 .14 4.9 r Front Range Dps r Min. Temperature − Geo. Dist. − Tree cover − Geo. Dist. − Tree cover − Impervious surface − Geo. Dist. − Tree cover − Geo. Dist. − Tree cover − Impervious surface − ΔBIC Note: Response variables were individual‐based genetic distances, i.e., proportion of shared alleles (Dps) and relatedness (r). Pairwise comparisons of individuals were controlled as a random effect. Fixed effects, after removing correlated variables, were the geographic (Euclidean) distance model (geo. dist.), percent impervious surface cover, percent tree canopy cover, vegetation density, river and stream riparian corridors, roads, and minimum temperature of the coldest month. Standardized beta coefficients were used to assess the direction of effect of each landscape variable on gene flow. Models reported are within the top 5 BIC units. Landscape factors are in order of standardized beta coefficients (largest to smallest). �4936 | TRUMBO et al. and strong evidence of inbreeding compared to Colorado pumas Range (Allendorf et al., 2013; Haasl & Payseur, 2011). As the human (Ernest et al., 2003, 2014; Johnson et al., 2010; Riley et al., 2014). population continues to expand, future urbanization could result in However, a smaller effective population size and higher among‐in- more fragmented populations and reductions in genetic diversity, as dividual genetic distances in the recently urbanized Front Range has been detected in other more urbanized landscapes like southern suggest habitat fragmentation has already impacted this popula- California and Florida (Ernest et al., 2003, 2014; Johnson et al., 2010; tion and could cause further reductions of genetic diversity as ur- Riley et al., 2014). banization continues to expand in Colorado (Theobald, 2005; US Despite similar geographic extents and levels of genetic diver- Census Bureau, 2017). Extensive losses of genetic diversity could sity in the Western Slope and Front Range, mean genetic distances eventually cause puma populations to decline, which in turn could among individuals were higher in the urban Front Range (Table 2), have important cascading effects into lower trophic levels, such as suggesting that fragmentation due to urbanization may be limiting overgrazing of vegetation by ungulate herbivores (Markovchick‐ puma dispersal and gene flow. In addition, a larger effective popu- Nicholls et al., 2008). lation size (Ne) of pumas was detected on the rural Western Slope (Ne = 81.2) compared to the urban Front Range (Ne = 48.1; Table 2), 4.1 | Population genomics and structure with the caveat that some assumptions of this estimator are violated in both regions (e.g., closed populations with no immigration, non‐ Using our genomic dataset of over 12,000 SNPs, we were able to overlapping generations). The effect of non‐overlapping generations distinguish the Western Slope and Front Range regions (i.e., K = 2; on Ne is difficult to predict (Waples et al., 2016), and this assump- Figures 1 and 2). Minimum temperature of the coldest month was tion is expected to be violated similarly in both the Western Slope also negatively associated with gene flow in one of the top landscape and Front Range populations. Immigration, however, is expected genomic models on the Western Slope (Table 5), suggesting there to downwardly bias Ne by creating linkage disequilibrium through may be restricted gene flow through high elevation, alpine tundra a multilocus Wahlund effect (Wahlund, 1928; Waples & England, habitats (McMahon et al., 2001). However, potential immigrants 2011). Thus, it is possible that the Front Range may be showing a and admixed individuals were identified moving in both directions lower Ne due to having more immigrants from outside populations (Figure 2), particularly using Admixture ancestry analysis, and DAPC than the Western Slope. This is possible, and perhaps likely, given the group assignments may be overfit given the method's approach to presence of more admixed individuals in the Front Range (Figure 2), minimize within population distances and maximize between popu- which could indicate more potential immigrants into this region. On lation distances (Jombart et al., 2015). In addition, overall genetic the other hand, if immigration rates are similar for both regions, the differentiation between the two populations was low (pairwise relatively smaller Front Range Ne may be due to: (a) urbanization and FST = 0.02; Table 2). Because our sample archive consisted of op- fragmentation impacting and limiting population size, and/or (b) spe- portunistically collected samples, our analyses were restricted to cies range limit theory (Abundant Center Hypothesis) predicting that populations in two distinct regions, whereas pumas occur through- smaller population sizes are likely to occur at the edge of the geo- out the southern Rocky Mountains in Colorado. Therefore, potential graphic range relative to core areas (Brown, 1984; Sagarin & Gaines, immigrants and admixed individuals are not necessarily moving be- 2002). These potential underlying factors are not mutually exclusive tween our specific Western Slope and Front Range study areas, but and may both be acting together. However, the lack of difference in may originate from other unsampled populations that share genetic most genetic diversity measures, in addition to slightly lower allelic ancestry with our two study regions. Nevertheless, results from our richness in the Front Range, which is the most sensitive metric to study suggest pumas maybe somewhat limited in dispersing across recent bottlenecks (Allendorf et al., 2013), suggests lower effective the high elevation peaks of the Continental Divide, and future stud- population size on the Front Range may be more consistent with re- ies should attempt to sample more intensively across the entire re- cent urbanization impacts than historical range boundary effects. gion to further investigate this trend. We identified similar levels of genetic diversity on the rural Western Slope compared to the more urbanized Front Range, al- 4.2 | Landscape genomics though Hexp, Ar, and FIS were slightly higher on the Western Slope With regard to general landscape genomics methodology, we found (Table 2). This suggests urbanization is not yet having a major impact MRDM to be a much more conservative approach that adds fewer on the genetic diversity of pumas in Colorado. One potential expla- explanatory variables to the models than MLPE (Tables 4 and 5). nation is that urbanization in the Front Range is primarily occurring Therefore only the strongest landscape genomic relationships were on the eastern edge of the region, possibly creating a relatively im- identified using MRDM, consisting of isolation by geographic dis- permeable urban boundary on the eastern border, but not isolating tance in both regions, as well as tree cover in the Western Slope. pumas in fragments or limiting their connectivity to wildland habitat Conversely, MLPE results in more complex models with more ex- to the west (Figure 1; Blecha et al., 2018; Lewis et al., 2015). Another planatory variables and higher r2 values (genetic variation explained) possibility is that many of the SNPs we sampled may not have high than MRDM (Tables 4 and 5). The different genetic distance meas- enough mutation rates to show a strong genomic signature of the ures we used (D PS and r) showed largely consistent relationships relatively recent effects of rapid urbanization occurring in the Front with landscape variables, but still provided a few different insights, �| TRUMBO et al. particularly using MLPE (Tables 4 and 5). Overall r2 values were 2 2 somewhat low (r = .03–.07 for MRDM, r = .12–.17 for MLPE), but this is expected for a large carnivore with extreme long distance dispersal abilities (e.g., Balkenhol et al., 2016; Short Bull et al., 2011). 4937 part by higher traffic mortality in the more urbanized region (Beier, 1995; Crooks, 2002). It is important to note that unsampled landscape variables for which we have no data, as well as correlated landscape variables On the rural Western Slope, tree canopy cover was most import- removed from the final models (Table S4), may also be contribut- ant for enhancing gene flow, suggesting pumas prefer to disperse ing to landscape genomic patterns. For example, we don't have data through forests rather than more open shrub and grassland habi- on snowpack or prey (e.g., mule deer) abundance, and elevation and tats in this landscape (Table 5). In addition, riparian areas were an annual precipitation were removed from final models because they important predictor in one top model (Table 5), further supporting were strongly correlated with several other landscape variables, in- the importance of tree cover for enhancing gene flow in this region. cluding tree cover (Table S4). Thus the unexpected negative rela- Forests and riparian areas provide more tree cover for concealment tionship of puma gene flow with tree cover on the Front Range may and ambush predation (Hornocker & Negri, 2009; Logan & Sweanor, also be due to pumas utilizing lower elevation areas with less snow- 2001; Warren, Wallin, Beausoleil, & Warheit, 2016). Use of open pack and higher prey abundance for dispersal and hunting during areas may also increase susceptibility to mortality by hunters and the late fall, winter, and early spring months. Moreover Row et al. ranchers (Newby et al., 2013), which are both more prevalent in the (2017) found that MLPE models were able to distinguish the true un- rural Western Slope than the more urbanized Front Range. In ad- derlying dispersal models in approximately 75% of their simulations dition, non‐forested and non‐riparian areas on the Western Slope when multiple, correlated landscape predictor variables influenced are dry, high elevation desert habitats (i.e., the Colorado Plateau dispersal. This highlights the inherent difficulty of identifying the un- ecoregion; McMahon et al., 2001), which may provide less prey and derlying landscape drivers of gene flow in empirical systems, while water resources, and thus be poorer habitats for hunting and disper- also underscoring the importance of choosing landscape genomic sal (Dickson, Roemer, McRae, & Rundall, 2013; McRae et al., 2005; hypotheses judiciously and interpreting the resulting landscape ge- Sweanor, Logan, & Hornocker, 2000). nomic associations carefully. In the more urbanized Front Range, impervious surface cover In conclusion, our findings are consistent with prior comparative restricted gene flow (Table 5). This suggests urbanization is limiting landscape genetic studies that have revealed varying effects of land- gene flow, despite high levels of genetic diversity (Table 2). Similarly, scape factors on movement and gene flow across different portions Lewis et al. (2015) found pumas were less likely to be detected in of a species' geographic range (e.g., Short Bull et al., 2011; Trumbo, habitats with residential development, even low‐density exurban Spear, Baumsteiger, & Storfer, 2013; Vandergast et al., 2007). We developments, which are increasingly encroaching into the foothills found that in the rural Western Slope with high hunting pressure, of the Front Range region. Genetic studies on pumas from more forests and riparian areas with high tree canopy cover are most im- urbanized and fragmented populations in southern California and portant for conserving puma genetic connectivity. In contrast, in the Florida have detected strong inbreeding and isolation associated more urbanized Front Range, non‐forested habitats such as shru- with urbanization (Ernest et al., 2003, 2014; Johnson et al., 2010; bland and grassland habitats are utilized more for dispersal and gene Riley et al., 2014). Our study detected more subtle impacts of urban- flow, effective population sizes are smaller, genetic distances among ization in a less fragmented landscape, within mountainous wildland individuals are higher, and gene flow is being restricted by urbaniza- habitats adjacent to a major metropolitan center, which experiences tion (Tables 2, 4, and 5). Next generation sequencing techniques can high levels of human outdoor recreation activities such as hiking and provide dense, genome‐scale SNP data sets of thousands of puta- skiing (Figure 1). In addition, in contrast with the rural Western Slope tively neutral markers, which gives researchers increased power to and contrary to our initial hypotheses, forest cover was negatively detect the often subtle effects of landscape factors, such as urban- associated with gene flow on the Front Range (Table 5). This pattern ization, on gene flow (Allendorf et al., 2013; Lowe & Allendorf, 2010; suggests pumas are more willing to disperse through open shrub and Luikart et al., 2003). This is particularly important for wide‐ranging grassland habitats in this region. The reasons for this are unclear, species with broad geographic distributions, as landscape effects on but pumas living in the more developed Front Range may be more gene flow occur at broader geographic scales and may be weaker acclimated to human activities and thus more willing to travel out- and more difficult to detect compared to more dispersal‐limited side of forested habitats, demonstrating that pumas have a range of species with smaller home ranges (Balkenhol et al., 2016; Epps et adaptable behaviors and will use and move through different types al., 2007; Holderegger et al., 2006). Indeed, prior work on pumas of habitat (Blecha et al., 2018; Dickson, Jenness, & Beier, 2005). using 16 microsatellites found no population structure across the Pumas may also be hunting more urban mesopredators, domestic, southern Rocky Mountains of Colorado and northern New Mexico and agricultural animals in these open habitats on the more devel- (McRae et al., 2005). Our results demonstrate that large SNP data oped Front Range, which was shown in a prior study using stable sets can allow researchers to identify impacts of urbanization on isotope analysis of Front Range puma diets (Moss, Alldredge, & Pauli, gene flow, effective population sizes, and patterns of population 2016). There is also less hunting of pumas in the Front Range com- genetic structure of wide‐ranging species, even before fragmen- pared to the rural Western Slope, so pumas may be less wary of open tation is extensive enough to cause substantial declines in genetic areas, although this effect would be expected to be counteracted in diversity. Maintaining genetic connectivity in these umbrella species �4938 | TRUMBO et al. can have outsized benefits towards conserving biodiversity, as preserving broad swathes of contiguous habitats that are necessary for their persistence also benefits many other species with smaller home ranges and narrower habitat requirements (Sergio et al., 2008, 2006; Thorne et al., 2006). AC K N OW L E D G E M E N T S Funding was provided by the National Science Foundation, Ecology of Infectious Disease Program (NSF‐EID 1413925 and 723676). Samples were collected by Colorado Parks and Wildlife. We also thank Michael Antolin, Kelly Pierce, and Jill Gerberich at Colorado State University for assistance in the laboratory. AU T H O R C O N T R I B U T I O N S D.R.T. performed laboratory work, analyzed landscape and population genomic data, and wrote the manuscript; P.S., R.B.G., C.P.K, S.K., and N.F.J. performed laboratory work and analyzed landscape and population genomic data; K.L., and M.A. directed fieldwork and collected field data; M.E.C., S.C., H.B.E., K.C., S.V., and W.C.F. conceived of study questions and directed research; and all authors contributed input to draft and final versions of the manuscript. DATA AVA I L A B I L I T Y S TAT E M E N T ddRADseq data used in genomic analyses are on Dryad (https​://doi. org/10.5061/dryad.12jm6​3xsr). ORCID Daryl R. Trumbo https://orcid.org/0000-0002-8438-2856 Christopher P. Kozakiewicz https://orcid. org/0000-0002-4868-9252 Nicholas M. Fountain‐Jones https://orcid. org/0000-0001-9248-8493 Holly B. Ernest W. Chris Funk https://orcid.org/0000-0002-0205-8818 https://orcid.org/0000-0002-6466-3618 REFERENCES Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model‐based estimation of ancestry in unrelated individuals. Genome Research, 19, 1655–1664. https​://doi.org/10.1101/gr.094052.109 Allendorf, F. W., Luikart, G., & Aitken, S. N. (2013). Conservation and the genetics of populations. Chichester, UK: Wiley‐Blackwell. Balkenhol, N., Cushman, S. A., Storfer, A. T., & Waits, L. P. (2016). Landscape genetics: Concepts, methods, applications. Chichester, UK: Wiley‐Blackwell. Beier, P. (1995). Dispersal of juvenile cougars in fragmented habitat. Journal of Wildlife Management, 59, 228–237. https​://doi. org/10.2307/3808935 Blecha, K. A., Boone, R. B., & Alldredge, M. W. (2018). Hunger mediates apex predator's risk avoidance response in wildland–urban interface. Journal of Animal Ecology, 87, 609–622. https​://doi. org/10.1111/1365-2656.12801​ Bowcock, A. M., Ruiz‐Linares, A., Tomfohrde, J., Minch, E., Kidd, J. R., & Cavalli‐Sforza, L. L. (1994). High resolution of human evolutionary trees with polymorphic microsatellites. Nature, 368, 455–457. https​ ://doi.org/10.1038/368455a0 Brown, J. H. (1984). On the relationship between abundance and distribution of species. American Naturalist, 124, 255–279. https​://doi. org/10.1086/284267 Carver, S., Bevins, S. N., Lappin, M. R., Boydston, E. E., Lyren, L. M., Alldredge, M., … Vandewoude, S. (2016). Pathogen exposure varies widely among sympatric populations of wild and domestic felids across the United States. Ecological Applications, 26, 367–381. Catchen, J. M., Hohenlohe, P. A., Bassham, S., Amores, A., & Cresko, W. A. (2013). Stacks: An analysis tool set for population genomics. Molecular Ecology, 22, 3124–3140. https​://doi.org/10.1111/ mec.12354​ Clarke, R. T., Rothery, P., & Raybould, A. F. (2002). Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural Biological and Environmental Statistics, 7, 361–372. https​://doi.org/10.1198/10857​1102320 Cohen, J. E. (2003). Human population: The next half century. Science, 302, 1172–1175. https​://doi.org/10.1126/scien​ce.1088665 Crooks, K. R. (2002). Relative sensitivities of mammalian carnivores to habitat fragmentation. Conservation Biology, 16, 488–502. https​:// doi.org/10.1046/j.1523-1739.2002.00386.x Crooks, K. R., Burdett, C. L., Theobald, D. M., King, S. R. B., Di Marco, M., Rondini, C., & Boitani, L. (2017). Quantification of habitat fragmentation reveals extinction risk in terrestrial mammals. Proceedings of the National Academy of Sciences, 114, 7635–7640. https​://doi. org/10.1073/pnas.17057​69114​ Cushman, S. A., Gutzweiler, K., Evans, J. S., & McGarigal, K. (2010). The gradient paradigm: A conceptual and analytical framework for landscape ecology. In S. A. Cushman, & F. Huettman (Eds.), Spatial com‐ plexity, informatics, and wildlife conservation. New York, NY: Springer. Dickson, B. G., Jenness, J. S., & Beier, P. (2005). Influence of vegetation, topography, and roads on cougar movement in southern California. Journal of Wildlife Management, 69, 264–276. https​://doi.org/10.2193 /0022-541X(2005)069<0264:IOVTA​R>2.0.CO;2 Dickson, B. G., Roemer, G. W., McRae, B. H., & Rundall, J. M. (2013). Models of regional habitat quality and connectivity for pumas (Puma concolor) in the southwestern United States. PLoS ONE, 8, 1–11. https​ ://doi.org/10.1371/journ​al.pone.0081898 Do, C., Wapels, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., & Ovenden, J. R. (2014). NeEstimator v2; re‐implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources, 14, 209–214. Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., … Lautenbach, S. (2012). Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36, 27–46. https​://doi.org/10.1111/j.1600-0587.2012.07348.x Epps, C. W., Wehausen, J. D., Bleich, V. C., Torres, S. G., & Brashares, J. S. (2007). Optimizing dispersal and corridor models using landscape genetics. Journal of Applied Ecology, 44, 714–724. https​://doi. org/10.1111/j.1365-2664.2007.01325.x Ernest, H. B., Boyce, W. M., Bleich, V. C., May, B., Stiver, S. J., & Torres, S. G. (2003). Genetic structure of mountain lion (Puma concolor) populations in California. Genetics, 4, 353–366. Ernest, H. B., Vickers, T. W., Morrison, S. A., Buchalski, M. R., & Boyce, W. M. (2014). Fractured genetic connectivity threatens a southern California puma (Puma concolor) population. PLoS ONE, 9, 1–12. https​ ://doi.org/10.1371/journ​al.pone.0107985 Fountain‐Jones, N. M., Craft, M. E., Funk, W. C., Kozakiewicz, C., Trumbo, D. R., Boydston, E. E., … Carver, S. (2017). Urban landscapes can change virus gene flow and evolution in a fragmentation‐sensitive �TRUMBO et al. carnivore. Molecular Ecology, 26, 6487–6498. https​://doi.org/10.1111/ mec.14375​ Haasl, R. F., & Payseur, B. J. (2011). Multi‐locus inference of population structure: A comparison between single nucleotide polymorphisms and microsatellites. Heredity, 106, 158–171. https​://doi.org/10.1038/ hdy.2010.21 Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., & Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 1965–1978. https​:// doi.org/10.1002/joc.1276 Hilty, J. A., & Merenlender, A. M. (2004). Use of riparian corridors and vineyards by mammalian predators in northern California. Conservation Biology, 18, 126–135. https​://doi. org/10.1111/j.1523-1739.2004.00225.x Holderegger, R., Kamm, U., & Gugerli, F. (2006). Adaptive vs. neutral genetic diversity: Implications for landscape genetics. Landscape Ecology, 21, 797–807. https​://doi.org/10.1007/s10980-005-5245-9 Homer, C., Huang, C., Yang, L., Wylie, B. K., & Coan, M. (2004). Development of a 2001 national land cover database for the United States. Photogrammetric Engineering and Remote Sensing, 70, 829– 840. https​://doi.org/10.14358/​PERS.70.7.829 Hornocker, M., & Negri, S. (2009). Cougar: Ecology and conservation. Chicago, IL: University of Chicago Press. Johnson, W. E., Onorato, D. P., Roelke, M. E., Land, E. D., Cunningham, M., Belden, R. C., … O'Brien, S. J. (2010). Genetic restoration of the Florida panther. Science, 329, 1641–1645. https​://doi.org/10.1126/ scien​ce.1192891 Jombart, T. (2008). adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics, 24, 1403–1405. https​://doi. org/10.1093/bioin​forma​tics/btn129 Jombart, T., Devillard, S., & Balloux, F. (2015). Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genetics, 11(94), 1–15. Kalinowski, S. T. (2005). hp‐rare 1.0: A computer program for performing rarefaction on measures of allelic richness. Molecular Ecology Notes, 5, 187–189. https​://doi.org/10.1111/j.1471-8286.2004.00845.x Legendre, P., Lapointe, F. J., & Casgrain, P. (1994). Modeling brain evolution from behavior: A permutational regression approach. Evolution, 48, 1487–1499. https​://doi.org/10.1111/j.1558-5646.1994.tb021​91.x Lewis, J. S., Logan, K. A., Alldredge, M. W., Bailey, L. L., Vandewoude, S., & Crooks, K. R. (2015). The effects of urbanization on population density, occupancy, and detection probability of wild felids. Ecological Applications, 25, 1880–1895. https​://doi.org/10.1890/14-1664.1 Logan, K. A., & Sweanor, L. L. (2001). Desert puma: Evolutionary ecology and conservation of an enduring carnivore. Washington, DC: Island Press. Lowe, W. H., & Allendorf, F. W. (2010). What can genetics tell us about population connectivity? Molecular Ecology, 19, 3038–3051. https​:// doi.org/10.1111/j.1365-294X.2010.04688.x Luikart, G., England, P. R., Tallmon, D., Jordan, S., & Taberlet, P. (2003). The power and promise of population genomics: From genotyping to genome typing. Nature Reviews Genetics, 4, 981–994. https​://doi. org/10.1038/nrg1226 Luu, K., Bazin, E., & Blum, M. G. B. (2016). PCAdapt: An r package to perform genome scans for selection based on principal component analysis. Molecular Ecology Resources, 17, 67–77. Magle, S. B., Hunt, V. M., Vernon, M., & Crooks, K. R. (2012). Urban wildlife research: Past, present, and future. Biological Conservation, 155, 23–32. https​://doi.org/10.1016/j.biocon.2012.06.018 Maletzke, B., Kertson, B., Swanson, M., Koehler, G., Beausoleil, R., Wielgus, R., & Cooley, H. (2017). Cougar response to a gradient of human development. Ecosphere, 8, 1–14. https​://doi.org/10.1002/ ecs2.1828 Markovchick‐Nicholls, L., Regan, H. M., Deutschman, D. H., Widyanata, A., Martin, B., Noreke, L., & Hunt, T. A. (2008). Relationships | 4939 between human disturbance and wildlife land use in urban habitat fragments. Conservation Biology, 22, 99–109. https​://doi. org/10.1111/j.1523-1739.2007.00846.x Mastretta‐Yanes, A., Arrigo, N., Alvaraz, N., Jorgensen, T. H., Piñero, D., & Emerson, B. C. (2015). Restriction site‐associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Molecular Ecology, 15, 28–41. McKinney, M. L. (2002). Urbanization, biodiversity, and conservation. BioScience, 52, 883–890. https​://doi. org/10.1641/0006-3568(2002)052[0883:UBAC]2.0.CO;2 McMahon, G., Gregonis, S. M., Waltman, S. W., Omernik, J. M., Thorson, T. D., Freeouf, J. A., … Keys, J. E. (2001). Developing a spatial framework of common ecological regions for the conterminous United States. Environmental Management, 28, 293–316. https​://doi. org/10.1007/s0026​702429 McRae, B. H. (2006). Isolation by resistance. Evolution, 60, 1551–1561. https​://doi.org/10.1111/j.0014-3820.2006.tb005​0 0.x McRae, B. H., Beier, P., Dewald, L. E., Huynh, L. Y., & Keim, P. (2005). Habitat barriers limit gene flow and illuminate historical events in a wide‐ranging carnivore, the American puma. Molecular Ecology, 14, 1965–1977. https​://doi.org/10.1111/j.1365-294x.2005.02571.x Moss, W. E., Alldredge, M. W., Logan, K. A., & Pauli, J. N. (2016). Human expansion precipitates niche expansion for an opportunistic apex predator (Puma concolor). Scientific Reports, 6, 1–5. https​://doi. org/10.1038/srep3​9639 Moss, W. E., Alldredge, M. W., & Pauli, J. N. (2016). Quantifying risk and resource use for a large carnivore in an expanding urban‐wildland interface. Journal of Applied Ecology, 53, 371–378. https​://doi. org/10.1111/1365-2664.12563​ Naiman, R. J., Decamps, H., & Polluck, M. (1993). The role of riparian corridors in maintaining regional biodiversity. Ecological Applications, 3, 209–212. https​://doi.org/10.2307/1941822 Newby, J. R., Scott Mills, L., Ruth, T. K., Pletscher, D. H., Mitchell, M. S., Quigley, H. B., … DeSimone, R. (2013). Human‐caused mortality influences spatial population dynamics: Pumas in landscapes with varying mortality risks. Biological Conservation, 159, 230–239. https​://doi.org/10.1016/j.biocon.2012.10.018 Paris, J. R., Stevens, J. R., & Catchen, J. M. (2017). Lost in parameter space: A road map for Stacks. Methods in Ecology and Evolution, 8, 1360–1373. Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., & Hoekstra, H. E. (2012). Double digest RADseq: An inexpensive method for denovo SNP discovery and genotyping in model and nonmodel species. PLoS ONE, 7, 1–11. https​://doi.org/10.1371/journ​al.pone.0037135 Radeloff, V. C., Hamer, R. B., Stewart, S. I., Fried, J. S., Holcomb, S. S., & McKeefry, J. F. (2005). The wildland‐urban interface in the United States. Ecological Applications, 15, 799–805. https​://doi. org/10.1890/04-1413 Richards, S. A. (2015). Likelihood and model selection. In G. A. Fox, S. Negrete‐Yankelevich, & V. J. Sosa (Eds.), Ecological statistics: Contemporary theory and application (pp. 58–80). Oxford, UK: Oxford University Press. Riley, S. P., Pollinger, J. P., Sauvajot, R. M., York, E. C., Bromley, C., Fuller, T. K., & Wayne, R. K. (2006). A southern California freeway is a physical and social barrier to gene flow in carnivores. Molecular Ecology, 15, 1733–1741. https​://doi.org/10.1111/j.1365-294X.2006.02907.x Riley, S. P., Serieys, L. E. K., Pollinger, J. P., Sikich, J. A., Dalbeck, L., Wayne, R. K., & Ernest, H. B. (2014). Individual behavior dominates the dynamics of a mountain lion population in a highly fragmented urban landscape. Current Biology, 24, 1989–1994. Row, J. R., Knick, S. T., Oyler‐McCance, S. J., Lougheed, S. C., & Fedy, B. C. (2017). Developing approaches for linear mixed modeling in landscape genetics through landscape‐directed dispersal simulations. Ecology and Evolution, 7, 3751–3761. https​://doi.org/10.1002/ ece3.2825 �4940 | Sagarin, R. D., & Gaines, S. D. (2002). The abundant centre distribution: To what extent is it a biogeographical rule? Ecology Letters, 5, 137– 147. https​://doi.org/10.1046/j.1461-0248.2002.00297.x Sergio, F., Caro, T., Brown, D., Clucas, B., Hunter, J., Ketchum, J., … Hiraldo, F. (2008). Top predators as conservation tools: Ecological rationale, assumptions, and efficacy. Annual Review of Ecology, Evolution, and Systematics, 39, 1–19. https​://doi.org/10.1146/annur​ ev.ecols​ys.39.110707.173545 Sergio, F., Newton, I., Marchesi, L., & Pedrini, P. (2006). Ecologically justified charisma: Preservation of top predators delivers biodiversity conservation. Journal of Applied Ecology, 43, 1049–1055. https​://doi. org/10.1111/j.1365-2664.2006.01218.x Shirk, A. J., Landguth, E. L., & Cushman, S. A. (2017). A comparison of regression methods for model selection in individual‐based landscape genetic analysis. Molecular Ecology Resources, 18, 55–67. https​://doi. org/10.1111/1755-0998.12709​ Short Bull, R. A., Cushman, S. A., Mace, R., Chilton, T., Kendall, K. C., Landguth, E. L., … Luikart, G. (2011). Why replication is important in landscape genetics: American black bear in the Rocky Mountains. Molecular Ecology, 6, 1092–1107. https​://doi. org/10.1111/j.1365-294X.2010.04944.x Smouse, P. E., & Peakall, R. (1999). Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity, 82, 561–573. https​://doi.org/10.1038/sj.hdy.6885180 Solberg, K. H., Bellemain, E., Drageset, O. M., Taberlet, P., & Swenson, J. E. (2006). An evaluation of field and non‐invasive genetic methods to estimate brown bear (Ursus arctos) population size. Biological Conservation, 128, 158–168. https​://doi.org/10.1016/j. biocon.2005.09.025 Sweanor, L. L., Logan, K. A., & Hornocker, M. G. (2000). Cougar dispersal patterns, metapopulation dynamics, and conservation. Conservation Biology, 14, 798–808. https​://doi. org/10.1046/j.1523-1739.2000.99079.x Theobald, D. M. (2005). Landscape patterns of exurban growth in the USA from 1980 to 2020. Ecological Society, 10, 1–34. Thorne, J. H., Cameron, D., & Quinn, J. F. (2006). A conservation design for the central coast of California and the evaluation of mountain lion as an umbrella species. Natural Areas Journal, 26, 137– 148. https​://doi.org/10.3375/0885-8608(2006)26[137:ACDFT​ C]2.0.CO;2 Trumbo, D. R., Spear, S. F., Baumsteiger, J., & Storfer, A. (2013). Rangewide landscape genetics of an endemic Pacific northwestern salamander. Molecular Ecology, 22, 1250–1266. https​://doi.org/10.1111/ mec.12168​ U.S. Census Bureau (2017). Cumulative estimates of resident population change for the United States, regions, states, and Puerto Rico and re‐ gion and state rankings: April 1, 2010 to July 1, 2017, Table 2 (NST‐ EST2017‐02). Source: U.S. Census Bureau, Population Division; Release date: December 2017. Retrieved from www.census.gov/ progr​ams-surve​ys/popest.html van Strien, M. J., Keller, D., & Holderegger, R. (2012). A new analytical approach to landscape genetic modelling: Least‐cost transect analysis TRUMBO et al. and linear mixed models. Molecular Ecology, 21, 4010–4023. https​:// doi.org/10.1111/j.1365-294X.2012.05687.x Vandergast, A. G., Bohanak, A. J., Weissman, D. B., & Fisher, R. N. (2007). Understanding the genetic effects of recent habitat fragmentation in the context of evolutionary history: Phylogeography and landscape genetics of a southern California endemic Jerusalem cricket (Orthoptera: Stenopelmatidae: Stenopelmatus). Molecular Ecology, 16, 977–992. https​://doi.org/10.1111/j.1365-294X.2006.03216.x Wahlund, S. (1928). Zuzammensetzung von populationen und korrelation‐serscheiunungen von standpunkt der vererbungslehre aus betrachtet. Hereditas, 11, 65–106. Waples, R. S. (2016). Making sense of genetic estimates of effective population size. Molecular Ecology, 25, 4689–4691. https​://doi. org/10.1111/mec.13814​ Waples, R. S., & England, P. R. (2011). Estimating contemporary effective population size on the basis of linkage disequilibrium in the face of migration. Genetics, 189, 633–644. https​://doi.org/10.1534/genet​ ics.111.132233 Waples, R. K., Larson, W. A., & Waples, R. S. (2016). Estimating contemporary effective population size in non‐model species using linkage disequilibrium across thousands of loci. Heredity, 117, 233–240. https​ ://doi.org/10.1038/hdy.2016.60 Warren, D. L., Glor, R. E., & Turelli, M. (2010). ENMTools: A toolbox for comparative studies of environmental niche models. Ecography, 33, 607–611. https​://doi.org/10.1111/j.1600-0587.2009.06142.x Warren, M. J., Wallin, D. O., Beausoleil, R. A., & Warheit, K. I. (2016). Forest cover mediates genetic connectivity of northwestern cougars. Conservation Genetics, 17, 1011–1024. https​://doi.org/10.1007/ s10592-016-0840-7 Wilmers, C. C., Wang, Y., Nickel, B., Houghtaling, P., Shakeri, Y., Allen, M. L., … Williams, T. (2013). Scale dependent behavioral responses to human development by a large predator, the puma. PLoS ONE, 8, 1–11. https​://doi.org/10.1371/journ​al.pone.0060590 Wright, S. (1942). Isolation by distance. Genetics, 28, 114–138. Zeller, K. A., Vickers, T. W., Ernest, H. B., & Boyce, W. M. (2017). Multi‐ level, multi‐scale resource selection functions and resistance surfaces for conservation planning: Pumas as a case study. PLoS ONE, 12, 1–20. https​://doi.org/10.1371/journ​al.pone.0179570 S U P P O R T I N G I N FO R M AT I O N Additional supporting information may be found online in the Supporting Information section.   How to cite this article: Trumbo DR, Salerno PE, Logan KA, et al. Urbanization impacts apex predator gene flow but not genetic diversity across an urban‐rural divide. Mol Ecol. 2019;28:4926–4940. https​://doi.org/10.1111/mec.15261​ � https://cpw.cvlcollections.org/files/original/f6dd1182c74ad3f028a1ba06910b7477.pdf e28ad57d4017ea3b92b438d5061d229a PDF Text Text (a) Western Slope (b) Front Range � https://cpw.cvlcollections.org/files/original/0e11354ebf558be3d8ab1e2c83c2b5e5.pdf 771b73b5b090882c728568e12dd53ec3 PDF Text Text 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Supporting Information Manuscript Title: Urbanization impacts apex predator gene flow but not genetic diversity across an urban-rural divide Authors: Trumbo DR1, Salerno PE1, Logan K2, Alldredge M3, Gagne RB4, Kozakiewicz CP5, Kraberger S4, Fountain-Jones N6, Craft ME6, Carver S5, Ernest HB7, Crooks K8, VandeWoude S4, Funk WC1,9 Department of Biology, Colorado State University, Fort Collins, CO 80523 USA Colorado Parks and Wildlife, Montrose, CO 81401 USA 3 Colorado Parks and Wildlife, Fort Collins, CO 80526 USA 4 Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, CO 80523 USA 5 Department of Biological Sciences, University of Tasmania, Hobart, TAS 7005 Australia 6 Department of Veterinary Population Medicine, University of Minnesota, Saint Paul, MN 55108 USA 7 Department of Veterinary Sciences, University of Wyoming, Laramie, WY 82070 USA 8 Department of Fish, Wildlife, and Conservation Biology, Colorado State University, Fort Collins, CO 80523 USA 9 Graduate Degree Program in Ecology, Colorado State University, Fort Collins, CO 80523 USA 1 2 �18 19 20 21 22 23 24 25 Table S1: Library replicate analysis of error rates from different Stacks parameter settings for minimum number of identical raw reads required to create a stack (-m), number of mismatches allowed between loci when processing a single individual (-M), number of mismatches allowed between loci when building the catalog (-n), and maximum number of stacks at a single de novo locus (max_locus_stacks). Locus error rate was the number of loci present in only one of the samples of a replicate pair divided by the total number of loci, allele error rate was the number of allele mismatches between replicate pairs divided by the number of loci, and SNP error rate was the proportion of SNP mismatches between replicate pairs. Final parameter settings are highlighted in bold. -m 3 3 3 3 3 3 3 3 3 5 7 Stacks parameter settings -M -n -max_locus_stacks 2 2 3 2 2 4 2 2 5 2 4 3 2 6 3 2 8 3 4 2 3 6 2 3 8 2 3 2 2 3 2 2 3 Error rates Total number of SNPs Locus error Allele error SNP error 0.217 0.129 0.0758 29,448 0.217 0.129 0.0761 29.434 0.212 0.129 0.0761 29,412 0.217 0.128 0.0740 29,814 0.216 0.128 0.0732 28,984 0.216 0.128 0.0730 28,507 0.217 0.128 0.0772 30,658 0.216 0.127 0.0763 29,636 0.217 0.127 0.0756 28,796 0.217 0.129 0.0758 23,062 0.219 0.065 0.0582 18,541 2 �26 27 Table S2: Landscape resistance transformations for the Western Slope and Front Range. Landscape Resistances Western Slope Min Max Median Mean Std Dev Elevation (m) 1 2910 2009 1974.4 Tree canopy cover (%) 1 101 81 Impervious surface (%) 1 101 Minimum temp. coldest month (°C) 1 Annual precipitation (mm) Front Range Transformations Min Max Median Mean Std Dev 529.3 1 2874 1984 1973.6 629.3 Convert to integer, subtract minimum value, add 1, invert raster. 71.1 31.4 1 101 69 66.0 33.6 Add 1, invert raster. 1 1.5 4.1 1 101 1 5.0 13.5 Add 1. 108 39 44.5 28.1 1 117 45 44.0 28.0 Multiply by 10, add minimum value, add 1, invert raster. 1 930 680 654.7 171.5 1 648 555 510.6 121.9 Subtract minimum value, add 1, invert raster. Enhanced vegetation index 1 10814 4374 4573.1 1858.3 1 11102 3717 4175.2 1800.9 Add minimum value, add 1, invert raster. Topographic roughness 1 27925 11 50.5 129.3 1 20068 20043 20010.7 101.8 Convert to integer, add 1, invert raster. Landcover 1 10 1 3.2 3.0 1 10 1 4.2 3.7 Forested habitats = 1, natural, nonforested habitats (e.g., grasslands, shrublands) = 5, human developed, barren, and agricultural areas = 10. Roads 1 10 1 1.7 2.5 1 10 1 2.7 3.6 Roads = 10, non-roads = 1. Riparian 1 10 10 9.4 2.3 1 10 10 9.4 2.3 Riparian = 1, non-riparian = 10. Raw Landscape Data Western Slope Min Max Median Mean Std Dev 1453.5 4362.9 2354 2418.0 Tree canopy cover (%) 0 100 20 Impervious surface (%) 0 100 Minimum temp. coldest month (°C) -20.2 Annual precipitation (mm) Enhanced vegetation index Elevation (m) Front Range Hypotheses** Min Max Median Mean Std Dev 552.5 1474.5 4347.1 2365 2374.9 629.3 29.9 31.4 0 100 32 35.0 33.6 0 0.5 4.1 0 100 0 4.0 13.5 Higher impervious surface cover has higher resistance to gene flow. -9.5 -13.3 -13.9 2.8 -19.9 -8.3 -12.7 -12.6 2.8 Higher minimum temperature has lower resistance to gene flow. 208 1137 458 483.3 171.5 359 1006 452 496.4 121.9 Higher precipitation has lower resistance to gene flow. -1806 8955 4634 4434.9 1858.3 -1969 9132 5416 4957.8 1800.9 Higher vegetation density has lower resistance to gene flow. 3 Higher elevation has lower resistance to gene flow. Higher tree cover has lower resistance to gene flow. �Topographic roughness 28 Higher topographic roughness has lower resistance to gene flow. Forested habitats have lowest resistance; natural, non-forested habitats (e.g., grasslands, shrublands) have intermediate resistance; and human developed, barren, and agricultural areas have highest resistance to gene flow. 0 27924.6 11 53.1 129.6 0 20067.0 25 56.2 100.8 Landcover* N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Roads* N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Roads have high resistance to gene flow. Riparian* N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Riparian corridors have low resistance to gene flow. *Categorical variables; no min, max, median, mean, or std dev ** Ecological justifications for each hypothesis are described in Table 1. 4 �29 30 31 32 33 Table S3: Correlations (Pearson’s r) between environmental raster surfaces used in landscape genomic analyses for the Western Slope and Front Range regions of Colorado. Pearson’s correlations > 0.7 are in bold. Western Slope Land cover Land cover 1 Impervious Roads Riparian Tree cover Veg. density Min. temp. Ann. precip. Topo. rough. Elevation Front Range Land cover Land cover 1 Impervious Roads Riparian Tree cover Veg. density Min. temp. Ann. precip. Topo. rough. Elevation Impervious Roads Riparian -0.18 1 -0.12 0.31 1 0.02 0.01 0.02 1 Tree cover 0.64 -0.08 -0.08 0.01 1 Veg. density 0.30 -0.01 0.04 0.04 0.32 Min. temp. -0.15 0.09 0.06 0.05 -0.29 Ann. precip. 0.11 -0.10 -0.09 -0.06 0.25 Topo. rough. -0.04 -0.04 -0.09 -0.02 0.08 1 -0.01 -0.16 -0.19 -0.10 1 -0.81 -0.27 -0.88 1 0.27 0.98 1 0.28 Elevation 0.14 -0.10 -0.09 -0.08 0.28 1 Impervious Roads Riparian -0.45 1 -0.32 0.47 1 -0.03 0.01 0.07 1 Tree cover 0.79 -0.27 -0.23 -0.04 1 Veg. density 0.15 -0.01 0.09 0.05 0.13 Min. temp. -0.37 0.33 0.30 0.09 -0.39 Ann. precip. 0.15 -0.21 -0.24 -0.08 0.19 Topo. rough. 0.05 -0.12 -0.15 -0.04 0.05 1 0.55 -0.71 -0.21 -0.60 1 -0.86 -0.23 -0.95 1 0.25 0.91 1 0.26 Elevation 0.40 -0.32 -0.31 -0.11 0.41 1 5 �34 35 36 Table S4: Correlations (Pearson’s r) between Circuitscape environmental resistances used in landscape genomic analyses for the Western Slope and Front Range regions of Colorado. Pearson’s correlations > 0.7 are in bold. Western Slope Null Geo. dist. Land cover Impervious Roads Riparian Tree cover Veg. density Min. temp. Ann precip. Topo. rough. Elevation Front Range 1 37 Impervious Roads Riparian 0.83 0.63 1 0.44 0.31 0.62 1 0.67 0.45 0.50 0.28 1 Tree cover 0.68 0.80 0.60 0.27 0.45 1 Veg. density 0.76 0.68 0.63 0.38 0.53 0.72 Min. temp. 0.61 0.35 0.49 0.28 0.54 0.19 Ann. precip. 0.95 0.62 0.81 0.43 0.58 0.73 Topo. rough. 0.99 0.60 0.83 0.44 0.67 0.68 1 0.37 0.80 0.76 0.78 1 0.38 0.61 0.38 1 0.95 0.99 1 0.95 Elevation 0.95 0.61 0.81 0.43 0.57 0.72 1 Null Geo. dist. Land cover Impervious Roads Riparian Tree cover Veg. density Min. temp. Ann. precip. Topo. rough. Elevation Land cover 0.60 1 1 Land cover 0.30 1 Impervious Roads Riparian 0.10 0.87 1 0.43 0.71 0.66 1 0.65 -0.04 -0.16 0.17 1 Tree cover 0.62 0.75 0.47 0.60 0.19 1 Veg. density 0.91 0.32 0.14 0.43 0.50 0.63 Min. temp. 0.50 -0.19 -0.22 0.15 0.61 0.04 Ann. precip. 0.98 0.33 0.13 0.42 0.58 0.65 Topo. rough. 0.99 0.30 0.10 0.43 0.65 0.63 1 0.47 0.89 0.90 0.76 1 0.35 0.50 0.08 1 0.98 0.94 1 0.89 Elevation 0.89 0.47 0.26 0.44 0.44 0.73 1 6 �38 39 40 41 Table S5: Effective population sizes (Ne) based on subsetting individuals into different time periods. Western Slope Ngen Ne (95% C.I.) Front Range Ngen Ne (95% C.I.) 2005-Dec.2010 34 indiv. 40.3 (38.7-41.9) 2007-May2010 25 indiv. 34.4 (33.0-35.8) Jan.2011-2014 42 indiv. 85.1 (81.8-88.4) June2010-2013 29 indiv. 53.4 (51.3-55.5) 7 �42 43 44 Figure S1: Principle Components Analyses (PCAs) and Admixture plots of (a) 76 Western Slope pumas and (b) 54 Front Range pumas, analyzed separately within each region. 45 8 � Dublin Core The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/. Title A name given to the resource Journal Articles Description An account of the resource CPW peer-reviewed journal publications Text A resource consisting primarily of words for reading. Examples include books, letters, dissertations, poems, newspapers, articles, archives of mailing lists. Note that facsimiles or images of texts are still of the genre Text. Dublin Core The Dublin Core metadata element set is common to all Omeka records, including items, files, and collections. For more information see, http://dublincore.org/documents/dces/. Title A name given to the resource Urbanization impacts apex predator gene flow but not genetic diversity across an urban‐rural divide Rights Information about rights held in and over the resource <a href="http://rightsstatements.org/vocab/InC-NC/1.0/" target="_blank" rel="noreferrer noopener">In Copyright - Non-Commercial Use Permitted</a> Date Created Date of creation of the resource. 2019-10-06 Subject The topic of the resource Effective population size Gene flow Landscape genomics <em>Puma concolor</em> Urbanization Genetic diversity Description An account of the resource <span>Apex predators are important indicators of intact natural ecosystems. They are also sensitive to urbanization because they require broad home ranges and extensive contiguous habitat to support their prey base. Pumas (</span><i>Puma concolor</i><span>) can persist near human developed areas, but urbanization may be detrimental to their movement ecology, population structure, and genetic diversity. To investigate potential effects of urbanization in population connectivity of pumas, we performed a landscape genomics study of 130 pumas on the rural Western Slope and more urbanized Front Range of Colorado, USA. Over 12,000 single nucleotide polymorphisms (SNPs) were genotyped using double-digest, restriction site-associated DNA sequencing (ddRADseq). We investigated patterns of gene flow and genetic diversity, and tested for correlations between key landscape variables and genetic distance to assess the effects of urbanization and other landscape factors on gene flow. Levels of genetic diversity were similar for the Western Slope and Front Range, but effective population sizes were smaller, genetic distances were higher, and there was more admixture in the more urbanized Front Range. Forest cover was strongly positively associated with puma gene flow on the Western Slope, while impervious surfaces restricted gene flow and more open, natural habitats enhanced gene flow on the Front Range. Landscape genomic analyses revealed differences in puma movement and gene flow patterns in rural versus urban settings. Our results highlight the utility of dense, genome-scale markers to document subtle impacts of urbanization on a wide-ranging carnivore living near a large urban center.</span> Creator An entity primarily responsible for making the resource Trumbo, Daryl R. Salerno, Patricia E. Logan, Kenneth A. Alldredge, Mathew W. Gagne, Roderick B. Kozakiewicz, Christopher P. Kraberger, Simona Fountain‐Jones, Nicholas M. Craft, Meggan E. Carver, Scott Ernest, Holly B. Crooks, Kevin R. VandeWoude, Sue Funk, W. Chris Language A language of the resource English Is Part Of A related resource in which the described resource is physically or logically included. Molecular Ecology Format The file format, physical medium, or dimensions of the resource application/pdf Extent The size or duration of the resource. 15 pages Bibliographic Citation A bibliographic reference for the resource. Recommended practice is to include sufficient bibliographic detail to identify the resource as unambiguously as possible. Trumbo, D. R., P. E. Salerno, K. A. Logan, M. W. Alldredge, R. B. Gagne, C. P. Kozakiewicz, S. Kraberger, N. M. Fountain‐Jones, M. E. Craft, S. Carver, H. B. Ernest, K. R. Crooks, S. VandeWoude, and W. C. Funk. 2019. Urbanization impacts apex predator gene flow but not genetic diversity across an urban‐rural divide. Molecular Ecology 28:4926–4940. <a href="https://doi.org/10.1111/mec.15261" target="_blank" rel="noreferrer noopener">https://doi.org/10.1111/mec.15261</a> Type The nature or genre of the resource Article