Fine-scale spatial genetic structure in a locally abundant native bunchgrass (Achnatherum thurberianum) including distinct lineages revealed within seed transfer zones

Abstract

Analyses of the factors shaping genetic variation in widespread plant species are important for understanding the evolutionary history and local adaptation and have applied significance for guiding conservation and restoration decisions. Thurber's needlegrass (Achnatherum thurberianum) is a widespread, locally abundant grass that inhabits heterogeneous arid environments of western North America and is of restoration significance. It is a common component of shrubland steppe communities in the Great Basin Desert, where drought, fire, and invasive grasses have degraded natural communities. Using a reduced representation sequencing approach, we generated SNP data at 5677 loci across 246 individuals from 17 A. thurberianum populations spanning five previously delineated seed zones from the western Great Basin. Analyses revealed a pronounced population genetic structure, with individuals forming consistent geographical clusters across a variety of population genetic analyses and spatial scales. Low levels of genetic diversity within populations, as well as high population estimates of linkage disequilibrium and relatedness, were consistent with self-fertilization as a contributor to population differentiation. Variance partitioning and partial redundancy analysis (pRDA) indicated local adaptation to environment as additionally influencing the spatial distribution of genetic variation. The environmental variables driving these results were similar to those implicated in recent genecological work which inferred local adaptation for seed zone delineation. Our analyses also revealed a complex evolutionary history of A. thurberianum in the Great Basin, where previously delineated seed zones contain distantly related populations. Our results indicate evolutionary history, mating system, and differentiation across distinct geographic and environmental scales have shaped genetic variation in A. thurberianum and illustrate how numerous aspects of population genetic variation might require consideration for restoration planning.

INTRODUCTION

Identifying the factors that drive patterns of genetic variation among plant populations is important for understanding ecological and evolutionary processes and has applied significance for conservation and ecological restoration (Hedrick, 2005; Holderegger & Wagner, 2008; Sommer et al., 2013; Sork et al., 1999). The spatial distribution of genetic variation reflects evolutionary processes, including drift, migration, and selection, which shape standing variation and the evolutionary potential of populations. Therefore, quantifying spatial genetic structure and the factors shaping it can help assess the degree of population connectivity, the scale of and potential for local adaptation to environmental variation, and, consequently, the persistence of plant populations faced with environmental change (Bauert et al., 1998; Booy et al., 2000; Manel et al., 2003). Such analyses can also be used to guide conservation and restoration decisions using biologically meaningful information (Carvalho et al., 2021; Ottewell et al., 2016). During the last decade, high throughput sequencing approaches have substantially improved our ability to quantify spatial genetic variation and infer its causes across populations of ecologically significant non-model organisms (Andrews et al., 2016; Breed et al., 2019; Hohenlohe et al., 2021).

For plant species with large distributions spanning heterogeneous environments, genetic variation can be shaped by numerous factors, including geological, historical, and environmental variation, as well as life-history strategies (Hamrick & Godt, 1996; Holderegger et al., 2010). Across large geographic scales, genetic differentiation among populations can be expected as gene flow decays with increasing geographic distance and across geological barriers, commonly resulting in isolation by distance (IBD; Gavrilets et al., 2000; Hoskin et al., 2005; Wright, 1943). Environmental and ecological factors may also shape geographic genetic structure (Alvarez et al., 2009; Mosca et al., 2018; Paz et al., 2015; Storfer et al., 2010). Environmental variation can directly influence genetic differentiation by causing local adaptation and indirectly by generating isolation by environment (IBE; Shafer & Wolf, 2013; Wang & Bradburd, 2014), where gene flow is reduced across environmental gradients or selection against migrants occurs (Kawecki & Ebert, 2004). Thus, strong population genetic differentiation can occur across regions experiencing different ecological and environmental conditions (Orsini et al., 2013; Ortego et al., 2012; Wang et al., 2013; Wang & Bradburd, 2014).

Mating system can also influence patterns of population genetic structure in plants (Duminil et al., 2007; Gamba & Muchhala, 2020; Williams et al., 2001), due to variation in the frequency with which offspring are produced asexually, through self-fertilization, or via sexual outcrossing (Holsinger, 2000). Plant species capable of self-fertilization often exhibit high inbreeding, reduced gene flow and reduced within-population genetic diversity, all of which can lead to relatively pronounced patterns of population differentiation (Hamrick & Godt, 1996; Huang et al., 2021; Volis et al., 2010; Wagner et al., 2012). Thus, selfing plant species often show stronger patterns of population structure at local scales and lower genetic diversity, while outcrossing species tend to maintain higher genetic diversity and occur in less structured populations (Hamrick & Godt, 1996). Evolutionary potential can also be influenced in selfing species experiencing spatially variable selection, although this will depend on the rate of outcrossing, the strength of drift, and the mechanism of (seed or pollen) and potential for gene flow (Linhart & Grant, 1996; Trickovic & Glémin, 2022).

In complex environments, isolation and convergent evolution can result in similar but independently derived phenotypes in populations with very different evolutionary histories (Massatti et al., 2018; St. Clair et al., 2013). While common in widely distributed species, these patterns may confound restoration efforts, leading to situations where practitioners may be choosing between seed sources that either possess phenotypes that are likely to be adaptive in a given site but are distantly related to plants that used to occur there, or are more closely related to former inhabitants, but with sub-optimal phenotypes that may not survive as well in restoration sites (Leger et al., 2021). Historically, evolutionary history has been commonly considered in conservation planning for rare species (Verdú et al., 2012), while an emphasis on adaptive phenotypes has been the focus for delineating seed transfer zones for restoration of widely distributed species (Baughman et al., 2019; Pedlar et al., 2021). However, due to advances in DNA sequencing technology, it has become increasingly possible to consider evolutionary history as well as both genomic and phenotypic evidence for local adaptation in widely distributed species (Breed et al., 2019; Parchman et al., 2012). Here, we consider how evolutionary history and environmental variation shape landscape genetic structure in a widespread grass species of restoration significance in the Great Basin Desert, for which seed transfer zones have been previously inferred based on phenotypic evidence for local adaptation (Johnson et al., 2017).

