Scale-dependent influence of the sagebrush community on genetic connectivity of the sagebrush obligate Gunnison sage-grouse.

Habitat fragmentation and degradation impacts an organism's ability to navigate the landscape, ultimately resulting in decreased gene flow and increased extinction risk. Understanding how landscape composition impacts gene flow (i.e., connectivity) and interacts with scale is essential to conservation decision-making. We used a landscape genetics approach implementing a recently developed statistical model based on the generalized Wishart probability distribution to identify the primary landscape features affecting gene flow and estimate the degree to which each component influences connectivity for Gunnison sage-grouse (Centrocercus minimus). We were interested in two spatial scales: among distinct populations rangewide and among leks (i.e., breeding grounds) within the largest population, Gunnison Basin. Populations and leks are nested within a landscape fragmented by rough terrain and anthropogenic features, although requisite sagebrush habitat is more contiguous within populations. Our best fit models for each scale confirm the importance of sagebrush habitat in connectivity, although the important sagebrush characteristics differ. For Gunnison Basin, taller shrubs and higher quality nesting habitat were the primary drivers of connectivity, while more sagebrush cover and less conifer cover facilitated connectivity rangewide. Our findings support previous assumptions that Gunnison sage-grouse range contraction is largely the result of habitat loss and degradation. Importantly, we report direct estimates of resistance for landscape components that can be used to create resistance surfaces for prioritization of specific locations for conservation or management (i.e., habitat preservation, restoration, or development) or as we demonstrated, can be combined with simulation techniques to predict impacts to connectivity from potential management actions.

INTRODUCTION

Landscape composition influences the way organisms navigate and distribute themselves (Fischer & Lindenmayer, 2007; Heinrichs et al., 2016; Sala et al., 2000) in space. When a species’ habitat is fragmented or degraded conservation concern increases. The maintenance of connectivity, or movement of individuals among populations, reduces the risk of genetically structured populations, loss of genetic diversity to drift, inbreeding depression, and extinction risk (Burkey, 1989; Hoffmann et al., 2017; Soule et al., 1992). Despite the link to known genetic consequences, much of what we understand about connectivity is derived from direct observations of movement (Epps et al., 2005; Fahrig & Merriam, 1985; Frankham, 2003) and may misrepresent the effect on gene flow (Gubili et al., 2017; Levin, 1981; Spear et al., 2010). Direct observation of individual movement can provide insight about how organisms use or move through the landscape within a limited time‐period. However, not all movements result in gene flow; a process which can occur over many generations and among populations not directly connected by dispersal events (Spear et al., 2010). Similarly, ecological processes shape populations differently across a spatiotemporal continuum (Johnson, 1980; Zeller et al., 2017). Yet few investigations of demographic or genetic connectivity have evaluated the influence of landscape components at more than one scale (Bissonette, 1997; Cushman & Landguth, 2010). Explicit characterization of the scale‐dependent relationship between gene flow and the landscape would result in a better understanding of the potential impact of anthropogenic or climate‐related landscape alteration on evolutionary capacity and provide more information for conservation and management decision‐making (Holderegger & Wagner, 2008; Manel et al., 2003; Storfer et al., 2010).

The early 2000s has seen numerous applications of landscape genetic connectivity analyses for wildlife species, attempting to gain insight on how landscape composition affects gene flow (e.g Gerlach & Musolf, 2000; Row et al., 2015). However, many analytical approaches have significant limitations; namely a reliance on expert opinion‐based parameterization of resistance, a lack of formal measures of uncertainty, and low statistical power (Bohonak, 1999; Graves et al., 2013; Keller et al., 2015; Koen et al., 2012; Landguth et al., 2010; Shirk et al., 2015; Zeller et al., 2016). Additionally, relatively few landscape genetic connectivity approaches consider more than a single scale (though see Galpern et al., 2012; Keller et al., 2013; Kozakiewicz et al., 2019). A recently developed statistical model of genetic connectivity extends the previous use of circuit theory as analogous to gene flow across a landscape (e.g McRae, 2006; Spear et al., 2010) and improves upon previous methods in multiple ways. Improvements include eliminating the use of expert opinion by directly estimating degree of impediment or facilitation for each considered landscape component, incorporating an estimate of uncertainty, providing a statistically supported framework for evaluating multiple variables in a single model, and allowing for model selection (Hanks & Hooten, 2013; Peterson et al., 2019). To date, this model has not been widely implemented.

One of the main goals of landscape genetics is to develop conservation planning tools, a goal frequently hindered by an inability to compare the relative influence of different actions (Keller et al., 2015). This formal statistical model for genetic connectivity lends itself to directly evaluating conservation actions. The best fit modelled relationship can be used to infer a resistance surface that reflects changes due to different potential conservation actions, proposed land‐use change, or climate predictions and can be used in simulations to evaluate the potential impact to evolutionary capacity and to compare the consequences of alternative scenarios, analogous to other simulation‐based conservation planning tools for wildlife (e.g Bennett et al., 2009; McRae et al., 2008). We utilized the statistical model for connectivity developed by Hanks and Hooten (2013) to evaluate the scale of influence for each landscape component and formed competing hypotheses of connectivity for two biologically relevant spatial scales. We subsequently extended the simulation approach described in Peterson et al. (2019) to predict impacts from hypothetical management actions using Gunnison sage‐grouse (Centrocercus minimus) as an example of the utility of landscape genetic analyses for conservation and management decision‐making.

Land‐use change in the sagebrush steppe of western North America has been extensive (Bock & Webb, 1984; Braun, 1998; Knick et al., 2013). Over 50% of the sagebrush (Artemisia sp.) distribution has been removed or degraded (Connelly et al., 2004; Knick et al., 2013; Schroeder et al., 2004), imperiling the resident sagebrush obligate wildlife species (Anderson & Inouye, 2001). The Gunnison sage‐grouse is one such sagebrush obligate species (Patterson, 1952; Young et al., 2000) impacted by the alteration of western landscapes (Oyler‐McCance et al., 2001; Primack, 1993; Theobald et al., 1996). The species historically occupied the naturally rough terrain of the southeastern extent of sagebrush steppe as population centres connected by dispersal through intervening sagebrush corridors (Braun et al., 2014). Anthropogenic fragmentation and removal of sagebrush habitat is assumed to have led to significant range contraction around the population centres, low migration among populations, genetic differentiation, low genetic diversity (Oyler‐McCance et al., 2005, 2015), and significant conservation concern (United States Fish and Wildlife Service, 2014). Currently the species occupies just 8% of the historical distribution (Braun et al., 2014; Schroeder et al., 2004) as a single, relatively large population in the Gunnison Basin of Colorado, USA (~85%–90% of remaining birds) and six smaller satellite populations (United States Fish and Wildlife Service, 2014). Individuals have relatively high fidelity for breeding grounds (i.e., leks) within populations, although movement among leks occurs relatively frequently (Aldridge et al., 2012; Oyler‐McCance et al., 2005). The natural hierarchical structure of the species distribution (e.g., discrete breeding grounds nested within discrete populations) increases the likelihood of scale‐dependent effects of the landscape on gene flow and reinforces the importance of considering scale when investigating gene flow patterns. Understanding the influence of landscape features on gene flow among and within populations can provide insights about conservation and management given that the human population is expected to continue growing for at least the next three decades (Colorado Water Conservation, 2009) and landscape composition to continue changing (Theobald et al., 1996; United States Fish and Wildlife Service, 2014).

