Quantifying functional connectivity: The role of breeding habitat, abundance, and landscape features on range‐wide gene flow in sage‐grouse

Abstract Functional connectivity, quantified using landscape genetics, can inform conservation through the identification of factors linking genetic structure to landscape mechanisms. We used breeding habitat metrics, landscape attributes, and indices of grouse abundance, to compare fit between structural connectivity and genetic differentiation within five long‐established Sage‐Grouse Management Zones (MZ) I‐V using microsatellite genotypes from 6,844 greater sage‐grouse (Centrocercus urophasianus) collected across their 10.7 million‐km2 range. We estimated structural connectivity using a circuit theory‐based approach where we built resistance surfaces using thresholds dividing the landscape into “habitat” and “nonhabitat” and nodes were clusters of sage‐grouse leks (where feather samples were collected using noninvasive techniques). As hypothesized, MZ‐specific habitat metrics were the best predictors of differentiation. To our surprise, inclusion of grouse abundance‐corrected indices did not greatly improve model fit in most MZs. Functional connectivity of breeding habitat was reduced when probability of lek occurrence dropped below 0.25 (MZs I, IV) and 0.5 (II), thresholds lower than those previously identified as required for the formation of breeding leks, which suggests that individuals are willing to travel through undesirable habitat. The individual MZ landscape results suggested terrain roughness and steepness shaped functional connectivity across all MZs. Across respective MZs, sagebrush availability (<10%–30%; II, IV, V), tree canopy cover (>10%; I, II, IV), and cultivation (>25%; I, II, IV, V) each reduced movement beyond their respective thresholds. Model validations confirmed variation in predictive ability across MZs with top resistance surfaces better predicting gene flow than geographic distance alone, especially in cases of low and high differentiation among lek groups. The resultant resistance maps we produced spatially depict the strength and redundancy of range‐wide gene flow and can help direct conservation actions to maintain and restore functional connectivity for sage‐grouse.


