Landscape effects on the contemporary genetic structure of Ruffed Grouse (Bonasa umbellus) populations

Abstract The amount of dispersal that occurs among populations can be limited by landscape heterogeneity, which is often due to both natural processes and anthropogenic activity leading to habitat loss or fragmentation. Understanding how populations are structured and mapping existing dispersal corridors among populations is imperative to both determining contemporary forces mediating population connectivity, and informing proper management of species with fragmented populations. Furthermore, the contemporary processes mediating gene flow across heterogeneous landscapes on a large scale are understudied, particularly with respect to widespread species. This study focuses on a widespread game bird, the Ruffed Grouse (Bonasa umbellus), for which we analyzed samples from the western extent of the range. Using three types of genetic markers, we uncovered multiple factors acting in concert that are responsible for mediating contemporary population connectivity in this species. Multiple genetically distinct groups were detected; microsatellite markers revealed six groups, and a mitochondrial marker revealed four. Many populations of Ruffed Grouse are genetically isolated, likely by macrogeographic barriers. Furthermore, the addition of landscape genetic methods not only corroborated genetic structure results, but also uncovered compelling evidence that dispersal resistance created by areas of unsuitable habitat is the most important factor mediating population connectivity among the sampled populations. This research has important implications for both our study species and other inhabitants of the early successional forest habitat preferred by Ruffed Grouse. Moreover, it adds to a growing body of evidence that isolation by resistance is more prevalent in shaping population structure of widespread species than previously thought.


| INTRODUC TI ON
Dispersal of organisms across the landscape ultimately determines gene flow among populations, and therefore population connectivity (Slatkin, 1985). Gene flow is essential in species perseverance as it maintains the genetic diversity necessary for populations to respond to changing ecological conditions (Frankham, 2005;Reed & Frankham, 2003). Using genetics, one can detect extrinsic factors restricting dispersal, such as mountain ranges (Funk et al., 2005;Worley et al., 2004), bodies of water (Díaz-Muñoz, 2012), or anthropogenic disturbance (Cegelski, Waits, & Anderson, 2003;Epps et al., 2005). For many species, the structure of the landscape is an important factor shaping contemporary distributions. Unsuitable habitat is a potential barrier to gene flow, but it is not necessarily an impermeable barrier. Habitat often varies in its degree of suitability (Cushman, McKelvey, Hayden, & Schwartz, 2006), resulting in a complex matrix of habitat types with varying dispersal costs or resistance to individuals moving across the landscape. With current landscape genetic methods, it is possible to identify areas of the landscape that are impeding or facilitating connectivity, and also identify the environmental factors that underlie patterns of contemporary gene flow (Keyghobadi, Roland, & Strobek, 1999;Manel, Schwartz, Luikart, & Taberlet, 2003;Storfer et al., 2007).
Differences in landscape resistance and physical distance can both dictate patterns of gene flow (Ruiz-Gonzalez, Cushman, Madeira, Randi, & Gómez-Moliner, 2015). When landscape heterogeneity exists between populations, suitable dispersal routes become more complex. For example, if one of two possible dispersal routes requires movement through habitat that is unsuitable for the study species, then it is likely to present more resistance to dispersing individuals than a route through suitable habitat even when geographic distances are the same. For this reason, landscape heterogeneity within a species' range means patterns of isolation by resistance (IBR) are more likely to occur (Fontaine et al., 2007;McRae & Beier, 2007;Ruiz-Gonzalez et al., 2015). Physical distance between populations can also act as a barrier by creating clinal genetic variation (Ruiz-Gonzalez et al., 2015), or isolation by distance (IBD), as species dispersal is as a function of geographic distance.
IBD and IBR are not, however, mutually exclusive and sometimes a combination of the two best explains genetic structuring (Metzger, Espindola, Waits, & Sullivan, 2015;Piertney, MacColl, Bacon, & Dallas, 1998). Species that are widespread and relatively continuously distributed are expected to exhibit either panmixia or clinal patterns of genetic structure explained by IBD, particularly when comparing populations at a large scale (Alcaide et al., 2009;Purdue, Smith, & Patton, 2000;Ralston & Kirchman, 2012). A few studies have emerged where widespread, continuously distributed species exhibit unexpected patterns of IBR (Pease et al., 2009;Pilot et al., 2006), but the extent to which species with broad geographic ranges exhibit IBD or IBR is unclear (Basto et al., 2016;Frankham, Ballou, & Briscoe, 2010). Understanding the roles that distance and resistance play in structuring populations is dependent on studying both species with limited distributions and broadly ranging species that are not experiencing obvious breaks in population connectivity.
The Ruffed Grouse (Bonasa umbellus; Figure 1) is widely distributed across North America (Figure 2), is resident throughout its broad distribution, and has relatively low dispersal distances for an avian species (approx. 2-4 km;Yoder, 2004). Furthermore, they inhabit early successional forest, and are closely tied to Quaking Aspen (Populus tremuloides), which is an integral part of their diet (Rusch, Destefano, Reynolds, & Lauten, 2000;Svoboda & Gullion, 1972;Zimmerman & Gutiérrez, 2008). Thus, the presence of suitable mixed forest habitat is important for survival, and required for successful dispersal events.
Although Ruffed Grouse have been well studied with respect to ecology and population dynamics (Atwater & Schnell, 1989;Gullion, 1984;Rusch et al., 2000;Zimmerman & Gutierrez, 2007), there is no published information on their population genetics. The Ruffed Grouse is one of the most extensively managed game birds due to heavy hunting pressure throughout most of its range (Rusch et al., 2000). It is also considered an indicator species in the management of early successional forest habitats (USDA Forest Service, 2006). Therefore, information on how macrogeographic barriers and habitat limit gene flow has important implications for managing not only Ruffed Grouse, but also other early successional forest species (e.g.,

philadelphia], Golden-winged Warbler [Vermivora chrysoptera]).
The aims of this study were to quantify population structure by assessing gene flow across a large section of the western extent F I G U R E 1 A Ruffed Grouse in the mixed aspen forest of Crowsnest Pass, Alberta, Canada. Copyright: Ashley Jensen (2016) of the species range and to identify geographic barriers and other landscape features that may be restricting or facilitating gene flow.
We chose to focus on the western extent of the range because this is where macrogeographic barriers are most likely to be influencing population structure, based on studies in a range of other species (Adams, & Burg, 2015b;Pulgarín-R & Burg, 2012;Vonhof et al., 2015). Although widespread species with a continuous distribution are expected to show evidence of IBD, we predicted that Ruffed Grouse populations would exhibit patterns of IBR due to the heterogeneous distribution of suitable habitat throughout their range, their low dispersal ability, and their preference for early successional forest habitat. In addition, we predicted that Ruffed Grouse populations would show significant population genetic structuring, and of the extrinsic factors that may be affecting gene flow, both mountains and swaths of unsuitable habitat would be the most likely geographic features to act as barriers.