Previous work directly observing sage‐grouse (Centrocercus spp.) habitat selection, or movement through the landscape (i.e., telemetry), provide some insight into the features likely to influence gene flow. For example, sage‐grouse require sagebrush as a food source and for concealment from predators (Barnett & Crawford, 1994; Braun et al., 2005; Patterson, 1952; Sveum et al., 1998), while they generally avoid conifer cover (Baruch‐Mordo et al., 2013; Hagen et al., 2011) and anthropogenic features (Aldridge et al., 2008, 2012; Blickley et al., 2012; Copeland et al., 2013). Additionally, a developed understory is important for sage‐grouse brood rearing (Aldridge & Brigham, 2002; Barnett & Crawford, 1994) and nest concealment (Aldridge & Brigham, 2002; Webb et al., 2012). Gunnison sage‐grouse populations, and leks within populations, are located within a variable range of climatic and topographic features, both of which are known to impact habitat use, movement, and persistence (Aldridge et al., 2008; Blomberg et al., 2012; Harju et al., 2013; Knick et al., 2013). Direct observation provides insight regarding short‐term impacts to movement, but the degree to which these insights apply to long‐term evolutionary impacts via gene flow is unknown (Levin, 1981; Roffler et al., 2016).

We fit a formal statistical model for genetic connectivity (Hanks & Hooten, 2013) to data for two biologically‐relevant, hierarchical, spatial scales for the Gunnison sage‐grouse, a species of conservation concern in an increasingly fragmented and degraded habitat. Our main objective was to identify the primary drivers of coarse‐scale interpopulation and fine‐scale intrapopulation genetic connectivity and directly estimate the degree of impediment or facilitation for each identified landscape component. We also wanted to illustrate the conservation applications of the identified relationships between genetic connectivity and landscape composition. We used our best fit models and simulations to demonstrate how hypothetical management actions might impact gene flow in future generations. We generally expected more contiguous sagebrush habitat to facilitate gene flow and rough terrain and anthropogenic modification to impede gene flow. An improved understanding of how the landscape influences genetic connectivity and how management actions might enable the achievement of specific targets for gene flow could ultimately result in better conservation outcomes.

MATERIALS AND METHODS

Study area

We defined our two study areas by buffering known lek locations with Gunnison sage‐grouse movement distances. We created the rangewide extent by applying a 54.68‐km buffer around each known active lek to account for maximum observed dispersal distance (2 × 17.34 km; D. J. Saher, unpublished data) and an additional 20‐km buffer to account for edge effects in our calculation of moving window averages for pixels (see Spatial Covariates section below and Appendix S1). The resulting rangewide extent covered ~56,818 km2 (Figure 1) of frequent changes in cover type, including anthropogenic features, land conversions (e.g., agriculture), and mountains (Gunnison sage‐grouse Rangewide Steering, 2005). We evaluated rangewide connectivity among populations as a function of the intervening landscape within the described study area. The Gunnison Basin extent covered ~6282 km2 and was previously created for development of seasonal resource selection models (Aldridge et al., 2012). Within the Gunnison Basin, we evaluated connectivity among leks within native sagebrush‐steppe (dominated by big sagebrush, Artemisia tridentata ssp.) interrupted by anthropogenic features and variable topography (Gunnison sage‐grouse Rangewide Steering, 2005).

FIGURE 1

Historical (grey) and current (yellow) distribution of Gunnison sage‐grouse in the southwestern United States. Populations are labeled with respective names. Black polygon designates the rangewide study area and the hatched polygon delineates the Gunnison Basin study area. The historical and current distribution maps were developed by Schroeder et al. (2004); the two northernmost portions of the historical distribution correspond to an unknown species of sage‐grouse and are not verified by Colorado Parks and Wildlife (Gunnison sage‐grouse Rangewide Steering, 2005)

Genetic samples

Our rangewide data set was composed of 254 blood samples collected from birds captured using spotlight trapping methods (Giesen et al., 1992; Wakkinen et al., 1992) between 1996 and 2004 as part of a previous study (Oyler‐McCance et al., 2005), avoiding confounding effects of translocations from the Gunnison Basin to the satellite populations that began in 2005 (United States Fish and Wildlife Service, 2014). Samples were distributed across populations as follows: Cimarron = 4, Crawford = 21, Dove Creek = 43, Gunnison Basin = 116, Piñon Mesa = 19, San Miguel = 51 (Figure 1). We excluded the recently re‐established Poncha Pass population (Nehring & Apa, 2000). Our Gunnison Basin data set was composed of 624 unique samples representing 49 of the 70 leks considered active between 2006 and 2014 (92% of samples collected after 2010), also as part of a previous study (Zimmerman et al., 2019). In brief, samples were primarily noninvasively collected feathers, but also included mortalities, and blood samples collected from captured individuals. Leks not represented in our data could not be accessed or no feathers were found when searched. The number of samples collected at a lek ranged from 1 to 57 (See Appendix Table S2.1 for all lek sample sizes). The Gunnison Basin samples used for each extent are distinct.

Details on DNA extraction, microsatellite characterization, and duplicate sample elimination were previously published (Zimmerman et al., 2019). Samples were genotyped with 22 microsatellite loci (Caizergues et al., 2003; Fike et al., 2015; Oyler‐McCance & St. John, 2010; Piertney & Höglund, 2001; Segelbacher et al., 2000; Taylor et al., 2003), using polymerase chain reaction (PCR) with components and concentrations as described in Oyler‐McCance and Fike (2011) and thermal profiles as originally published. Duplicate noninvasive samples (Gunnison Basin only) were identified and removed using a combination of the R package allelematch (Galpern et al., 2012) and the stand‐alone program Dropout (McKelvey & Schwartz, 2005). We calculated an individual‐based pairwise absolute genetic distance matrix (D) based on a Manhattan distance in R (R Core Team, 2018) as in Hanks and Hooten (2013), for the rangewide and Gunnison Basin data sets. Genetic distance matrices were used as the dependent variable in our connectivity modelling as described in Hanks and Hooten (2013). Individuals sampled within the same population (rangewide analysis) or from the same lek (Gunnison Basin analysis) were treated as repeated samples at a location in our connectivity model.

Spatial covariates

