Environmental variation shapes genetic variation in Bouteloua gracilis: Implications for restoration management of natural populations and cultivated varieties in the southwestern United States

Abstract With the increasing frequency of large‐scale restoration efforts, the need to understand the adaptive genetic structure of natural plant populations and their relation to heavily utilized cultivars is critical. Bouteloua gracilis (blue grama) is a wind‐dispersed, perennial grass consisting of several cytotypes (2n = 2×–6×) with a widespread distribution in western North America. The species is locally dominant and used regularly in restoration treatments. Using amplified fragment length polymorphism (AFLP) and cpDNA analyses, we assessed the genetic variability and adaptive genetic structure of blue grama within and among 44 sampling sites that are representative of the species’ environmental and habitat diversity in the southwestern United States. Five cultivars were also included to investigate genetic diversity and differentiation in natural versus cultivated populations. Three main findings resulted from this study: (a) Ninety‐four polymorphic AFLP markers distinguished two population clusters defined largely by samples on and off the Colorado Plateau; (b) substructure of samples on the Colorado Plateau was indicated by genetic divergence between boundary and interior regions, and was supported by cytotype distribution and cpDNA analysis; and (c) six AFLP markers were identified as “outliers,” consistent with being under selection. These loci were significantly correlated to mean annual temperature, mean annual precipitation, precipitation of driest quarter, and precipitation of wettest quarter in natural populations, but not in cultivated samples. Marker × environment relationships were found to be largely influenced by cytotype and cultivar development. Our results demonstrate that blue grama is genetically variable, and exhibits genetic structure, which is shaped, in part, by environmental variability across the Colorado Plateau. Information from our study can be used to guide the selection of seed source populations for commercial development and long‐term conservation management of B. gracilis, which could include genetic assessments of diversity and the adaptive potential of both natural and cultivated populations for wildland restoration.