| Samples
Fieldwork was conducted from mid-April through May 2016, during the peak activity period of the male Ruffed Grouse's drumming display (Rusch et al., 2000). Birds were located aurally by drumming activity, and the location of each male's drumming log was marked with a handheld GPS unit. We caught birds with a lift net (Fischer, 1974), a carbon dioxide-powered net gun, or mirror traps (Gullion, 1965), which were placed on males' drumming logs adjacent to the drumming stage. Brachial venipuncture was used to collect a blood sample, which was stored in 99% ethanol. For this study, we collected 26 Ruffed Grouse samples at two Alberta locations: Buck Lake (52.91 N,115.01 W),and Crowsnest Pass (49.35 N,114.40 W).
An additional 49 samples were taken from birds collected at the same two sites between 2010 and 2015 for unrelated anatomical studies (Corfield, Harada, & Iwaniuk, 2013;Corfield, Krilow, Vande Ligt, & Iwaniuk, 2013;Krilow & Iwaniuk, 2015). All procedures ad- F I G U R E 2 Map showing the current range of Ruffed Grouse (Bonasa umbellus), and sampling sites for this study. Sampling sites with the same label were pooled for analyses due to close proximity or lack of sufficient sampling at one or more of these sites. Site abbreviations available in Table 1. The data for the range distribution were taken from Birds of North America Online and were projected and overlaid onto a digital elevation map of North America in ArcGIS® v10.2. Digital elevation map courtesy of ESRI®

| DNA extraction and amplification
Genomic DNA was isolated from each blood sample with a modified chelex extraction method (Walsh, Metzger, & Higuchi, 1991).
Samples were screened at intron 6 of Aldolase B on the Z-chromosome (Cheviron & Brumfield, 2009), a 394 bp portion of the mitochondrial control region (domain I and II), and an intron of the nuclear SLC45a2 gene (involved in the melanin production pathway; Gunnarsson et al., 2007). We chose these markers because they are presumed to be neutrally evolving, and this diversity of marker types represents different modes of inheritance and different rates of mutation. Samples were amplified with polymerase chain reaction (PCR) in a 25 µl reaction containing Green GoTaq® Flexi buffer (Promega), 1.5 mM MgCl 2 , 0.08 mM dNTP, 0.4 µM of each primer and 0.5 U GoTaq® Flexi DNA polymerase (Promega) for Aldolase B.
These reaction mixes were the same for control region and SLC45a2, except for MgCl 2 , for which 2.0 mM was used. Amplification of Aldolase B consisted of one cycle at 95ºC for 120 s, 62°C for 45 s, and 72°C for 60 s; 37 cycles of 94°C for 30 s, 62°C for 45 s, and 72°C for 60 s; followed by a final cycle of 72°C for 300 s. For the control region and SLC45a2 primers, the annealing temperature was decreased to 54°C. Successfully amplified samples (Aldolase B = 28, and control region = 57, SLC45a2 = 80) were sequenced at Genome Quebec (Montréal, QC, Canada).
The Aldolase B sequences contained a 7 bp indel. The frequency of the indel differed among populations, so a set of three primers was designed to screen for this indel. The forward primer was placed upstream from the indel, while the other two primers were designed to bind to the insertion and deletion regions, respectively. An M13 tag was added to the 5′ end of the reverse primer for the insertion to increase the size difference between the fragments. Resulting PCR products were 161 bp for the insertion and 118 bp for the deletion.
All 351 samples were screened on a 3% agarose gel.