The Great Basin is the most extensive cold desert in North America, with an area of about 540,000 km², and harbours significant environmental heterogeneity and biological diversity, including high levels of genetic diversity and local adaptation across extensive environmental gradients created by the complex basin and range topography (Baughman et al., 2019; Faske et al., 2021; Pilliod et al., 2017; West, 1983). However, in recent decades, the combined effects of altered fire regimes, invasive annual grasses, and human land use have led to widespread degradation and fragmentation of habitats in the Great Basin (Balch et al., 2013; Pilliod & Arkle, 2013). Cheatgrass (Bromus tectorum) and other non-native annual species have transformed native shrublands into invasive-dominated grasslands (Knapp, 1996; Nagy et al., 2021; Parkinson et al., 2013). These factors, along with climate change, are contributing to native plant species declines, especially native grasses, many of which are decreasing in abundance (Chambers & Wisdom, 2009; Svejcar et al., 2017). Overall, restoration efforts are increasing in response to global initiatives (Dudley et al., 2020; Stange et al., 2021), especially in drylands such as the Great Basin (Pilliod et al., 2017; Shackelford et al., 2021). Despite the widely recognized importance of considering spatial genetic variation and local adaptation for restoration planning (Breed et al., 2019; Hufford & Mazer, 2003; Knapp & Rice, 1994; McKay et al., 2005) and the growing number of plants in this region with phenotype-based seed transfer zones (TRM Seed Zone Applications: https://www.fs.fed.us/wwetac/threat-map/TRMSeedZoneMapper.php), we lack population genomic perspectives for most Great Basin native plants of restoration significance (but see Massatti et al., 2018).

Among the Great Basin grass species of restoration interest is Thurber's needlegrass, Achnatherum thurberianum (Piper) Barkworth, a widespread perennial bunchgrass species and an essential component of many sagebrush communities (Johnson et al., 2017). Achnatherum (Poaceae, subfamily Pooideae, tribe Stipeae) consists of large perennial grasses that grow in temperate grassland and savannah habitats (Soreng et al., 2015, 2017), many of which are thought to self-fertilize (Durka et al., 2013; Jones & Nielson, 1989; Kraehmer, 2019). Efforts to infer species relationships in Achnatherum have been difficult (Cialdella et al., 2010; Jacobs et al., 2007; Peterson et al., 2019; Romaschenko et al., 2010, 2012), and accordingly, the taxonomy of the genus has been reviewed several times (Hamasha et al., 2012). A. thurberianum has been previously classified as Stipa thurberiana (Piper) (Circ. Div. Agrostol. U.S.D.A., 1900), A. thurberianum (Piper) (Barkworth, 1993), and the more recent although not yet widely in use Eriocoma thurberiana (Piper) (Peterson et al., 2019; see the Missouri Botanical Garden's taxonomic database; http://www.tropicos.org). Achnatherum thurberianum is one of the easier species of its genus to identify, and while it is known to hybridize with Indian ricegrass (A. hymenoides) (creating A. × bloomeri), these hybrids are readily identifiable (Barkworth et al., 2007). Recent genecological work on A. thurberianum phenotypes across the Great Basin illustrated local adaptation in response to temperature and precipitation variation, which led to the formation of seed transfer zones for this species (Johnson et al., 2017). Specifically, in that study, populations from warmer and drier regions exhibited earlier flowering time and narrower leaves than those from cooler wetter regions (Johnson et al., 2017). Nonetheless, we lack knowledge of the spatial distribution of genetic variation for A. thurberianum in the Great Basin, which is important because a complex evolutionary history could have unintended consequences for seed sourcing if transfer zones are delineated based on only climate, proximity or phenotypic evidence for local adaptation (Massatti et al., 2018). More generally, analyses of evolutionary history and how landscape genetic variation is shaped by environmental, geographic and life history variation could provide a critical context for restoration planning, including refinements in the application of phenotype-based seed transfer zones.

Here, we used reduced representation sequencing to characterize spatial genetic variation across 17 A. thurberianum populations distributed across the southwestern Great Basin, in an area representing 5 of the 12 seed zones proposed by Johnson et al. (2017). We examined phylogeographic patterns and spatial genetic structure across populations, considered how geographical and environmental factors influenced them and asked to what degree existing seed zones reflected evolutionary history. We also quantified levels of genetic diversity, linkage disequilibrium, inbreeding, and relatedness within populations to evaluate the extent to which mating system might influence standing variation and fine-scale differentiation across the sampled populations. We predicted pronounced regional and fine-scale genetic structure, influenced by both mating system and local adaptation to environmental variation in this widespread grass. Further, we predicted we may find evidence of multiple lineages occupying individual seed zones, given the complexity of the basin and range topography and the potential for convergent evolution during the process of local adaptation.

MATERIALS AND METHODS Plant material

We collected leaf material from 17 localities in the western Great Basin during the Fall of 2017 from a range of 20 to 39 plants per location (Figure 1a, Table 1). We sampled these locations because they hosted multiple native species that could potentially serve as restoration seed sources for this region; collections of other species from these locations are being used for additional restoration genetic studies (e.g., Agneray et al., 2022; Faske et al., 2021). Although A. thurberianum can hybridize with A. hymenoides, hybrids are easily distinguishable (A. thurberianum has elongated ligules and pilose awns and the hybrids have a deciduous awn and relatively short ligules; Barkworth et al., 2007) and were not encountered in our sampling locations. Means for representative environmental variables for each population were obtained as described below and can be found in Table S1. In addition, these localities spanned five of the seed zones proposed by Johnson et al. (2017). Because of the complex topography in this region, note that proximate populations are not necessarily within the same seed zone (Figure 1a).

FIGURE 1. Population genetic structure of 17 populations of Achnatherum thurberianum based on 5677 SNPs. (a) Map of the sampled locations with each population code. Each population is coloured consistently in panels a, b, and c, and is represented by one of five shapes corresponding to the seed zone of Johnson et al. (2017) containing each population (the seed zone details can be found in Table 1). (b) The first and second principal components (PCs) resulting from PCA on the genotypic data plotted for each individual. (c) Each individual plotted for the first two axes from a Uniform Manifold Approximation and Projection (UMAP) clustering analysis based on the genotypic data. (d) ADMIXTURE plot representing estimated ancestry coefficients for each individual, following results of an analysis with K set to eight ancestral populations. The inferred ancestral populations (K) are additionally labelled by one of five shapes corresponding to seed zones of Johnson et al. (2017).

TABLE 1 Population genetic summary statistics for each sampled Achnatherum thurberianum population.

Population name, state	Population abbr.	N	H e	H o	F IS	π	D	${\overline{r}}_d$	I A	Seed zone
Austin Hwy, NV	AH	14	0.029	0.007	0.719*	0.039	0.487	0.283	122.350	8
Austin Summit, NV	AS	15	0.028	0.010	0.649*	0.027	0.519	0.226	93.120	5
Bald Mountain, NV	BM	18	0.025	0.008	0.671*	0.025	1.030	0.186	58.549	4
Buena Vista, OR	BV	15	0.066	0.029	0.509*	0.067	0.339	0.019	21.443	8
Dayton Hill, NV	DH	14	0.044	0.011	0.728	0.045	0.577	0.137	87.699	8
East Walker, CA	EW	15	0.028	0.011	0.704	0.028	0.307	0.314	148.414	11
Finger Rock, NV	FR	14	0.026	0.010	0.675*	0.026	0.487	0.271	106.050	7
Grey's Butte, NV	GB	14	0.060	0.024	0.522*	0.061	0.311	0.090	89.900	5
Hwy 140, NV	HO	15	0.052	0.024	0.484*	0.052	0.252	0.025	22.540	8
Jones Canyon, NV	JC	15	0.059	0.014	0.740*	0.058	0.752	0.102	86.020	11
Long Valley, CA	LV	14	0.075	0.012	0.816*	0.074	0.512	0.236	270.389	8
Modoc, CA	MD	15	0.058	0.016	0.715*	0.057	0.483	0.047	43.020	NA
Peavine Low, NV	PL	15	0.049	0.032	0.368*	0.054	0.787	0.180	121.440	5
Patagonia, NV	PT	12	0.035	0.024	0.344*	0.041	0.183	0.133	81.633	5
Smith Creek, NV	SC	15	0.030	0.010	0.611	0.029	0.536	0.329	143.590	8
Spanish Springs, CA	SS	14	0.060	0.024	0.551	0.067	0.273	0.029	29.786	5
Virginia Mountains, NV	VM	12	0.049	0.017	0.603	0.054	0.401	0.084	64.681	8

Note: For each sampling location (population), the table shows the number of individuals (N), expected (H_e) and observed heterozygosity (H_o), inbreeding coefficient (F_IS), nucleotide diversity (π), Tajima's D (D), standardized index of association (${\overline{r}}_d$), and index of association (I_A). Significant F_IS results based on bootstrap coefficient intervals are depicted with an asterisk. The seed zone for each population corresponds to those proposed by Johnson et al. (2017). The MD population has no seed zone information as it was located outside of the Johnson et al. (2017) projected area.

Library preparation, sequencing and variant calling

DNA was extracted from dried tissue of 20 plants per population using Qiagen DNeasy Plant Mini Kits and quantified with a Qiagen QIAxpert microfluidic analyzer (Qiagen Inc., Valencia, CA, USA). We constructed reduced representation libraries for Illumina sequencing using a ddRADseq method (Parchman et al., 2012; Peterson et al., 2012). The genomic DNA was digested with two restriction enzymes, EcoRI and MseI, and custom oligos with Illumina base adaptors and unique barcodes were ligated to the digested fragments (ranging from eight to 10 base pairs in length). Ligated fragments were amplified by PCR using a high-fidelity proofreading polymerase (iProof High-Fidelity DNA Polymerase, BioRad Inc.) and subsequently pooled into a single library. Libraries were size selected for fragments between 350 and 450 bp in length with the Pippin Prep System (Sage Sciences) at the University of Texas Genome Sequencing and Analysis Facility (UTGSAF). Sequence data were generated for the full set of individuals using a partial lane of sequencing with S2 chemistry on the Illumina NovaSeq platform at the UTGSAF.

We used the tapioca pipeline (https://github.com/ncgr/tapioca) to identify and discard Illumina primer/adapter sequences and potential biological sequence contaminants (e.g., PhiX, E. coli). We then demultiplexed the reads using a custom Perl script that corrects one or two base sequencing errors in barcoded regions, parses reads according to their associated barcode sequence, and trims restriction site-associated bases.

Filtered reads were clustered for variant identification and filtering with the software Stacks v 1.46 (Catchen et al., 2013). We followed a de novo assembly approach, using the “denovo_map.pl” module, which allows genotype inference by identifying SNP loci without a reference genome. The parameters were set as follows: the minimum number of identical reads required to call an allele was set to 3 (m = 3), the maximum number of mismatches between loci for individuals was set to 2 (M = 2), and the maximum number of mismatches among loci when comparing across individuals was set to 2 (n = 2). These parameters were selected through an optimization process following recommendations from Mastretta-Yanes et al. (2015) and Paris et al. (2017). In brief, we set the optimal m among values ranging from 2 to 7 (for M and n = 2) and the optimal M value among values ranging from 2 to 6 for the m optimal value (for n = M). Then, we used the “populations” module in Stacks v 1.46 (Catchen et al., 2011, 2013; Rochette et al., 2019) to extract loci that were present in at least 80% of the individuals (--r = 0.80) with a maximum observed heterozygosity of 0.65 (--max_obs_het = 0.65) and to generate and export the SNP data in vcf format for further analyses. We used vcftools v 4.2 (Danecek et al., 2011) to estimate allele frequencies, the mean depth per individual, the mean depth per site, and the proportion of missing data per site for the vcf outputs. We explored these statistics in R in order to decide the optimal m, M, and n parameters. Additionally, we filtered the obtained pool of loci using vcftools v 4.2. We allowed a maximum of 20% missing data (--max-missing 0.8), a minimum minor allele frequency of 0.03 (--maf 0.03), and specified a thin value of 5 (--thin 5), which specifies that no two sites are within the specified distance from one another. Also, we only included sites with quality scores above 10 (--minQ 10). After removing individuals that did not pass strict quality screening, the SNP dataset consisted of 246 individuals from a range of 12–18 plants per location. Trimmed fastq files for each individual are available at Sequence Read Archive (SRA; see Data Availability Statement).

Patterns of population genetic diversity and differentiation

Genetic diversity estimates were calculated using the vcf file with the R package hierfstat v 0.5–7 (Goudet, 2005). We used the “basic.stats” function to estimate mean observed heterozygosity (H_o), mean expected heterozygosity (H_e), and individual inbreeding coefficients (F_IS) within each population. Pairwise Nei's F_ST (Nei, 1987) and pairwise genetic distances (Nei's D) were estimated for all pairs of populations with the “genet.dist” function. The confidence intervals over F_IS and F_ST values were estimated using 1000 bootstraps with the “boot.ppfis” function. Additionally, we estimated nucleotide diversity (π) and Tajima's D for each population using vcftools v 4.2.

We further quantified genetic structure within and among populations with an analysis of molecular variance (AMOVA) using the R package poppr v 2.9.2 (Kamvar et al., 2014). We tested whether most genetic variances was observed among individuals within populations (i.e., no population structure) or between populations (i.e., population structure). The significance of AMOVA results was tested with the function “randtest” from the R package adegenet v 1.3–1 (Jombart, 2008) using 999 simulations. A linkage disequilibrium (LD) analysis was conducted based on the index of association (I_A; Brown et al., 1980) and the standardized index of association (${\overline{r}}_d$) over all loci to further consider the mode of reproduction within populations. Linkage disequilibrium is expected to be more pronounced in populations engaging in self-fertilization or asexual reproduction in comparison to those reproducing through outcrossing. We used the package poppr v 2.9.2 and performed the analysis using 999 permutations. To additionally consider how self-fertilization may influence population genetic structure, we estimated relatedness among individuals within each population using two different approaches that estimate kinship coefficients ranging from 0 (no relationship) to 0.50 (self). We applied the “snpgdsIBDMoM” and “snpgdsIBDMLE” functions in SNPRelate v1.30.1 (Zheng & Zheng, 2013) to estimate kinship coefficients using the method of moments (MoM) (Purcell et al., 2007) and maximum likelihood estimation (MLE) (Choi et al., 2009; Milligan, 2003), respectively. The MLE method does not assume linkage equilibrium or any particular population structure and can be used to estimate relatedness even when populations contain inbred individuals (Zheng et al., 2012; Zheng & Zheng, 2013). It should be noted that in contrast to the MLE method, the MoM method assumes marker independence, which could introduce bias in relatedness estimates when linkage disequilibrium is present (Min et al., 2022; Speed et al., 2012).

Spatial genetic structure

We tested whether populations exhibited isolation by distance (IBD; Wright, 1943) and/or isolation by environment (IBE; Wang & Bradburd, 2014) by comparing the population pairwise matrices of genetic distances (Nei's D; see above) with geographic and environmental distances (see environmental data details below) using Mantel tests. We used the “mantel” function from the R package vegan v 2.5–7 (Oksanen et al., 2013), with Spearman correlation and 9999 permutations. We estimated the geographic distances among populations as Haversine distances using the “distm” function of the geosphere v 1.5–14 R package (Hijmans et al., 2017) and the environmental distances (obtained as described below) as Euclidean distances, using the “dist” function from the R package stats v 3.3.1 (R Core Team, 2013).

Patterns of population differentiation were further assessed using model-free and model-based inference of genetic variation among individuals. First, we inferred population structure and individual ancestries without a priori information on sample origin using ADMIXTURE v 1.3.0 (Alexander & Lange, 2011). We used PLINK v 1.07 (Purcell et al., 2007) to convert the vcf file into unlinked SNPs (i.e., LD-pruned SNPs) and then ran ADMIXTURE with K values ranging from 2 to 10. The optimal value of K was estimated by evaluating cross-validation errors. We visualized ancestry variation for each model using individual based ancestry bar plots, and additionally used the R package LEA v 3.8.0 (Frichot & François, 2015) to plot admixture proportions in population-specific pie charts projected onto maps of the sampled distribution. Patterns of genetic variation across all sampled individuals were further summarized by principal component analyses (PCA; Patterson et al., 2006) using the “prcomp” function from the R package stats v 3.3.1 (R Core Team, 2013). We conducted additional PCAs on select pairs of neighbouring populations (DH and EW, VM and PL, and LV and SS) to illustrate variation among and within populations across smaller spatial scales. We also performed uniform manifold approximation and projection analyses (UMAP; Leland et al., 2018; McInnes et al., 2018) using the “umap” function from the R package umap v 0.2.7.0 (Konopka, 2020). UMAP has recently been shown to excel at detecting fine-scale spatial genetic structure of populations (Diaz-Papkovich et al., 2019, 2021). We ran UMAP with the minimum distance between points in low-dimensional space (MD) ranging from 0.1 to 0.99 and the number of approximate nearest neighbours used to construct the initial high-dimensional graph (NN) ranging from 2 to 16. Based on an assessment of clustering results across this range of parameter space (Figure S2a–h), we present results here with min_dist = 0.25 and n_neighbours = 16.

Finally, we analysed and visualized population structure using a phylogenomic approach. First, we converted the vcf file to fasta format using vcf2phylip v 2.0 (Ortiz, 2019). We trimmed the fasta alignment to exclude unreliably aligned positions with trimAl v 1.2 using the “gappyout” method (Capella-Gutiérrez et al., 2009). Then, we ran IQ-TREE v 1.6.10 (Nguyen et al., 2015) using the “Model Finder Plus” parameter (-m MFP) to determine the best substitution model (choosing the model that minimizes the BIC score), the ascertainment bias correction method (ASC; Lewis, 2001), and the ultrafast bootstrap option with 1000 bootstrap replicates (-bb 1000). We visualized the tree in Figtree v 1.4.4 (Rambaut, 2018) and indicated the corresponding subset of five seed zones from Johnson et al. (2017) for each sampled individual. We then plotted the tree linked with the population's geographical coordinates after simplifying the tree to retain one sample tip per population using the drop.tip function in the R APE package v 5.5 (Paradis et al., 2004). We used the “phylo.to.map” function in the R package phytools v 0.7–80 to plot a map (Revell & Revell, 2014), using the dropped tree and the geographical coordinates and choosing the “state” database, again indicating seed zone for each population. To quantify the extent to which seed zones reflect evolutionary history, we tested for phylogenetic signal by estimating Pagel's λ (Pagel, 1999) for the distribution of seed zones across sampled populations. We used the “fitDiscrete” function from the R package geiger v 2.0.7 (Harmon et al., 2015). To test the significance of our results, we estimated the log-likelihood if λ = 0 and if λ = 1 using the function “rescale” and conducted a likelihood ratio test.

Influence of environmental variation on spatial genetic variation

We conducted genetic-environment association (GEA) analyses using partial redundancy analysis (pRDA) to identify environmental variables that covary with genetic differentiation among populations. Climate environmental variables for each site were obtained from the PRISM database (https://prism.nacse.org) using the “get_prism_normals” function from the R prism library (Hart et al., 2015) with a data range from 1981 to 2010 and an 800 m × 800 m resolution. Following Faske et al. (2021), we converted monthly normals to estimates of potential evapotranspiration, actual evapotranspiration, soil water balance and climatic water deficit, which have been shown to effectively predict aspects of spatial and distributional variation across sagebrush steppe plant communities (Barga et al., 2018). Elevation data was acquired from the R library elevatr v 0.2.0 (Hollister & Shah, 2017). We also included several climate variables that previously predicted genecological variation and were used for seed zone delineation for A. thurberianum (Johnson et al., 2017; see Table S2 for details). Before any analyses, we examined multicollinearity among the pool of variables using the “pairs.panels” function from the pyshc v 2.1.9 R package (Revelle & Revelle, 2015), based on Pearson's |r| ≤ 0.60, to select the most orthogonal subset of variables possible. We removed all the environmental variables that were highly correlated with other variables, thus reducing the data set from 46 to 8 variables (Table 2). We then applied pRDA variance partitioning to decompose the contribution of climate, population structure, and geography in explaining genetic variation. We used three sets of variables: (1) climate environmental variables (Table 2); (2) two proxies of genetic structure (population scores along the two previously estimated PC axes); and (3) each population's coordinates (longitude and latitude). As a response variable, we used the individual-based genotypes (coded as the count of one allele, i.e., 0/1/2). We used the “rda” function from the vegan v 2.5–7 R package (Oksanen et al., 2013) for pRDA. Following Capblancq and Forester (2021), we first tested the significance of the full RDA model (with all the variables included). Subsequently, explanatory variables were added one by one, using the “ordiR2step” function of the vegan v 2.5–7 R package, with the following stopping criteria: variable significance of p < 0.01, 1000 permutations, and the adjusted R² of the global model. Then, we performed three different pRDA models: first, a model accounting for environmental variables only (conditioning the model by geography and population structure); second, a model accounting for population structure (conditioning the model by geography and environmental variables); third, a model accounting for geography (conditioning the model by population structure and environmental variables). We then compared the amount of variance explained by each pRDA to the variance of the full model (including all explanatory variables) to estimate the independent contribution of each set of variables together with any confounding effects induced by collinearity.

TABLE 2 ANOVA results show the variance explained by each environmental variable in the partial RDA, the F value, and the associated p value.

Environmental variables	Description	Df	Variance	F	p
MaxCWD	Magnitude of climatic water deficit	1	21.001	28.743	0.0009
AfRad	Folded aspect converted to radians	1	17.772	24.323	0.0009
prcpsum	PRISM total precipitation from June-August	1	16.801	22.994	0.0009
WsAETspr	Difference between water supply and AET during the spring	1	15.934	21.807	0.0009
mxtmpwin	PRISM average max temp from December-February (°C)	1	15.715	21.507	0.0009
FallAET	Difference between actual evapotranspiration summer low and fall peak	1	14.417	19.731	0.0009
SlopRad	Slope converted to radians	1	13.412	18.356	0.0009
AETgdd	Cumulative annual actual evapotranspiration during the growing season	1	11.450	15.671	0.0009
	Residuals	233	170.245

Given results from the above analyses, we conducted pRDA on the genotypic and environmental data to infer the influence of specific environmental variables on spatial genetic structure and detect the genetic signal of local adaptation and its environmental causes. We used a pRDA conditioning by population structure (PC1 and PC2) and geography (latitude and longitude) to assess whether the degree of adaptive genetic variation among individuals is explained by a particular set of environmental variables. The significance of models and RDA axes and the proportion of variation explained by each environmental variable were tested with an analysis of variance (ANOVA) and permutation (n = 999), using the “anova.cca” function of the vegan v 2.5–7 R package. Also, we used RDA to identify outlier loci potentially under selection using loadings of SNPs from the first three constrained ordination axes. We used a stringent outlier filtering threshold of 3.5 standard deviations (p < 0.0005) (Forester et al., 2018). We checked for duplicate candidate loci associated with more than one RDA axis and used Pearson's correlation (r) to identify the strongest predictor. Given the limited marker density of our SNP data, these analyses were aimed at detecting the environmental drivers of local adaptation, and less at describing the specific genetic signature of local adaptation.

RESULTS Patterns of genetic diversity and differentiation

After bioinformatic and quality filtering, we retained 5677 SNPs in a subset of 246 individuals with a mean coverage depth per individual of 21.83. Population genetic statistics were obtained for each population (Table 1). Observed heterozygosity (H_o) values were lower than the expected (H_e) in all cases, indicating reduced heterozygosity in all sampled populations. As follows, inbreeding coefficient (F_IS) values were positive for all populations, consistent with self-fertilization. F_IS values were variable across populations but were significantly positive for all populations except DH, EW, SC, SS and VM, based on bootstrap coefficient intervals. Nucleotide diversity (π) estimates were low, congruent with previous results. Moreover, all populations had positive values of Tajima's D consistent with population size contraction. The index of association (I_A) and the standardized index of association (${\overline{r}}_d$) were different from zero and significant in all cases (p < 0.001), indicating elevated levels of linkage disequilibrium consistent with self-fertilization influencing variation within populations. Pairwise relatedness estimates within populations generated with the MoM and MLE methods differed but were strongly correlated (r = 0.92; Table S3). The MoM method suggested that 82% of individuals within populations were unrelated and 14.3% had relatedness estimates between 0.4 and 0.5. For the MLE method, 65% of individuals were inferred as unrelated while 14.5% had relatedness estimates between 0.4 and 0.5 (Table S3, Figure 4). Both methods suggested considerable variability in the degree of relatedness within sample locations, with more evidence for self-fertilization in some populations (EW, FR, AH, SC; Figure 4) and far less in others (SS, BV, HO; Figure 4). Nonetheless, these results suggested the presence of both highly related and unrelated individuals, potentially reflecting both self-fertilization and outcrossing. The presence of linkage disequilibrium in our data may have affected relatedness estimation with the MoM method (Min et al., 2022; Speed et al., 2012), and could be a reason this method suggested fewer related individuals within populations. As the MLE method can be more robust to linkage disequilibrium, results of this analysis may more realistically reflect the consequences of self-fertilization in populations of A. thurberianum (Zheng, 2013; Zheng et al., 2012).

Spatial genetic structure

AMOVA analyses indicated that 79.66% (p < 0.01) of the observed genetic variance was explained by variation between populations, consistent with strong population differentiation, with the remaining 20.33% (p < 0.01) reflecting variation among individuals within populations (Table 3). Mantel tests indicated a positive association between geographic distance and genetic distance (IBD: Mantel r: 0.257, p: 0.026; Figure 3a) and between environmental and geographic distances (Mantel r: 0.238, p: 0.025; Figure 3c). No significant association was found between environmental and genetic distances based on simple population pairwise analyses (IBE: Mantel statistic r: −0.001, p: 0.437; Figure 3b).

TABLE 3 Molecular analysis of variance (AMOVA) results for the 17 populations of Achnatherum thurberianum. Here, we provide the degrees of freedom (df), sum of squares (SS), mean squares (MS), total variance (sigma) and the percentage of variance explained by each source of variation (Variance %).

Source of variation	Df	SS	MS	Sigma	Variance %
Between populations	16	203377.59	12711.09	863.60	79.66
Within populations	229	50468.42	220.38	220.38	20.33
Total	245	253846.01	1036.10	1083.99	100

Population pairwise F_ST values were significant in all pairwise comparisons (mean: 0.18, range: 0.02–0.31), even for populations that were highly spatially proximate, indicating pronounced genetic differentiation between populations (Table S4). PCA revealed three strongly separated population genetic clusters (Figure 1b), but also suggested a high degree of identifiability of individuals from most populations. The first two principal components accounted for 16.75 and 13.27% of the variation. One cluster grouped the eastern populations (AH, AS, BM, FR, SC) and the HO population, the second cluster grouped the western populations (BV, GB, JC, LV, MD, PL, PT, SS and VM), and lastly, the third cluster grouped EW and DH. PCAs on select pairs of populations illustrated variation among individuals within and between pairs of neighbouring populations (Figure S1). UMAP analyses revealed a striking fine-scale population genetic structure in which all individuals for each population clustered tightly together (Figure 1c). We note, however, that distances among UMAP clusters do not reliably represent levels of genetic differentiation among them; thus, caution is suggested when inferring evolutionary distance from UMAP results (Diaz-Papkovich et al., 2021). UMAP analyses across the ranges of the minimum distance (MD) and the number of nearest neighbours (NN) parameters produced generally consistent clustering patterns (Figure S2a–h).

ADMIXTURE identified K = 8 as the most likely number of clusters among the 17 populations sampled, based on cross-validation error values (Figures 1d, S3). Ancestry coefficient estimates for this model generally produced similar patterns to ordination results from above (as did models for smaller values of K; Figures 2, S4). Eastern populations formed an ancestry cluster (AH, AS, BM, FR, and SC), while western populations split into two different ancestry clusters (GB, JC, MD and SS, in one, and PL, PT and VM, in another). Lastly, the BV, DH, EW, HO and LV populations were assigned to additional clusters, reflecting relatively stronger differentiation of these populations in relation to those within the larger ancestral clusters above. Notably, populations from individual seed zones of Johnson et al. (2017) were commonly assigned to multiple ancestral groups in the ADMIXTURE results (e.g., populations in the seed zone represented by circles belong to three different ancestry groups; Figure 1d). ADMIXTURE models for K ranging from 2 to 5 (Figures 2, S4a–d) additionally illustrated regionally structured ancestry variation across a continuum of spatial and evolutionary scales.

FIGURE 2. Pie charts depicting the mean ancestry coefficients for each population embedded in maps of the study region, with K ranging from 2 to 5 ancestral populations (a–d). The individual ancestry coefficients for each population studied, along with the pie charts showing the mean admixture proportions ranging from K = 2 to K = 5, are represented in Figure S4.

All main clades yielded bootstrap branch support higher than 99 in the maximum likelihood tree, which produced a topology illustrating pronounced population structure. Similar to the patterns of clustering in PCA and ADMIXTURE analyses above, the topology resolved four main clades: (1) the eastern populations (AH, AS, BM, FR and SC) with HO, (2) the BV population, (3) the EW and DH populations, and lastly, (4) the western populations (GB, JC, MD, SS, LV, PL, PT and VM) (Figure 5a,b). Individuals from the same populations predominantly clustered together in the maximum likelihood tree, further illustrating population differentiation at fine spatial scales. Consistent with ancestry-based analyses above, there was no evidence of phylogenetic signal for seed zones (λ = 0.000); populations from the same seed zones often appeared in multiple distantly related clades (e.g. DH and EW or PL and VM; Figure 5a,b).

Genetic-environment association analyses

The pool of environmental variables was reduced from 46 to 8 after removing highly correlated variables (based on Pearson's |r| ≤ 0.60). The specific environmental variables used for seed zone delineation in Johnson et al. (2017) were not included in the analyses due to multicollinearity, but some of them were correlated with the eight variables selected from our analyses (see Table S2). Results from the pRDA suggested that specific environmental variables may influence spatial patterns of genetic variation. In particular, the climatic variables explained 23% of the total genetic variation (39% of the variance explained by the full model), suggesting an association between genetic variation and environmental gradients (IBE) (Table 4). The environmental variables with the greatest explained variance in the pRDA were the magnitude of climate water deficit (MaxCWD; 21.01%, p < 0.001) and the folded aspect (AfRad; 17.80%, p < 0.001) (Table 2, Figure 6). The first two constrained axes were significant (p < 0.001), explaining 24.63% and 23.44% of the total variation (Figure 6). We identified 10 outlier loci associated with environmental variation across the second and third RDA axes (Table S2). Four of these loci were associated with the difference between water supply and actual evapotranspiration during the spring (WsAETspr), two with the maximum temperature during the winter season (mxtmpwin), and the four remaining were associated with the cumulative annual actual evapotranspiration during the growing season (AETgdd), the difference between actual evapotranspiration summer low and fall peak (FallAET), the magnitude of climatic water deficit (MaxCWD) and the slope (SlopRad). One of the variables, the maximum temperature during the winter season (mxtmpwin; Table 2), was highly correlated with the mean average temperature (MAT; Pearson's |r| = 0.82), which was among the variables most strongly predicting genecological variation in Johnson et al. (2017).

TABLE 4 The influence of climate, genetic structure, and geography on genetic variation decomposed with partial RDA. The proportion of explainable variance represents the total constrained variation explained by the full model.

Partial RDA models	Variance	R 2	p value	Proportion of explainable variance	Proportion of total variance
Full model: F ~ clim. + geog. + struct.	276.21	0.59	0.001	1	0.59
Pure climate: F ~ clim. | (geog. + struct.)	110.23	0.23	0.001	0.39	0.23
Pure structure: F ~ struct. | (clim. + geog.)	49.93	0.10	0.001	0.17	0.10
Pure geography: F ~ geog. | (clim. + struct.)	19.53	0.04	0.001	0.06	0.03
Confounded climate/structure/geography	96.52			0.33	0.19
Total unexplained	189.92
Total variance	466.13

DISCUSSION

Understanding the distribution of genetic variation in native plants is crucial not only for understanding the origin and maintenance of diversity but also for conserving and restoring populations. Our analyses of population genetic variation in Thurber's needlegrass (A. thurberianum), a widespread bunchgrass in the Great Basin, illustrated strong regional differentiation as well as remarkably fine-scale genetic structure among populations. These patterns, along with reduced heterozygosity, elevated linkage disequilibrium, and evidence of high relatedness estimates in some populations, are consistent with self-fertilization reducing genetic diversity and contributing to population structure. Despite low genetic diversity within the sampled populations, our analyses indicated that environmental variation has shaped spatial genetic structure and influenced local adaptation. Our results suggested the potential for local adaptation driven by climate water deficit, aspect and summer precipitation, as well as strongly correlated variables that were excluded from analyses. Notably, our results illustrated previously unidentified differences in evolutionary history within seed zones proposed for A. thurberianum based on phenotype-environment associations (Johnson et al., 2017). Altogether, our findings suggest that numerous factors have shaped the genetic structure of A. thurberianum across fine geographic and environmental scales and provide a landscape genomic perspective that could allow managers to consider both phenotypic variation and evolutionary history when making restoration seed source decisions.

Strongly congruent results across multiple analyses (Figures 1, 2, 5) illustrated spatial genomic structure across both broad geographic regions and among proximate populations at finer scales than might otherwise be expected for a primarily wind-dispersed species (Linder et al., 2018). As discussed in more detail below, our analyses consistently illustrated pronounced differentiation among eastern and western groups of populations (Figures 1, 2, 5) including differentiation of two populations at the southwestern limit of sampling. Interestingly, the HO population clustered with the eastern populations, despite being geographically distant and in closer proximity to populations in the northwestern portion of the sampled area (e.g. JC; see Figure 1a,b). ADMIXTURE models for K ranging from 2 to 5 (Figures 2, S4a–d) illustrated patterns of subsequently finer regional differentiation across a likely continuum of evolutionary and spatial scales (Figures 2, S4). The broad patterns of regional differentiation could be explained by historical habitat contractions during the multiple Pleistocene glacial cycles that shaped the Great Basin landscape (Beck & Jones, 1997), in addition to environmental and topographical complexity of the region. During the Last Glacial Maximum and throughout the last deglaciation, the Great Basin region was marked by the formation of multiple extensive lacustrine systems (Lyle et al., 2012; Tchakerian & Lancaster, 2002), which would have reduced the availability of suitable habitat for A. thurberianum. More specifically, a large portion of our study area was historically occupied by Lake Lahontan (38.75–40.75°N, 117.5–120.5°W; Russell, 1885; Matsubara & Howard, 2009), which may have played a central role in shaping regional patterns of marked differentiation in A. thurberianum. Consistent with such history, positive values of Tajima's D indicated past population contraction (Table 1; Tajima, 1989). Temperate grass species of Europe (Blanco-Pastor et al., 2019, 2021; Hensen et al., 2010), Africa (Mairal et al., 2021), and North America (Avendaño-González et al., 2019; Palacio-Mejía et al., 2021) demonstrate comparable regional patterns of differentiation and gene flow also influenced by Quaternary glacial cycles.

In addition to the larger-scale geographic patterns, genetic differentiation across smaller scales was evident in phylogenetic analyses, where individuals clustered largely by population (Figures 5a, S5, Table S5), and in the UMAP analyses, where individuals formed population-specific clusters (Figure 1c). UMAP can excel at detecting fine-scale population structure (Diaz-Papkovich et al., 2019, 2021) and here provided remarkable resolution, with all A. thurberianum individuals having 100% identifiability to their population of origin. While parameter settings can influence the UMAP depiction of clustering (MD and NN; Diaz-Papkovich et al., 2021), our analyses across a range of two parameters produced largely congruent results with minor variation in the strength of clustering (Figure S2). These results suggest that the degree of differentiation among populations seen here could possibly be used to retroactively identify the constituents of historically seeded populations with a high degree of certainty. While gene flow among spatially proximate populations can be high in some wind-dispersed grasses (Stritt et al., 2022; Vogel et al., 2009), pronounced spatial genetic structure has also been reported in Eurasian species of Stipeae (Durka et al., 2013; Wagner et al., 2011) and could be additionally driven by environmental variation in A. thurberianum. Patterns of population differentiation and identifiability across both large and small geographic scales indicate that genetic variation in A. thurberianum has been shaped by a combination of historical isolation, local adaptation to environment and life history variation.

Species that employ self-fertilization can be prone to reduced within-population diversity and pronounced among-population differentiation due to reductions in effective population size and effective migration (Hamrick & Godt, 1996; Honnay & Jacquemyn, 2007; Huang et al., 2021; Volis et al., 2010). The populations studied here were characterized by low levels of heterozygosity and diversity (mean H_e: 0.04; mean H_o: 0.02; mean π: 0.05), strongly positive inbreeding coefficient (F_IS) estimates, and a high degree of linkage disequilibrium (${\overline{r}}_d$ and I_A) (Table 1), all of which are consistent with a mating system that includes self-fertilization.

The distribution of relatedness estimates varied across sampled locations and methods, but generally suggested the presence of both highly-related and unrelated individuals (Figure 4). This could suggest the occurrence of both self-fertilization and outcrossing, which would be consistent with the observed pattern of IBD across populations (Figure 3). While studies directly assessing the mating system of A. thurberianum are lacking, closely related species from the Poaceae family have been generally described as self-fertilizing (Arnesen et al., 2017; Jones & Nielson, 1989; Marques et al., 2017; Stritt et al., 2022), and grass species employing self-fertilization have been commonly documented to have low genetic diversity and pronounced spatial genetic structure (Dell'Acqua et al., 2014; Guo et al., 2017; Marques et al., 2017). Indeed, past studies of related species in Eurasia (from the tribe Stipeae) based on traditional molecular markers (i.e., SSRs, AFLP, single mtDNA loci) indicated pronounced population differentiation and low diversity consistent with self-fertilization and demographic processes shaping population differentiation at small spatial scales (Durka et al., 2013; Wagner et al., 2012).

FIGURE 3. Mantel test results and relationships between genetic (Nei's D) and geographic distances (km) (a), genetic (Nei's D) and environmental distances (b), and geographic (km) and environmental distances (c).

Although mating system and low genetic diversity have allowed genetic drift to influence spatial genetic structure, environmental variation also appears to have contributed to genetic variation through local adaptation in A. thurberianum. Studies on Eurasian species from the Stipeae tribe have similarly indicated the potential for climate to shape local adaptation and population genetic structure (Gao et al., 2018; Hamasha et al., 2013; Schubert et al., 2019). The environmental variables influencing genetic variation in A. thurberianum are known to predict genetic and phenotypic variation across populations of other plant species (Barga et al., 2018; Dilts et al., 2015; Faske et al., 2021) and have been commonly implicated as underlying local adaptation in genecological studies of Great Basin plant species (Baughman et al., 2019; Johnson et al., 2017). Variance partitioning with partial RDA illustrated that neutral genetic structure, geography and environmental variation together explained a substantial proportion of genetic variance (Table 4). However, the individual contribution of environmental variation was substantially stronger than that of both genetic structure and geography (Table 4). Variation in climate variables (Table 2), especially climate water deficit, aspect, and summer precipitation, explained a substantial proportion of genetic variation and influenced spatial genetic structure (Table 2, Figure 6). The pronounced environmental heterogeneity of the Great Basin (Figure 3c) could commonly give rise to IBE, where environmental transitions act as barriers to gene flow or where local adaptation to different environments leads to low migrant fitness (Shafer & Wolf, 2013; Wang & Bradburd, 2014). As we lack extensive phenotypic data for the specific populations analysed here, we cannot directly assess the degree to which phenotypic and genetic data implicate similar environmental variables in shaping local adaptation. Nonetheless, genomic evidence for environment contributing to spatial genetic structure and local adaptation here is strongly consistent with results from a common garden study of phenotypic variation across 66 Great Basin A. thurberianum populations that inferred local adaptation as a basis for delineating seed zones (Johnson et al., 2017).

Self-fertilization has long been thought of as a limit to adaptation as it reduces heterozygosity, reduces recombination, and can allow deleterious mutations to accumulate (Charlesworth, 2006; Igic & Busch, 2013; Stebbins, 1957; Wright et al., 2013). However, self-fertilization could also indirectly facilitate local adaptation by limiting gene flow among populations inhabiting different environments (Hereford, 2010; Linhart & Grant, 1996; Trickovic & Glémin, 2022). These contrasting evolutionary pressures may be one reason why past reviews and meta-analyses have not detected clear evidence for mating system predicting the degree of local adaptation in phenotypic studies of plants (Hereford, 2010; Leimu & Fischer, 2008). Despite low diversity, high linkage disequilibrium and evidence for highly related individuals within our sampling locations, variation among individuals and evidence for an abundance of unrelated individuals within some populations (Figure 4) suggest the possibility that evolutionary potential could be maintained in A. thurberianum by occasional outcrossing. More generally, despite the potential for self-fertilization to limit evolutionary potential, empirical evidence for local adaptation in species with strongly reduced diversity has not been uncommon (e.g. Dlugosch & Parker, 2008; Gasca-Pineda et al., 2020; Verhoeven et al., 2008), and theoretical work has identified conditions under which selfing may favour local adaptation across heterogeneous environments (e.g. Hodgins & Yeaman, 2019; Trickovic & Glémin, 2022).

FIGURE 4. Relatedness among individuals varied across sampling locations as estimated using two different methods and illustrated in dot plots. These analyses provide information about the relationship of any pair of individuals by assessing their kinship coefficient, which ranges from 0 (no relationship) to 0.50 (self). For each plot, populations are sorted along the X-axis from highest to lowest mean relatedness estimates, and dots are jittered to avoid overlap. (a) Relatedness estimation based on maximum likelihood (MLE; Milligan, 2003; Choi et al., 2009). (b) Relatedness estimation based on the method of moments (MoM; Zheng & Zheng, 2013).

Given the influence of environmental variables on genetic variation observed here and in Johnson et al. (2017), along with the topographically and environmentally heterogeneous nature of the Great Basin, we suspected that proposed seed zones might span populations from distantly related clades. Indeed, geographically differentiated groups of populations commonly crossed multiple seed zones, and specific seed zones encompassed populations from distantly related clades (Figures 1, 5). It is worth noting that our sampling focused on an area that was not well represented in the sampling design for the phenotype-based seed zone development (Johnson et al., 2017), which might indicate that caution is warranted when designating seed zones outside areas of extensive sampling. Additionally, our sampling covered only some of the seed zones inferred by Johnson et al. (2017), and more thorough population genetic sampling would improve understanding of the discordance between evolutionary history and seed zones. Nonetheless, convergent evolution to parallel environmental conditions across divergent lineages is well documented across diverse groups of plants (e.g. Rellstab et al., 2020; Xu et al., 2020), and this possibility should be considered for seed zone design in the Great Basin, especially if particular geographic regions show consistent barriers to gene flow for multiple species. Indeed, other plant species from the Great Basin have been found to have similar patterns of population genetic structure in the vicinity of historic Lake Lahontan (e.g. Faske et al., 2021). Convergent evolution to environment could have consequences for restoration if outbreeding depression (reviewed in Edmands, 2007) or genetic incompatibilities (Etterson et al., 2007) occur in seed mixes containing distantly related populations. As discussed above, the strong east/west division observed in our sampling area may be an impetus to keep seeds from these regions segregated during ecological restoration, even if they share similar phenotype-based seed zones.

FIGURE 5. Population phylogenetic differentiation of the 17 A. thurberianum populations. (a) Maximum likelihood topology inferred with IQ-TREE. The scale bar represents the expected number of nucleotide substitutions per site. The Operational Taxonomic Units (OTUs) correspond to individual plants. Ultrafast bootstrap support values (UFBS) are indicated at the nodes; bootstrap values higher than 95 are depicted with an asterisk. (b) Projected phylogeny onto the geographic map showing each population's location. The tree was pruned to one individual per population to improve visualization. In both panels, colours correspond to populations and the different shapes at the tips correspond to the individual seed zones of Johnson et al. (2017) (Table 1).

FIGURE 6. Population genetic structure was predicted by environmental variation which predicted allele frequency shifts at loci potentially under selection. Redundancy analysis (RDA) plot depicting the association of environmental variables with genotypic variation for each individual, with colours representing populations as in the other figures. The direction and length of arrows correspond to the loadings of each variable on the two RDA axes.

In addition to its relevance for seed zones, a general understanding of genetic diversity, population differentiation, and local adaptation in A. thurberianum could have utility for guiding provenance and seed-sourcing decisions (Breed et al., 2019; Mijangos et al., 2015). Consistent results across multiple analyses (Figures 1, 2, 5) suggested a general pattern of historical divergence among three differentiated groups, which are then further sub-structured across finer geographic and environmental gradients. Although such structure would suggest relevant scales and processes for guiding restoration in A. thurberianum, it is important to also consider that our work focused on a subset of its Great Basin distribution. Populations exhibited very low levels of standing variation, which presumably contributed to population differentiation across fine spatial scales. Low diversity can be a concern when sourcing seeds for restoration, as genetic diversity can be a proxy for evolutionary potential (Ellstrand & Elam, 1993; Willi et al., 2022), and some have proposed mixing populations of low-diversity species as a way to increase the chances of long-term persistence of restored populations (Bischoff et al., 2010; Bucharova et al., 2019; St. Clair et al., 2020). Mating system could thus impact another decision increasingly common in restoration, which is whether to combine collections from multiple populations, either to deliberately increase diversity or as a practical decision when there are simply not enough seeds available from one population (St. Clair et al., 2020). Concern about the potential for outbreeding depression when combining populations (Hufford et al., 2012; Templeton, 1986) could be reduced for self-fertilizing species (Johnson et al., 2017). Despite low diversity, our analyses indicated that environmental variation has shaped spatial genetic structure and influenced local adaptation across A. thurberianum populations of the Great Basin. The consistency between the inference of local adaptation in our study and phenotypic analyses of genecological experiments (Baughman et al., 2019; Johnson et al., 2017) in this species highlights the utility of population genomic analyses for characterizing the environmental variables contributing to local adaptation.

ACKNOWLEDGEMENTS

This project was funded by grant 2017-67019-26336 from the USDA National Institute of Food and Agriculture Sustainable Agroecosystems program and Cooperative Agreements L16AC00318 and L19AC00013 from the Bureau of Land Management through the Great Basin CESU to EAL and TLP. Illumina S2 NovaSeq sequencing was performed by the Genomic Sequencing and Analysis Facility at UT Austin, Center for Biomedical Research Support (RRID#: SCR_021713). We benefitted from intellectual, editorial and analytical discussion during the completion of this research from Trevor Faske, Joshua Jahner, Kathryn Uckele and Joshua Hallas. We would also like to thank Owen Baughman, Shannon Swim, Meagan O'Farrell, Sage Ellis and Lana Sheta for assistance with collecting material in the field, and Lana Sheta and Shelby Burdo for assistance with DNA extractions.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The trimmed vcf file and scripts used for analyses can be found at the Dryad Digital Repository: https://doi.org/10.5061/dryad.pvmcvdnpn. The raw data from this project were submitted to NCBI Sequence Read Archive (SRA) and can be found by the BioProject ID PRJNA849003. The individual fastq files for each population can be found under the following accession numbers: AH (SRR19646741), AS (SRR19646740), BM (SRR19646732), BV (SRR19646731), DH (SRR19646730), EW (SRR19646729), FR (SRR19646728), GB (SRR19646727), HO (SRR19646726), JC (SRR19646725), LV (SRR19646739), MD (SRR19646738), PL (SRR19646737), PT (SRR19646736), SC (SRR19646735), SS (SRR19646734) and VM (SRR19646733).