| INTRODUC TI ON
It is becoming increasingly important to understand how environmental factors shape genetic variation in plants, especially in the southwestern United States where climate change is expected to have major impacts (Archer & Predick, 2008;Gremer, Bradford, Munson, & Duniway, 2015). Such information is critical for species that are commonly used for restoration, such as grasses, which are often dominant components of ecosystems on the Colorado Plateau.
Given their broad distribution, local adaptation is likely to play a key role in how grasses will respond to climate change (Alberto et al., 2013;Davis & Shaw, 2001;Jump & Penuelas, 2005) and may be predictive of the success of cultivated seeds for restoration purposes (Hufford & Mazer, 2003;Linhart, 1995;Mijangos, Pacioni, Spencer, & Craig, 2015).
Due to anthropogenic activity, drastic and widespread changes to ecosystems have and will continue to occur (International Panel on Climate Change, 2014; Parmesan & Yohe, 2003;Vitousek, Mooney, Lubchenco, & Melillo, 1997). Many of these changes are also accompanied by uncharacteristic disturbance regimes (Dale et al., 2001;Flannigan, Stocks, & Wotton, 2000). Because local adaptation is so widespread among plant species, a changing climate is likely to influence all species, even those whose distributions span large environmental gradients (Davis & Shaw, 2001;Reusch & Wood, 2007).
It is unclear, however, whether species distributed across these gradients will persist through natural dispersal alone (Jump & Penuelas, 2005), perhaps warranting assisted migration in the near future.
Information on genetic-based adaptation in widespread grasses is particularly critical given the increased emphasis on restoration of wildlands, which are currently under accelerated change due to anthropogenic forces such as climate change, fire, and habitat disturbance. In the southwestern United States, for example, increased fire frequency in grasslands drives increased demand for seed. The average acres burned per year by wildland fires has nearly tripled since the 1980s, with an average of 7.2 million acres burned annually in the United States between 2000 and 2007, with a particularly alarming increase in the western United States (National Interagency Fire Center, 2015).
Although its efficacy is debated, reseeding is extensively used to reduce soil erosion, suppress non-native plants, and restore desirable plant communities following high-intensity wildfires (Beyers, 2004;Burton & Burton, 2002). Cost associated with reseeding burned areas has also substantially increased in recent years: The Forest Service Burned Area Emergency Response (BAER) seeding expenditure alone has reached an average of $3.3 million per year, which is nearly triple the average annual expense of the previous 30 years (Peppin, Fule, Sieg, Beyers, & Hunter, 2010). Such high costs could be more defensible by thoroughly understanding the effectiveness of restoration activities. One way to address such concerns is to investigate the degree to which source populations used for restoration purposes, including both natural and cultivated sources, are locally adapted to their environment. A first step toward investigating adaptive potential in these populations can be achieved by assessing genome-wide variability and determining the extent to which this variability is influenced by environmental variation.
A key issue for conservation management concerns the degree to which cultivated varieties maintain similar levels of genetic variability and structure as wild source populations. Although seed sources that are locally adapted have consistently been recognized for their increased success for restoration purposes (Langlet, 1971;Leimu & Fischer, 2008;Lesica & Allendorf, 1999), cultivated varieties are presumably less genetically varied than their wild counterparts.
Because of this, there is growing concern that cultivated varieties widely used in large-scale restoration can be futile or even detrimental to meeting land management objectives. Many cultivated varieties, for example, may have adapted to constant human care with frequent selection for large biomass and high seed yield. These potential cultivation pressures may lead to the loss of traits that allow for survival in a variable wildland climate (Coyne & Lande, 1985).
Cultivar varieties may outcompete or hybridize with locally adapted plants, potentially causing native populations to lose important genetic traits that enable them to respond to natural fluctuations in their local environments (Schröder & Prasse, 2013a, 2013b. Thus, understanding the extent of genetic variability in widespread species and cultivated varieties and how it is structured across broad environmental gradients can provide important clues for conservation management. Here, we assess genetic variability, structure, and the potential for local adaptation in the native grass species, Bouteloua gracilis, through the analysis of allelic variation within and among populations across a broad environmental gradient. Allelic variance was coupled with spatial and environmental variables (Manel et al., 2010;Schoville et al., 2012) to determine relationships between potentially adaptive "outlier" loci and source environment for individual populations (Beaumont & Balding, 2004;Riebler, Held, & Stephan, 2008).
To better understand how environment, cytotype, and cultivation shape genetic variation in B. gracilis, we addressed the following questions: 1. How genetically diverse is B. gracilis in the Colorado Plateau region and how is genetic variation structured across populations? Are phylogenetic relationships of these populations similarly structured? 2. Does population genetic structure in B. gracilis correlate with key climatic variables (e.g., temperature and precipitation) and how do these relate to adaptation to specific environments?
blue grama, Bouteloua gracilis, cultivars, population genetics, restoration, seed source 3. Does genetic structure and variability covary with cytotype? If so, is it linked to differentiation of functional traits as suggested for other species (Chao et al., 2014;Khazaei, Monneveux, Hongbo, & Mohammady, 2010;Manzaneda et al., 2011;Yang, Huang, Quin, Zhao, & Zhou, 2013), Madlung (2013)? 4. Are cultivated varieties of B. gracilis genetically distinct from and less genetically diverse than natural populations? We consider this to be a critical question because locally adapted seed has consistently been recognized for its increased success for restoration (Langlet, 1971;Leimu & Fischer, 2008;Lesica & Allendorf, 1999), while many cultivars have demonstrated susceptibility to agricultural pressures that lead to loss of traits that allow for survival in a variable climate (Schröder & Prasse, 2013a, 2013b.

| Study site
The Colorado Plateau is approximately 362,600 square kilometers in the Four Corners region of Arizona, Utah, Colorado, and New Mexico. It is characterized by diverse climates that range from Sonoran Desert to Alpine, with elevations ranging from 914 to 4,267 m. The region is dominated by semiarid conditions and includes the watersheds of the Colorado River and its tributaries including the Green, San Juan, and Little Colorado Rivers (Foos, 1999).
Annual precipitation is broadly varied on the Colorado Plateau.
Though the average annual precipitation across the region is 254 mm (Foos, 1999), the minimum average rainfall is 136 mm, while areas above 2,440 m in elevation receive over 508 mm and mountaintops over 3,353 m can receive nearly 1 m (U.S. Bureau of Land Management, 2014). Seasonal monsoon activity is also widespread across the region, with a significant north-south and east-west trend in precipitation seasonality (PS), with the monsoon being stronger in the southern and eastern portions of the Plateau (Adams & Comrie, 1997). Temperatures are also highly variable, with southern and lower elevation temperatures ranging from −7°C in the winter to 35°C in the summer, and mid and upper elevation temperatures ranging from −20°C in the winter to 15-25°C in the summer.
The region also undergoes multidecadal drought cycles that are associated with the Pacific Decadal Oscillation (USBLM, 2014).

| Study species
Bouteloua gracilis is a highly valued species for conservation and restoration (Herbel, Steger, & Gould, 1974) because of its broad range, adaptability, ease of establishment, and year-round forage value for livestock and wildlife (Morris, Booth, Payne, & Stitt, 1950). It is a densely tufted, C4 perennial grass that is widespread from Alberta and Manitoba, Canada, south through the Rocky Mountains and Great Plains, Midwest United States, and northwestern México (Cronquist, Holmgren, Holmgren, Reveal, & Holmgren, 1977). Blue grama inhabits a wide array of habitat types including sagebrush, salt-desert shrub, oak woodland, ponderosa pine, pinyon-juniper, prairies, and montane (Anderson, 2003) and is regularly a dominant species in the ponderosa pine and pinyon-juniper ecosystems of the Colorado Plateau.
Fragmented populations typically have few prospects for long-distance seed dispersal and may exhibit restricted gene flow (Anderson, 2003), though gene flow could be expected to remain high even in patchy environments due to long-distance wind dispersal of pollen (Llorens et al., 2016). Bouteloua gracilis is also evidently an autopolyploid across the study region, with diploid, tetraploid, and mixed cytotype sites distributed across the Colorado Plateau (Butterfield & Wood, 2015).

| Collections
Forty-four collection sites (Table 1, Figure 1) were selected as representative of the species distribution both geographically and climatically. To establish a representative distribution to capture climate variability, the species' geographic range was coupled with data obtained from WorldClim (Hijmans, Cameron, & Parra, 2005), including annual mean temperature, mean diurnal range, isothermality, temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, temperature annual range, mean temperature of the wettest/driest and warmest/coldest quarter, annual precipitation, precipitation of the wettest/driest month, PS, and precipitation of the wettest/driest and warmest/coldest quarter. This information was used to differentiate 25 climate zones within the range of B. gracilis across the Colorado Plateau and adjacent regions (using spatial-climatic stratification methods described in Doherty, Butterfield, & Wood, 2017). Twentythree of 25 zones were sampled. Approximately two sample sites within each zone were selected, and the resulting sample site distribution was geographically varied. Sample sites are identified with the zone as the prefix (0, 1, 2, etc.), followed by an alphabetic identifier (A, B, C, etc.) to differentiate multiple sample sites within a single climate zone. Sites that begin with "0" do not indicate occupation of similar climate zones, but rather a region that fell outside of delineated climate-model zones.
In addition to the 44 sampling locations of native populations, five  Note. Population ID is defined by Bouteloua gracilis-specific Colorado Plateau climate zone number (Doherty et al., 2017) and subsequent assignment of A, B, C, or D to differentiate between separate sampling occurrences within a single climate zone. Population IDs beginning with "0" are an exception and do not indicate a shared climate zone but rather indicate a site that fell outside of the expected B. gracilis range based on climate data. Latitude and longitude are formatted in decimal degrees and are projected using WGS 84. Ploidy is indicated as the proportion of sampled individuals within each collection site that is diploid (2×), tetraploid (4×), or hexaploid (6×). Ploidy fields are grayed out if no information is available for any given ploidy race at a site.
Mexico, and selected for increased seedling vigor and drought tolerance (USDA, 1982). Alma was released in 1992, traces to Lovington, Hatchita, and an experimental composite (Texas, Kansas), and was selected for seedling vigor and wide adaptability (USDA, 1992).
Samples for this study were derived from both natural sites and collections that were obtained from natural sites, but maintained locally either in a common garden or greenhouse for approximately 1 year. Thus, all samples were originally sourced from the wild, with the exception of the cultivars, which were obtained as germplasm (seed) and subsequently germinated to obtain leaf tissue. One hundred and thirteen samples were obtained from 17 natural sites.
One hundred and ten samples were obtained from a common garden representing 20 natural sites. Twelve samples from two natural sites were obtained from living specimens growing at the Northern Arizona University greenhouse in October of 2015. Forty-four samples from six natural sites were grown from seed in December of 2015. Eighty-one samples from five cultivated varieties were also grown from seed (~16 samples/cultivar). All materials were harvested while still green, desiccated upon collection in science-grade silica and stored at room temperature.
As indicated above, some samples were obtained from spec-

| Cytotype distribution and ploidy measurements
Cytotypic distribution, including diploids, tetraploids, and a few hexaploid sites, has been identified across the Colorado Plateau for B. gracilis (Table 1)

| AFLP analysis
All extracted samples were analyzed using AFLP genotyping using   (Table 2). These selective primer combinations relied on template produced from the PCR preselection amplification using primer combination M02/E01.
Electropherograms were analyzed using GeneMarker (SoftGenetics) software. Panels were manually created, and loci were automatically scored and manually confirmed as present (1) or absent (0) for each sample. GeneMarker parameters included smoothing, inclusion of fragment lengths of 50-500 base pairs, a peak height threshold of 50, local and global detection percentages set to off, and the stutter peak filter set to off. Percent error rates were calculated by comparing reproducibility and consistency among replicates as a whole and at each locus.

| Chloroplast DNA sequencing
To provide a framework for examining patterns detected by the AFLP analysis and to investigate colonization history, chloroplast DNA (cpDNA) sequences were acquired from select samples using the Electropherograms were imported into Staden (Bonfield, Smith, & Staden, 1995) to generate and assess quality of consensus sequences and subsequently aligned in BioEdit (Hall, 2004)

| Statistical analysis
Genetic distance calculations, allelic frequency measures, and analysis of molecular variance (AMOVA) were performed in GenAlEx (Peakall & Smouse, 2006). Structure (Pritchard, Stephens, & Donnelly, 2000;ver. 2.3) was used to identify statistically significant population structure through the identification of population clusters or K CLUMPP (Jakobsson & Rosenberg, 2007) was used to correct for label switching and multimodality of the data. Lastly, the CLUMPP output was graphically illustrated using Excel. The analysis was performed including and excluding outlier loci (see BayeScan below), with negligible differences detected.
The program BayeScan 2.01 (Foll & Gaggiotti, 2008) was used to detect outlier loci. BayeScan can directly estimate the probability that each locus is under selection by identifying population and marker-specific components of F ST values using a logistic regression.
A burn-in of 50,000 was used with chains of 1,000,000 iterations.
The posterior probability was set to 0.76 and above, and the false discovery rate was set to 5%. The remaining parameters were de- To assess environmental and outlier loci frequency relationships, general linear models and correlation analyses (Table 3) were conducted in R and verified using the Analysis Toolpak in Excel. Because unequal sample sizes were an issue in population and cytotype comparison, Fisher-Z statistics (Meng, Rosenthal, & Rubin, 1992) were used to test the comparison of correlations derived from independent samples. The influence of geographic distance on outlier loci was examined using Mantel and partial-Mantel tests (Legendre & Legendre, 1998;Mantel, 1967).
Molecular evolutionary genetic analysis, or MEGA version 5.1 (Tamura et al., 2011), was used for generating phylogenetic trees, including distance-based neighbor-joining and UPGMA trees to examine population genetic clustering as indicated by AFLP loci. Phylogenetic trees were generated for cpDNA sequences using maximum-likelihood, maximum parsimony, and neighbor-joining methods. For all three trees generated from the cpDNA data, 1,119 positions were identified including 24 variable sites. Positions containing gaps or missing data were omitted. Bootstrap values were calculated with 1,000 iterations (Felsenstein, 1985). For the cpDNA neighbor-joining tree, evolutionary distances were computed using the maximum composite likelihood (MCL) method (Tamura, Nei, & Kumar, 2004) and are in the units of the number of base substitutions per site. The maximum parsimony tree was generated with the subtree-pruning-regrafting algorithm (Takahashi & Nei, 2000) with search level 1, in which the initial trees were obtained by the random addition of sequences (10 replicates). The maximum-likelihood method, using cpDNA data, was based on the Tamura-Nei model (Tamura & Nei, 1993). Initial trees for the heuristic search were obtained automatically by applying Neighbor-Joining and BioNJ algorithms to a matrix of pairwise distances estimated using the MCL approach and then selecting the topology with the superior log-likelihood value. Because all three methods produced similar results, only the neighbor-joining tree is presented. When all 44 natural populations were analyzed together, the AMOVA showed that 17% of genetic variation was attributed to differentiation between the two regions: on and off the Colorado Plateau (Table 4). Approximately 13% was due to differences between populations within each region, and 70% was attributed to   Note. All AMOVA analyses were based on Nei's genetic distance calculated from AFLP data. Sampling sites were standardized at n = 5 individuals, though AMOVA runs with unequal and/or larger sample sizes did not result in significantly different results. AFLP: amplified fragment length polymorphism; AMOVA: analysis of molecular variance.

| Genetic and environmental correlations
Six loci (described above) were identified by BayeScan as outliers and showed significant correlations with one or more of the six nonredundant environmental variables (Table 3, Figure 5, Supporting Information Table S1). Genetic correlation was most strongly asso- All cultivars were originally sourced from off the Colorado Plateau and genetically grouped with the natural off-Plateau sampling sites (Figures 2 and 3). However, unlike the natural populations, no significant correlation between outlier loci and any of the environmental variables was evident among the cultivars ( Figure 5, Supporting Information Table S1). These results suggest that the cultivars may be adapted to the agricultural environment in which they have been cultivated. Alternatively, cultivars may have lost adaptive trait variation through intense selection during the cultivation process (Knapp & Rice, 1994;Schröder & Prasse, 2013a, 2013b. The different cytotypes demonstrated the most pronounced variance in response to environment ( Figure 5, Supporting Information Table S1). Tetraploids were found to occupy climates that were, on average, 50% warmer (MAT of 9°C and 6°C, respectively; p = 0.066) and nearly 40% drier than their diploid counterparts (PDQ of 48 and 78 mm, respectively; p = 0.023). These climatic differences likely drive adaptive divergence between the ploidy races, which is suggested in five genetic × environmental

| Genetic variation and structure in B. gracilis
As  -Santacruz et al., 2004;Phan et al., 2004). While these studies demonstrated higher variation among individuals within populations (88% and 98%, respectively), both studies covered relatively small sample areas, and used different primer combinations, which may account for the higher within-population variation observed in these studies.
The AFLP data clearly define two distinct groups: one on and the other off the Colorado Plateau. Extensive gene flow throughout the species' range is supported by this finding, with a distinct interruption of gene flow at the boundary, presumably in large part due to the elevational and geographic barriers presented by the Plateau itself. Although tillering is common, gene flow in B. gracilis is also facilitated through seed and pollen dispersal (Aguilera & Lauenroth, 1993;Shaw & Cooper, 1973). Pollen has been shown to have a short window of viability, with dispersal up to 1,800 m within 30 min (Copeland & Hardin, 1970), and seed dispersal through physical attachment to ungulates has been found to be critical to the dispersal of the related species, Bouteloua curtipendula (Poaceae) (Laughlin, 2003). These studies, in addition to the study presented here, suggest that gene flow via long-distance pollen dispersal can greatly contribute to consistent gene flow across continuous habitat.
While genetic structure is observed in populations on and off the Colorado Plateau, the AFLP analysis indicated additional substructure, implying a general divergence of sample sites occupying the south boundary/interior of the Colorado Plateau and sample sites occupying the north/east/west boundary of the Colorado Plateau ( Figure 3).
The divergence of the boundary and interior sites across the Plateau was mirrored in the resulting cpDNA phylogenies, where an interior Colorado Plateau clade was identified with moderate bootstrap support ( Figure 4). Further strengthening the distinction between the boundary and the interior of the Colorado Plateau is the distribution of cytotypes (Table 1) While distinct climate patterns, such as divergence in MAP (National Research Council, 2007;PRISM Climate Group, 2004) and MAT (PRISM Climate group, 2004), drive the differentiation of ecosystems on the Colorado Plateau, with drylands dominant in the interior and pinyon-juniper and montane forests dominant at and above the boundary (National Park Service, 2016), other biotic and abiotic factors may further shape the genetic variation observed in blue grama between these two regions. For example, the interior and the boundary regions of the Colorado Plateau also have markedly different topologies, soils, fire return intervals, and drought intensity patterns (Fule, Covington, & Moore, 1997;Miller & Tausch, 2001;Reynolds, Belnap, Reheis, Lamothe, & Luiszer, 2001;The National Drought Mitigation Center, 2017;Webb, Griffiths, & Rudd, 2007).
In addition to biotic and abiotic variables, migration history is also important to take into consideration when evaluating factors that may influence the genetic structure and variability of a spe- Range (Brown & Gersmehl, 1985). Colonization history and isolation by distance (IBD) may also play a role in determining how genetic variation is shaped and ultimately the detection of loci that may be under divergent selection (Lotterhos & Whitlock, 2014). Although our study did not directly test for the influence of these factors, we recognize the importance of doing so, especially within the framework of a more exhaustive sampling design that captures the maximum amount of environmental variation across a species distribution, while minimizing the potential effects of IBD and colonization history. Incorporating this approach can be challenging (Nadeau, Meirmans, Aitken, Ritland, & Isabel, 2016), but may yield additional insight into how environmental variation shapes genetic variation in widespread species, including grasses (e.g., Gray et al., 2014). Given that blue grama occupies a range well beyond our study area, future investigation of genetic × environmental correlations would benefit from more extensive sampling and assessment of divergent loci using recently developed programs that can take into account demographic history (e.g., Bayenv2; Günther & Coop, 2013).

| Correlation with environmental variables
As expected [2], B. gracilis demonstrates significant genetic covariation with climate, particularly MAT and PDQ (Table 3, Figure 5, Supporting Information Table S1). In addressing question [3], we also found evidence that genetic structure and variability covary with cytotype. While the total genetic divergence between diploids and tetraploids was less than the overall divergence between populations on and off the Colorado Plateau, the influence of cytotype was considerable when assessed in terms of the genetic response to climate (Figure 5), suggesting that genetic variation may be linked to functional trait differentiation. The AFLP data suggest that diploids (prevalent along the boundary) are more sensitive to climate gradients than are tetraploids (prevalent within the interior of the Colorado Plateau) and typically occupy cooler regions that receive more precipitation.
The occupation of blue grama diploids and tetraploids in distinct environments suggests an adaptive advantage of tetraploids across the interior of the Colorado Plateau, where the climate is markedly harsher. In a study by Manzaneda et al. (2011), cytotypes of Brachypodium distachyon (Poaceae) were distributed along an aridity gradient and tetraploids were found to have greater water use efficiency than diploids under water-restricted growing conditions.
In B. gracilis, temperature was shown to significantly affect physiologically driven traits, including Net CO 2 assimilation and Rubisco activity in plants collected from high and low elevation sites in the Rocky Mountains (Pitterman & Sage, 2000). This study suggests local adaptation of these photosynthetic-related traits, which may be functionally divergent due, in part, to genetic differentiation, and the potential influence of cytotype.
In a similar study, tetraploid populations of Allium przewalkski-  (Wu et al., 2010). This finding supports the claim that the tetraploids may have an adaptive advantage in more arid climates.
Interestingly, a common garden study of B. gracilis did not detect any effect of cytotype on measured functional traits (Butterfield & Wood, 2015). The current study, however, detected substantial differences between diploid and tetraploid populations regarding gene-environment relationships ( Figure 5). This discrepancy may be explained by the location of the common garden used in the 2015 study (Flagstaff, AZ), where the MAT is 7.9°C (potentially favorable to tetraploids) and the PDQ is 63 mm (potentially favorable to diploids). Additionally, a diploid site, several tetraploid sites, and mixed cytotype sites all occur within a 20-mile radius, lending more evidence to the suggestion that the Flagstaff, AZ location, is particularly neutral in regard to climatic pressures on the cytotypic varieties grown in the garden.
In partial support of our question on the genetic distinctness of cultivated B. gracilis [4], we found that cultivars currently available on the market were significantly different than samples derived from sites across the Colorado Plateau. Although, contrary to our prediction, the cultivars strongly reflected genetic patterns present in natural sampling sites located off of the Colorado Plateau (Figures 2   and 3), they did not demonstrate the gene × environment correlation observed in their natural counterparts ( Figure 5). The observed lack of correlation may, in part, be explained by the agricultural environment in which the cultivars are grown, whereby selective pressure is relaxed for genes (or marker regions) most responsive to environmental variation. In this case, the cultivation environment might include increased climatic stability, warmer growing conditions, irrigation, and the removal of inter-and intraspecific competition. However, we recognize that our sampling of cultivated varieties is not exhaustive and that additional surveys are needed to confirm the patterns of adaptive variation observed here. Furthermore, since natural population analogs of the cultivars were not included in this study, additional research is needed to conclusively determine the extent of cultivation effects on environmental adaptation.

| Implications for restoration
Since B. gracilis is predominantly used in restoration across the Colorado Plateau, findings from this study suggest that increased awareness of seed transfer zones and associated utility will likely be important for restoration success. For example, the distinction between the boundary and the interior of the Colorado Plateau, which was supported collectively by the AFLP data, the cpDNA data, and the cytotypic distribution of this species, suggests that the boundary and interior populations should be considered different ecotypes, consistent with a hypothesis of adaptation within different ecoregions (Hufford & Mazer, 2003). Cytotypes, due to their boundary/ interior pattern, may also play a role in assisting with selection of the right seed for restoration and should be a consideration when delineating these broad Colorado Plateau ecoregions. Ecoregions could be further subdivided by environmental variables that appear to shape genetic variation in this species (particularly MAT, MAP, PDQ, PCQ), to develop seed transfer zones for B. gracilis across the Colorado Plateau.
The selection of sample sites for this study was based on species-specific climate modeling (Figure 1) presented by Doherty et al. (2017), which also suggests that B. gracilis may be adequately represented across the entire southwest United States with as few as five climate centers. This is reflected, in part, in the genetic analysis, where genetic variation due to differences among individuals within sampling sites across the Colorado Plateau was fairly high (PhiPT = 0.15). This level of within-population genetic diversity has been documented in other wind-pollinated species, even within populations that occupy fragmented habitat (Gray et al., 2014;Llorens et al., 2016). This same model could be bolstered with genetic and cytotypic patterns across environmental gradients to help delineate the most appropriate seed transfer zones.
Unfortunately, matching seed transfer zones is not always a possibility when high-volume seed is required to restore large-scale disturbances such as wildfire. In addition to the development of robust seed transfer zones, this study also shows that the agricultural development of ecotypes specific to the Colorado Plateau is warranted, with special emphasis on a boundary and interior variety.
Methods to produce "facilitated adaptability" in germ lines can be employed from the crossing of individuals collected from diverse natural populations across broad ecoregions to develop new cultivated seed sources (Burton & Burton, 2002).
While ecotype cultivar development of B. gracilis has proven difficult in some regions (Carr & Rea, 2013), multisourced ecotype development of the species has been shown to be successful as reported in studies conducted by Fu et al. (2003) and Phan et al. (2004). Ecotypes that were generated with a greater number of individuals (99 vs. 25) and included more source populations (11 vs. 8) had significantly higher genetic diversity (Phan et al., 2004).
Further, multisourced ecotype varieties were found to be more genetically varied than the more commonly used single-sourced cultivars (Fu et al., 2003). While the ecotypes only underwent direct human-mediated selection for large seed in the source population (G 0 ), subtle though significant shifts in marker frequency in cultivated ecotypes was evident after only a single generation of breeding (G 1 ). The shift was most pronounced in the ecotype developed from fewer local sources with the fixation of some polymorphic markers arising within only a single generation (Phan et al., 2004). While overall genetic diversity remained quite high after yet another generation (G 2 ), these studies highlight that any population in cultivation, regardless of a multisource origin, is susceptible to the loss of genetic diversity. In this regard, marker identification, as was used in this study and others, could be developed to be an affordable and rapid test for recurrent genetic monitoring of cultivated varieties. If such ill-effects are recognized early, corrective measures could be more efficiently and effectively applied (Burton & Burton, 2002).

| Future directions
Additional common garden experiments would be useful in further investigating trends and resolving uncertainties presented by this study. Specifically, a Colorado Plateau garden at a field site colonized by diploids (boundary) and a field site colonized by tetraploids (interior) could substantiate differences observed between cytotypes.
A diploid and tetraploid site at off-Plateau locations could further bring to light the differences observed in this study between populations on and off the Colorado Plateau. The inclusion of currently available cultivars across sites could test their performance across different environments in direct comparison with native genotypes. Seed viability and germination as well as mortality would be important measurements to include to more fully investigate restoration success of different seed sources. A common garden study of this nature could strengthen efforts to develop seed transfer guidelines as well as assist in population selection for cultivated ecotype development.
Bouteloua gracilis plays a crucial and expanding role in the restoration of the Colorado Plateau. This study confirms that this species is genetically variable and structurally distinct on and off the Colorado Plateau and is genetically divergent across this area due to environmental and cytotypic variation. Because genetic, cytotypic, and functional trait data have all been acquired for this species, along with methodologies for spatially delineating this information (Butterfield & Wood, 2015;Doherty et al., 2017), the establishment of well-defined seed transfer zones is readily attainable. Additionally, seed currently available on the market for this species has been sourced from off the Colorado Plateau and has been demonstrated in this study to exhibit extensive genetic departure from and lack of covariation with environmental variability relative to natural populations on the Plateau. This study underscores the need for the development of new cultivated varieties to provide suitable seed for large-scale restoration that encompasses the variation observed across diverse environments on the Colorado Plateau.

ACK N OWLED G M ENTS
This research was conducted in the Environmental Genetics and Genomics Facility (https://in.nau.edu/enggen/) at Northern Arizona University and was funded by the Colorado Plateau Native Plant Program, jointly supported by the United States Geological Survey and the Bureau of Land Management. Additional research support was provided by Troy Wood who provided collection material for the study and assistance with obtaining financial support from the USGS and BLM to GJA.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
KT performed all of the molecular genetic experiments and data analysis and wrote the first draft of the manuscript. GJA oversaw the laboratory work and analyses and edited the final version of the manuscript.

DATA ACCE SS I B I LIT Y
The AFLP data generated for this study will be made available for public use on http://data/dryad.org. The regression summary of outlier loci and environmental variables is included as a Supporting Information Table S1.