| Genetic diversity analyses
All sequences were checked, manually aligned and assessed for variation using MEGAv6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013 locus had a significant probability of null alleles for a few sampling sites, so analyses were performed with and without this marker to determine if the potential presence of null alleles was biasing the data. SGCA5 was retained in the final analyses because its exclusion did not cause noticeable variation in the results. SGCA5 had more missing data than the other markers (>25% for some sampling sites) and had to be excluded from F-statistic calculations. Of the 351 genotyped samples, 324 were used for analyses after removing samples that amplified at fewer than six loci. Samples collected in the same area on the same day (i.e., hunter-donated or museum collection samples harvested on the same day) were checked for shared ancestry that would indicate multiple individuals from the same family group; none were found. For analyses that required a priori population assignments, sampling sites within 100 km from each other were grouped as a single "population". All sampling sites in Washington and all sites in Minnesota were grouped together, respectively, due to low sample sizes at some sites within each state (n ≤ 5).

Individuals were sexed prior to compiling final genotypes for
Aldolase B to determine if each individual was hemizygous (females) or homozygous (males). The allele frequencies were then tested for significant pairwise population differentiation using Fisher's exact test (Fisher, 1922).
The data were analyzed with STRUCTURE v2.3.4 (Pritchard, Stephens, & Donnelly, 2000) using correlated allele frequencies in the admixture model, and sampling locations as locpriors. The locpriors option allows sampling location information to be input into the model, but will not create population structuring where there is none. Ten independent runs were performed with 50,000 MCMC repetitions and a 10,000 burn-in period for K values varying from 1 to 10. After these initial runs, values from each K for both LnPr(X|K) and delta K (ΔK; Evanno, Regnaut, & Goudet, 2005) were averaged in STRUCTURE HARVESTER v0.6.94 (Earl & vonHoldt, 2012) to determine the most likely value of K or number of genetic clusters. For the optimal K, any clusters that included more than one population were run through the program independently using the same settings to test for additional substructure. TESS v2.3 (Chen, Durand, Forbes, & François, 2007) was implemented for K values from 2 to 10 using 100,000 sweeps and 50,000 burn-in, and Ψ (value determining how much geographic coordinates influence clustering) was set to 0.6. K was selected based on the runs with the highest posterior probability and highest deviance information criterion (DIC). As with STRUCTURE, once K was determined, any clusters including more than one population were analyzed independently to test for additional substructure. GENELAND v4.0.6 (Guillot, Mortier, & Estoup, 2005) was used to evaluate the optimal value of K using a correlated alleles model, 500,000 iterations, thinning of 200, a burn-in of 500, and uncertainty of spatial coordinates set to 10 km. Default settings were used for the maximum rate of the Poisson process, and the maximum number of nuclei in the Poisson-Voronoi tessellation. The most likely value of K was determined by examining the posterior probabilities averaged over multiple runs (10 runs, K = 1-10), and choosing the K value with the highest average posterior probability. Ten additional runs were conducted at this fixed K value.

| Principal coordinates analysis
To examine genetic structure from a multivariate perspective, we ran a principal coordinates analysis (PCoA) using GenAlEx v6.5 (Peakall & Smouse, 2012). Because it does not make any assumptions about the input data (e.g., Hardy-Weinberg Equilibrium), PCoA is well suited for genetic data (Jombart, Pontier, & Dufour, 2009). Furthermore, patterns revealed by multivariate analyses of genetic data are increasingly being used to further validate Bayesian clustering patterns (Basto et al., 2016). The PCoA analysis was conducted on the matrix of F ′ ST values for the microsatellite data, and the three axes containing the most variation were retained. A three-dimensional plot was made in R using the 3D Scatter Plot package (R Core Team, 2016) to visualize the first three principal coordinates.