We identified spatial covariates likely to influence the ability of sage‐grouse (Centrocercus spp.) to navigate the landscape based on existing literature. Thus we considered anything that affects dispersal, resource selection, survival, fecundity, occupancy, avoidance, or behaviour. The landscape components identified included the following: habitat structure (e.g., shrub height) and amount of sagebrush cover (Aldridge & Boyce, 2007; Aldridge et al., 2008, 2012; Baruch‐Mordo et al., 2013; Doherty et al., 2010; Fedy et al., 2014; Harju et al., 2013; Knick et al., 2013; Oyler‐McCance et al., 2001; Rice et al., 2017), amount of conifer cover (Baruch‐Mordo et al., 2013; Commons et al., 1999; Fedy et al., 2014; Hagen et al., 2011; Knick et al., 2013; Severson et al., 2017), conifer configuration or aggregation (Baruch‐Mordo et al., 2013), normalized difference vegetation index (NDVI) as a proxy for phenology (Aldridge & Boyce, 2007; Aldridge et al., 2012), soil wetness as described by compound topographic index or CTI (Aldridge & Boyce, 2007; Carpenter et al., 2010), seasonal habitat selection models for Gunnison sage‐grouse (Aldridge et al., 2012), agricultural cover (Aldridge & Boyce, 2008; Beck & Maxfield, 2003; Bush et al., 2011; Fedy et al., 2014; Knick et al., 2013), proportion of anthropogenic development (Aldridge et al., 2012; Knick et al., 2013; Rice et al., 2017), distance to development (Aldridge & Boyce, 2007; Aldridge et al., 2012), human population density (Aldridge et al., 2008), distance to human population density (Aldridge & Boyce, 2007; Aldridge et al., 2008), road density (Aldridge & Boyce, 2007; Aldridge et al., 2012; Fedy et al., 2014; Knick et al., 2013), distance to roads (Aldridge et al., 2012; Rice et al., 2017), oil and gas well density (Aldridge & Boyce, 2007; Blickley et al., 2012; Copeland et al., 2013; Green et al., 2017; Smith et al., 2014; Tack et al., 2011), distance to oil and gas wells (Aldridge & Boyce, 2007; Blickley et al., 2012; Copeland et al., 2013; Tack et al., 2011), slope (Harju et al., 2013; Knick et al., 2013), topographic roughness as described by the terrain ruggedness index or TRI (Fedy et al., 2014; Harju et al., 2013; Knick et al., 2013), annual precipitation (Blomberg et al., 2012; Fedy et al., 2014), mean maximum temperature (Blomberg et al., 2012), growing degree days, and a dryness index (Aldridge & Boyce, 2008; Table 1). Wind energy infrastructure can also impact movement and survival of sage‐grouse (LeBeau et al., 2013, 2017). Presently, there are no major existing wind energy developments within the Gunnison sage‐grouse range (United States Fish and Wildlife Service, 2014). We considered all identified landscape components influencers of connectivity.

TABLE 1

Covariates used to model connectivity

Abv.	Covariate a	Source	Predicted impact
Habitat composition and configuration
AS	Proportion of all sagebrush cover	Landfire	+
BS	Proportion of big sagebrush cover	Landfire	+
CON	Proportion of conifer cover	Landfire	−
LS	Proportion of low sagebrush cover	Landfire	+
CC	Conifer configuration	Landfire derived	−
PBS	Percentage big sagebrush cover	NLCD shrubland product	+
PSB	Percentage all sagebrush cover	NLCD shrubland product	+
SBHT	Shrub height	NLCD shrubland product	+
Habitat selection
N	Nest habitat RSF	Aldridge et al. (2012)	+
Climate and phenology
MMT	Mean maximum temp.	PRISM	−
MAR	Mean annual rainfall	PRISM	+
DRI	Dryness index	Daymet derived	−
GDD	Growing degree days	Daymet derived	+
FNDVI	Fall NDVI	MODIS	+
SNDVI	Spring NDVI	MODIS	+
Terrain morphology
CTI	Compound Topo. Index	EROS	+
S	Slope	DEM derived	−
TRI	Terrain Ruggedness Index	DEM derived	−
Anthropogenic change
DI14	Dist. To Class 1–4 Roads	NAIP derived	+
DI12	Dist. To Class 1 & 2 Roads	NAIP derived	+
DI47	Dist. To Class 4–7 Roads	NAIP derived	+
DIA	Dist. To All Roads	NAIP derived	+
DIB	Dist. To BLM Roads	NAIP derived	+
DID	Dist. To Development	NLCD	+
DILD	Dist. To Light Duty Roads	aerial imagery/TIGER derived	+
DIOG	Dist. To Oil & Gas Wells	States of CO & UT	+
DIP	Dist. To Pop. Dens.	LandScan	+
DOG	Dens. Of Oil & Gas Wells	States of CO & UT	−
DAG	Proportion Of Agriculture	NLCD	−
DA	Dens. Of All Roads	NAIP/TIGER dervied	−
DB	Dens. Of BLM Roads	NAIP derived	−
DD	Proportion Of Development	NLCD	−
DP	Dens. Of Population	LandScan	−
D14	Dens. Of Roads Class 1–4	NAIP/TIGER derived	−
D12	Dens. Of Roads Class 1 & 2	NAIP/TIGER derived	−
D2	Dens. Of Roads Class 2	TIGER	−
D47	Dens. Of Roads Class 4–7	NAIP derived	−
DLD	Dens. Of Light Duty Roads	aerial imagery/TIGER derived	−

The data source (source), abbreviation (Abv.), and predicted impact (positive [+] or negative [−]) to connectivity associated with each covariate are listed below. Covariates are organized into general categories (italic). The different sagebrush layers were based on different data types: Landfire derived layers are based on a binary landcover classification and the NLCD shrubland product layers are based on percentage of cover.

Calculation of all covariates and details on data source are described in Appendix S1.

Like other wildlife species, sage‐grouse respond to different landscape components and characteristics in different ways and at different scales (Aldridge et al., 2012; Wiens, 1989; Wiens & Milne, 1989). Thus, we used a circular moving window to summarize the landscape variables in a neighbourhood larger than a single pixel to create versions of each landscape variable representing different scales. We used a circular moving window centred on each pixel for all landscape variables calculating the mean within a radius of 564, 1000, 3000, 6400, 10,000, 15,000, and 20,000‐m. Each considered radius had some support in the sage‐grouse literature for actual movement distances up to 6400‐m (e.g Aldridge et al., 2012; Fedy et al., 2012; Fremgen et al., 2016; Taylor et al., 2013). Additional radii were considered given potential differences in scale of direct observations versus gene flow. All spatial data processing was performed in ArcMap 10.1 (Environmental Systems Research, 2012). Details on covariate calculation and spatial data processing are included in Appendix S1. We predicted high sagebrush cover, high‐quality habitat (as defined by resource selection models), more rain and growing degree days, and high CTI would facilitate gene flow. Conversely, we predicted high conifer cover and/or aggregation, high temperatures, drier climate, steeper slopes, rough terrain, and anthropogenic modification would impede gene flow (Table 1).

Model fitting

We fit a formal statistical model based on circuit theory (Hanks & Hooten, 2013) to our genetic data to directly estimate the degree of resistance or facilitation to gene flow imposed by landscape covariates in a Bayesian framework. Our model was specified as follows: where is the absolute genetic distance matrix for individuals (we used the Manhattan distance), is the the covariance matrix of the observed nodes (i.e., populations/leks) accounting for variability () in repeated measures (i.e., multiple individuals) at a node. In the above model, must be negated to use the mathematical relationships for estimation of edge weights as analogous to electrical conductance. The graph structure of the nodes connected by edges (i.e., effective dispersal or gene flow) is used to define the correlation structure of the measured genetic distance () among nodes (Equation 1). We obtained the covariance matrix by calculating edge weights () for each pair of nodes for a spatial covariate (i.e., raster layer) for the entire spatial extent of interest, and then found the inverse of the covariance matrix incorporating the influence of intervening nodes on the correlation structure. In calculating edge weights (), and are locations of two nodes, is the vector of values of spatial covariates at node , is the vector of values of spatial covariates at node , and is the Euclidean distance between nodes and . The edge weight () is the conductance between nodes and is proportional to the transition rate from one location () to another (). The exponential link function (Equation 2) ensures intuitive interpretation of coefficients, where negative values denote degree of impediment and positive values degree of facilitation. If and are not first order neighbours () the edge weight is zero; nodes not directly connected, however, may still be connected through intervening nodes. The covariance matrix of the observed nodes, , is calculated as , where is a matrix relating each sample to a node, is the identity matrix, and is a spatial nugget parameter representing the variability in repeated observations (i.e., individual birds) collected from the same location (i.e., node in the graph). We used functions in the rwc R package (Hanks, 2018) for data processing, the generalized Wishart probability distribution, to relate sample locations to raster grids, and to calculate covariance structure. Before fitting models to our data, we confirmed our model fitting approach could recover the true relationships with simulated data (see Appendix S1 for validation of our approach). Variability associated with coefficient estimates for each landscape covariate was typically large and credible intervals (CRIs) frequently included zero (Appendix S1 Table S2.2–2.5); however our modelling approach appropriately identified the true model and the mean of the coefficients approximated the true resistance values.

