Signatures of adaptive divergence among populations of an avian species of conservation concern

Abstract Understanding the genetic underpinning of adaptive divergence among populations is a key goal of evolutionary biology and conservation. Gunnison sage‐grouse (Centrocercus minimus) is a sagebrush obligate species with a constricted range consisting of seven discrete populations, each with distinctly different habitat and climatic conditions. Though geographically close, populations have low levels of natural gene flow resulting in relatively high levels of differentiation. Here, we use 15,033 SNP loci in genomic outlier analyses, genotype–environment association analyses, and gene ontology enrichment tests to examine patterns of putatively adaptive genetic differentiation in an avian species of conservation concern. We found 411 loci within 5 kbp of 289 putative genes associated with biological functions or pathways that were overrepresented in the assemblage of outlier SNPs. The identified gene set was enriched for cytochrome P450 gene family members (CYP4V2, CYP2R1, CYP2C23B, CYP4B1) and could impact metabolism of plant secondary metabolites, a critical challenge for sagebrush obligates. Additionally, the gene set was also enriched with members potentially involved in antiviral response (DEAD box helicase gene family and SETX). Our results provide a first look at local adaption for isolated populations of a single species and suggest adaptive divergence in multiple metabolic and biochemical pathways may be occurring. This information can be useful in managing this species of conservation concern, for example, to identify unique populations to conserve, avoid translocation or release of individuals that may swamp locally adapted genetic diversity, or guide habitat restoration efforts.