| Species distribution modeling
We constructed a species distribution model (SDM) using 53,145 Ruffed Grouse occurrences from the Global Biodiversity Information Facility (GBIF; http://data.gbif.org/). Observations from nonscientific institutions that were not reviewed or moderated were removed, and we further excluded any occurrences prior to 1980 to ensure accuracy of georeferencing. Environmental data were obtained from the WORLDCLIM dataset (v1.4, http://www.worldclim. org/) at a resolution of 2.5 min, and we used the 19 variables in the BIOCLIM layers (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) for the most recent time period . The MODIS-based Global Land Cover Climatology layer (Broxton et al., 2014)  Data were prepared for ecological niche modeling using the SDMtoolbox v1.1c (Brown, 2014) for ArcGIS. All duplicate occurrence records were removed and rarified at a resolution of 30 km 2 to reduce sampling bias toward human settlements and roads; 2,421 occurrences were retained. Due to the similarity of some of the climatic variables used in the layers, we tested for layer autocorrelation at the spatial scale of the North American continent. For pairs of layers that were highly correlated (R > 0.70), one of each pair was removed from the model until no correlated pairs remained, so as not to bias the SDM. The remaining 10 BIOCLIM layers (1, 2, 3, 4, 8, 9, 12, 14, 15, and 18) were used along with the MODIS land cover layer, and rarefied occurrence data to create the SDM. The Gaussian kernel density tool in SDMtoolbox was used to create a bias layer that was added to the model to aid in further accounting for anthropomorphic bias (Phillips et al., 2009).

| Dispersal route analyses
To evaluate whether the intervening landscape matrix leads to population differentiation by influencing dispersal routes and dispersal costs (i.e., IBR), we conducted least cost path (LCP) and least cost corridor (LCC) analyses using SDMtoolbox v1.1c (Brown, 2014) in ArcGIS v10.2 (ESRI®). The SDM was inverted to create a friction layer, and coordinates for each sampling site were entered in decimal degrees. LCPs and LCCs were calculated between each pair of populations using the friction values. LCPs are calculated as the dispersal paths of least resistance between each pair of sampling sites based on the resistance surface. To calculate LCCs, the LCPs were weighted by resistance values based on the friction layer and then categorized using a "percentage of LCP" method with cutoffs for inclusion into high-, mid-, and low-classes set at 5%, 2%, and 1% of the LCP value, respectively. The weighted and categorized LCPs were then summed to create a LCC dispersal network.

| Isolation by distance and resistance
We tested for IBD with a Mantel test in GenAlEx v6.5 (Peakall & Smouse, 2012) using pairwise F ′ ST (Meirmans & Hedrick, 2011) and geographic distance. We calculated geographic distance as a straight line between populations, where lines were kept within the boundaries of the species' geographic distribution. This analysis was also conducted using least cost path distance instead of Euclidean distance. To explicitly test for IBR, we used a similar analysis to IBD.

ST
) and resistance values were assessed for correlations using a paired Mantel test implemented in GenAlEx v6.5 (Peakall & Smouse, 2012). The matrix of resistance values was created from the LCC map by weighting the distance of each LCC by the resistance values along the corridor. After calculating the Mantel tests for all population pairs, the same tests were performed on subsets of populations to examine patterns at a regional scale (if the number of sampling sites permitted).
Because Mantel tests can be prone to Type I error (Legendre & Fortin, 2010;Legendre, Fortin, & Borcard, 2015), the matrices described above were also analyzed using distance-based Moran's eigenvector map analysis in the MEMGENE package (Galpern, Peres-Neto, Polfus, & Manseau, 2014) for R. This method finds Moran's eigenvectors in the spatial data (input as geographic, least cost, or resistance distance) using principal coordinates analysis. Then redundancy analysis is used to select a reduced set of vectors based on their contribution as predictors of the response variable (genetic distance  (Fisher, 1922)  to reveal complex patterns (Griffith & Peres-Neto, 2006;Legendre & Fortin, 2010;Richardson, Brady, Wang, & Spear, 2016). However, due to the constraints on matrix size that can be analyzed in MEMGENE (P. Galpern, pers. comm.), only the full dataset including all 15 populations could be analyzed using both this package and Mantel tests.

| Genetic diversity
The    (Table 2). Most populations had one to three private alleles, but BL and BV both had five. Also notable was the high frequency (0.14) of a single private allele for the PR population.

| Genetic structure
Pairwise Ф ST values for the CR locus ranged from −0.057 for the CP:EA comparison to 0.660 for the AK:WA comparison (Table 3).  Table 4 values. Samples from Alberta (

| Bayesian clustering analyses
Plots of delta K (ΔK) and mean log likelihood (LnPr(X|K); Supporting Information Figure S1) from the initial STRUCTURE analyses showed K = 5. The five clusters were as follows: Alaska + Yukon, Washington, the Great Lakes, Crowsnest Pass, and all remaining sites in Alberta ( Figure 6a). To investigate additional population structure, we analyzed each cluster independently. The AK-YT cluster showed evidence of additional structure (K = 2; Figure 6b). Splitting the AK-YT cluster in two creates a total of six clusters from the STRUCTURE analyses which is concordant with the pairwise F ′ ST values (Table 4). Both DIC and log likelihood plots of the Bayesian clustering analysis performed in TESS showed the most likely number of clusters as K = 4 with potential substructure (Supporting Information Figure S2). The DIC plot was bimodal with a second peak at K = 7 (Supporting Information Figure S2a); however, when examined, the Q plots for K = 7 showed clear oversplitting of clusters. We therefore concluded that once hierarchical analysis was performed to reveal substructure, the true number of clusters was K = 6 (Supporting Information Figure S3), which is concordant with both STRUCTURE and F ′ ST results. The GENELAND analysis showed evidence of K = 7 clusters at the highest frequency over the MCMC chain, which was in agreement with the highest value for the averaged posterior probabilities of the initial set of runs. Five of the seven groupings identified by GENELAND corroborated the clusters inferred by STRUCTURE and TESS: AK, YT, WA, CP and Great Lakes. In addition to those five, GENELAND split GP from the remaining Alberta sampling sites ( Figure 6c; Supporting Information Figure S4).

| Principal coordinates analysis
The PCoA using F ′ ST values showed distinct genetic groupings, with the first and second axes accounting for 35.6% and 19.7% of the variation, respectively, and the third axis explaining 16.0% F I G U R E 5 The species divergence map made using the Landscape Genetics toolbox (Vandergast, Perry, Lugo, & Hathaway, 2011) in ArcGIS®. Pairwise cluster together as they do in all other analyses, and the GP and PR sites clustered together. Although GP and PR were separated from other Alberta sampling sites, they were in much closer proximity to these Alberta sites than to the other sampling sites. The groupings of the PCoA confirm groupings identified by TESS and STRUCTURE and provide evidence of GP population divergence from the rest of Alberta as indicated by the GENELAND results ( Figure 6b). F I G U R E 6 Ruffed Grouse population structure as inferred by hierarchical runs in STRUCTURE v2.3.4 (Pritchard et al., 2000) with eight microsatellite loci. At (a) K = 5, and further substructure was apparent when samples from Alaska and Yukon were analyzed independently, which resulted in (b) K = 2. No further substructure was found when the remaining clusters were analyzed independently. A (c) map of multiple Bayesian clustering programs where GENELAND clusters have been color-coded, mapped in geographic space and clipped to limits of the species' range, with the gradient of colors representing clines created by the contour lines of the posterior probability maps in GENELAND (Supporting Information Figure S3). Circles represent sampling sites, and circle color corresponds to the STRUCTURE and TESS consensus cluster assignments. There was only one instance of discordance between the programs: additional cluster (GP) assigned by GENELAND

| Isolation by distance
The results of the Mantel test for IBD showed a significant correlation between Euclidean distance and pairwise F ′ ST values when all sampling sites were compared (R 2 = 0.378; p = 0.01; Table 5; Figure 9). MEMGENE results were similar to those of the Mantel test, although a lower proportion of genetic distance was explained by Euclidean distance (R adj 2 = 0.211; p = 0.02; Table 5). When subsets of data were tested for IBD using Mantel tests, only the analysis containing populations east of the Rockies (Alberta and Great Lakes clusters) provided evidence that geographic distance is a significant predictor of pairwise genetic distance (R 2 = 0.567; p = 0.02; Table 5

| Isolation by resistance
When a Mantel test was performed to test the correlation between LCP distance and genetic distance (F ′ ST ), the correlation was F I G U R E 7 A three-dimensional plot of the first three axes of the PCoA as calculated in GenAlEx v6.5 (Peakall & Smouse, 2012). Populations are labeled, and principal components are labeled on their respective axes including the amount of variation captured by each in R (R Core Team 2016) considerably higher (R 2 = 0.649, p = 0.01; Table 5) than the value calculated for IBD (R 2 = 0.370, p = 0.01). The correlation between resistance distance and genetic distance was higher yet (R 2 = 0.674, p = 0.001). MEMGENE analysis using least cost distances provided an adjusted R 2 value (R adj 2 = 0.189; p = 0.01; Table 5) similar to that of the IBD analysis. However, using resistance distances, both the MEMEGENE and Mantel tests were similar and accounted for considerably more of the variation in the genetic data than the other two distance measures (R adj 2 = 0.487; p = 0.005; Table 5). When Mantel tests were performed on subsets of populations, the IBR model explained more of the variation in genetic data than either the IBD or LCP models for all subsets (Table 5). The only group that F I G U R E 8 (a) The Species distribution model (SDM) created with SDM toolbox (Brown, 2014) for ArcGIS® and MaxEnt (Phillips et al., 2006). Occurrences from Global Biodiversity Information Facility (GBIF; http://data.gbif.org/) and environmental layers (climate and vegetation data) were input into the model. The SDM shows areas where the environmental conditions are suitable for the Ruffed Grouse to occur (i.e., ecological niche). The scale depicted is cumulative and represents the percent likelihood of habitat suitability for Ruffed Grouse based on the model variables. Using a resistance layer created from inverse of the SDM, (b) least cost corridors (LCC) were calculated among the 15 sampled populations of Ruffed Grouse. The LCC provides more information than least cost paths (Supporting Information Figure  S4) and shows the most likely dispersal routes among populations as corridors instead of paths. It also provides dispersal costs along these corridors coded by color; red representing areas where there is low resistance (i.e., low dispersal cost), and blue representing areas of high resistance did not yield a significant correlation were the Alberta populations (excluding CP). No other groups had enough sampling sites to perform within region comparisons.

| D ISCUSS I ON
Using multi-locus genetic data and environmental variables, we found significant genetic differentiation and limited connectivity among western populations of Ruffed Grouse. Macrogeographic barriers, tracts of unsuitable habitat, and the species' preference for aspen-dominated mixed forest are likely playing important roles in creating genetically structured populations.

| Contemporary population genetic structure and macrogeographic barriers
Data from multiple neutral genetic markers show structuring of Ruffed Grouse populations across their range at multiple spatial scales. Aside from the most highly differentiated populations, AK and WA, at least four other distinct genetic groups exist: Yukon, southwest Alberta (CP), a large group including most of central/ northern Alberta, and one group near the Great Lakes (Figure 6b).
Although these results support our hypothesis of genetic differentiation due to low dispersal ability, the degree of differentiation is somewhat unexpected.
A number of landscape features co-occur with the boundaries of genetic clusters for Ruffed Grouse across western North America.
The Columbia River basin (southeast of WA sites) and northern ex- Notes. To examine multiple spatial scales, an overall correlation was run for all sampling sites (both Mantel tests and dbMEM), as well as subsets of the different sampling sites (Mantel tests only). The correlation value for each comparison is reported (R 2 or adjusted R 2 ), along with the significance level of each test.
TA B L E 5 Results of Mantel tests and distance-based Moran's eigenvector map analysis (dbMEM) for three different models; isolation by distance (IBD), isolation by distance using least cost path distance (LCP), and isolation by resistance (IBR) The Washington population is genetically distinct, as supported by the microsatellite data (Table 4), Z-linked marker data (Table 4) (Broxton et al., 2014;Pater, Bryce, Thorson, Kagan, & Chapell, 1998), and generally lower elevation passes than the northern Rockies (Franklin & Dyrness, 1973). Similar patterns of differentiation occur in a widespread generalist passerine, the Black-capped F I G U R E 9 Plots showing Mantel tests of (a) isolation by distance (IBD) comparing genetic distance as measured by F ′ ST and geographic distance between populations, (b) isolation by least cost distance using least cost paths (LCP), and (c) isolation by resistance (IBR) using resistance-weighted corridor distances from the LCC map. Correlation values (R 2 ) and p-values of each Mantel test are given for each respective plot as calculated in GenAlEx v6.5 (Peakall & Smouse, 2012). Plots shown are comparisons including all 15 sampled populations Chickadee (Poecile atricapillus), where the habitat composition of a mountain range corresponds to restricted gene flow (Adams, & Burg, 2015b). The Black-capped Chickadee also has similar patterns of isolation in Alaska, northwest British Columbia, and on either side of the southern Rockies (Adams, & Burg, 2015b;Hindley, Graham, & Burg, 2018).
The increased number of sampling sites in Alberta allowed us to assess genetic structure on a finer scale. While most of the Alberta populations are not differentiated from one another, evidence from the nuclear loci show that the Crowsnest Pass population is isolated from all other populations (Figures 4, 5 and 6b). Because CP does not show significant differentiation at the control region ( Figure 3), it is unlikely that divergence of this population reflects historical isolation. Instead, the differentiation of the CP population likely arose due to contemporary barriers to gene flow. In some species, southwest Alberta populations are divergent from individuals sampled throughout the rest of Alberta, and instead group with either British Columbia populations (Adams, & Burg, 2015a) or Intermountain West (i.e., Montana, Idaho, Wyoming) populations (Dohms, Graham, & Burg, 2017;Pulgarín-R & Burg, 2012) implying that the geography of the Rocky Mountains may affect the genetic structure of species differently depending on their life history.
The presence of unsuitable habitat may also be restricting gene flow, particularly for the CP population. Only a narrow swath of suitable Ruffed Grouse habitat presently connects southwestern Alberta and the rest of the province; most of the southeast part of the province is open grassland, which this species is reluctant to disperse through (Yoder, 2004), and the Rocky Mountains run along the western edge of the province. The habitat in the Rocky Mountains consists mainly of contiguous coniferous forest, with suitable mixed forest habitat occurring mostly on low elevation slopes and valleys (Broxton et al., 2014;Natural Regions Committee, 2006). Although Ruffed Grouse are more likely to disperse through coniferous forests than grasslands, their short dispersal distances (approx. 2-4 km;Yoder, 2004) suggest that dispersal through vast expanses of coniferous forest are likely to be infrequent. Because the CP population is in close spatial proximity to some of the other populations sampled in Alberta, geographic distance is unlikely to be a causal factor for population differentiation and this is corroborated by IBD analysis (Table 5) While there is not complete consensus across our analyses for the GP cluster, there is certainly evidence of differentiation of this population, which could be due to its proximity to the Peace River valley. The Peace River is the only river to cut a continuous valley through the entire width of the Rocky Mountain range (Feinstein, 2010). It is possible that genotype frequencies at GP are subject to an influx of genes from British Columbia through the Peace River valley corridor, a low elevation valley through the Rocky Mountains. This is supported by the STRUCTURE results, in which PR and GP show some admixture with the AK/YT cluster ( Figure 6a). The Peace River valley may be a contact zone for Ruffed Grouse populations on either side of the Rocky Mountains.
There is evidence that this important corridor facilitates connectivity for multiple species that are otherwise geographically isolated by the Rockies, particularly those reliant on mixed forest or shrubby habitats, and those that would have difficulty dispersing through coniferous forests (Irwin, Brelsford, Toews, MacDonald, & Phinney, 2009;Seneviratne, Davidson, Martin, & Irwin, 2016;Toews, Brelsford, & Irwin, 2011). Furthermore, the permeability of a mountain barrier may range from porous (Vonhof et al., 2015) to impermeable (Irwin, Irwin, & Smith, 2011)

| Landscape genetics: isolation by distance or resistance?
Aside from mountain ranges, the presence of unsuitable habitat is the most prevalent potential barrier between our sampling sites.
Due to the seemingly high degree of habitat heterogeneity across the landscape, we incorporated environmental variables into our analyses to test their effects on the genetic structuring present, and to help further test the presence of putative geographic barriers discussed previously. LCC analysis revealed that dispersal cost varies across the landscape and confirms our earlier hypothesis that mountain ranges are likely to impede dispersal among populations.
Mountain ranges in Alaska, as well as the Rockies have markedly high levels of resistance to Ruffed Grouse dispersal, with the exception of two corridors through the mountains: one through the Peace River valley, and one through the Intermountain West, connecting south-central Alberta and northeast Washington (Figure 8b). These corridors are lower elevation areas with slightly milder climate and a higher proportion of mixed forest than the surrounding mountain slopes (Broxton et al., 2014;Hijmans et al., 2005). No direct dispersal routes could be identified among Yukon, Alaska, and Washington populations, and the only dispersal route connecting Alaska to the other sampled populations has high resistance. It should, however, be noted that the lack of sampling from within British Columbia may have prevented the identification of dispersal routes among these population. That said, the LCCs occurred mainly in areas with tracts of mesic, mixed forest, which implies that variation in climate and forest type across the landscape may be important in creating population structuring.
Across all populations, IBR explained significantly more of the genetic differentiation than IBD at most spatial scales (Table 5); the only exception being the comparisons within Alberta. The LCC between Alberta and the Great Lakes had low resistance, stretching across the parkland/boreal forest in a direct path (Figure 8b). Because much of the landscape between these two regions presents low dispersal resistance, it is not unusual that IBR only moderately outperformed IBD. This implies that the genetic differences between Alberta and Great Lakes populations are explained by a combination of physical distance and dispersal cost through intervening habitat at this large spatial scale, although this should be verified further by sampling individuals in between these two populations (e.g., Saskatchewan, Manitoba). In contrast, the genetic differentiation among western populations (AK, YT, WA) was better explained by IBR than by IBD alone (Table 5). This is concordant with the LCP/LCC maps, which do not show any direct dispersal routes between the AK, YT, and WA populations. Furthermore, the routes that were detected have moderate to high resistance. The patterns lend support to the idea that the higher heterogeneity of habitat types west of the Rocky Mountains are restricting dispersal, but sampling of BC populations is required to rule out the effects of small sample sizes.
Due to large gaps among our sampling sites outside of Alberta, it is difficult to be certain whether genetic boundaries between populations are gradual genetic clines or genetic breaks. IBD manifests as smooth, clinal gradients between genetic clusters (Mims et al., 2016), whereas distinct boundaries among genetic clusters are more likely to occur in populations mediated by IBR (Coulon et al., 2006). The high pairwise F ′ ST values and relatively steep genetic cluster boundaries (Table 4; Figure 5; Supporting Information Figure   S3) provide further evidence that Ruffed Grouse populations are distinct genetic clusters mediated by IBR. Furthermore, the patterns of IBR among Alberta populations and those west of the Rockies point to the Rocky Mountain range as a barrier. The conifer-dominated habitat characterizing high-elevation mountains has high dispersal costs and would explain patterns of reduced connectivity in Ruffed Grouse found in the SDM and LCC (Figure 8a,b).
Although IBR performed consistently better than IBD, both of these models only explained a small to moderate proportion of the genetic differentiation between CP and adjacent populations (Table 5). CP might therefore have additional factors affecting connectivity with other Ruffed Grouse populations. The environmental variables used in the LCC analysis explained much of the differentiation present among the populations sampled for this study (Table 5), but it is possible that additional, unsampled environmental factors are contributing to the genetic break at the CP population. In some species, genetic structure may not be evident at a broad spatial scale and is only detected when populations are assessed at a finer scale (Adams, & Burg, 2015a). Sampling in other areas such as British Columbia, Saskatchewan, and Manitoba could also provide additional insight into Ruffed Grouse population structure. Therefore, a smaller scale analysis with more environmental variables and higher resolution sampling may aid in identifying the factors underlying the genetic differentiation of the CP population. Furthermore, differences in the vegetation communities of the Intermountain West (including southwestern Alberta) compared to the boreal region (including central-northern Alberta) may be important. The shift in Populus species (e.g., P. trichocarpa, P. angustifolia, P. balsamifera) between these regions is an example of how the environmental conditions in these areas differ and may be contributing to genetic structure through local adaptation. These shifts in vegetation communities correspond to patterns of genetic variation in other avian species in the Intermountain West and central Alberta (Adams, & Burg, 2015a;Dohms et al., 2017).

| CON CLUS IONS
Our study is a first look into the population genetics of Ruffed Grouse, and we found that contemporary populations in the western extent of the range were highly genetically structured, with the strongest genetic breaks co-occurring with high-elevation mountain ranges and habitat. Resistance modeling revealed genetic structure in Ruffed Grouse is primarily influenced by the heterogeneous habitat mosaic of the western North American landscape. Dispersal routes appear to be restricted to areas with mixed forest habitat supporting our hypothesis of IBR in Ruffed Grouse, despite its widespread and relatively continuous distribution. In doing so, this study represents one of relatively few contemporary population genetic studies focusing on broadly distributed organisms, and one of even fewer revealing patterns of IBR in widespread, fairly continuously distributed organisms that may be expected to exhibit IBD (Ball, Finnegan, Manseau, & Wilson, 2010;Pease et al., 2009;Pilot et al., 2006). As such, we should no longer assume that dispersal distance and distribution are the only factors driving dispersal patterns, and that IBR is only likely to affect species at small scales. This adds to the growing body of work highlighting the importance of evaluating the role of environmental variables in population genetic structuring (McRae & Beier, 2007;Vergara et al., 2015). It also underscores the need for more landscape genetic studies focusing on broadly distributed taxa because they may be experiencing genetic isolation regardless of their relatively ubiquitous distributions.

ACK N OWLED G M ENTS
We wish to thank all of the hunters who assisted us in obtaining samples, and the Alberta Conservation Association for their assistance in reaching out to hunters across Alberta. We also wish to thank the Royal Alberta Museum, Burke Museum, Field Museum of Natural History, University of Alaska Museum, Emily Herdman with Alberta Environment and Parks, and Mike Suitor with Environment Yukon for providing us with additional samples.
We also thank Barb and George Iwaniuk, Karen Dean, and Bob and Joyce Mackenzie-Grieve for logistical support during fieldwork.

Funding was provided by an Alberta Conservation Association
Biodiversity Grant to AMJ, Natural Science and Engineering Council of Canada (NSERC) Discovery grants to ANI and TMB and an NSERC Discovery Accelerator Supplement and Canada Foundation for Innovation grant to ANI.

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

AUTH O R CO NTR I B UTI O N S
AMJ, TMB, and ANI designed the study; AMJ, NPO, and ANI conducted fieldwork and obtained samples; AMJ performed all laboratory work and analyses with help from TMB; AMJ, TMB, and ANI interpreted the results; AMJ, TMB, and ANI prepared the manuscript. All authors revised, read, and approved the final manuscript.