| INTRODUC TI ON
Functional connectivity describes how landscapes influence the movement of individuals between habitat patches (Tischendorf & Fahrig, 2000). Its conservation is fundamental to the protection of diverse and viable populations able to persist and adapt to changing conditions (Fahrig & Merriam, 1985;Lamy et al., 2013;Morrissey & de Kerckhove, 2009). The distribution of habitat on the landscape (i.e., structural connectivity) often does not correlate directly with functional connectivity, because it does not consider the behavioral responses of individuals to habitat structure, nor their willingness to disperse through undesirable habitats. Landscape genetic approaches compare gene flow (i.e., functional connectivity) to connectivity measures and yield direct insights into the relative importance of landscape features and their configuration (Manel, Schwartz, Luikart, & Taberlet, 2003;Row et al., 2015).
Landscape features influencing genetic structure can range from natural barriers, such as mountains or lakes, to more recent barriers such as roads and development. Thus, insight from landscape genetics can help to guide management and conservation efforts by identifying specific features that reduce or facilitate gene flow and by identifying locations where mitigation of impedance is required (Epps et al., 2005;Roever, van Aarde, & Leggett, 2013). Despite the potential implications, the output from functional connectivity analyses is often omitted from conservation plans because they are not carried out in ways that are conducive to management objectives (Keller, Holderegger, van Strien, & Bolliger, 2015). For example, both the spatial extent of the study area and the resolution of landscape variables can influence the estimated importance of landscape features to functional connectivity (Anderson et al., 2010), which makes extrapolation across scales and comparison among studies challenging. Thus, conducting research at scales too large or too small for planning tools will hinder implementation (Keller et al., 2015). Furthermore, investigating both landscape variables that can be actively managed and those that are likely to impact connectivity, but cannot be actively managed (e.g., elevation), will likely lead to biologically meaningful results that can also be incorporated into conservation plans.
Although functional connectivity is likely influenced by a multitude of landscape features, it can prove challenging to disentangle their relative effects (Row, Knick, Oyler-McCance, Lougheed, & Fedy, 2017). One way of modeling the biological importance of multiple landscape features is through the use of habitat indices.
Habitat indices are typically derived from occurrence or radiotelemetry data and are used to predict the probability of occupancy across a landscape based on the suitability of habitat (Fedy et al., 2012;Phillips, Anderson, & Schapire, 2006) or the probability of selection for habitat types (Gubili et al., 2017;Shafer et al., 2012). Using habitat indices in landscape genetic models can assist with establishing a direct link between model results and management actions, and can alleviate some of the issues related to testing a multitude of landscape variables (i.e., type I error). If structural connectivity quantified using habitat indices proves to be a strong predictor of functional connectivity (Row, Blouin-Demers, & Lougheed, 2010;Row et al., 2015;Wang, Yang, Bridgman, & Lin, 2008), identifying regions of importance to maintaining connectivity and predicting impacts from changes to habitat suitability become more clear (Roever et al., 2013). However, in some cases, the relationship between habitat indices and functional connectivity can be weak (Roffler et al., 2016). This lack of relationship is likely the result of a mismatch between dispersal and the habitats modeled in the indices (Ribe, Morganti, Hulse, & Shull, 1998;Spear, Balkenhol, Fortin, McRae, & Scribner, 2010).
Habitat models used in landscape genetic analysis generally do not include, nor model, variation in local population abundance.
However, local population abundance can influence dispersal (Matthysen, 2005;Pflüger & Balkenhol, 2014;Strevens & Bonsall, 2011) and regional differences in population abundance can impact genetic drift and spatial genetic structure (Row, Wilson, & Murray, 2016). Therefore, the inclusion of effective population sizes in landscape genetic models can improve the fit between genetic differentiation and landscape variables hypothesized to affect genetic differentiation (Weckworth et al., 2013). For example, even when the probability of occupancy is high, low abundance could lead to increased genetic differentiation. Therefore, including abundance estimates in habitat indices should improve the prediction of functional connectivity. However, due to the difficulty in obtaining abundance information at large spatial scales, abundance is rarely incorporated into landscape genetic models.
In this study, we establish the drivers of functional connectivity for the greater sage-grouse (Centrocercus urophasianus; hereafter sage-grouse) across the species' range in North America. Sagegrouse are distributed across eleven US states and two Canadian provinces in western North America. However, there has been a 44% range contraction since European settlement and substantial declines in population size since the 1960s Schroeder et al., 2004). These decreases in range and population size have largely been attributed to habitat loss and fragmentation and have resulted in the species being petitioned for listing under the U.S. Endangered Species Act multiple times. Although populations are declining, there is a wealth of information available on the species' habitat utilization (Doherty, Naugle, & Walker, 2010;Doherty, Naugle, Walker, & Graham, 2008;Fedy et al., 2014) and range-wide estimates of abundance through counts of males at breeding leks K E Y W O R D S dispersal, gene flow, genetic differentiation, genetic diversity, habitat selection models, isolation by resistance, landscape resistance (Doherty, Evans, Coates, Juliusson, & Fedy, 2016; where they congregate to attract and compete for access to females. Combining this habitat utilization and abundance information with our range-wide genetic data makes sage-grouse an ideal candidate for exploration of the link between habitat, abundance, and functional connectivity. The relationship between functional connectivity and landscape or habitat variables is often quantified using regression-based approaches in which pairwise genetic differentiation between groups or individuals is compared to pairwise landscape-based resistance (or landscape cost) distances derived from a variety of single or combined landscape features McRae, 2006;Van Strien, Keller, & Holderegger, 2012). Here, we used this approach to address three main questions. 1) Do habitat indices provide a better fit with genetic data than individual or combined landscape predictors? 2) Does the incorporation of abundance estimates improve the fit between genetic and habitat-based resistance distances? 3) How well do the best fitting resistance models predict genetic differentiation across landscapes when compared to distance alone? Because the abundance and distribution of very low-quality habitat can be disproportionately important for functional connectivity (Row et al., 2010(Row et al., , 2015, identifying these thresholds (Keller et al., 2015;Méndez, Vögeli, Tella, & Godoy, 2014) where functional connectivity is reduced could help focus conservation efforts on improving the habitat quality of landscapes to reach the necessary thresholds. Thus, we also determined whether there were consistent thresholds in habitat or landscape predictors which, when exceeded, resulted in disruption of functional connectivity for sage-grouse.
Overall, our results should provide guidelines for large-scale evaluations of functional connectivity and assist with management goalms of maintaining or restoring functional connectivity.

| Study area and sample collection
Our study area includes nearly the entire range of sage-grouse in North America, with the exception of a few small, geographically disjunct portions of the range which include individuals within the Canadian provinces of Alberta and Saskatchewan and sage-grouse MZs VI and VII (Figure 1, Figure S1). The sage-grouse range largely coincides with the distribution of sagebrush (Artemisia spp.) in western North America, an ecosystem that has become increasingly fragmented due to habitat conversion for agriculture, housing development, oil and gas development, wildfire, exotic grasses and invasive conifer encroachment (Knick et al., 2003). Within the range, we stratified our analysis based on long-standing Sage-Grouse Management Zones (here after MZ; Figure 1), which were established using differences in environmental attributes that influenced vegetation communities in each zone and were uninfluenced F I G U R E 1 Locations (black circles) of 267 breeding lek clusters of greater sage-grouse samples (6,844 samples in total) collected from 2005 to 2014 and used to compare genetic differentiation and landscape resistance. Full distribution of all samples can be found in Figure S1. Shaded gray area represents the current distribution of greater sage-grouse, and dotted lines represent the extent of seven management zones for the species. Samples from management zones I-V were included in this analysis  (Stiver et al., 2006).

MZ
This stratification ensured relevance to previous studies, ongoing management initiatives, and consistency with the development of existing habitat and abundance layers (Doherty et al., 2016).
We used 6,844 spatially referenced feather and tissue samples collected across the sage-grouse range ( Figure 1, tracted from feather and blood using previously described methods (Cross, Naugle, Carlson, & Schwartz, 2016;Row et al., 2015). Both FORT and NGC amplified 15 microsatellite loci in eight multiplex polymerase chain reactions (PCR) and genotypes were determined using electrophoresis (Cross et al., 2016;Row et al., 2015). Full genetic methods can also be found in Appendix S1.
Feather DNA samples can have low-quality and quantity DNA.
Therefore, to ensure correct genotypes from feather samples, each sample was PCR amplified at least twice across the 15 microsatellite loci to screen for allele dropout, stutter artifacts, and false alleles.
Alleles for each locus were coded as missing if they did not match across at least two independent runs. Samples with missing genotypes for more than five loci were removed. Genotypes were then screened to ensure consistency between allele length and length of the microsatellite repeat motif. We used program DROPOUT v2.3 (McKelvey & Schwartz, 2005) and package ALLELEMATCH v2.5 (Galpern, Manseau, Hettinga, Smith, & Wilson, 2012) in R (R Core Team, 2016) to screen for genotyping error and to identify and remove multiple captures of the same individual.
To combine the genotype data sets from both labs, we first genotyped the same 70 individuals. Each laboratory's genotypes for these individuals were compared, and shifts in allele size were implemented to synchronize allele calls for all samples where necessary.
Following the combination of samples, an additional ALLELEMATCH analysis was performed on the complete, combined data. Finally, we quantified the power of our microsatellite locus panel to discern individuals using probability identity (P ID ; Evett & Weir, 1998) which calculates the probability that two individuals drawn at random from the population could have the same genotype across all loci.

| Genetic differentiation within management zones
Sage-grouse have clustered distributions (Doherty et al., 2016), particularly during the breeding season when individuals attend centralized breeding leks where the majority of our samples were collected. Additionally, sage-grouse can have substantial withinand interseasonal movement distances (Cross, Naugle, Carlson, & Schwartz, 2017;Fedy et al., 2012), so we expected little differentiation between spatially proximate leks and individual samples (Row et al., 2015). Thus, we used a clustered (i.e., "group-based") approach to evaluate genetic differentiation as this likely best represents the ecology of the species at the spatial scale of this study and the sampling scheme (i.e., multiple samples from each lek). We derived lek clusters by grouping spatially proximate sampling locations using hierarchical clustering (hclust function with complete method) in R (R Core Team, 2016). The algorithm iteratively joined each sampling location based on distance using the Lance-Williams dissimilarity formula (Legendre & Legendre, 2012). Subsequently, we used a cut distance to differentiate clusters separated by a distance >25 km, which represented the maximum average summer to winter movement distance for any population in Wyoming (Fedy et al., 2012). We used a minimum sample size of eight individuals per genetic cluster for all subsequent analyses.
TA B L E 1 Number of Greater Sage-Grouse samples, breeding leks and grouped lek clusters used to establish patterns of genetic diversity and population structure in each of five long-established management zones

| Landscape resistance within management zones
We quantified pairwise effective resistance among lek clusters using circuit theory. Landscape surfaces were coded such that each pixel was assigned a value representing resistance to gene flow.
Subsequently, we used the gdistance package in R to quantify pairwise resistance between two groups as the expected correlate to the amount of gene flow between two groups (i.e., higher pairwise We generated pairwise effective resistances from landscape resistance surfaces using three published models representing accepted abiotic and biotic variables influential in sage-grouse biology (Doherty et al., 2016; see below for specific predictors). For each resistance surface, we tested three threshold values above or below which the landscape is hypothesized to be resistant to grouse movements (Keller et al., 2015). In all cases, we used the thresholds to derive a binary resistance surface representing habitat and nonhabitat. For variables with a known negative association with sagegrouse (e.g., human disturbance, tillage), raster pixel values below the thresholds were set as habitat and assigned a resistance of 1 with nonhabitat cells receiving a higher value. When the landscape variable represented a positive association with sage-grouse (e.g., percentage sagebrush cover, habitat utilization), values above the threshold were classified as habitat (resistance of 1) and we assigned higher values to nonhabitat (below the threshold). By determining the threshold values that best describe functional connectivity, we can make specific predictions as to when landscape degradation will lead to decreased connectivity. Below, we describe each resistance surface and their associated resistance and threshold values.
All raster surfaces were resampled using bilinear interpolation to TA B L E 2 Landscape variables used as resistance surfaces and associated thresholds used to delineate "habitat" and "nonhabitat" and the top selected resistance values. Additional details on habitat layers can be found in Doherty et al. (2016) Description Abbrev

| Development of habitat predictors
Our first set of resistance surfaces were built on a breeding habitat utilization model (BH) developed from landscape characteristics surrounding active breeding leks. The dominant variables varied for each MZ, but primarily included positive relationships with percent aerial coverage of sagebrush and negative relationships with increasing tree canopy cover and anthropogenic variables such as cultivation and human disturbance (Table 2). Details on the development of the BH models and top variables for each MZ are presented in Doherty et al. (2016). Raster values associated with the BH model represent the predicted probability that the raster pixel contains sufficient breeding habitat to support sage-grouse lek formation.
All active leks within the sage-grouse range occurred in areas with a predicted probability (p) ≥ .65, so we used this as one of three thresholds, classifying all raster cell values with p ≥ .65 as habitat and anything below as nonhabitat and, as such, resistant to dispersal (Table 2). We expected that not all nonlek habitat (e.g., p < .65) would reduce dispersal. Thus, we considered two lower probability thresholds (p ≥ .5 and p ≥ .25) which split the distribution and in which all cells below the thresholds in p were assigned as nonhabitat to represent increased resistance to movement.
We also derived a set of resistance layers based on a kernel index

| Development of landscape predictors
We used eight landscape variables to develop landscape resistance surfaces independent of breeding habitat and abundance (Table 2).
We used a subset of the variables that proved most important in Doherty et al. (2016) and which have previously been shown to influence sage-grouse habitat. Landscape predictors included features that could be managed through protection (e.g., percent sagebrush cover) and restoration (e.g., canopy cover, tilled agriculture) and those that would be impossible to manage but have existed as barriers or which have restricted movement for a longer period of time (e.g., terrain topology, measured by steepness and roughness).
Environmental conditions can also limit dispersal (Row et al., 2012).
Thus, we also included two relevant environmental predictors: Annual Drought Index and degree days above 5°C (Doherty et al., 2016). In both cases, thresholds were primarily derived based on the distribution of values and values that split the distribution of values into evenly spaced bins to minimize correlation between resistance distances within a predictor (Table 2).

| Within-variable resistance thresholds
We identified the top resistance surfaces using maximum-likelihood population-effects models (MLPE), which is a mixed modeling approach that accounts for nonindependence in pairwise datasets (Clarke, Rothery, & Raybould, 2002;Van Strien et al., 2012). For each of the individual habitat, abundance, and landscape predictors, we used two analyses for each MZ to determine (i) the thresholds and resistance values that produced resistance distances (RD) that were best correlated with genetic differentiation, and (ii) whether that correlation was significantly greater than its correlation with distance.
We determined the top model with Akaike's information criteria (AIC; Akaike, 1973) and estimated relative support using ΔAIC values compared with a model including only distance (i.e., the null model).
We also determined the significance of standardized coefficients for the underlying model (GenDist ~ RD1) and a model that also includes geographic distance (GenDist ~ RD1 + RD distance ). When the resistance coefficient confidence intervals were positive in both models, we considered the variable to have had a significant effect on gene flow over-and-above distance alone (Row et al., 2017).

| Relative importance of landscape and habitat predictors
Within each MZ, we compared all of the models with significant resistance variables against each other using ΔAIC in a betweenvariable analysis. There is some difficulty in identifying the correct multivariate models when resistance values are at different scales and have different distributions, but these results can be improved by calculating resistances from single landscapes that summarize multiple variables (Row et al., 2017). Therefore, in addition to comparing the univariate landscape predictors, we developed a combined landscape surface using all of the significant landscape predictors summed together. Cells with increased resistance for more than one landscape variable had an additive resistance value. We calculated pairwise resistances using this new combined additive surface and compared the fit with genetic differentiation to the fit for values derived from individual landscape variables and to the retained habitat predictors for each MZ.

| Individual-based genetic differentiation
Although using a population-based evaluation is most appropriate for our species and sampling regime, using an individual-based approach in landscape genetics may influence the resulting perceived influence of landscape features (Prunier et al., 2013). Therefore, as a first validation of our results, we repeated the between-variable analysis using an individual-based approach and compared the results. We calculated individual-based genetic differentiation using

ST
We validated our top selected resistance surface for each MZ using a Monte Carlo cross-validation approach within each MZ.
For each of 100 iterations, we split the data set into a training and testing dataset. Using the training dataset (80% of the data), we developed an MLPE model predicting genetic differentiation from the best resistance surface (GenDist ~ RD1) and distance alone (GenDist ~ RD distance ). Subsequently, we predicted genetic differentiation based on the model, and we performed a crossvalidation using the testing dataset (the remaining 20% of the data) by determining the absolute difference between the predicted and actual genetic differentiation for both models using 100 iterations. We grouped differences into predictions ranging across five increasing categories of genetic differentiation to evaluate the ability of the model to predict genetic differentiation across a range of values. lek cluster to all others sequentially and then averaged.

| Genetic differentiation within management zones
The number of lek clusters within MZs ranged from 15 to 81, and the mean number of individuals sampled per lek cluster was rela-  Figure S2).

| Landscape resistance surfaces
Effects of changing thresholds and resistance values on pairwise resistance varied among the MZs (Table S1). Most of this variation was tied to the abundance of the landscape classified as habitat, which varied between each MZ (Table S1). When a given landscape variable was rare on the landscape (e.g., tillage in MZ III and MZ V), changing the threshold or resistance had little effect on the pairwise resistance, resulting in high correlations in resistance across the thresholds. In most cases, changes in thresholds had a greater effect on resulting pairwise resistances than changes in resistance values (Table S1).
The importance of each landscape variable to functional con-  Anthropogenic landscape disturbances also had consistent negative influences on functional connectivity. Human disturbance, tillage, or both human disturbance and tillage had significant effects on functional connectivity in all but MZ III (Figure 2).
In MZ I, II, and IV, the thresholds for tillage were all 25% (Figure 2), suggesting disrupted functional connectivity on landscapes where row crops and small grain fields composed a greater percentage of the landscape than this threshold. In MZ V, the selected threshold was only 5%, but the overall percentage of the tilled agricultural in the landscape was also low (~2% on average; Table S1), and pairwise resistances were highly correlated (.99) with geographic distance. In MZ III, tilled agricultural similarly represented a small percentage of the landscape.

peared to have little effect on functional connectivity for most
MZs. Annual drought index (adi) was not significant for any MZs, and degree days above 5 (dd5) was only significant for MZ II ( Figure 2).

| Breeding habitat utilization and abundance resistance surfaces
Similar to the landscape variables, the effects of changing thresholds and resistance values for habitat surfaces on pairwise resistances varied among the MZs, but was much greater for changes in thresholds overall (Table S1). In three of the five MZs, breeding

| Relative importance of landscape and habitat predictors
In four of the five MZs, resistance values derived from BH or BPI provided the best fit with genetic differentiation (Appendix S2). In MZ Habitat variables TA B L E 3 Top models predicting pairwise genetic differentiation from pairwise resistance between lek groups of sage-grouse samples.

F I G U R E 2 Top univariate coefficients for each landscape and habitat variable (see details in
Resistance surface codes are a combination of the base predictor (

TA B L E 3 (Continued)
F I G U R E 3 Absolute difference between predicted and actual pairwise genetic differentiation for groups of greater sage-grouse samples with increasing levels of differentiation. Predictions for resistances from the top resistance model (see Table 3) and distance alone are Undiff MZ III, BPI provided a better fit over distance. In MZ IV, BH was the top surface over BPI-the best fit in the population-based analysiswhich was the second-best fitting model.

ST
In all zones, the top resistance model was better able to predict G ′

| Range-wide patterns of gene flow
Average gene flow between clusters identified connectivity between populations at local scales and at some regional scales within MZs (Figure 4). At the range-wide scale, connectivity was more variable with large areas of low modeled gene flow separating some On the California/Nevada border, two weak and diffuse paths connect the Bistate population to others to the north and east. Given F I G U R E 4 Range-wide gene flow map describing low (dark blue) to high (yellow) connectivity between greater sage-grouse lek groups used in landscape genetic analysis. For display flow was split into percentiles with bin labels representing the upper bound that many of these paths span between management zones, more research is required in these regions to test the extent to which they represent areas of active gene flow.

| D ISCUSS I ON
Our landscape genetic analysis of gene flow and causal resistance surfaces represents one of the largest single evaluations of functional connectivity for a wild terrestrial vertebrate. Our primary advancements herein include accounting for regional variation in factors driving functional connectivity, deriving thresholds in functional connectivity for a species at risk, and depicting those responses as spatially explicit resistance surfaces for use in range-wide conservation evaluations. Through noninvasive sampling techniques, the feasibility of developing broad-scale range-wide datasets is increasing.
Furthermore, ever-improving molecular genetic technologies are making it possible to generate large amounts of genotypic information from individuals in these datasets. Our findings can serve as a roadmap for other large-scale evaluations in landscape genetics aimed at identifying connectivity thresholds and validating established resistance surfaces within regions of importance to wildlife management.
As hypothesized, MZ-specific habitat metrics were the best predictors of differentiation. However, to our surprise, contributions of bird abundance-corrected indices were equivocal to model fit.
Collectively, these findings suggest that maintenance of breeding habitat above critical thresholds so as not to reduce genetic exchange is fundamental to conservation of connected populations.
Our inferences also add to a body of literature suggesting that structural connectivity estimated from habitat indices can be a strong predictor of functional connectivity for many species (Hagerty, Nussear, Esque, & Tracy, 2011;Row et al., 2010Row et al., , 2015Shafer et al., 2012;Wang et al., 2008), but we stress the importance of including model validations to have a clear understanding of the predictive capabilities of landscape genetic models.
Herein, BH models in each zone were developed independently and the relative importance of the variables underlying each breeding habitat model varied. Broadly speaking, the variables with the greatest importance were a positive association with sagebrush cover and a negative association with increased tree canopy cover for most zones. When resistance surfaces from these landscape variables were compared independently to genetic differentiation, they were similarly found have a greater fit than distance alone, demonstrating their overall importance. Interestingly, the distribution of Annual Drought Index (MZ II, MZ IV, MZ V) and degree days >5°C (MZ II, MZ III, MZ IV) were prominent variables in some of the habitat models, but when analyzed independently they had little explanatory power in our assessment of functional connectivity. Thus, their importance may be connected or correlated with other variables.
Using our threshold approaches, we found that landscapes with a probability of occurrence for breeding leks <0.25 or 0.5 reduced functional connectivity. These values are likely below the threshold for persistence, as all known active breeding leks are present in regions with values >0.65 (Doherty et al., 2016). The lower value for functional connectivity is not surprising, as individuals are often willing to disperse through undesirable habitat (Tischendorf & Fahrig, 2000), but has important conservation implications on how to manage landscapes to preserve functional connectivity for sage-grouse.
For example, although the habitat may be degraded below the 0.65 threshold for breeding lek formation, it will still be important to maintain habitat above the lower thresholds identified here in order to maintain functional connectivity.
It is important to point out that landscape genetic results can vary at different spatial scales (Anderson et al., 2010) and habitat configuration can influence both dispersal (D'Eon, Glenn, Parfitt, & Fortin, 2002) and habitat selection (Wisdom, Meinke, Knick, & Schroeder, 2011). Furthermore, different threshold values had different overall effects on pairwise resistance and on model fit for each management zone. Thus, in planning for conservation, conducting analyses at different spatial extents and testing the effect of the different thresholds identified here (0.5, 0.25) and the resulting habitat configuration will provide valuable insights.
How individuals respond to habitat structure will vary between species. Thus, comparing other species that overlap our study area, especially other sagebrush obligate species, could provide insight into the generality of the results found here. We are not aware of other studies that have utilized a threshold approach. However, both Wang et al. (2008) and Row et al. (2010) found that the best fitting resistance models assigned very low-quality habitat disproportionately higher resistances or set these habitats as absolute barriers.
These results suggest that thresholds in habitat indices may be an effective approach for other species. Yet, in some cases, continuous habitat surfaces have performed better than discrete values (Hagerty et al., 2011), or individual landscape predictors have performed better than combined habitat indices altogether (Roffler et al., 2016;Wasserman, Cushman, Schwartz, & Wallin, 2010). A weak association between habitat and functional connectivity is likely when there are large differences between daily use and dispersal habitat or when one or a few landscape components are the primary drivers of functional connectivity. It is also hard to compare between studies when different approaches are used to derive resistance surfaces. For example, both Wasserman et al. (2010) and Roffler et al. (2016) only tested habitat resistance surfaces with a direct linear relationship between habitat and landscape resistance.
In contrast, discrete values (tested here and by Wang et al., 2008 and by Row et al., 2010) or the exponential relationships between habitat and resistance (Row et al., 2015) often appear to provide better model fit.
We predicted that the inclusion of abundance information would improve the fit between habitat and genetic differentiation. Local population size has long been linked to patterns of dispersal and can influence source-sink dynamics (Matthysen, 2005;Ozgul, Oli, Armitage, Blumstein, & Van Vuren, 2009). Furthermore, population size will influence the effects of genetic drift and have impacts on spatial patterns of genetic differentiation (Hutchison & Templeton, 1999). However, contrary to our predictions including abundance information into the resistance surfaces did not improve the model fit for most MZs. In our analysis, inclusion of abundance indices acted to decrease resistance around population centers and increase the resistance in areas with low abundance, but our approach could not test the effect of population size on the attractiveness of individual locations to dispersers. Thus, using models where local abundance can modify the number of dispersers into, or out of, a population (Murphy, Dezzani, Pilliod, & Storfer, 2010;Row, Oyler-McCance, & Fedy, 2016), might provide more insight into the effects of abundance on differentiation. We only included abundance estimates of males during the breeding season, but it is possible that abundance during other seasons or estimates of effective population size may be of greater importance to functional connectivity.

| Resistance thresholds and spatial variation in landscape predictors
Terrain roughness shaped functional connectivity across all MZs, and sagebrush availability (+) or tree canopy cover (−), and extent of cultivation (−) were influential in all but MZ III, the most highly fragmented zone. Time lags exist from when a barrier to dispersal first arises on the landscape and when the influence of that barrier on functional connectivity can be detected . For sage-grouse, the topography of the terrain can have a strong influence on habitat use and distribution (Davis, Reese, Gardner, & Bird, 2015;Fedy et al., 2014), and these imposed restrictions to movement or barriers have likely remained static for millennia. Indeed, in other large-scale landscape genetic evaluations of sage-grouse, areas with higher landscape ruggedness, a measure of sharp changes in elevation, restricted gene flow (Row et al., 2015). Similarly, in our individual landscape analysis, we found increases in steepness or roughness reduced functional connectivity in all of the five MZs, likely resulting from avoidance of these areas by dispersing sagegrouse. Although these are natural landscape features and will not be the target of conservation initiatives, it is clear that they impact dispersal and so it is important to consider and control for terrain in future evaluations.
In contrast to terrain topography, the current distribution of sagebrush has been modified through development and land conversion (Davies et al., 2011;Welch, 2005). Thus, many resistance features have likely been created or modified since human settlement. Sage-grouse rely on sagebrush for both food and shelter and sagebrush is a strong predictor of sage-grouse habitat across seasons and spatial scales (Connelly, Knick, Schroeder, & Stiver, 2004;Doherty et al., 2016;Fedy et al., 2014), and it influences functional connectivity (Row et al., 2015). Similarly, we found that areas with low sagebrush cover impeded gene flow in three of the five MZs. In two cases, we found that sagebrush cover <30% impacted dispersal, with a 10% threshold in the third, suggesting that sagebrush cover <10%-30% reduces gene flow. As with the habitat thresholds, sagegrouse are capable of dispersal through habitats in which they are unlikely to persist (suggested persistence thresholds are in the range of 40%-65% sagebrush cover; Aldridge, Nielsen, & Boyce, 2008;Wisdom et al., 2011;Knick, Hanser, & Preston, 2013).
Given the importance of sagebrush to sage-grouse, it is surprising that the distribution of sagebrush cover was not a significant predictor of genetic differentiation in two of the MZs. Variation in the importance of a predictor can be related to its abundance and distribution, but this does not seem to be the case here, as mean values were similar across all zones. In MZ I, where individuals have been shown to move far distances (Cross et al., 2017;Newton et al., 2017), isolation by distance appeared to drive differentiation. None of the landscape predictors were very strong in MZ III, a zonal boundary spanning multiple populations and/or habitat-population relationships. Perhaps this should be expected as MZ III is set in basin and range topography with this natural fragmentation exacerbated by conifer encroachment and land-use change (Chambers et al., 2017).
In both of these zones, sagebrush was an important component of the derived habitat indices suggesting that its distribution is important in combinations with other landscape variables.
A long history of fire control has enabled encroaching conifer woodlands to degrade sagebrush habitats into areas with higher amounts of tree canopy cover. Sage-grouse avoid canopy cover at low levels (<4%, Miller, Naugle, Maestas, Hagen, & Hall, 2017) or stay and suffer demographic impacts . In conifer removal areas, females readily nested in restored sites Severson et al., 2017) and were more successful in raising their broods (Sandford et al., 2017). We build on this knowledge to add that connectivity among population centers is reduced when conifer expansion exceeds a 10% threshold in canopy cover. Connectivity is relevant to management because conifer-encroached habitats stimulate faster yet riskier movements, especially in juveniles, that may make sage-grouse more vulnerable to visually acute predators with demonstrated fitness consequences . Future restoration planning with the goal of improving genetic connectivity can use our range-wide resistance surfaces to select areas where the greatest benefit may be found.
Cultivation is known to reduce breeding populations (Doherty et al., 2016;Smith et al., 2016). Findings here further suggest that cultivation reduces gene flow when >25% of the landscape is converted to cropland in three of five MZs tested. In eastern Montana, where cultivation is most prevalent, 70% of the best sage-grouse habitat is privately owned. Therefore, activities that keep sagebrush habitats intact (such as large working cattle ranches) as opposed to those that do not (such as agriculture) should help maintain connectivity between population strongholds. Cultivation land use was not prevalent in the other two zones (MZ III, V) so we had low power to determine its effects on functional connectivity therein.

| Resistance model validation
Landscape genetic studies typically quantify correlations between patterns of genetic differentiation or gene flow with landscape resistance (or cost) distances and provide insight into the relative importance of landscape variables influencing functional connectivity (Manel et al., 2003). Here, we not only evaluated these common objectives, but also tested the ability of the top linear mixed models to predict genetic differentiation among populations based on landscape resistance. Predictive ability of our resistance maps varied quite dramatically among MZs and also with the level of differentiation between the populations considered. In most cases, improvement in predictive ability of our models when compared to the null distance model was greatest when pairwise comparisons had little (MZ III) or much (MZ IV, MZ V, MZ I) differentiation among populations and in one case where differentiation among populations varied (MZ II). Where our models did not perform better, results indicated that distance was the predominate driver of differentiation and that our models did not add much predictive improvement.
For example, in MZ IV, plots of genetic differentiation and distance ( Figure S2) reveal that highly differentiated populations are not explained well by distance alone due to the high value of residuals; the resistance models for this zone have much better predictive power for these highly differentiated groups. Conversely, in MZ III, the residuals are larger for populations with low differentiation, which corresponds to the improved performance of the resistance models.
Comparing performance of our models to other studies is not possible as we are unaware of others who have validated the predictive power of their landscape resistance models. The potential power of landscape genetics for informing management will be greatly improved if the model results can be used to predict genetic differentiation between regions where genetic data are lacking or to predict how changes to the landscape are likely to impact functional connectivity. Predictions backed by any confidence require more validation than is current practice, and we present a novel framework for providing this necessary model validation.

| Implications for sage-grouse management
Currently, sage-grouse conservation is largely focused on implementing beneficial conservation measures within population strongholds (e.g., Priority Areas of Conservation; U.S. Fish and Wildlife Service 2013) and around known sage-grouse leks as they represent areas of importance for breeding and early brood rearing habitat (Doherty et al., 2016;Fedy et al., 2014). However, our analyses and resulting resistance surfaces point to several measures that can be taken to help improve and maintain functional connectivity for sagegrouse. First, although population strongholds likely have much higher suitability values, maintaining areas outside of these regions above habitat thresholds of 0.5, or potentially 0.25 in some management zones, should help maintain connectivity between these existing protection areas. Secondly, our models could help identify landscapes where targeted conservation would maximize conservation return on investment. For example, a 100-year history of fire suppression has enabled conifer expansion into sagebrush habitats, reducing lek attendance, breeding habitat quality and survival of sage-grouse Miller et al., 2017). We found that functional connectivity appeared to drop at around 10% and, thus, a conifer removal strategy incorporating known dispersal pathways from our current flow map (Figure 4) may help to maintain connectivity between population strongholds.
In addition to the thresholds we identified, the resistance surfaces and gene flow maps we generated help identify areas within which to prioritize management actions. Resistance maps identify areas that are above and below threshold values that obstruct gene flow and direct conservation actions within these areas where it is possible to maintain or improve habitat above or below a targeted threshold. Furthermore, because the nodes in our analysis represent clusters of active sage-grouse leks, the modeled gene flow should reflect movement from these high density areas and, as such, can be used to help locate and protect dispersal corridors. Additional gene flow maps can be produced among management areas of particular interest to managers and used to target conservation initiatives that will maintain connectivity among population strongholds.
It was clear from our cross-validation that the predictive ability of our resistance models varied with the levels of genetic differentiation and among management zones. Even when our results strongly suggested an improvement in model fit when compared to the null distance model, the overall predictive ability of our models was at times marginal or poor depending on the amount of genetic differentiation among populations. Without our cross-validation to provide an estimate of predictive ability, conservation initiatives could direct actions that will not have the desired improvement on connectivity.
Overall, our cross-validated approach used in developing our threshold resistance surfaces for sage-grouse should initiate a new era of spatial analyses which emphasizes the value of functional connectivity and the identification of habitats supporting it.

ACK N OWLED G EM ENTS
This research would not have been possible without the many state, federal, and non-governmental organization (NGO) biologists and technicians who have spent lifetimes working tirelessly in the field and behind desks for the conservation of the greater sagegrouse. We specifically acknowledge field biologists and techni- Ecology Laboratory. This draft manuscript is distributed solely for purposes of scientific peer review. Its content is deliberative and predecisional, so it must not be disclosed or released by reviewers.
Because the manuscript has not yet been approved for publication by the U.S. Geological Survey (USGS), it does not represent any official USGS finding or policy. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service.

DATA A R C H I V I N G S TAT E M E N T
Cluster-based genotypic data used for the analysis in this article are available on the USGS ScienceBase website: https://www.sciencebase.gov/ (https://doi.org/10.5066/f7rb73v0).