We first fit single covariate models for each scale (i.e., moving window radius) as both a linear and quadratic relationship. Each model fit included two independent chains of 50,000 Markov chain Monte Carlo (MCMC) iterations with the first 5000 iterations discarded as a burn‐in period, using a random walk Metropolis‐Hastings sampler. We evaluated convergence of independent chains through visual inspection of trace and density plots (not included here) and with the Gelman‐Rubin diagnostic (Gelman & Rubin, 1992). We also fit a null hypothesis model based on an undifferentiated landscape (i.e., an intercept‐only model) testing the isolation by distance null hypothesis (Wright, 1943) for comparison. We next fit multicovariate models using all combinations of covariates. Different forms of a single covariate were highly correlated (Pearson's r > 0.8). Consequently, we retained the form with the lowest deviance information criterion (DIC; Spiegelhalter et al., 2002) for combination in multicovariate models. Convergence was problematic when a model included moderately correlated covariates (Pearson's r > |0.50|), we therefore chose variables based on DIC rank and excluded lower ranked variables correlated at this level. We did not eliminate variables for combination with coefficient CRIs that included zero because our simulation‐based evaluation of model fitting showed that coefficient CRIs frequently included zero but the mean coefficient value closely approximated the simulated true coefficient. Each model was fit and checked for convergence as described above, with the exception that our multicovariate models required longer runs for convergence resulting in 150,000 MCMC iterations, discarding the first 50,000 iterations as a burn‐in period. The top model was chosen through DIC rank. All model fitting was accomplished in R.

We used a 5‐fold cross‐validation procedure to evaluate the generalizability of our models to unobserved data (Conn et al., 2018; Hooten & Hobbs, 2015). We split our data into five approximately equal partitions. We fit the model five times, each time withholding a different partition of individuals, using the withheld data to calculate the log predictive density (lpd) for each iteration. A cross‐validation (CV) score based on the lpd for all five folds of the data was calculated as in Hobbs and Hooten (2015). We performed CV for the top 10 DIC ranked models for each spatial extent (rangewide and Gunnison Basin). The highest CV score corresponds to the model with the best predictive ability in the evaluated model set. We also evaluated the ability of the resistance surface created from our top ranked model to give rise to the observed data. We simulated genetic data from resistance surfaces created from the top ranked modelled relationship and an undifferentiated surface (isolation by distance or null hypothesis) for both the rangewide and Gunnison Basin extents using the PopGenReport R package (Adamack & Gruber, 2014) as described in Peterson et al. (2019). We simulated 4900 individuals (49 leks × 100 individuals) for Gunnison Basin and 5750 individuals for the rangewide analysis (6 subpopulations with 50, 200, 200, 4775, 175, and 350 individuals simulated representing true population estimates for Cimarron, Crawford, Dove Creek, Gunnison Basin, Piñon Mesa, and San Miguel, respectively; United States Fish and Wildlife Service, 2013). The true population sizes are highly variable; thus it was important to modify some functions in the R package to accommodate variable population size in our simulations (Appendix S3). For simulation parameters at both extents, we assumed a 1.6:1 female to male sex ratio (Gunnison sage‐grouse Rangewide Steering, 2005), 22 loci with 11 alleles per locus (median alleles per locus in observed data), four offspring per reproductive event (based on the product of average clutch = 8 [Young et al., 2020], average hatchability = 82.5% [Young et al., 2020], average nest success = 46.7% [Davis, 2012; Stanley et al., 2019; Young et al., 1994] and rounded up to the nearest whole number), 0.01 as the percentage of individuals moving the maximum dispersal distance (proportion moved > the 99th percentile of all max distances in Aldridge et al., 2012), an empirical estimate of migrants per year (Wilson & Rannala, 2003), and simulated data for 1000 and 500 years for the rangewide and Gunnison Basin extents respectively. We then randomly subsampled an even number of individuals from each population (n = 50) or lek (n = 15) from the simulated data to approximate the dimensions of the observed data. We calculated a genetic distance matrix based on the individual‐based Manhattan distance from the subsampled simulated data and fit the top 10 and bottom 10 ranked multicovariate models as described above. If the gene flow pattern observed in our original data was primarily shaped by the landscape variables we identified in our top model, we expected the model fitting process using simulated data from our top‐ranked model to result in a DIC rank similar to the observed data, and DIC rank of models fit to data simulated from the undifferentiated surface to result in the opposite DIC rank order. We evaluated rank similarity with a Spearman's rank correlation coefficient (ρ).

Conservation simulation

To illustrate the conservation potential of our modelled relationships we simulated management actions. We modified resistance surfaces based on the top ranked Gunnison Basin and rangewide modelled relationships. Rangewide we simulated actions via a modified surface (based on our best fit rangewide model—see Section 3) by maximizing the observed percentage of sagebrush cover (increasing to 26% in a 3000‐m window) in an artificial 708‐km2 habitat corridor connecting the Gunnison Basin population to the Piñon Mesa population by way of Crawford (Figure 2b), an unmodified surface but with an artificial new population of 200 individuals (i.e., reintroduction) between Crawford and Piñon Mesa (Figure 2c), and the modified surface including the habitat corridor and new population combined (Figure 2d). For the Gunnison Basin, we artificially increased shrub height by 50% (to ~33 cm) in a 3000‐m window over a 7.2‐km2 patch between two leks (Figure 3b), with the expectation that increased shrub height would result in greater cover of taller shrubs, and thus, increase gene flow between the two leks (based on our best fit Gunnison Basin model—see Section 3). We also modified the same area to decrease shrub height by 50% (Figure 3c). We used the above‐described simulation procedure to simulate from the modelled resistance surface and the modified versions using the PopGenReport R package with 100 independent repetitions for 1000 and 150 years post landscape modification for the population and lek comparisons, respectively. We compared the change in Manhattan genetic distance () between the targeted leks or populations over time for the modified and unmodified resistance surfaces. This exercise was intended as an example of how our modelling approach can facilitate decision making and does not make any suggestion regarding management actions.

FIGURE 2

The sagebrush cover surrounding three populations before (a) and after (b) a hypothetical management action increasing the intervening sagebrush cover by 26% in a 3000‐m window over 708‐km2 (white = less, green = more), (c) establishing a new population (♦) in existing sagebrush habitat, and the combination of a sagebrush corridor and a new population (d). The average genetic distance (y‐axis; Manhattan genetic distance) between Gunnison Basin and Piñon Mesa (e), Gunnison Basin and Crawford (f), and Piñon Mesa and Crawford (g) from simulations of gene flow based on a resistance surface developed with the original sagebrush cover raster (●), new sagebrush corridor (), and/or establishment of a new population (,) and the relationships described by the top model. Vertical lines in panels (e–f) are 95% confidence intervals for the 100 independent simulations at each time step (1–1000 years)

FIGURE 3

