BACKGROUND
population genetic studies are essential for understanding how behavioral, ecological, and evolutionary processes shape the geographic distribution of genetic variation.
recent dna sequencing advances now allow variation to be readily assayed at thousands of loci for large numbers of individuals and populations in non-model organisms  <cit> , potentially illuminating the relationships among ecological, genetic, and geographic variation at previously unattainable scales .
these approaches have value for uncovering previously unrecognized patterns of genetic structure, generating more precise estimates of demographic and evolutionary parameters, and providing increased opportunity for understanding adaptation .
such approaches also have special promise for analyses of genetic variation related to contemporary ecological processes and population demography in species of ecological or conservation significance  <cit> .
variation in behavioral and life history traits may have important consequences for the geographic structuring of genetic variation across populations .
for example, a number of avian species have lek-based mating systems in which a small number of males contribute disproportionately to reproduction, dispersal is restricted by strong site fidelity, and leks  remain at the same geographic locations for long periods of time .
strong site fidelity could reduce gene flow among leks and give rise to spatial genetic structure even in the absence of geographic barriers.
because a subset of males typically contribute disproportionately to reproduction on leks, reproductive skew  can also lower effective population size  and increase genetic drift  <cit> .
these factors, in addition to the potential for kin selection in lek-based mating systems  <cit> , could enhance the effects of drift in structured populations, give rise to fine-scale spatial genetic structure, and increase extinction probability of isolated populations .
greater sage-grouse  are an iconic western north american bird species with a lek-based mating system  <cit> .
male sage-grouse display for females at leks residing in open patches of sagebrush , where a subset of dominant males can account for a large proportion of the successful matings .
degradation of sagebrush ecosystems across north america has substantially reduced the geographic range and abundance of sage-grouse over the last century.
population declines ranged from  <dig> to 90 % across the distribution and local extirpations have occurred in many regions  <cit> .
habitat degradation and population declines  <cit>  have raised concerns about population persistence , though sage-grouse were recently found unwarranted for listing under the u.s. endangered species act  <cit> .
low movement rates among leks due to high site fidelity , in addition to reproductive skew, could interact with habitat variation to limit gene flow, enhance genetic drift, and shape fine-scale genetic structure in sage-grouse populations, as has been found in several lek-breeding species .
previous studies utilizing small data sets of traditional molecular markers  have reported variable patterns of genetic structure across sage-grouse populations, leaving uncertainty about the scale of geographic genetic structure and obscuring the expected link between demographic and population genetic parameters across the range  <cit> .
the detection of genetic differentiation across broad geographic scales resulted in the recognition of three distinct populations: gunnison sage-grouse , greater sage-grouse, and the bi-state population of greater sage-grouse  .
however, most studies at finer geographic scales have not identified genetic structure among geographically proximate leks  <cit> , leaving uncertainty about the potential for site fidelity and reproductive skew to influence genetic structure  <cit> .
in contrast, a recent study of geographically proximate leks in northwestern colorado documented among-lek differentiation and evidence for isolation by distance  <cit> .
although these studies suggest movements among leks may limit differentiation in some parts of the range, higher resolution data sets should improve our understanding of the geographic scale of genetic structure and the factors associated with it.
here we use a genotyping-by-sequencing  approach to quantify genetic variation at  <dig>  single-nucleotide polymorphisms  across  <dig> sage-grouse from ten leks located throughout a small geographic region in central nevada .
very strong breeding-site fidelity has been documented for these specific leks across multiple years of monitoring, where  <dig> of  <dig>  among-year observations involved males returning to the same leks  <cit> .
we focus on five groups of leks with low rates of inter-group movement  to quantify genetic structure among leks and lek complexes and the extent to which geography and habitat suitability predict genetic differentiation.
our results indicate that geography, habitat, and potentially behaviors associated with a lek-based mating system may contribute to genetic differentiation and highlight the utility of large snp data sets for characterizing genetic structure at fine geographic scales.fig.
<dig> locations of leks within eureka county, nv, usa.
the study area varies across the landscape in both  topology and  habitat resistance.
the habitat resistance surface was derived from the inverse of a sage-grouse nest selection model  <cit>  and included the following variables : 1) nest site elevation ; 2) slope ; 3) slope * elevation ; 4) amount of habitat classified as sagebrush with 1000 m ; 5) distance from nearest lek ; 6) amount of habitat classified as sagebrush * distance to lek ; 7) amount of habitat classified as pinyon-juniper woodlands ; 8) amount of habitat classified as pinyon-juniper woodlands <dig> ; 9) distance from nearest water source ; and 10) distance from nearest water source <dig> .
warm colors are indicative of areas dominated by sagebrush or other shrubs, while cool colors are indicative of playas, pinyon-juniper woodlands, exotic grasslands, and pivot agricultural fields.
the legend is ordered from highest to lowest latitude, with the number of individual sage-grouse sampled at each lek listed parenthetically.
leks within the same lek complex share the same color
methods
sage-grouse sample collection
during the spring breeding season, we sampled ten leks separated by an average distance of  <dig>  km in eureka county, nevada  with known histories of individual movements and sufficiently large sample sizes  <cit> .
although 4– <dig> additional leks occur within this immediate area, the sampled leks have been the subject of past behavioral and ecological studies, remained constantly active from  <dig> to  <dig>  and were evenly distributed across the region  <cit> .
a priori, we grouped leks into five lek complexes based on geographic proximity , patterns of male movements documented across multiple years of monitoring, and potential habitat barriers  <cit> .
observed movements among leks were highly uncommon , but typically occurred within the same lek complex  <cit> .
two state highways, unpaved secondary roads, and power transmission lines are present in the study area, but sage-grouse regularly crossed these infrastructures .
there are no major geographic barriers to dispersal in the study area, but it is possible that individuals are less likely to move through low-quality habitat.
sage-grouse here make regular seasonal movements from spring breeding leks to higher elevation summer/fall habitat : individuals from the pine valley, kobeh valley, and pony express lek complexes utilize the roberts creek range, while individuals from the cortez and north cortez lek complexes typically occupy the cortez range  <cit> .
although birds from different leks and lek complexes share summer habitat and can have large yearly and summer movement distances , nearly all birds return to their specific lekking grounds each spring  <cit> .
thus, the sage-grouse in our study area exhibit strong lek-site fidelity.
we captured  <dig> male sage-grouse on these leks during the spring over a 5 year period .
we only included males in our study because we were unable to collect a sufficient number of female samples.
approximately  <dig>  ml of blood was drawn from the basilic vein of each captured bird.
protocols for the capture and handling of sage-grouse were approved by the university of nevada, reno institutional animal care and use committee .
genotyping-by-sequencing library preparation
dna was extracted from blood using qiagen dneasy blood and tissue kits  and quantified using a nanodrop spectrophotometer .
reduced-representation libraries were constructed using a gbs approach  <cit> .
genomic dna was first cut with two restriction enzymes, ecori and msei, and a unique dna barcode was ligated to digested fragments from each individual.
each ecori adaptor consisted of bases matching the endonuclease cut site, a unique eight to ten base pair dna barcode sequence, and an illumina adaptor; the msei adaptors consisted of the endonuclease cut site and the opposite illumina adaptor.
uniquely barcoded ligation products from all individuals were pooled and pcr-amplified using standard illumina primers.
libraries were size-selected for a region between  <dig> and  <dig> bases using a bluepippin quantitative electrophoresis unit .
we generated single end, 100 bp reads on one lane of an illumina hiseq  <dig> platform at the university of texas genomic sequencing and analysis facility .
assembly and variant calling
illumina reads were first filtered to remove potential contaminant dna  and low quality or aberrant reads.
we used a perl script to recognize barcodes assigned to each individual bird, to correct single-base errors in barcode sequences, and to remove sequences containing portions of the illumina adaptors.
after removing barcodes and cut site bases, retained sequences were  <dig>   <dig>  or 89 bp in length.
a reference of gbs sequenced regions was created with de novo assembly of a random subset of  <dig> million reads using the seqman ngen software .
this used a minimum match percentage of  <dig>  a gap penalty of  <dig>  and only retained contigs containing a minimum of  <dig> reads .
we generated the reference by removing contigs that were over-assembled or that were shorter than  <dig> or longer than  <dig> bases.
reads from each bird were subsequently aligned to the reference using the aln and samse algorithms in bwa .
the number of reads representing alternative nucleotide states  at each variant site was quantified using samtools and bcftools  <cit> .
we only called snps when 90 % of the individuals had at least one read at the locus.
genotype likelihoods were calculated with bcftools, stored in variant call format, and converted to a composite genotype likelihood format for downstream analysis.
we excluded variable sites with more than one alternate allele and loci with minor allele frequencies less than 5 %.
for contigs containing more than one snp, we randomly selected a single snp to increase the independence of loci used in subsequent analyses.
assembly related files and genotype matrices are available at dryad , and additional parameter settings for these analyses are available from the authors by request.
we used a hierarchical bayesian model that incorporates uncertainty in sequencing coverage across loci and individuals to estimate allele frequencies and genotype probabilities simultaneously for all individuals based on estimated genotype likelihoods  <cit> .
this model treats allele frequencies as priors, and coincidentally estimates allele frequencies and genotype probabilities while incorporating uncertainty arising from variation in sequence coverage.
thus, the estimates incorporate information on the probability of sampling each of the two alleles at a locus from the population, as well as coverage.
we obtained posterior estimates of genotype probabilities by running  <dig>  mcmc steps after a  <dig> step burn-in and thinning every other step.
these genotype probabilities were then converted to composite genotype values, where an individual’s genotype ranged from  <dig> to  <dig> at a locus .
the relationship between geography and genetic structure
to summarize genotypic variation across individuals and to explore the relationship between geography and genetic structure, we utilized three complementary ordination methods  <cit> .
for all three approaches, we implemented permutational multivariate analysis of variance  in the vegan package  <cit>  in r  <cit>  as a post hoc test of differentiation among leks and lek complexes based on euclidean distances of the first two ordination axes.
we first conducted principal component analysis  on genotype probabilities using the prcomp function in r. to quantify the relationship between individual genotypic variation and geography, we conducted procrustes analyses on principal components   <cit>  using the vegan package in r. we rotated a genomic matrix of pc <dig> and pc <dig> onto a geographic matrix of latitudes and longitudes and used  <dig> permutations to test significance of procrustes correlations.
to test the hypothesis that sage-grouse within a lek complex are more genetically similar to one another than expected, we calculated the euclidean genetic distance among all individuals using the first two pcs.
for each subset, the observed genetic distance of individual pairs sharing the same lek complex was compared to a null distribution of genetic distances constructed from  <dig> permutations of lek complex-sharing status, which was randomly assigned among pairs.
the same analyses were conducted at the lek level.
we conducted discriminant analysis of pcs   <cit>  using the adegenet package  <cit>  in r. unlike pca, dapc maximizes the ratio between the variance explained between groups and the total variance explained  <cit>  and is often used to group individuals into clusters.
dapc is useful for gbs datasets because the analysis effectively reduces the dimensionality of large datasets while delineating population genetic structure in a somewhat analogous manner to more computationally intensive bayesian clustering methods   <cit> .
to reduce over- or under-discrimination, we selected the number of retained pcs by predicting the maximum α-score with the optim.a.score function , following  <cit> .
we first conducted dapc without a priori group assignments by inferring the most likely number of genetic clusters  using the find.clusters function in the adegenet package.
this function utilizes k-means clustering to calculate a bayesian information criterion  value for each potential value of k  and delineates individual group assignments for dapc  <cit> .
additionally, we conducted dapc using a priori group assignments at both the lek  and lek complex  scales, which were based on patterns of male movement across multiple years of monitoring  <cit> .
we additionally used redundancy analysis   <cit>  to assess the relationship between geography and genetic structure.
rda is a constrained ordination analogous to multivariate regression in which a predictor dataset  is used to maximize the variance partitioned in a response dataset   <cit> .
rda has been previously employed to compare the relative strengths of geography and other potential barriers to gene flow  as predictors of population genetic structure .
we performed rda with latitude and longitude as predictors of our genotype probability matrix using the vegan package in r. the correlation between lek latitude and longitude was− <dig>  , which is smaller than the maximum correlation of rda predictor variables  suggested by  <cit> .
genetic structure and habitat resistance
to further summarize genetic differentiation among leks, we calculated pairwise fst  <cit>  based on allele frequencies within leks across all loci, and assessed significance using a permutation approach.
as an additional metric of differentiation, we calculated nei’s genetic distance  and used these values to generate neighbor-joining trees with the ape package  <cit>  in r.
to assess if environmental features contribute to observed differentiation among leks, we used a landscape genetic approach that compared genetic distance  with the summed environmental resistance  between each pair of leks using a least-cost paths framework  <cit> .
we derived a habitat resistance surface by taking the inverse of female sage-grouse nest site habitat selection model developed for this system  <cit>  .
additionally, we developed resistance surfaces including anthropogenic features  as barriers at different thresholds of movement permeability.
movement barrier resistance surfaces were derived from a 100 m buffer around all roads and power lines within the study area.
as outlined by  <cit> , each pixel within 100 m of a barrier was assigned either the maximum environmental resistance value from the habitat model  or assigned a resistance value four orders of magnitude greater than the maximum resistance value .
for each resistance surface, we used the landscape genetics toolbox  <cit>  in arcmap  <dig>   to produce a pairwise matrix of least-cost distance paths among leks.
we assessed the association between genetic distance and least-cost path for each resistance surface with mantel tests calculated in the ecodist package  <cit>  in r.
the effect of genomic sampling on genetic structure estimates
to determine the number of snps required to capture geographic genetic structure and to examine the robustness of our results to the number of loci included in our analyses, we performed power analyses for pca and dapc.
we constructed ten randomly sampled datasets of  <dig> different sizes, ranging from  <dig> to  <dig>  snps .
for each dataset, we conducted pca and procrustes analysis as described above.
additionally, we employed dapc on each subset to determine how data set size influenced the proportions of individuals correctly assigned to lek and lek complex.
we also evaluated the possible influence of loci with exceptionally high fst  on genetic structure parameter estimates by generating genotype sets with highly differentiated loci removed.
we generated locus-specific pairwise fst estimates  <cit>  for all  <dig> pairwise combinations of leks and for all ten pairwise combinations of lek complexes.
we removed loci with fst estimates above the 98th and 97th quantiles of the genome-wide fst distribution for any pairwise analysis at the lek and lek complex levels, and conducted pca and procrustes analyses on each of the four trimmed data sets.
RESULTS
after quality screening, we retained  <dig> , <dig> dna sequences from  <dig> individuals.
the initial de novo assembly placed  <dig> , <dig> reads into a set of  <dig>  high-coverage contigs; the consensus sequences of these contigs were then used as a reference for the genomic regions represented by our gbs data.
the final assembly contained 14 ×  <dig> reads  with an average of  <dig>  reads assembled per individual.
after identifying variant sites where >90 % of individuals had at least one read, we retained  <dig>  snps with minor allele frequencies >5 % and with an average coverage of  <dig> x per locus per individual .
analyses below are based on genotype probabilities estimated with the hierarchical bayesian model of  <cit> .
the relationship between geography and genetic structure
the first two pcs explained  <dig>  and  <dig>  % of the total genotypic variation  and clustered individuals from leks and lek complexes together to varying degrees.
although there was considerable overlap of adjacent leks, those separated by >30 km showed almost no overlap in pc space and lek complexes formed largely non-overlapping groups .
pc <dig> and pc <dig> scores were statistically different among both individual leks  and lek complexes .
there was a hierarchical pattern of genetic structure, where lek complexes exhibited stronger and more consistent differentiation, while neighboring leks exhibited subtle if any differentiation.
additionally, the mean genetic distance among individuals sharing the same lek complex was significantly smaller than a null distribution of genetic distances , suggesting that individuals within lek complexes have elevated relatedness .
the procrustes correlation between genomic and geographic matrices was  <dig>   and latitude explained  <dig>  % of the variation in pc <dig> .fig.
<dig> sage-grouse population genomic structure mirrors the fine-scale geography of leks.
a mean pca scores for each lek are denoted by symbols, with segments connecting means and individual points.
the proportion of variation explained by each pc is listed on each axis.
b pca scores were rotated onto a standardized geographic matrix to maximize the procrustes correlation  between genomic and geographic structure.
c a regression with latitude predicting pc <dig>  the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color
for the dapc without a priori group assignments, two pcs were retained based on  <dig> replicate α-score tests .
the k =  <dig> model fit the data best; bic steadily increased from k =  <dig> to k =  <dig> .
the first discriminant function  clearly separated sage-grouse individuals into two clusters ; one group mostly comprised individuals from the north cortez and cortez lek complexes, while the other group was composed of individuals from the pine valley, pony express, and kobeh valley lek complexes .
for the dapc analyses using a priori lek assignments,  <dig> pcs were retained for the lek dapc  and  <dig> pcs were retained for the lek complex dapc .
the lek complex dapc  performed much better at clustering individuals than the lek dapc ; in both, leks within a complex largely overlapped.
df <dig> separated the north cortez and cortez lek complexes from the pine valley, pony express, and kobeh valley lek complexes , similar to the results from the dapc without a priori group assignments .
although df <dig> largely separated the north cortez and cortez lek complexes, some overlap remained between the pine valley, pony express, kobeh valley lek complexes .
despite this, individual scores for df <dig> and df <dig> were significantly different for both individual leks  and lek complexes .
dapc assigned only  <dig>  % of individuals to the correct lek, but  <dig>  % of individuals were assigned to the correct lek complex .fig.
<dig> discriminant analysis of principal components  clusters sage-grouse individuals by a priori lek and lek complex designations.
scatterplots of df <dig> and df <dig> are displayed for both the  lek complex and  lek analyses.
insets display the eigenvalues for each df from dapc.
for the lek complex dapc, density distributions are displayed to aid in the visualization of how lek complexes cluster along  df <dig> and  df <dig>   individual assignment probabilities for each lek complex are displayed in a bar graph, with a priori lek and lek complex designations displayed above.
the legends are ordered from highest to lowest latitude and leks within the same lek complex share the same color
the first two rda axes  explained  <dig>  and  <dig>  % of the variance in total genotypic variation, respectively .
as in pca and dapc, there was a significant effect of lek  and lek complex  designations.fig.
<dig> the rda plot illustrates the strong relationship between geography and population genetic structure for sage-grouse in central nevada.
mean rda scores are plotted for each individual and vectors depict the direction of the two predictor variables  relative to the rda axes.
the proportion of variance explained by each rda axis is listed on each axis label.
the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color
genetic structure and habitat resistance
estimates of d were highly correlated with fst  and grouped leks within the same lek complex , consistent with multivariate analyses .
although the pairwise fst estimates were quite small, all but two were significantly larger than expected .
mantel tests indicated that geographic distance most strongly predicted variation in d among leks ; however, inclusion of the habitat resistance variable to the model explained an increased proportion of the variation in d among leks .
we did not find any support for an influence of anthropogenic barriers on genetic distance .fig.
<dig> a graphical summary of genetic distances among sage-grouse leks.
a a heat map displays pairwise genetic distance  and pairwise f
st  among leks ordered from highest to lowest latitude.
all but two pairwise f
st estimates were significantly larger than expected .
b a neighbor-joining tree made from d depicts the genetic relationships among individual leks and coincides with the results from pca  and dapc .
the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color
r
2
p
distance + habitata
distance + permeable barriersb
distance + low-permeable barriersb
athe habitat suitability resistance surface was derived from the inverse of a sage-grouse nest selection model
bmovement barrier resistance surfaces were derived from a 100 m buffer around all roads and power lines within the study area.
each pixel within 100 m of a barrier was assigned either the maximum environmental resistance value observed in the habitat suitability model  or assigned an environmental resistance value four orders of magnitude greater than the maximum environmental resistance value
fig.
<dig> results from a mantel test with pairwise geographic distance among leks predicting pairwise genetic distance  are consistent with a pattern of isolation by distance
the effect of genomic sampling on genetic structure estimates
the magnitude and precision of procrustes correlations between genetic and geographic matrices increased strongly in subsets containing more snps, ultimately plateauing around ~ <dig> snps .
a similar pattern was found for dapc, with the proportion of correctly assigned individuals reaching an asymptote at ~ <dig> snps at both the lek and lek complex scales .
thus, data sets containing fewer than ~ <dig> snps generated parameter estimates with reduced precision , and failed to quantify the geographic scale of genetic structure.
however, data sets containing more than ~ <dig> snps generated similar parameter estimates, indicating that the pattern of genetic structure is not driven by small numbers of exceptional loci or those residing in specific genomic regions.fig.
<dig> to determine the effect of increased snp sampling on the precision and information content of population genetic analyses for this system,  procrustes analysis and  dapc were performed for  <dig> subsets of the total number of snps .
ten replicates were analyzed for each subset.
note that approximately  <dig> loci are needed to ensure precise estimates of genetic structure .
inset depicts a male sage-grouse displaying at a lek in central nevada
we also examined how the results displayed in fig.
<dig>  were affected by removing loci with exceptionally high fst estimates.
after removing loci with fst estimates above the 97th and 98th quantile for all pairwise lek comparisons, subsets of  <dig>  and  <dig>  loci were retained, while subsets of  <dig>  and  <dig>  loci were retained after the lek complex comparisons.
results from pca, procrustes, and regression analyses based on the four “outlier” trimmed genotype subsets were indistinguishable from those for the full data set .
discussion
our results indicate that lek geography and habitat features influence sage-grouse population genetic variation at a fine geographic scale across the individual, lek, and lek complex levels.
consistent with isolation by lek distance, geography explained a large proportion of individual genetic variation in ordination space , predetermined delineations of lek complexes clustered together , and lek-level estimates of genetic distance were strongly related to geographic distance .
although statistically significant, levels of genetic differentiation among neighboring leks were low, as indicated by the low percentage of variation explained by the first two axes of pca, dapc, and rda and by the small estimates of genetic differentiation among leks .
such low levels of differentiation are not surprising and are consistent with past studies based on smaller datasets not detecting genetic differentiation across geographically proximate leks .
sage-grouse from the different lek complexes formed distinct and statistically identifiable clusters in ordination space.
the most pronounced genetic structure divided northern lek complexes  from more southern complexes  .
the dapc analysis without a priori information inferred k =  <dig> as the best model and assigned individuals almost entirely to these two groups .
this pattern is consistent with sage-grouse from the cortez and north cortez complexes predominantly sharing summer brood rearing habitat in the cortez range, while individuals from other complexes typically cohabitate in the roberts creek range  <cit> .
nonetheless, individuals are highly faithful to their specific leks during subsequent spring breeding seasons  <cit> .
additionally these results are in accord with other studies reporting that the pattern and strength of genetic structure can depend on which stage in the breeding cycle samples are collected  .
within these groups, a priori designated lek complexes formed largely non-overlapping clusters in ordination analyses , with ~80 % of individuals assigned to the correct lek complex in dapc .
differentiation among neighboring leks was much less pronounced; only ~49 % of individuals were assigned to the correct lek in dapc, and individuals from neighboring leks overlapped substantially in ordination space.
even with limited gene flow among neighboring leks, our results illustrate a hierarchical pattern of fine-scale spatial genetic structure consistent with isolation by effective lek and lek complex distance .
the strong site fidelity previously documented in this system  <cit>  could reduce gene flow among leks at some geographic scales.
capture-mark-recapture data from these leks indicate that inter-lek movements by males are very rare  even though individuals have overlapping summer distributions  <cit>  and can move large distances across the landscape within a year .
by using geographic distances among leks in our population-level genetic analyses we more directly accounted for the effect of lek location on genetic structure.
nearly all of the movements we did observe were among neighboring leks, which could contribute to limited differentiation within complexes .
our landscape genetic analyses based on habitat resistance surfaces indicate that habitat suitability also influences genetic distance among leks and could reinforce the effect of site fidelity on hierarchical genetic structure at the lek and lek complex levels .
similar to  <cit> , we did not find evidence of an influence of anthropogenic barriers on genetic distance among leks, even though certain structures within this system  were associated with reduced population growth  <cit> .
these results contrast with another study that documented significant impacts of anthropogenic barriers on population genetic structure in sage-grouse from the state of washington  <cit> .
reproductive skew caused by a subset of males contributing disproportionately to reproduction  <cit>  could also potentially influence population differentiation by reducing effective population sizes and reducing the impact of migration on genetic differentiation  <cit> .
empirical evidence consistent with reproductive skew and/or kin selection has been found for other bird species with lek-based or cooperative mating systems ; however, the strength of kin selection can be context-dependent  <cit> .
previous analyses of sage-grouse have not found evidence for increased relatedness within leks .
nevertheless, sage-grouse that shared lek complexes  in our study were more genetically similar than would be expected by chance , a pattern one would predict if reproductive skew were influencing geographic genetic variation.
although it is difficult to directly account for the effects of reproductive skew and kin selection on levels of population genetic structure, our results suggest that the lek-breeding system could potentially influence the hierarchical pattern of genetic differentiation we observe.
indeed, social and reproductive behaviors are commonly associated with genetic structure in vertebrates and insects .
as habitat fragmentation and reproductive skew could enhance differentiation and lead to loss of variation, such mating systems could pose special challenges for protecting genetic diversity in declining species.
our analyses provide a finer view of geographic genetic structure than most previous studies of sage-grouse , although some studies of lek-mating birds have detected a subtle influence of lek geography .
in most previous studies, genetic differentiation among sage-grouse populations has been documented when leks have been isolated by large expanses of low-quality habitat .
studies at smaller geographic scales have typically recovered less fine-scale structure, but have found evidence for isolation by distance  and, in one case, genetic differentiation among groups of leks found in neighboring management zones  <cit> .
it is possible that the strength of genetic structure in our study was elevated as a result of either excluding females  or including leks that occupy a part of the range where site fidelity, relatedness, and genetic differentiation interact at different levels.
however, our ability to detect such fine-scale genetic structure was likely influenced by the much larger number of loci employed in our study .
not surprisingly, we found that the level of genetic structure detected and the precision of population genetic parameter estimates initially increased as the number of loci analyzed rose .
larger numbers of loci  more successfully captured the relationship between genetic variation and geography  and increased the fraction of individuals assigned to the correct lek and lek complex in discriminant analyses .
however, parameter estimates and their precision were consistent in randomly built data sets with > <dig> loci .
in addition, parameter estimates in our analyses of the full data set of snps were indistinguishable from those for data sets with high fst outlier loci removed.
this indicates that the overall patterns in our results were not strongly influenced by loci with exceptional patterns of differentiation or residing in particular genomic locations.
our results are consistent with many recent studies that illustrate how the information content of population genetic analyses can depend on the extent of genomic sampling .
for instance, a recent study of soay sheep found that the precision and accuracy of heritability estimates grew with increasing marker numbers , and parameter estimates including this many markers were comparable to results calculated using pedigrees  <cit> .
CONCLUSIONS
along with other recent studies detecting previously unrecognized structure, our results demonstrate the promise of population genomic data sets for quantifying fine-scale, ecologically relevant population genetic variation  <cit> .
the increase in resolution such approaches could bring is highly relevant to the field of conservation biology  <cit>  because the discovery of cryptic genetic structure in species of conservation concern has the potential to inform management strategies .
indeed, the fine-scale geographic genetic structure of numerous species may be shaped in part by behavioral processes, including habitat selection , site fidelity , and reproductive skew .
future investigations utilizing genome-level data sets and covering wider ecological and geographic contexts should help to tease apart the influence of habitat features, behavioral ecology, and gene flow in shaping fine-scale genetic variation in natural populations.
abbreviations
bic, bayesian information criterion; dapc, discriminant analysis of principal components; df, discriminant function; gbs, genotyping-by-sequencing; pc, principal component; pca, principal component analysis; permanova, permutational multivariate analysis of variance; rda, redundancy analysis; snp, single nucleotide polymorphism
additional file
additional file 1: table s <dig>  geographic coordinates and number of individuals included in final analyses  for each lek and lek complex.
table s <dig>  mean observed f
st estimates among leks was compared to a null distribution of f
st values .
figure s <dig>  histograms display the distribution of  yearly and  summer movement distances for sage-grouse individuals in this system based on telemetry data.
figure s <dig>  histograms display the distribution of euclidean genetic distances for  <dig> permutations of randomly sampled subsets of sage grouse pairs.
these permutations were used to test if sage-grouse within a  lek complex and  lek were genetically more similar than expected from a null distribution.
figure s <dig>  the correlation between the genomic and geographic matrices obtained from the procrustes analysis  is compared to the distribution of correlations from  <dig> random permutations.
figure s <dig>  summary results from the dapc that did not use a priori group assignments are displayed.
bayesian information criterion  steadily increases as the number of clusters  increases in the model .
for each group, a density histogram is displayed to ease in the visualization of how individuals cluster along discriminant function one.
additionally, individual assignment probabilities for each lek complex are displayed in a bar graph, with a priori lek and lek complex designations displayed above .
figure s <dig>  results from pca using the entire snp dataset  were compared to subsets of data with the top 2 %  and top 3 % of fst  outliers removed.
any loci that were in the top  <dig> or 3 % of the fst distribution for any pairwise lek comparison were removed to create the subsets.
panels , , and  display scatterplots of pc <dig> and pc <dig>  with symbols representing lek means and segments drawn from the means to individual scores.
panels , , and  depict the procrustes correlation between the first two pcs and a standardized geographic matrix.
panels , , and  show the results of a regression with latitude predicting pc <dig>

 population genetic studies are essential for understanding how behavioral, ecological, and evolutionary processes shape the geographic distribution of genetic variation. recent dna sequencing advances now allow variation to be readily assayed at thousands of loci for large numbers of individuals and populations in non-model organisms, potentially illuminating the relationships among ecological, genetic, and geographic variation at previously unattainable scales. these approaches have value for uncovering previously unrecognized patterns of genetic structure, generating more precise estimates of demographic and evolutionary parameters, and providing increased opportunity for understanding adaptation. such approaches also have special promise for analyses of genetic variation related to contemporary ecological processes and population demography in species of ecological or conservation significance. variation in behavioral and life history traits may have important consequences for the geographic structuring of genetic variation across populations. for example, a number of avian species have lek-based mating systems in which a small number of males contribute disproportionately to reproduction, dispersal is restricted by strong site fidelity, and leks remain at the same geographic locations for long periods of time. strong site fidelity could reduce gene flow among leks and give rise to spatial genetic structure even in the absence of geographic barriers. because a subset of males typically contribute disproportionately to reproduction on leks, reproductive skew can also lower effective population size and increase genetic drift. these factors, in addition to the potential for kin selection in lek-based mating systems, could enhance the effects of drift in structured populations, give rise to fine-scale spatial genetic structure, and increase extinction probability of isolated populations. greater sage-grouse are an iconic western north american bird species with a lek-based mating system. male sage-grouse display for females at leks residing in open patches of sagebrush, where a subset of dominant males can account for a large proportion of the successful matings. degradation of sagebrush ecosystems across north america has substantially reduced the geographic range and abundance of sage-grouse over the last century. population declines ranged from to 90 % across the distribution and local extirpations have occurred in many regions. habitat degradation and population declines have raised concerns about population persistence, though sage-grouse were recently found unwarranted for listing under the u.s. endangered species act. low movement rates among leks due to high site fidelity, in addition to reproductive skew, could interact with habitat variation to limit gene flow, enhance genetic drift, and shape fine-scale genetic structure in sage-grouse populations, as has been found in several lek-breeding species. previous studies utilizing small data sets of traditional molecular markers have reported variable patterns of genetic structure across sage-grouse populations, leaving uncertainty about the scale of geographic genetic structure and obscuring the expected link between demographic and population genetic parameters across the range. the detection of genetic differentiation across broad geographic scales resulted in the recognition of three distinct populations: gunnison sage-grouse, greater sage-grouse, and the bi-state population of greater sage-grouse. however, most studies at finer geographic scales have not identified genetic structure among geographically proximate leks, leaving uncertainty about the potential for site fidelity and reproductive skew to influence genetic structure. in contrast, a recent study of geographically proximate leks in northwestern colorado documented among-lek differentiation and evidence for isolation by distance. although these studies suggest movements among leks may limit differentiation in some parts of the range, higher resolution data sets should improve our understanding of the geographic scale of genetic structure and the factors associated with it. here we use a genotyping-by-sequencing approach to quantify genetic variation at single-nucleotide polymorphisms across sage-grouse from ten leks located throughout a small geographic region in central nevada. very strong breeding-site fidelity has been documented for these specific leks across multiple years of monitoring, where of among-year observations involved males returning to the same leks. we focus on five groups of leks with low rates of inter-group movement to quantify genetic structure among leks and lek complexes and the extent to which geography and habitat suitability predict genetic differentiation. our results indicate that geography, habitat, and potentially behaviors associated with a lek-based mating system may contribute to genetic differentiation and highlight the utility of large snp data sets for characterizing genetic structure at fine geographic scales.fig. locations of leks within eureka county, nv, usa. the study area varies across the landscape in both topology and habitat resistance. the habitat resistance surface was derived from the inverse of a sage-grouse nest selection model and included the following variables: 1) nest site elevation; 2) slope; 3) slope * elevation; 4) amount of habitat classified as sagebrush with 1000 m; 5) distance from nearest lek; 6) amount of habitat classified as sagebrush * distance to lek; 7) amount of habitat classified as pinyon-juniper woodlands; 8) amount of habitat classified as pinyon-juniper woodlands; 9) distance from nearest water source; and 10) distance from nearest water source. warm colors are indicative of areas dominated by sagebrush or other shrubs, while cool colors are indicative of playas, pinyon-juniper woodlands, exotic grasslands, and pivot agricultural fields. the legend is ordered from highest to lowest latitude, with the number of individual sage-grouse sampled at each lek listed parenthetically. leks within the same lek complex share the same color methods sage-grouse sample collection during the spring breeding season, we sampled ten leks separated by an average distance of  km in eureka county, nevada with known histories of individual movements and sufficiently large sample sizes. although 4– additional leks occur within this immediate area, the sampled leks have been the subject of past behavioral and ecological studies, remained constantly active from to and were evenly distributed across the region. a priori, we grouped leks into five lek complexes based on geographic proximity, patterns of male movements documented across multiple years of monitoring, and potential habitat barriers. observed movements among leks were highly uncommon, but typically occurred within the same lek complex. two state highways, unpaved secondary roads, and power transmission lines are present in the study area, but sage-grouse regularly crossed these infrastructures. there are no major geographic barriers to dispersal in the study area, but it is possible that individuals are less likely to move through low-quality habitat. sage-grouse here make regular seasonal movements from spring breeding leks to higher elevation summer/fall habitat: individuals from the pine valley, kobeh valley, and pony express lek complexes utilize the roberts creek range, while individuals from the cortez and north cortez lek complexes typically occupy the cortez range. although birds from different leks and lek complexes share summer habitat and can have large yearly and summer movement distances, nearly all birds return to their specific lekking grounds each spring. thus, the sage-grouse in our study area exhibit strong lek-site fidelity. we captured male sage-grouse on these leks during the spring over a 5 year period. we only included males in our study because we were unable to collect a sufficient number of female samples. approximately  ml of blood was drawn from the basilic vein of each captured bird. protocols for the capture and handling of sage-grouse were approved by the university of nevada, reno institutional animal care and use committee. genotyping-by-sequencing library preparation dna was extracted from blood using qiagen dneasy blood and tissue kits and quantified using a nanodrop spectrophotometer. reduced-representation libraries were constructed using a gbs approach. genomic dna was first cut with two restriction enzymes, ecori and msei, and a unique dna barcode was ligated to digested fragments from each individual. each ecori adaptor consisted of bases matching the endonuclease cut site, a unique eight to ten base pair dna barcode sequence, and an illumina adaptor; the msei adaptors consisted of the endonuclease cut site and the opposite illumina adaptor. uniquely barcoded ligation products from all individuals were pooled and pcr-amplified using standard illumina primers. libraries were size-selected for a region between and bases using a bluepippin quantitative electrophoresis unit. we generated single end, 100 bp reads on one lane of an illumina hiseq platform at the university of texas genomic sequencing and analysis facility. assembly and variant calling illumina reads were first filtered to remove potential contaminant dna and low quality or aberrant reads. we used a perl script to recognize barcodes assigned to each individual bird, to correct single-base errors in barcode sequences, and to remove sequences containing portions of the illumina adaptors. after removing barcodes and cut site bases, retained sequences were or 89 bp in length. a reference of gbs sequenced regions was created with de novo assembly of a random subset of million reads using the seqman ngen software. this used a minimum match percentage of a gap penalty of and only retained contigs containing a minimum of reads. we generated the reference by removing contigs that were over-assembled or that were shorter than or longer than bases. reads from each bird were subsequently aligned to the reference using the aln and samse algorithms in bwa. the number of reads representing alternative nucleotide states at each variant site was quantified using samtools and bcftools. we only called snps when 90 % of the individuals had at least one read at the locus. genotype likelihoods were calculated with bcftools, stored in variant call format, and converted to a composite genotype likelihood format for downstream analysis. we excluded variable sites with more than one alternate allele and loci with minor allele frequencies less than 5 %. for contigs containing more than one snp, we randomly selected a single snp to increase the independence of loci used in subsequent analyses. assembly related files and genotype matrices are available at dryad, and additional parameter settings for these analyses are available from the authors by request. we used a hierarchical bayesian model that incorporates uncertainty in sequencing coverage across loci and individuals to estimate allele frequencies and genotype probabilities simultaneously for all individuals based on estimated genotype likelihoods. this model treats allele frequencies as priors, and coincidentally estimates allele frequencies and genotype probabilities while incorporating uncertainty arising from variation in sequence coverage. thus, the estimates incorporate information on the probability of sampling each of the two alleles at a locus from the population, as well as coverage. we obtained posterior estimates of genotype probabilities by running mcmc steps after a step burn-in and thinning every other step. these genotype probabilities were then converted to composite genotype values, where an individual’s genotype ranged from to at a locus. the relationship between geography and genetic structure to summarize genotypic variation across individuals and to explore the relationship between geography and genetic structure, we utilized three complementary ordination methods. for all three approaches, we implemented permutational multivariate analysis of variance in the vegan package in r as a post hoc test of differentiation among leks and lek complexes based on euclidean distances of the first two ordination axes. we first conducted principal component analysis on genotype probabilities using the prcomp function in r. to quantify the relationship between individual genotypic variation and geography, we conducted procrustes analyses on principal components using the vegan package in r. we rotated a genomic matrix of pc and pc onto a geographic matrix of latitudes and longitudes and used permutations to test significance of procrustes correlations. to test the hypothesis that sage-grouse within a lek complex are more genetically similar to one another than expected, we calculated the euclidean genetic distance among all individuals using the first two pcs. for each subset, the observed genetic distance of individual pairs sharing the same lek complex was compared to a null distribution of genetic distances constructed from permutations of lek complex-sharing status, which was randomly assigned among pairs. the same analyses were conducted at the lek level. we conducted discriminant analysis of pcs using the adegenet package in r. unlike pca, dapc maximizes the ratio between the variance explained between groups and the total variance explained and is often used to group individuals into clusters. dapc is useful for gbs datasets because the analysis effectively reduces the dimensionality of large datasets while delineating population genetic structure in a somewhat analogous manner to more computationally intensive bayesian clustering methods. to reduce over- or under-discrimination, we selected the number of retained pcs by predicting the maximum α-score with the optim.a.score function, following. we first conducted dapc without a priori group assignments by inferring the most likely number of genetic clusters using the find.clusters function in the adegenet package. this function utilizes k-means clustering to calculate a bayesian information criterion value for each potential value of k and delineates individual group assignments for dapc. additionally, we conducted dapc using a priori group assignments at both the lek and lek complex scales, which were based on patterns of male movement across multiple years of monitoring. we additionally used redundancy analysis to assess the relationship between geography and genetic structure. rda is a constrained ordination analogous to multivariate regression in which a predictor dataset is used to maximize the variance partitioned in a response dataset. rda has been previously employed to compare the relative strengths of geography and other potential barriers to gene flow as predictors of population genetic structure. we performed rda with latitude and longitude as predictors of our genotype probability matrix using the vegan package in r. the correlation between lek latitude and longitude was−, which is smaller than the maximum correlation of rda predictor variables suggested by. genetic structure and habitat resistance to further summarize genetic differentiation among leks, we calculated pairwise fst based on allele frequencies within leks across all loci, and assessed significance using a permutation approach. as an additional metric of differentiation, we calculated nei’s genetic distance and used these values to generate neighbor-joining trees with the ape package in r. to assess if environmental features contribute to observed differentiation among leks, we used a landscape genetic approach that compared genetic distance with the summed environmental resistance between each pair of leks using a least-cost paths framework. we derived a habitat resistance surface by taking the inverse of female sage-grouse nest site habitat selection model developed for this system. additionally, we developed resistance surfaces including anthropogenic features as barriers at different thresholds of movement permeability. movement barrier resistance surfaces were derived from a 100 m buffer around all roads and power lines within the study area. as outlined by, each pixel within 100 m of a barrier was assigned either the maximum environmental resistance value from the habitat model or assigned a resistance value four orders of magnitude greater than the maximum resistance value. for each resistance surface, we used the landscape genetics toolbox in arcmap to produce a pairwise matrix of least-cost distance paths among leks. we assessed the association between genetic distance and least-cost path for each resistance surface with mantel tests calculated in the ecodist package in r. the effect of genomic sampling on genetic structure estimates to determine the number of snps required to capture geographic genetic structure and to examine the robustness of our results to the number of loci included in our analyses, we performed power analyses for pca and dapc. we constructed ten randomly sampled datasets of different sizes, ranging from to snps. for each dataset, we conducted pca and procrustes analysis as described above. additionally, we employed dapc on each subset to determine how data set size influenced the proportions of individuals correctly assigned to lek and lek complex. we also evaluated the possible influence of loci with exceptionally high fst on genetic structure parameter estimates by generating genotype sets with highly differentiated loci removed. we generated locus-specific pairwise fst estimates for all pairwise combinations of leks and for all ten pairwise combinations of lek complexes. we removed loci with fst estimates above the 98th and 97th quantiles of the genome-wide fst distribution for any pairwise analysis at the lek and lek complex levels, and conducted pca and procrustes analyses on each of the four trimmed data sets.  after quality screening, we retained, dna sequences from individuals. the initial de novo assembly placed, reads into a set of high-coverage contigs; the consensus sequences of these contigs were then used as a reference for the genomic regions represented by our gbs data. the final assembly contained 14 ×  reads with an average of reads assembled per individual. after identifying variant sites where >90 % of individuals had at least one read, we retained snps with minor allele frequencies >5 % and with an average coverage of x per locus per individual. analyses below are based on genotype probabilities estimated with the hierarchical bayesian model of. the relationship between geography and genetic structure the first two pcs explained and  % of the total genotypic variation and clustered individuals from leks and lek complexes together to varying degrees. although there was considerable overlap of adjacent leks, those separated by >30 km showed almost no overlap in pc space and lek complexes formed largely non-overlapping groups. pc and pc scores were statistically different among both individual leks and lek complexes. there was a hierarchical pattern of genetic structure, where lek complexes exhibited stronger and more consistent differentiation, while neighboring leks exhibited subtle if any differentiation. additionally, the mean genetic distance among individuals sharing the same lek complex was significantly smaller than a null distribution of genetic distances, suggesting that individuals within lek complexes have elevated relatedness. the procrustes correlation between genomic and geographic matrices was and latitude explained  % of the variation in pc.fig. sage-grouse population genomic structure mirrors the fine-scale geography of leks. a mean pca scores for each lek are denoted by symbols, with segments connecting means and individual points. the proportion of variation explained by each pc is listed on each axis. b pca scores were rotated onto a standardized geographic matrix to maximize the procrustes correlation between genomic and geographic structure. c a regression with latitude predicting pc the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color for the dapc without a priori group assignments, two pcs were retained based on replicate α-score tests. the k =  model fit the data best; bic steadily increased from k =  to k = . the first discriminant function clearly separated sage-grouse individuals into two clusters; one group mostly comprised individuals from the north cortez and cortez lek complexes, while the other group was composed of individuals from the pine valley, pony express, and kobeh valley lek complexes. for the dapc analyses using a priori lek assignments, pcs were retained for the lek dapc and pcs were retained for the lek complex dapc. the lek complex dapc performed much better at clustering individuals than the lek dapc; in both, leks within a complex largely overlapped. df separated the north cortez and cortez lek complexes from the pine valley, pony express, and kobeh valley lek complexes, similar to the results from the dapc without a priori group assignments. although df largely separated the north cortez and cortez lek complexes, some overlap remained between the pine valley, pony express, kobeh valley lek complexes. despite this, individual scores for df and df were significantly different for both individual leks and lek complexes. dapc assigned only  % of individuals to the correct lek, but  % of individuals were assigned to the correct lek complex.fig. discriminant analysis of principal components clusters sage-grouse individuals by a priori lek and lek complex designations. scatterplots of df and df are displayed for both the lek complex and lek analyses. insets display the eigenvalues for each df from dapc. for the lek complex dapc, density distributions are displayed to aid in the visualization of how lek complexes cluster along df and df individual assignment probabilities for each lek complex are displayed in a bar graph, with a priori lek and lek complex designations displayed above. the legends are ordered from highest to lowest latitude and leks within the same lek complex share the same color the first two rda axes explained and  % of the variance in total genotypic variation, respectively. as in pca and dapc, there was a significant effect of lek and lek complex designations.fig. the rda plot illustrates the strong relationship between geography and population genetic structure for sage-grouse in central nevada. mean rda scores are plotted for each individual and vectors depict the direction of the two predictor variables relative to the rda axes. the proportion of variance explained by each rda axis is listed on each axis label. the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color genetic structure and habitat resistance estimates of d were highly correlated with fst and grouped leks within the same lek complex, consistent with multivariate analyses. although the pairwise fst estimates were quite small, all but two were significantly larger than expected. mantel tests indicated that geographic distance most strongly predicted variation in d among leks; however, inclusion of the habitat resistance variable to the model explained an increased proportion of the variation in d among leks. we did not find any support for an influence of anthropogenic barriers on genetic distance.fig. a graphical summary of genetic distances among sage-grouse leks. a a heat map displays pairwise genetic distance and pairwise f st among leks ordered from highest to lowest latitude. all but two pairwise f st estimates were significantly larger than expected. b a neighbor-joining tree made from d depicts the genetic relationships among individual leks and coincides with the results from pca and dapc. the legend is ordered from highest to lowest latitude and leks within the same lek complex share the same color r 2 p distance + habitata distance + permeable barriersb distance + low-permeable barriersb athe habitat suitability resistance surface was derived from the inverse of a sage-grouse nest selection model bmovement barrier resistance surfaces were derived from a 100 m buffer around all roads and power lines within the study area. each pixel within 100 m of a barrier was assigned either the maximum environmental resistance value observed in the habitat suitability model or assigned an environmental resistance value four orders of magnitude greater than the maximum environmental resistance value fig. results from a mantel test with pairwise geographic distance among leks predicting pairwise genetic distance are consistent with a pattern of isolation by distance the effect of genomic sampling on genetic structure estimates the magnitude and precision of procrustes correlations between genetic and geographic matrices increased strongly in subsets containing more snps, ultimately plateauing around ~ snps. a similar pattern was found for dapc, with the proportion of correctly assigned individuals reaching an asymptote at ~ snps at both the lek and lek complex scales. thus, data sets containing fewer than ~ snps generated parameter estimates with reduced precision, and failed to quantify the geographic scale of genetic structure. however, data sets containing more than ~ snps generated similar parameter estimates, indicating that the pattern of genetic structure is not driven by small numbers of exceptional loci or those residing in specific genomic regions.fig. to determine the effect of increased snp sampling on the precision and information content of population genetic analyses for this system, procrustes analysis and dapc were performed for subsets of the total number of snps. ten replicates were analyzed for each subset. note that approximately loci are needed to ensure precise estimates of genetic structure. inset depicts a male sage-grouse displaying at a lek in central nevada we also examined how the results displayed in fig.  were affected by removing loci with exceptionally high fst estimates. after removing loci with fst estimates above the 97th and 98th quantile for all pairwise lek comparisons, subsets of and loci were retained, while subsets of and loci were retained after the lek complex comparisons. results from pca, procrustes, and regression analyses based on the four “outlier” trimmed genotype subsets were indistinguishable from those for the full data set. discussion our results indicate that lek geography and habitat features influence sage-grouse population genetic variation at a fine geographic scale across the individual, lek, and lek complex levels. consistent with isolation by lek distance, geography explained a large proportion of individual genetic variation in ordination space, predetermined delineations of lek complexes clustered together, and lek-level estimates of genetic distance were strongly related to geographic distance. although statistically significant, levels of genetic differentiation among neighboring leks were low, as indicated by the low percentage of variation explained by the first two axes of pca, dapc, and rda and by the small estimates of genetic differentiation among leks. such low levels of differentiation are not surprising and are consistent with past studies based on smaller datasets not detecting genetic differentiation across geographically proximate leks. sage-grouse from the different lek complexes formed distinct and statistically identifiable clusters in ordination space. the most pronounced genetic structure divided northern lek complexes from more southern complexes. the dapc analysis without a priori information inferred k =  as the best model and assigned individuals almost entirely to these two groups. this pattern is consistent with sage-grouse from the cortez and north cortez complexes predominantly sharing summer brood rearing habitat in the cortez range, while individuals from other complexes typically cohabitate in the roberts creek range. nonetheless, individuals are highly faithful to their specific leks during subsequent spring breeding seasons. additionally these results are in accord with other studies reporting that the pattern and strength of genetic structure can depend on which stage in the breeding cycle samples are collected. within these groups, a priori designated lek complexes formed largely non-overlapping clusters in ordination analyses, with ~80 % of individuals assigned to the correct lek complex in dapc. differentiation among neighboring leks was much less pronounced; only ~49 % of individuals were assigned to the correct lek in dapc, and individuals from neighboring leks overlapped substantially in ordination space. even with limited gene flow among neighboring leks, our results illustrate a hierarchical pattern of fine-scale spatial genetic structure consistent with isolation by effective lek and lek complex distance. the strong site fidelity previously documented in this system could reduce gene flow among leks at some geographic scales. capture-mark-recapture data from these leks indicate that inter-lek movements by males are very rare even though individuals have overlapping summer distributions and can move large distances across the landscape within a year. by using geographic distances among leks in our population-level genetic analyses we more directly accounted for the effect of lek location on genetic structure. nearly all of the movements we did observe were among neighboring leks, which could contribute to limited differentiation within complexes. our landscape genetic analyses based on habitat resistance surfaces indicate that habitat suitability also influences genetic distance among leks and could reinforce the effect of site fidelity on hierarchical genetic structure at the lek and lek complex levels. similar to, we did not find evidence of an influence of anthropogenic barriers on genetic distance among leks, even though certain structures within this system were associated with reduced population growth. these results contrast with another study that documented significant impacts of anthropogenic barriers on population genetic structure in sage-grouse from the state of washington. reproductive skew caused by a subset of males contributing disproportionately to reproduction could also potentially influence population differentiation by reducing effective population sizes and reducing the impact of migration on genetic differentiation. empirical evidence consistent with reproductive skew and/or kin selection has been found for other bird species with lek-based or cooperative mating systems; however, the strength of kin selection can be context-dependent. previous analyses of sage-grouse have not found evidence for increased relatedness within leks. nevertheless, sage-grouse that shared lek complexes in our study were more genetically similar than would be expected by chance, a pattern one would predict if reproductive skew were influencing geographic genetic variation. although it is difficult to directly account for the effects of reproductive skew and kin selection on levels of population genetic structure, our results suggest that the lek-breeding system could potentially influence the hierarchical pattern of genetic differentiation we observe. indeed, social and reproductive behaviors are commonly associated with genetic structure in vertebrates and insects. as habitat fragmentation and reproductive skew could enhance differentiation and lead to loss of variation, such mating systems could pose special challenges for protecting genetic diversity in declining species. our analyses provide a finer view of geographic genetic structure than most previous studies of sage-grouse, although some studies of lek-mating birds have detected a subtle influence of lek geography. in most previous studies, genetic differentiation among sage-grouse populations has been documented when leks have been isolated by large expanses of low-quality habitat. studies at smaller geographic scales have typically recovered less fine-scale structure, but have found evidence for isolation by distance and, in one case, genetic differentiation among groups of leks found in neighboring management zones. it is possible that the strength of genetic structure in our study was elevated as a result of either excluding females or including leks that occupy a part of the range where site fidelity, relatedness, and genetic differentiation interact at different levels. however, our ability to detect such fine-scale genetic structure was likely influenced by the much larger number of loci employed in our study. not surprisingly, we found that the level of genetic structure detected and the precision of population genetic parameter estimates initially increased as the number of loci analyzed rose. larger numbers of loci more successfully captured the relationship between genetic variation and geography and increased the fraction of individuals assigned to the correct lek and lek complex in discriminant analyses. however, parameter estimates and their precision were consistent in randomly built data sets with > loci. in addition, parameter estimates in our analyses of the full data set of snps were indistinguishable from those for data sets with high fst outlier loci removed. this indicates that the overall patterns in our results were not strongly influenced by loci with exceptional patterns of differentiation or residing in particular genomic locations. our results are consistent with many recent studies that illustrate how the information content of population genetic analyses can depend on the extent of genomic sampling. for instance, a recent study of soay sheep found that the precision and accuracy of heritability estimates grew with increasing marker numbers, and parameter estimates including this many markers were comparable to results calculated using pedigrees.  along with other recent studies detecting previously unrecognized structure, our results demonstrate the promise of population genomic data sets for quantifying fine-scale, ecologically relevant population genetic variation. the increase in resolution such approaches could bring is highly relevant to the field of conservation biology because the discovery of cryptic genetic structure in species of conservation concern has the potential to inform management strategies. indeed, the fine-scale geographic genetic structure of numerous species may be shaped in part by behavioral processes, including habitat selection, site fidelity, and reproductive skew. future investigations utilizing genome-level data sets and covering wider ecological and geographic contexts should help to tease apart the influence of habitat features, behavioral ecology, and gene flow in shaping fine-scale genetic variation in natural populations. abbreviations bic, bayesian information criterion; dapc, discriminant analysis of principal components; df, discriminant function; gbs, genotyping-by-sequencing; pc, principal component; pca, principal component analysis; permanova, permutational multivariate analysis of variance; rda, redundancy analysis; snp, single nucleotide polymorphism additional file additional file 1: table s geographic coordinates and number of individuals included in final analyses for each lek and lek complex. table s mean observed f st estimates among leks was compared to a null distribution of f st values. figure s histograms display the distribution of yearly and summer movement distances for sage-grouse individuals in this system based on telemetry data. figure s histograms display the distribution of euclidean genetic distances for permutations of randomly sampled subsets of sage grouse pairs. these permutations were used to test if sage-grouse within a lek complex and lek were genetically more similar than expected from a null distribution. figure s the correlation between the genomic and geographic matrices obtained from the procrustes analysis is compared to the distribution of correlations from random permutations. figure s summary results from the dapc that did not use a priori group assignments are displayed. bayesian information criterion steadily increases as the number of clusters increases in the model. for each group, a density histogram is displayed to ease in the visualization of how individuals cluster along discriminant function one. additionally, individual assignment probabilities for each lek complex are displayed in a bar graph, with a priori lek and lek complex designations displayed above. figure s results from pca using the entire snp dataset were compared to subsets of data with the top 2 % and top 3 % of fst outliers removed. any loci that were in the top or 3 % of the fst distribution for any pairwise lek comparison were removed to create the subsets. panels,, and display scatterplots of pc and pc with symbols representing lek means and segments drawn from the means to individual scores. panels,, and depict the procrustes correlation between the first two pcs and a standardized geographic matrix. panels,, and show the results of a regression with latitude predicting pc 