Byrne, 2016). Additionally, adaptive variation could inform whether augmentation should be done at all (Benedict, Oyler-McCance, Braun, & Quinn, 2003), guide development of captive breeding programs (Williams & Hoffman, 2009), aid in monitoring and maintaining locally adapted variation in populations, or be used to identify evolutionarily significant units (ESUs; Funk, McKay, Hohenlohe, & Allendorf, 2012). While diversity at putatively neutral genetic markers has long been used to characterize populations, advances in DNA sequencing technology (Mardis, 2008;Metzker, 2010;Shendure & Ji, 2008) and methods to separate neutral and functional genetic variation (Allendorf, Hohenlohe, & Luikart, 2010) have facilitated a shift in focus to understanding the role genetic diversity plays in adaptation to local environments (Nielsen, 2005;Schweizer et al., 2016;Wenzel & Piertney, 2015;De Wit & Palumbi, 2013). Genomic methods can be particularly valuable for characterizing adaptive divergence in species where traditional approaches to evaluate local adaptation (i.e., reciprocal transplant experiments) are not feasible, such as with federally protected species (Funk et al., 2012).
The Gunnison sage-grouse (Centrocercus minimus) is a sagebrush (Artemisia spp.) obligate avian species persisting as seven isolated populations with low gene flow and high genetic differentiation (Oyler-McCance, St John, Taylor, Apa, & Quinn, 2005). A single population, the Gunnison Basin, supports the majority of the species (~85%-90% of ~5,000 individuals) with the remaining birds residing in smaller satellite populations (United States Fish & Wildlife Service, 2014). Historically, Gunnison sage-grouse occurred across ~46,521 km 2 of sagebrush habitat in Colorado, Utah, New Mexico, and Arizona (Schroeder et al., 2004). Land-use change in sagebrush habitat has reduced the species to just 8% of the historical range with birds remaining only in southwestern Colorado and southeastern Utah (Figure 1; Braun et al., 2014;Schroeder et al., 2004). In 2014, the species was listed as threatened under the Endangered Species Act (United States Fish & Wildlife Service, 2014). As a sagebrush obligate, Gunnison sage-grouse requires sagebrush cover for habitat during all life stages (Patterson, 1952;Wallestad & Eng, 1975), and as a source of forage, with up to 99% of winter diet consisting of F I G U R E 1 Historical (gray) and current (yellow) distribution of Gunnison sage-grouse in the southwestern United States. Populations labeled with respective names. Black rectangle designates the study area. The historic range map is as described by Braun et al. (2014); the two northernmost portions of the historic range correspond to an unknown species of sage-grouse and are not verified by Colorado Parks and Wildlife (Gunnison sage-grouse Rangewide Steering Committee, 2005). Sample locations are indicated on the map as point of variable size, scale by number of samples collected at the location sagebrush leaves (Braun, Britt, & Wallestad, 1977;Braun, Connelly, & Schroeder, 2005;Young, 1994). Differences in local population en- Although populations are close in proximity (33.34 to 203.72 km apart) relative to observed dispersal capabilities (up to 120-240 km for greater sage-grouse in mostly contiguous habitat; Cross, Naugle, Carlson, & Schwartz, 2017;Newton et al., 2017;Tack, Naugle, Carlson, & Fargey, 2011), genetic differentiation between populations is relatively high (Oyler-McCance et al., 2005), suggesting low levels of homogenizing gene flow which might otherwise limit local adaptation. Conversely, some gene flow can increase the local genetic variation in a population and therefore provide more opportunities for natural selection to result in local adaptation (Lenormand, 2002;Whiteley, Fitzpatrick, Funk, & Tallmon, 2015), which suggests observed low levels of gene may promote local adaptation in the different local habitat patches. The male-dominant polygynous mating system of sagegrouse skews mating success among males (Wiley, 1973;Young, Braun, Oyler-McCance, Hupp, & Quinn, 2000) and imposes strong sexual selection which could lead to rapid morphological and/or behavioral changes and further divergence among isolated groups (Ellsworth, Honeycutt, & Silvy, 1995;Oyler-McCance, St. John, & Quinn, 2010;Spaulding, 2007;Uy & Borgia, 2000). The skew in mating success decreases effective population size (Stiver, Apa, Remington, & Gibson, 2008). This mating skew, along with small population size, also indicates genetic drift could overwhelm the efficacy of selection for local adaptation.
Previous studies have found evidence for significant genetic divergence within some sage-grouse populations. Isolated populations of greater sage-grouse (C. urophasianus) are genetically distinct enough at neutral loci to warrant consideration for special protection (Benedict et al., 2003;Oh, Aldridge, Forbey, Dadabay, & Oyler-McCance, 2019). An evaluation of genetic variation at cytochrome P450 genes and additional candidate genes related to metabolism of plant secondary metabolites (PSMs) in greater sage-grouse identified evidence for positive selection, potentially pointing to local dietary adaptation . The cytochrome P450 superfamily of genes have broad roles in physiological and toxicological processes (Kubota et al., 2011). Importantly, some of the members of this gene family are involved in metabolism of PSMs (Miyazawa, Shindo, & Shimada, 2001;Skopec, Malenke, Halpert, & Denise Dearing, 2013), like the monoterpenes, sesquiterpene lactones, and phenolics found in sagebrush species (Kelsey, Stephens, & Shafizadeh, 1982

| Study system
Our study area encompassed the entire species range excluding the eastern most population, Poncha Pass ( Figure 1). The Poncha Pass population is thought to have been extirpated in the 1950s and reestablished with Gunnison Basin individuals beginning in the 1970s, persisting only due to ongoing translocation (Nehring & Braun, 2000). For these reasons, the Poncha Pass population was excluded from our analyses.

| Library preparation
We accomplished SNP identification using an adapted version of the ddRAD protocol as first described by Peterson, Weber, Kay, Fisher, and Hoekstra (2012

| Sequence data processing and genotyping
Raw sequencing reads were trimmed at a maximum error probability of 0.05 using CLC Genomics v. 9.5 (Qiagen), allowing at most two ambiguous bases. Reads were mapped to the Cmin_1.0 Gunnison sage-grouse draft genome assembly ; GenBank accession: SPOS00000000) with Bowtie2 (Langmead & Salzberg, 2012) using the "very-sensitive" and "end-to-end" parameter sets and filtered on a mapping quality of 20 (Phred-scaled) with the samtools/bcftools package, v. 1.3 (Li et al., 2009). Potential PCR duplicates were removed by processing unique molecular identifiers with the UMI-tools package (Smith, Heger, & Sudbery, 2017), using the "unique" identifier detection algorithm.
The samtools/bcftools package was used to merge alignments and identify variant sites in the reference genome. Base composition at sites was computed with the mpileup function using the recommended map-quality adjustment ("-C" set to 50) and base-alignment qualities recalculated from the combined data. Indels were called for the purposes of filtering nearby SNP sites (within 3 bp) that could be affected by local misalignment, but were not otherwise used. Genotype likelihoods were estimated with the bcftools call function using the multiallelic model, although only biallelic loci were retained. SNP loci were further filtered by requiring a minimum coverage of 960× across all individuals (based on an average of 15X per sequenced individual) and called genotypes for at least 50 of the total 64 individuals. Loci potentially located on sex chromosomes were removed using both coverage and homology information: SNPs on scaffolds putatively homologous with the sex chromosomes of Gallus gallus  were excluded, as were SNPs with unequal coverage in males and females (i.e., if the ratio of male to female mean coverage was outside the range 0.9-1.2). We excluded potentially sex-linked loci because the proportion of each sex sampled in each population was variable and we wanted to reduce the likelihood of false positives for adaptive divergence due to sampling bias at sex-linked loci. Sites with low-frequency minor alleles (below 5%) were also excluded. In addition to filtering variant sites based on these locus-based criteria, individual genotype calls were removed if coverage was less than 10X for an individual at a given location, regardless of whether a genotype was called by bcftools. After excluding four of the 64 sequenced individuals due to low coverage overall, the final data set included 15,033 loci across 35 "pseudo-chromosomes" (chromosome scaffolds inferred from synteny with chicken) for 60 individuals (four Cimarron, 12 Crawford, 12 Dove Creek, 12 Gunnison Basin, 10 Pinon Mesa, and 10 San Miguel). Because sample size in Cimarron was low, the power to detect unique outliers in this population was also low. However, inclusion of the Cimarron samples would still help estimate population structure and identify global outliers.

| Outlier locus analyses
We identified outlier SNPs using the BayPass core model (Gautier, 2015) and pcadapt (Luu, Bazin, & Blum, 2017) to balance the tradeoff between power and false-positive rates observed when using parametric (as in BayPass) versus nonparametric (as in pcadapt) outlier analysis approaches. Both approaches perform well under high population structure demographic scenarios, similar to that of Gunnison sage-grouse (Gautier, 2015;Luu et al., 2017). The core BayPass model expands on the approach implemented in BAYENV (Coop, Witonsky, Di Rienzo, & Pritchard, 2010;Gunther & Coop, 2013) by providing greater computational efficiency, flexibility, and a formal procedure for calculating outlier thresholds. As with BAYENV, this method incorporates a scaled covariance matrix accounting for background population structure which can confound analyses for adaptive variation (Meirmans, 2012). The population covariance matrix was directly estimated with the core model, the Inverse-Wishart prior set to 1, and both hyperpriors for beta (a.pi, b.pi) set to 1. Five thousand MCMC iterations were performed after discarding a 5,000-iteration burn-in and thinning by a factor of 25.
Twenty pilot runs of 1,000 MCMC iterations were performed in order to adjust the parameters in the proposal distribution of the Metropolis-Hastings algorithm so that an acceptance rate between 0.25 and 0.40 was achieved. The adjustment parameter (set to 1.25) was used in the pilot runs to adjust the range of possible values from the proposal distributions if the acceptance rate fell outside the desired window. Inference on the core model was through estimates of the XtX statistic. Allele counts were simulated with the simulate.
baypass R function and the population covariance matrix to generate a pseudo-observed data set from the core model. Outliers were loci with XtX values exceeding the 99th quantile of the XtX distribution that resulted from the simulated pseudo-observed data set (false discovery rate (FDR) of 0.01). We verified that the scaled covariance matrix of population allele frequencies estimated from the simulated data was close to the matrix estimated from our data (FMD distance = 0.19, see Gautier, 2015). For pcadapt, we used the pcadapt R package (Luu et al., 2017) with K = 5 based on visual inspection of a scree plot ( Figure S3.1). The first 5 principal components generally separate individuals into loose populations ( Figure S3.2). We considered loci with a q-value < 0.01 as outliers.
As an independent evaluation of the SNP data and the ability of BayPass to control for population structure, we compared pairwise F ST (as in Weir & Cockerham, 1984) as calculated with the R package diveRsity (Keenan, McGinnity, Cross, Crozier, & Prodohl, 2013) for the neutral SNPs (the SNP data with all outliers removed) to values previously obtained using a microsatellite genotype data set created from the whole available sample set. If SNPs showed the same pattern of differentiation as microsatellite loci, we were confident in the ability of the SNP data and BayPass to estimate population structure ( Figure S4).

| Genotype-environment association analyses
Environmental covariates tested for association with SNPs were A total of 72 covariates were initially considered (Table S1.1). We reduced this set to eight variables with a Pearson correlation coefficient <|0.70| for our analyses. This reduced list included spring and fall precipitation, spring maximum temperature, winter vapor pressure deficit (i.e., evapotranspiration), compound topographic index (CTI; a wetness index), green-up rate, big sagebrush cover, and a dryness index. We also used the loadings of the first three principal components (PCs; PC1 = 37.59%, PC2 = 29.53%, PC3 = 18.85% of the variance) of the eight minimally correlated variables as a covariate in attempt to incorporate multiple covariates in a single model.
The principal components analysis was performed with the prcomp function in R (see Appendix S1 in Supporting Information for full details on covariates).
Correlation of environmental covariates with SNP genotypes was evaluated using the standard covariate model in BayPass. The model was implemented as in the core model, though including the population covariance matrix estimated with the core model and the addition of regression coefficients which had a uniform prior bounded between −0.3 and 0.3. Covariates with an empirical Bayesian p-value (eBPmc) greater than 4 were considered significantly associated. Gautier (2015) recommends an eBPmc of 3 as a threshold for candidacy; however, we used 4 as a threshold to further control for false discovery. We also used a partial redundancy analysis (RDA), an approach based on a combination of multivariate linear regression and principal components analysis (PCA) that can identify SNPs weakly associated with environmental covariates (Forester, Lasky, Wagner, & Urban, 2018;Rellstab, Gugerli, Eckert, Hancock, & Holderegger, 2015). RDA is a nonparametric approach to identifying candidate adaptive loci which has high power and low false-positive rates, yet has not been evaluated under a high population structure demographic scenario as observed in Gunnison sage-grouse. Further, evaluation of a partial RDA (accounting for demographic structure) on a low-structure system resulted in reduced power and increased false-positive rates (Forester et al., 2018). However, results of a recent application of a partial RDA to a high-structure system (global Brauer, Hammer, & Beheregaray, 2016) were consistent with an independent transcriptomics study (Brauer, Unmack, & Beheregaray, 2017), suggesting the partial RDA performs well in high-structure systems. We therefore report the results of a partial RDA as a complement to the BayPass standard covariate model. For the partial RDA, the genetic data were recoded to a 0,1,2 format, where individuals heterozygous at a locus were designated as 1, and homozygotes were designated as 0 and 2 for the reference and alternate alleles, respectively. We accomplished the partial RDA by calculating a genetic distance matrix (vegdist function using the "bray" method in the vegan R package; Oksanen et al., 2017) and identifying the significant genetic axes of a principal coordinates analysis on the genetic distance matrix (pcoa function in the ape R package; Paradis & Schlier, 2018); significance of genetic axes was based on the broken-stick criterion (Legendre & Legendre, 2012). We then created spatial predictors by calculating Moran eigenvector maps (MEM; based on a Gabriel graph neighborhood and inverse distance weights) in the R packages spdep (Bivand & Piras, 2015) and adespatial (Dray et al., 2019).
We used the first 4 genetic axes (as determined by the broken-  Figure S5 for corresponding scree plot). Loci in the tails of the distribution of the SNP loadings on each axis were considered outliers. In attempt to keep false positives low, we used a two-tailed p-value of 0.0027 (based on 3 standard deviations) as a cut off for candidacy.

| Linkage disequilibrium and gene ontology enrichment analyses
We phased our SNPs and estimated linkage disequilibrium (LD) decay to determine at what distance candidate loci could be considered physically linked to a putative gene region on the reference genome based on the position of SNPs on pseudo-chromosomes. To phase our SNPs we used BEAGLE 5.0, setting N E to 1,000 as recommended to indicate our data are from a small and inbred population (Browning & Browning, 2007). With the phased SNPs, we calculated LD in vcftools (-hap-r2 command) at multiple distances, from SNPs 10 bp to 1Mbp apart. We considered SNPs at the distance where LD as measured by r 2~0 .10 to be physically linked.
We then further investigated the relationships between putative gene products using a gene ontology (GO) enrichment analysis. We used Gowinda v1.12 (Kofler & Schlötterer, 2012) (Wickham, 2009). For comparison, we also included plots of the first three principal components for analyses on all SNPs and putatively neutral SNPs.

| Population genetic structure check
Pairwise multilocus F ST values as calculated from only putatively neutral SNPs corresponded well by rank with previous microsatellite estimates, although the latter appear to have consistently lower means ( Figure S4). The population covariance matrix inferred from neutral SNPs also confirms our previous understanding of population structure. We have included a heat-plot of the correlation matrix derived from the allele frequency covariance matrix estimated in the BayPass program to illustrate how populations are related (Figure 2).

| Genome scans for adaptive divergence and association with environmental variables
The BayPass core model identified 76 outlier loci located on 13 of the 35 pseudo-chromosomes that had SNPs, and pcadapt identified 157 outlier loci located on 17 pseudo-chromosomes (  Figure S6). Overlap of loci identified with each analysis varied; no two analyses identified identical lists (Table S7).

| Gene ontology enrichment analyses
Eight unique genes from gene sets associated with 51 GO terms (FDR 0.05) were identified in a gene ontology enrichment analysis of outliers identified by the BayPass core model and 15 unique genes from gene sets associated with 161 GO terms (FDR 0.05) were identified with outliers from pcadapt (  Table 2 for summarization of all tests).
Global LD was estimated to be 0.02, and dropped to r 2 < 0.1 at ~ 350 kbp ( Figure S8), a distance much greater than the default distance we used to align candidate SNPS to gene regions (5 kbp). There are likely many genes within the 350 kbp LD blocks, any of which could be the true target of selection. Restricting the distance between candidate SNPs and gene regions increased our confidence that the association between candidate SNP and gene was nonrandom although it ensures that many potential targets of selection would not be included.  (59) by SnpEff (Table 3 and Table S10). Additionally, three SNPs in three putative genes were indicated as nonsynonymous variants and moderate or high impacts to putative gene function ( 1.42), ANK sequence repeat (cluster 2 enrichment score: 1.14), cytochrome P450, E-class, group 1 (cluster 3 enrichment score: 1.08), and short sequence motif: DEAD box (cluster 4 enrichment score: 0.90) ( Table 4; see Table S11 for a complete list).
Individuals generally clustered by population when candidate loci were used in a PCA (Figure 4c)   Identification of signals of adaptive divergence in Gunnison sagegrouse populations also provides more evidence of natural selection occurring in unexpected situations. First, effective population TA B L E 3 Summary of the outlier loci from Gunnison sage-grouse populations in enriched families or proteins. Findings from pcadapt, the core BayPass model ("core"), the standard BayPass model including principal component 1 ("PC1"), spring maximum temperature ("Spring Tmax."), green-up rate, dryness index, and RDA with associated predictor variable ("RDA: Spring PPT," "RDA: Fall PPT.," "RDA: CTI"). The gene code is listed in the left-hand column (see Table S9 for a list of the corresponding full gene names and Table S10 for all putative adaptive genes) followed by the pseudo-chromosome number where it is located ("Chromosome"), the total number of SNPs identified as outliers in each gene region ("# SNPs"), indication of significance at Large effective population sizes are less influenced by genetic drift and therefore natural selection is expected to be more efficient (Frankham, 1996 (Rousset, 1997;Slatkin, 1993). At the microgeographic scale (when geographic distances between populations are within the known physical dispersal range of an organism), high gene flow is expected to impede local adaptation (Slatkin, 1987), although some argue microgeographic local adaptation is more common than previously appreciated (Richardson, Brady, Wang, & Spear, 2016). Though not a perfect system to eval-  (Karlin & McGregor, 1972;Levene, 1953;Urban et al., 2017), a physical limit to dispersal (Fischer & Lindenmayer, 2007;Slatkin, 1987), mating signal divergence (e.g., Langin et al., 2015;Langin et al., 2017), or that any type of assortative mating may be facilitating natural selection.

| Population-level divergence
Across all candidate loci three Gunnison sage-grouse popula-  Table S11 for complete list of functional annotation clusters and corresponding terms. See Table S9

| Ecological importance of identified signals of selection
One of the top enrichment clusters included terms indicating detoxification (oxidoreductase activity) as a biological process potentially underlying adaptive divergence and identified the same four cytochrome P450 family genes in many of the gene sets (Table 4) (Frye, Connelly, Musil, & Forbey, 2013;Kelsey et al., 1982) and divergence at genes involved in PSM metabolism may reflect local adaptation to consuming different species or subspecies of sagebrush.
Because sage-grouse have mechanisms to mitigate inhibitory action of PSM on digestive enzymes (Kohl, Connelly, Dearing, & Forbey, 2016), these genes could potentially be responsible for proteins or enzymes that aid in sagebrush digestion. The candidate SNPs associated with all cytochrome P450 gene regions were identified with one or more of the environmental association analyses: CYP4V2 with green-up rate, CYP2R1 with the dryness index and fall precipitation, CYP4B1 with fall precipitation, and CYP2C23B with CTI, respectively (Table 3).
The fourth enrichment cluster indentified by DAVID (Table 4) contained gene sets dominated by members of the DEAD box helicase gene family, generally known to function within multiprotein  (Linder & Jankowsky, 2011). Of particular interest to our findings are the members of this gene family known to play a role in detecting viral RNA in the cytoplasm of chicken (Schoggins et al., 2011;Zhang et al., 2016). SETX (Table S10), a gene associated with candidate adaptive loci, has been implicated in response to viral pathogens as well, including West Nile virus in chicken (WNV; Miller et al., 2015). Signals of diversifying selection associated with putative genes involved in antiviral activity could indicate the populations may have had different exposure histories which may result in differing abilities to respond to viral pathogens. Though it has yet to affect Gunnison sage-grouse specifically, WNV has impacted susceptible greater sage-grouse populations (Naugle et al., 2004) and the virus has been reported in other species within the Gunnison sage-grouse range at varying levels (see Table S13 for information on reported WNV incidence in populations), suggesting a potential for exposure.
The reduced representation approach used here allowed us to break the entire genome down into smaller pieces and obtain higher confidence genotypes for more individuals than we would have been able to obtain with whole-genome resequencing. However, this resulted in a low density of SNPs (~16 SNPs/Mb), and many regions of the genome were not sampled (Tiffin & Ross-Ibarra, 2014).
Additionally, our use of a threshold for linkage between SNPs and gene regions (5 kbp) was much lower than LD blocks (~350 kbp) that contain multiple gene regions. Consequently, candidate SNPs are likely linked to more than one gene region, any of which could be the target of selection. Therefore, it is likely there are more regions of the genome under adaptive divergence and more processes involved.

| Conservation implications
We have identified signals of adaptive divergence associated with potentially ecologically important genes and groups of genes which may underlie adaptive divergence among populations of Gunnison sage-grouse. Populations with different functional genetic variants could potentially impact management and conservation decisions (Savolainen, Lascoux, & Merilä, 2013). Theoretically, gene flow can have either a positive or negative impact on local adaptation of populations (Slatkin, 1987;Wright, 1931). If populations are locally adapted, increasing gene flow could risk outbreeding depression (Edmands, 2007), especially if populations are small. This has been exemplified in populations of streamside salamanders with and without predators where gene flow constrained the evolution of effective antipredator behaviors (Storfer & Sih, 1998 (Frankham, 2015;Miller, Poissant, Hogg, & Coltman, 2012), the new genetic variants may allow the population to respond to the local environmental conditions and potentially occupy new niche space (Aitken & Whitlock, 2013;Lenormand, 2002). However, there is still much to understand about the relative contributions of gene flow and natural selection to local adaptation (Kawecki & Ebert, 2004 (Mahalovich & McArthur, 2004), which may result in planting a species or subspecies to which the local population is maladapted. Although matching the local sagebrush type during restoration could be important, efforts to do so could be complicated because seed sources for different sagebrush species or subspecies are not always available and factors involved in establishment of seedlings are just starting to be understood (Brabec, Germino, & Richardson, 2016).
Captive-rearing of sage-grouse has been attempted in recent years (Apa & Wiechman, 2015). Knowledge of adaptive differences could guide selection of targeted populations for release of captivereared birds. In the case of sagebrush digestion or disease response, releasing individuals with maladapted genotypes could not only result in wasted effort and resources, but may even lead to further reduction of average population fitness.
In conserving species with fragmented ranges and declining populations, restoration of gene flow between isolated groups is a common objective. Our findings suggest increasing gene flow between Gunnison sage-grouse may require careful consideration of local adaptation. On the other hand, locally adapted variation might persist in the face of gene flow (Fitzpatrick, Gerberich, Kronenberger, Angeloni, & Funk, 2015) and the existence of adaptive environmental clines suggests gene flow via assisted migration can facilitate adaptive responses to climate change (Kelly & Phillips, 2016). We would be remiss to not acknowledge the potential for false positives in our analyses, however. Our outlier analysis methods generally control for demographic structure (i.e., incorporation of a kinship matrix or nonparametric approaches), though observed differentiation could still be a result of background selection, or linkage of neutrally evolving sites to sites under purifying selection (Shafer et al., 2015).

| CON CLUS ION
Our results are consistent with the hypothesis of adaptive divergence among populations of Gunnison sage-grouse for potentially ecologically important metabolic phenotypes. This study takes the first step in understanding and characterizing local adaptation within populations of Gunnison sage-grouse. The correlative approach we used assumes high-frequency alleles in a population correspond to a higher fitness phenotype locally. This relationship could be confirmed or further probed through genomic methods that more directly evaluate fitness effects and function (Carneiro et al., 2014;Prasad et al., 2013). We used historical samples and publicly available geospatial data sets. More insight from our historical samples could be obtained by using the gene families as the subject of resequencing, or target enrichment, to identify functional variants supporting a putative role in adaptation and confirming signals of selection on a larger sample size (Jones & Good, 2016).
Many approaches used to draw more direct lines between the underlying genetic controls and phenotype, such as quantitative trait analysis (Kearsey, 1998), and gene expression and/or reciprocal transplant studies (Kawecki & Ebert, 2004), may be attractive options to provide a phenotype link, especially given that many loci of varying effect size underlie adaptive divergence (Rockman, 2012).
However, these strategies are unlikely feasible due to difficulty in generating large segregating populations in captivity and given federal protection of the species under the Endangered Species Act. Genome-wide association studies (GWAS), on the other hand, can also identify genetic regions underlying phenotypes and can be accomplished without the use of captive populations making it a much more likely approach for future studies investigating local adaptation in Gunnison sage-grouse. Nevertheless, our results have provided many avenues for future investigations of adaptation for this avian species of conservation concern.

ACK N OWLED G EM ENTS
We are grateful to those at Colorado Parks and Wildlife who helped in sample collection, especially Dr. Anthony D. Apa. We would like to thank those who helped in the review process: W. Chris Funk and three anonymous reviewers. We are also grateful to our funders: U.S. Geological Survey, U.S. Fish and Wildlife Service, and Colorado Parks and Wildlife. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare. contributed new reagents or analytical tools. S.J.Z wrote the manuscript and analyzed the data. All authors contributed to manuscript editing.

DATA ACCE SS I B I LIT Y
Genomic sequencing data for this study have been deposited in GenBank (biosample accession numbers: SAMN1084489-SAMN10844548; bioproject number: PRJNA517770). SNP genotypes for this study were deposited in the U.S. Geological Survey