The shrub height surrounding two leks within the Gunnison Basin before (a) and after (b) a hypothetical management action increasing the intervening shrub height by 50% in a 3000‐m window over a 7.2‐km2 patch (white = shorter shrub, green = taller shrub), and (c) reducing shrub height in the same patch. The change in average genetic distance (y‐axis; Manhattan genetic distance) between two selected leks from simulations of gene flow among active leks in the Gunnison Basin based on a resistance surface developed with the original shrub height raster (●), increased shrub height (), and decreased shrub height () and the relationships described by the top ranked model, mimicking management actions (d). Vertical lines in panel (d) are 95% confidence intervals for the 100 independent simulations at each time step (1–150 years)

RESULTS

Rangewide model fitting

Individual pairwise Manhattan genetic distances for the rangewide analysis ranged between 8.00 and 92.00. The best fit scale (i.e., moving window radius) for each variable suggests habitat influenced gene flow primarily at a fine‐moderate scale, where the best fit radius for proportions of conifer cover (CON) and low‐statured sagebrush cover (LS) was 1000‐m, percentage of sagebrush cover (PSB) was 3000‐m, and proportion of big sagebrush cover (BS), percentage of big sagebrush cover (PBS), and shrub height (SBHT) was 6400‐m (Table 1; see Appendix Table S2.2 for complete results). Climate and phenology also influenced gene flow at the fine‐moderate scale, where mean annual rainfall (MAR), spring NDVI (SNDVI), and a dryness index (DRI) had a best fit radius of 1000‐m, mean maximum temperature (MMT) and fall NDVI (FNDVI) with 3000‐m, and growing degree days (GDD) with 6400‐m. Anthropogenic and topographic variables influenced gene flow at more coarse scales: compound topographic index (CTI), density of class 1 and 2 roads (D12, D2), proportion of agricultural cover (DAG) and development (DD), density of oil and gas wells (DOG), population density (DP), slope (S) ≥6400‐m. Most of the modelled relationships match our predicted impact on gene flow (Table 1). The few notable exceptions included negative coefficient estimates for MAR, GDD, and SNDVI, and positive coefficient estimates for DAG, DA, DD, DP, D14, D12, D2.

Of the 33 variables considered in single covariate models, nine were retained for combination in multicovariate models after considering DIC rank and correlation: proportion of conifer cover (CON), density of class 2 roads (D2), proportion of development (DD), distance to population density (DIP), density of light duty roads (DLD), density of oil and gas wells (DOG), proportion of low‐statured sagebrush cover (LS), mean maximum temperature (MMT), percent big sagebrush cover (PBS) (Figure 4; Appendix Table S2.2). We fit a total of 247 multicovariate models, the top 20 of which included a relatively large positive relationship with percent big sagebrush cover (βPBS = 2.3 [–1.9, 7.0] in top model) and a smaller but consistent negative relationship with proportion of conifer cover (βCON = –0.5 [–6.3, 3.9] in top model) indicating these 2 variables are consistently important in explaining observed variation in genetic distance (Table 2; Appendix Table S2.3). In addition to CON and PBS, the top ranked multicovariate model (DIC = −72,497.79; lpd = 3842.39) also included three variables with a small but positive coefficient estimate: D2 (βD2 = 0.5 [–4.3, 5.9]), DLD (βDLD = 0.3 [–5.0, 5.7]), DOG (βDOG = 0.2 [–5.9, 6.3]). The conductance surface created from the top modelled relationships indicated relatively low gene flow rangewide with predicted high gene flow areas in population centres, primarily in the Gunnison Basin (Figure 5a). In our evaluation of whether the conductance surface from the top modelled relationship could result in our observed genetic distances, we found distance between populations probably plays a role in addition to landscape resistance. The rank of models fit to data simulated from the top modelled resistance surface (inverse of the conductance surface) closely approximated the DIC rank of the models based on observed data (ρ = 0.90). The DIC rank of models fit to data simulated from an undifferentiated surface also approximated DIC rank of models fit to the observed data; however, the correlation coefficient was smaller (ρ = 0.83) (Appendix Table S2.6).

FIGURE 4

Rangewide single‐covariate model estimates (a) and associated deviance information criterion (DIC ) (b). Coefficient estimates and 95% credible intervals (CRI) for each covariate and associated DIC. Only the top ranked version (moving window radius, linear vs. quadratic) for each covariate was included above. Only a subset of the above covariates () were included in multicovariate model after considering correlation and DIC rank. Where appropriate, quadratic relationships were plotted just below the estimate for the linear term. Descriptions of covariates are provided in Table 1

TABLE 2

Rangewide top 10 Gunnison sage‐grouse multicovariate model parameter estimates and associated 95% credible intervals (CRI). Covariates associated with each coefficient (β) estimate are abbreviated for clarity

Model no.	βINT [95% CRI]	βCON [95% CRI]	βPBS [95% CRI]	βD2 [95% CRI]	βDLD [95% CRI]	βDOG [95% CRI]	βDD [95% CRI]	DIC	CV	DIC rank	CV rank
RW173	−0.1 [−1.7, 1.6]	−0.5 [−6.3, 3.9]	2.3 [−1.9, 7.0]	0.5 [−4.3, 5.9]	0.3 [−5.0, 5.7]	0.2 [−5.9, 6.3]		−72497.8	3842.39	1	1
RW119	0.0 [−1.7, 1.6]	−0.5 [−6.1, 3.9]	2.4 [−1.8, 7.2]	0.5 [−4.3, 5.9]	0.3 [−5.0, 5.9]			−72436.3	3842.36	2	2
RW145	0.0 [−1.3, 1.5]	−0.5 [−6.2, 3.8]	2.3 [−1.7, 7.1]	0.5 [−4.1, 5.5]		0.1 [−6.2, 6.0]		−72066.2	3842.33	8	3
RW148	0.0 [−1.4, 1.5]	−0.5 [−6.1, 4.0]	2.3 [−1.8, 7.1]	0.5 [−4.2, 5.8]			0.0 [−6.1, 6.1]	−72228.0	3842.31	6	4
RW186	0.0 [−1.4, 1.5]	−0.4 [−6.2, 4.0]	2.3 [−1.8, 7.0]	0.5 [−4.3, 5.7]		0.2 [−5.9, 6.4]	−0.0 [−6.2, 6.1]	−72037.6	3842.26	10	5
RW239	−0.1 [−1.8, 1.6]	−0.5 [−6.2, 3.9]	2.5 [−1.8, 7.3]	0.4 [−4.4, 5.6]	0.4 [−4.9, 5.9]	0.1 [−6.3, 6.1]	−0.1 [−6.3, 6.2]	−72348.3	3842.16	3	6
RW174	−0.1 [−1.7, 1.6]	−0.5 [−6.2, 4.0]	2.4 [−1.9, 7.1]	0.5 [−4.5, 5.7]	0.4 [−5.0, 6.0]		−0.1 [−6.4, 6.3]	−72269.1	3842.07	4	7
RW37	0.0 [−1.4, 1.4]	−0.4 [−6.2, 4.0]	2.3 [−1.7, 7.0]	0.5 [−4.1, 5.7]				−72064.0	3842.05	9	8
RW72	0.0 [−1.6, 1.6]	−0.6 [−5.9, 3.6]	2.4 [−1.8, 7.2]		0.5 [−4.7, 5.9]			−72144.1	3841.84	7	9
RW165	0.0 [−1.6, 1.6]	−0.6 [−6.0, 3.5]	2.5 [−1.6, 7.2]		0.4 [−4.9, 5.8]	0.1 [−6.0, 6.4]	0.0 [−6.3, 6.2]	−72229.5	3841.79	5	10
RW0	0.2 [−0.1, 0.5]							−67603.0		216