models incorporating a habitat resistance surface as a proxy for potentially reduced movement due to landscape features indicated that both geographic distance and habitat suitability predicted genetic structure, with no significant effect of man-made barriers to movement. lek geography predicted genetic differentiation even at a small geographic scale, which could be sharpened by strong site fidelity. relatedness was also higher among individuals within lek complexes, suggesting that reproductive skew, where few males participate in most of the successful matings, could also potentially contribute to genetic differentiation. lek complexes, and to a lesser extent individual leks, formed statistically identifiable clusters in ordination analyses, providing evidence for fine-scale geographic genetic differentiation. mating systems that reduce dispersal and lead to non-random mating might increase the potential for genetic structure to arise at fine geographic scales. finally, we illustrate how data sets containing fewer loci had less statistical precision and failed to detect the full degree of genetic structure. our results suggest that habitat features and lek site geography of sage-grouse shape fine scale genetic structure, and highlight how larger data sets can have increased precision and accuracy for quantifying ecologically relevant genetic structure over small geographic scales. we quantified population structure by analyzing variation at single-nucleotide polymorphisms in males from ten leks occurring in a small geographic region in central nevada. centrocercus urophasianusgenotyping-by-sequencinglek matingpopulation genomicspopulation structurereproductive skewuniversity of nevada, agricultural experiment stationnevada department of wildlifeu.s.
 mating systems that reduce dispersal and lead to non-random mating might increase the potential for genetic structure to arise at fine geographic scales. greater sage-grouse have a lek-based mating system and exhibit high site fidelity and skewed mating ratios. we quantified population structure by analyzing variation at single-nucleotide polymorphisms in males from ten leks occurring in a small geographic region in central nevada.  lek complexes, and to a lesser extent individual leks, formed statistically identifiable clusters in ordination analyses, providing evidence for fine-scale geographic genetic differentiation. lek geography predicted genetic differentiation even at a small geographic scale, which could be sharpened by strong site fidelity. relatedness was also higher among individuals within lek complexes, suggesting that reproductive skew, where few males participate in most of the successful matings, could also potentially contribute to genetic differentiation. models incorporating a habitat resistance surface as a proxy for potentially reduced movement due to landscape features indicated that both geographic distance and habitat suitability predicted genetic structure, with no significant effect of man-made barriers to movement. finally, we illustrate how data sets containing fewer loci had less statistical precision and failed to detect the full degree of genetic structure.  our results suggest that habitat features and lek site geography of sage-grouse shape fine scale genetic structure, and highlight how larger data sets can have increased precision and accuracy for quantifying ecologically relevant genetic structure over small geographic scales. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords centrocercus urophasianusgenotyping-by-sequencinglek matingpopulation genomicspopulation structurereproductive skewuniversity of nevada, agricultural experiment stationnevada department of wildlifeu.s. bureau of land management national fish and wildlife foundation nv energy corpuniversity of nevada renohttp://dx.doi.org/ /100000155division of environmental biology1145609division of graduate education 1447692weitzman chava l. issue-copyright-statement© the author 2016 