Each model has an associated deviance information criterion (DIC) estimate, DIC rank, cross‐validation score (CV), and CV rank. Estimates of τ are consistent across all models: 0.96 [0.92, 0.99].

CON, proportion of conifer cover; D2, density of class 2 roads; DD, proportion of development; DLD, density of light duty roads; DOG, density of oil and gas wells; INT, intercept; PBS, percent big sagebrush.

FIGURE 5

Posterior mean conductance (gene flow) for the top rangewide multicovariate model (a) and across the Gunnison Basin (b). Rangewide conductance (i.e., gene flow) increased with percent big sagebrush cover (PBS), density of class 2 roads (D2), light duty roads (DLD) and oil and gas wells (DOG), and decreased with proportion of conifer cover (CON). Within the Gunnison Basin, conductance increased with shrub height (SBHT), population density (DP) and high‐quality nesting habitat (N). For reference major roads, Gunnison sage‐grouse population centres (a) and the city of Gunnison (b; ) were included

Gunnison Basin model fitting

Individual pairwise Manhattan genetic distance for Gunnison Basin ranged between 5.12 and 87.13. The shrub component influenced gene flow across multiple scales with best fit radii ranging from 564‐m to 10,000‐m, while the conifer component had a best fit with a radius of 20,000‐m (Appendix Table S2.4). All but a single climate variable was best fit at a coarse scale (15,000–20,000‐m radius), the exception being MAR with a best fit to the 3000‐m radius. Phenology variables were best fit at a fine‐scale in the fall (FNDVI 3000‐m radius) and coarse‐scale in the spring (SNDVI 15,000‐m). Similar to shrub variables, anthropogenic variables influenced gene flow across nearly all possible radii. Topographic variables were all best fit when averaged over a 20,000‐m radius. As for the rangewide single covariate models, most of our predicted relationships for the effect of individual variables on gene flow were affirmed. Notably, we found negative coefficient estimates for LS, MAR, FNDVI, SNDVI, and CTI and positive coefficient estimates for MMT, DRI, GDD, S, and TRI when we predicted the opposite.

Of the 36 variables considered, seven were retained for combination after considering correlation and DIC rank: proportion of agriculture (DAG), density of class 1–4 roads (D14), proportion of low‐statured sagebrush cover (LS), population density (DP), N (nesting habitat quality), shrub height (SBHT), and terrain ruggedness index (TRI) (Figure 6; Appendix Table S2.4). We initially evaluated habitat quality variables for the nesting, brood rearing, and winter seasons; only the nesting habitat quality variable was previously published. Data for the three seasons were highly correlated and the two unpublished seasonal habitat quality variables were poor predictors and excluded for further consideration. A total of 120 multicovariate models were fit to the data. A single variable, shrub height (SBHT), was included in the top 62 DIC ranked models with a relatively large positive coefficient estimate (Appendix Table S2.5). The top ranked Gunnison Basin model (Table 3; DIC = –31,665.91, lpd = 9181.49) included shrub height (βSBHT = 2.95 [0.08, 6.22]), DP (βDP = 0.55 [–5.13, 5.83]), and nesting habitat quality (βN = 0.22 [–1.27, 1.72]). The conductance surface representing the top modelled Gunnison Basin relationship predicted high gene flow where there is a robust shrub community, with areas of low conductance near roads and significant development (Figure 5b). The coefficients for N, D14, and DP occasionally switch sign, possibly from multicollinearity. The coefficient for N is only negative when combined with DAG or D14 and a sagebrush cover variable (SBHT or LS) (Table3 and Appendix Table S2.5). The nesting habitat quality layer (N) was developed with some of the same spatial data sources we used here (Aldridge et al., 2012), notably the basin‐wide road data and Landfire data. The covariance structure is calculated through matrix decomposition of covariates into a precision matrix for the observed nodes (i.e., populations or leks). We checked variables for correlation based on the entire raster for the study area, thus leaving the possibility that the observed nodes and directly intervening landscape are more correlated than was apparent at the study area scale. For DP, the coefficient is often small and CRIs bound zero nearly evenly suggesting the true effect of DP is either very small or zero. The DIC rank of the top and bottom 10 models fit to data simulated from the resistance surface and an undifferentiated surface, clearly show the data simulated from the resistance surface resulted in a DIC rank more closely aligned to the DIC rank based on the observed data (ρ = 0.72) than the rank of models fit to data simulated from an undifferentiated surface (ρ = –0.12) and support the ability of our top ranked modelled relationships to have shaped the observed genetic distance patterns (Appendix Table S2.6).

FIGURE 6

Gunnison Basin single‐covariate model estimates (a) and associated deviance information criterion (DIC) (b). Coefficient estimates and 95% credible intervals (CRI) for each covariate and associated DIC. Only the top ranked version (moving window radius, linear versus quadratic) for each covariate was included above. Only a subset of the above covariates () were included in multicovariate model after considering correlation and DIC rank. Where appropriate, quadratic relationships were plotted just below the estimate for the linear term. Descriptions of covariates are provided in Table 1

TABLE 3

Gunnison Basin top 10 Gunnison sage‐grouse multicovariate model parameter estimates and associated 95% credible intervals (CRI). Covariates associated with each coefficient (β) estimate are abbreviated for clarity

Model no.	βINT [95% CRI]	βD14 [95% CRI]	βDP [95% CRI]	βN [95% CRI]	βSBHT [95% CRI]	βTRI [95% CRI]	DIC	CV	DIC rank	CV rank
GB51	−0.5 [−2.8, 1.5]		0.6 [−5.1, 5.8]	0.2 [−1.3, 1.7]	3.0 [0.1, 6.2]		−31665.9	9181.49	1	1
GB87	−0.2 [−2.6, 1.6]	−0.5 [−1.7, 0.5]	1.1 [−4.3, 6.2]	−0.4 [−2.4, 1.6]	3.2 [0.3, 6.5]		−31346.1	9181.19	8	2
GB26	−0.4 [−2.5, 1.4]			0.2 [−1.2, 1.7]	2.8 [0.1, 5.9]		−31572.4	9181.07	3	3
GB93	−0.8 [−3.2, 1.0]		0.4 [−5.3, 5.4]	0.2 [−1.4, 1.8]	2.6 [−0.2, 5.9]	1.9 [−0.8, 4.4]	−31402.7	9180.87	7	4
GB45	−0.2 [−2.6, 1.6]	−0.4 [−1.2, 0.5]	1.2 [−4.6, 6.3]		2.9 [0.4, 6.3]		−31564.4	9180.76	4	5
GB16	−0.5 [−2.7, 1.4]		0.4 [−5.2, 5.5]		3.0 [0.4, 6.3]		−31642.5	9180.41	2	6
GB6	−0.4 [−2.5, 1.5]				2.8 [0.2, 6.0]		−31523.0	9180.28	5	7
GB52	−0.8 [−3.1, 1.1]		0.1 [−5.6, 5.0]		2.6 [0.0, 6.0]	1.8 [−0.8, 4.3]	−31442.1	9179.94	6	8
GB27	−0.8 [−3.0, 1.0]				2.5 [0.1, 5.7]	1.9 [−0.7, 4.5]	−31343.8	9179.58	9	9
GB20	−0.1 [−2.3, 1.7]	−0.3 [−1.1, 0.5]			2.7 [0.2, 5.7]		−31314.8	9179.33	10	10
GB0	1.6 [1.3, 1.8]						−29155.2		114

Each model has an associated deviance information criterion (DIC) estimate, DIC rank, cross validation score (CV), and CV rank. Estimates of τ are consistent across all models: 0.94 [0.92, 0.97].

D14, density of road classes 1–4; DP, population density; INT, intercept; N, nesting habitat; SBHT, shrub height; TRI, terrain ruggedness index.

Conservation action simulation

Pairwise population comparison of genetic distance for simulation scenarios including an artificial sagebrush corridor were slightly lower, though not significantly, on average after 1000 years when compared to no action or a reintroduction alone (Figure 2e–g). Notably, differences in pairwise genetic distance were only observed after hundreds of years of simulation (~600 for comparisons between Gunnison Basin and Piñon Mesa [Figure 2e] or Crawford [Figure 2f] and ~200 for the comparison between Crawford and Piñon Mesa [Figure 2g]). Variation around the mean genetic distance across years was high and scenarios that included an artificial habitat corridor did not consistently result in lower genetic distance relative to no action. The observed noise in relative change to genetic distance resulting from hypothetical management actions may be a result of the significant role of distance among populations in moderating gene flow. Our simulation of hypothetical management actions in the Gunnison Basin resulted in the expected effect: increasing shrub height in the area between two leks resulted in higher gene flow than if the shrub height were to remain the same and decreasing shrub height increased genetic distance (Figure 3). Twenty years after the hypothetical increase in shrub height, genetic distance between the two leks was lower than if no action were taken (Figure 3d). After approximately 60 years post shrub height increase, a slight decline in genetic distance was apparent. When the shrub height was decreased, genetic distance started to increase relative to no action at approximately 50 years and continued to increase slightly through 150 years (Figure 3d).

DISCUSSION

We used a newly‐developed formal statistical model to identify the main drivers of gene flow and estimate the degree of impact to genetic connectivity for Gunnison sage‐grouse at two spatial scales. Our approach allowed us to identify the spatial scale at which each landscape feature was most likely to influence gene flow and characterize subtle differences in the primary drivers for long‐ and short‐distance effective dispersal events. Our results reinforce the known importance of sagebrush‐steppe habitat to the Gunnison sage‐grouse, but also showed that interpopulation connectivity was primarily influenced by the amount of sagebrush cover on the landscape while intrapopulation connectivity was primarily facilitated by the height of the shrub community, mostly composed of sagebrush species. These findings highlight that different considerations may be necessary for conserving different biological processes (inter‐ vs. intrapopulation connectivity).

The role of scale in gene flow

The landscape features influencing a species’ behavior or movement capabilities vary across spatial scales (Wiens, 1989) and may determine corridors of gene flow and spatial genetic structure. Identification of the appropriate scale of impact to gene flow is crucial to understanding the underlying environmental drivers (Cushman & Landguth, 2010) and providing conservation insight. Our approach suggested different types of landscape variables impacted connectivity at different scales when considered individually. Among leks in the Gunnison Basin, we found that most habitat and anthropogenic variables best predicted the data at smaller scales (564–6400‐m radius), while climate, topography, and other anthropogenic and habitat variables (D14, LS, CON) were more predictive at coarser scales (10,000–20,000‐m radius). Among populations rangewide all but six variables (CON, DA, DLD, DRI, LS, MAR) were best fit at moderate‐coarse scales (≥6400‐m radius). The differences in best‐fit scales for inter‐ and intrapopulation single covariate models support the role of hierarchical decisions in genetic connectivity, as previously described for Gunnison sage‐grouse habitat selection (Aldridge et al., 2012; Saher et al., 2021). Our finding of primarily fine‐scale drivers of intrapopulation gene flow within relatively contiguous sagebrush are in contrast to coarse‐scale landscape effects on sage‐grouse gene flow (Row et al., 2015) and large inter‐seasonal movements (Fedy et al., 2012; Taylor et al., 2012) within similarly contiguous habitat patches. While the two species of sage‐grouse have commonalities, Gunnison sage‐grouse distribution has been described as naturally patchy, including decidedly nonhabitat barriers such as mountain passes and deep canyons (Braun et al., 2014; Gunnison sage‐grouse Rangewide Steering Committee, 2005). The naturally patchy habitat of Gunnison sage‐grouse may force effective dispersal behavior to be dependent on successive short‐range movements within patches of favourable conditions.

The model best describing interpopulation gene flow included a positive relationship with percent big sagebrush cover and a smaller negative relationship with conifer cover (Table 2) suggesting interpopulation connectivity was primarily facilitated by relatively large, intact patches of sagebrush cover unobstructed by conifer. While small positive relationships with road and oil and gas well density (D2, DLD, DOG) were also included in the best fit model, this was probably due to the placement of these features within sagebrush habitat (e.g., low but positive correlation among sagebrush cover and these variables; Appendix Table S2.7) especially considering wildlife generally (Galpern, Manseau, & Wilson, 2012; Gerlach & Musolf, 2000), and sage‐grouse specifically (Aldridge et al., 2012; Bush et al., 2010; Green et al., 2017; Kirol et al., 2015; Saher et al., 2021), avoid anthropogenic features. Sagebrush cover within the topographically diverse distribution of Gunnison sage‐grouse coincides with areas that are less rugged and therefore more favorable sites for development. The positive correlation between sagebrush cover and some anthropogenic variables in addition to the positive coefficient estimates in our top model suggest that current densities of roads and energy extraction infrastructure do not inhibit gene flow completely, despite many documented instances of disruption to behaviour (e.g., Green et al., 2017), survival (e.gAldridge & Boyce, 2007; Holloran et al., 2010), and habitat quantity and quality (e.g Aldridge et al., 2012; Copeland et al., 2013). Much of the cited research on the negative impacts of anthropogenic development come from greater sage‐grouse whose distribution encompasses significantly more human alteration. It is possible that presence of higher densities of some anthropogenic modifications can increase gene flow precisely because they create unfavorable conditions, prompting individuals to disperse through the area despite the lack of available resources (Ribe et al., 1998; Spear et al., 2010). These relationships highlight the importance of carefully considering the biological mechanism underlying modelled relationships for conservation insight, as increased anthropogenic footprint will necessarily remove and degrade sagebrush habitat essential for this species to persist (Young et al., 2020). Generally, our findings are consistent with previous assumptions that land alteration that removes or degrades sagebrush cover has contributed to the formation of physically isolated and genetically distinct populations via reduced gene flow (Oyler‐McCance et al., 2005, 2015).

Among leks in the Gunnison Basin, our top model included a relatively large positive relationship between shrub height and gene flow, and smaller positive relationships with nesting habitat quality and human population density (Table 3). As in the rangewide model, the positive relationship of human population density with gene flow could be the result of one of the following: correlation with sagebrush habitat (Appendix Table S2.8) and a small anthropogenic footprint relative to locally available habitat with a small positive effect facilitating dispersal through unfavorable habitat, or a statistically zero effect given the large credible intervals evenly bounding zero. The inclusion of shrub height in the top 62 DIC ranked models (Appendix Table S2.5) further reinforces the primary importance of the shrub community in maintaining connectivity among leks. Protecting the remaining sagebrush habitat is considered essential to the survival of this sagebrush obligate species (Oyler‐McCance et al., 2001). We found that shrub height (not sagebrush cover) was important to intrapopulation connectivity, as was previously observed in greater sage‐grouse (Row et al., 2015). The multicovariate model including shrub height was a better fit to the data than when shrub height was replaced with sagebrush cover (DIC SBHT = −31,665.9 versus DIC PBS = −30,735.2), suggesting shrub height is capturing a pattern important to explaining gene flow that is not captured by sagebrush cover alone. However, sagebrush cover and shrub height are highly correlated in the region (Pearson's r > 0.9; Appendix Table S2.1) and any management efforts increasing the cover of tall sagebrush correspond to increased shrub height and should benefit genetic connectivity for Gunnison sage‐grouse. Additionally, we know that locally abundant and essential landscape components may not be identified as a primary influencer of gene flow if heterogeneity is low (Shirk et al., 2010; Short Bull et al., 2011). The Gunnison Basin has some of the most contiguous sagebrush habitat within the species’ range (Oyler‐McCance et al., 2001), however it is not without fragmentation. While the shrub‐steppe plant community is essential for the species, variables included in our top model suggest habitat structure and quality are important for inter‐lek connectivity within relatively contiguous habitat.

Conifer cover & gene flow

Conifer encroachment has displaced and degraded >170,000‐km2 of shrubland in western North America over the last 150 years (Miller et al., 2019; Romme et al., 2009). Thus, conifer removal has widely been suggested as a mechanism of habitat improvement for shrub‐steppe species (Miller et al., 2017), including Gunnison sage‐grouse (i.e., Doherty et al., 2018). Previous studies have shown negative population‐level impacts from conifer cover for both greater (Coates et al., 2016; Severson et al., 2017) and Gunnison sage‐grouse (Aldridge et al., 2012; Commons et al., 1999; Doherty et al., 2018; Saher et al., 2021). We expected high conifer cover to correspond to low gene flow. While we found the expected relationship in our top 20 rangewide models, conifer cover was not a significant influencer of gene flow within the Gunnison Basin. The nesting habitat quality layer included in our top Gunnison Basin model (N), however, included a negative relationship between habitat selection and conifer cover (Aldridge et al., 2012) suggesting increased conifer cover may play a role in influencing gene flow indirectly via nesting habitat quality. Although conifer cover has always been part of the naturally fragmented landscape encompassing the Gunnison sage‐grouse distribution, encroachment has contributed to the displacement of sagebrush cover in recent decades (Bukowski & Baker, 2013). Productive, early‐phase woodland sites, such as those formed by encroachment, may be attractive to sage‐grouse despite the potential negative effects on vital rates, and may act as an ecological trap (Coates et al., 2016). We could not evaluate intra‐lek connectivity for all satellite populations due to data limitations. Region specific impacts to gene flow from landscape components are not uncommon (Kozakiewicz et al., 2019; Row et al., 2018). Conifer encroachment and known presence of pinyon‐juniper woodlands in some satellite populations may be important to consider in the future given the known hazards.

Conservation applications

Development of conservation planning tools is a key goal of many landscape genetic analyses (Keller et al., 2015). The relationships identified in our findings could provide insight into areas of high connectivity that might be considered dispersal corridors, as well as areas that have low gene flow that might benefit from some type of habitat improvement or restoration. We identified an area in the Gunnison Basin where our top model predicted low gene flow and evaluated alternative scenarios for impact to genetic distance if shrub height increased, shrub height decreased, or if no action was taken (Figure 3a–c). The predicted effects of our hypothetical management actions indicated that if the shrub height (i.e., tall sagebrush cover in this area) across the selected 7.2 km2 area between the two leks was increased by 50% we should expect genetic divergence to continue for 20 years, remain approximately stable through 75 years and then begin to decrease (Figure 3d). Changes in genetic distance occur over multiple generations and lag behind changes in population dynamics, thus we expect genetic distance to continue increasing for a period of time post modification as we observed. If no action is taken, we expect genetic distance to increase and stabilize after 50 years. If shrub height was decreased by 50%, we expect genetic distance to increase for longer than if no action was taken for 150 years. We found a similar, although less pronounced, result when we simulated restoration of a sagebrush corridor between populations and the establishment of a new population; after ~200–600 years we might expect to see genetic divergence decrease among the targeted populations (Figure 2e–g). Notably, we did not see an effect from the establishment of a new population in the absence of the sagebrush corridor across the 1000 simulated years, variation in mean change across years was highly variable among different scenarios. Our example landscape alterations were arbitrarily chosen and depict a uniform increase in shrub height (between leks) or the establishment of dense sagebrush cover (among populations) over relatively large areas; conditions that may not be realistic for management actions but are useful for illustrative purposes. Predicted changes to the landscape from further anthropogenic alteration, climate change, or potential management actions or conservation could be used to evaluate a predicted response in connectivity. Similar approaches have been utilized for sage‐grouse to prioritize habitat for conservation or potential management actions based on population size (Heinrichs et al., 2017) and lek persistence (Doherty et al., 2018); we showed how those concepts could be extended to a gene flow scenario.

CONCLUSIONS

The ability to explicitly consider scale at every step of our landscape genetic analysis and estimate the degree of influence on gene flow are important steps to making conservation relevant inference (Keller et al., 2015). Similarly, the ability to compare relative impact of land‐use change or conservation action underscores the usefulness of the modelled relationships for decision making. Our comparison of inter‐ and intrapopulation genetic connectivity highlights the necessity of understanding scale‐dependent relationships between environment and gene flow. While sagebrush cover was important for Gunnison sage‐grouse regardless of scale, we found that large patches of dense sagebrush cover free from conifer cover was important for gene flow among populations, which contrasts with the importance of tall shrubs (primarily composed of sagebrush species) and high‐quality nesting habitat among leks within the Gunnison Basin (Table 3 and Appendix Table S2.6). Like the indirect effect of conifer cover on gene flow via inclusion in the nesting habitat quality layer (N) in our top Gunnison Basin model, sagebrush cover had a positive relationship with nesting habitat selection supporting indirect effects of sagebrush cover to gene flow (Aldridge et al., 2012). Altogether this suggests the structure of the sagebrush habitat (i.e., shrub height) might be important within a population, while longer effective dispersal events appear to be governed more by the presence or absence of sagebrush habitat in general. Additionally, we found that landscape features observed to influence behaviour or movement for wildlife species generally (e.g., roads, agriculture), and including for sage‐grouse, were not important direct influencers of Gunnison sage‐grouse genetic connectivity, underscoring the potential follies of relying on direct observations as an indicator of gene flow. Future analyses could explore the impacts of lag effects in genetic signal from landscape change (Epps & Keyghobadi, 2015), relaxing the assumption of symmetric gene flow (Hanks, 2017; Peterson et al., 2019), allowing the coefficient to vary across the landscape using splines (Hanks & Hooten, 2013), considering the influence of natural selection (Peterson et al., 2019) or ways to incorporate behavior or fitness measures (Spear et al., 2010) to further refine our understanding of drivers of gene flow. Understanding and characterizing the primary drivers of gene flow can facilitate prediction of the impact of multiple conservation relevant scenarios and decision making.

AUTHOR CONTRIBUTION

Shawna J. Zimmerman, Sara J. Oyler‐McCance, Cameron L. Aldridge, and Mevin B. Hooten designed the study; Shawna J. Zimmerman analysed the data and wrote the first manuscript draft. All authors contributed to the final version of the manuscript.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.