The demographic and ecological factors shaping diversification among rare Astragalus species

Evolutionary radiations are central to the origin and maintenance of biodiversity, yet we rarely understand how they are jointly shaped by demography and ecological opportunity. Astragalus is the largest plant genus in the world and is disproportionately comprised of rare species restricted to narrow geographic and ecological regions. Here, we explored the demographic and ecological mechanisms underlying patterns of diversification in a threatened Astragalus species complex endemic to a small desert region in the western United States.


| INTRODUC TI ON
Evolutionary radiations, or the rapid diversification of clades, are central to the formation of new species and their coexistence across space and time (Gittenberger, 1991;Simões et al., 2016). A central model to emerge after nearly 100 years of research attributes rapid diversification to adaptation across unoccupied ecological niches (Schluter, 2001;Stroud & Losos, 2016) that often occurs in the face of substantial gene flow (i.e., adaptive radiation; Simões et al., 2016).
However, a growing number of studies suggest that many radiations are primarily driven by demographic factors, including population fragmentation, bottlenecks and restricted gene flow (Kozak et al., 2006;Prunier & Holsinger, 2010;Simões et al., 2016) and, in some cases, by intrinsic genomic characteristics (e.g., polyploidization; Van de Peer et al., 2017). In order to fully understand the mechanisms underlying evolutionary radiations, we must jointly examine the influence of demographic history and ecological opportunity on species diversification and persistence (Foote et al., 2016;Gavrilets & Vose, 2005).
As the efficacy of natural selection is positively related to long-term population size (Kimura, 1983), historically large populations may be more likely to rapidly adapt to and exploit ecological opportunities.
Yet, evolutionary radiations are often associated with population bottlenecks linked to the colonization of a novel environment (Foxe et al., 2009;Särkinen et al., 2007;Wessel et al., 2013). In rare cases, genetic drift associated with founder events may actually trigger ecological diversification (Wessel et al., 2013, although see Barton & Charlesworth, 1984). However, if bottlenecks are severe, populations may lose substantial standing genetic variation (James, 1971), potentially limiting their ability to rapidly adapt to and persist in novel environments (Frankham et al., 1999). Gene flow between populations or species in early stages of colonization may counteract the effects of founder events by increasing adaptive genetic diversity (Dlugosch & Parker, 2008). Recent genomic studies of classic radiations (e.g., Darwin's finches [Lamichhaney et al., 2015] and East African cichlids [Meier et al., 2017]) suggest gene flow played a crucial role in seeding genetic variation that ultimately facilitated adaptive divergence. On the other hand, local adaptation is likely to be swamped by the influx of maladapted alleles if rates of gene flow exceed a certain threshold (Wright, 1932), which may prevent the establishment of stable populations. Understanding the dynamics of genetic drift and gene flow during diversification therefore remains a fundamental research goal (Lawton-Rauh et al., 2007).
Spatial and temporal environmental heterogeneity facilitate the generation and maintenance of biodiversity by increasing opportunities for populations to experience aberrant climatic conditions (Rainey & Travisano, 1998;Williams & Jackson, 2007) or geographic isolation (Hewitt, 2000;Knowles & Massatti, 2017).
The intermountain region of western North America (hereafter referred to as the Intermountain West) is a spatially complex landscape with a history of dramatic climate changes (e.g., repeated glaciations altering regional climatic gradients throughout the Pleistocene) and therefore is an important system for investigating how demographic and ecological conditions interact to shape evolutionary radiations (Egan & Crandall, 2008;Tidwell et al., 1972).
The Intermountain West hosts an abundance of endemic plant clades (e.g., Weber, 2003) that span sharp climate and soil gradients (Tidwell et al., 1972). Diversification of many of these plant clades has been linked to temporal environmental heterogeneity via glacial cycles (e.g., Egan & Crandall, 2008;Levsen et al., 2012). Moreover, ongoing environmental changes (e.g., drought, grazing, wildfires and the spread of invasive species) continue to threaten native species and communities (Litt & Pearson, 2011;Winkler et al., 2019;Zedler et al., 1983). Uncovering the precise mechanisms that shape biodiversity in this region therefore likely has broader implications for evolution and conservation during periods of extreme environmental change (Thorpe et al., 2015).
To test how demographic and ecological processes affect evolutionary diversification, we combined population and landscape genomic analyses to study a rare group of Astragalus species narrowly distributed across a dryland ecosystem in the Intermountain West (see Methods for study system details). Astragalus (Linnaeus 1753; family Fabaceae) is the most species-rich plant genus in the world, comprising ~3,000 species mostly distributed across temperate regions of the Northern hemisphere, including at least 160 species restricted to the Intermountain West (Barneby, 1989).
Interestingly, a disproportionate number of Astragalus species are rare and endemic to narrow geographic or ecological regions, with approximately one third of extant species considered either vulnerable, endangered or critically endangered (IUCN, 2020;Rundel et al., 2015). Despite widespread interest in Astragalus evolution and conservation (Sanderson & Wojciechowski, 1996;Scherson et al., 2008;Stebbins, 1981;Wojciechowski et al., 1999), the lack of detailed genomic studies of closely related species has hindered our understanding of the mechanisms underlying patterns of diversification, endemism and rarity in this group (Rundel et al., 2015; although see Záveská et al., 2019).
The broader Astragalean clade does not appear to possess defining ecological traits (i.e., key innovations), and thus, Sanderson and Wojciechowski (1996) hypothesized that diversification within this group may be better explained by demographic factors, such as limited gene flow and population fragmentation, as opposed to ecological opportunity. However, the apparent prevalence K E Y W O R D S conservation, environmental change, glaciation, natural selection, Pleistocene, population genetics, speciation of ecological specialization (particularly edaphic specialization) among Astragalus suggests that local adaptation may play a role in promoting the divergence and persistence of lineages (Rundel et al., 2015), although genetic evidence for local adaptation within Astragalus is lacking. In either case, environmental change may be predicted to play an important role in establishing the demographic or ecological conditions for divergence between Astragalus taxa. Here, we test predictions that rare and endemic Astragalus species (a) diversified during periods of major environmental change, (b) diversified largely in the absence of gene flow (according to Sanderson & Wojciechowski, 1996) and (c) are locally adapted across their narrow distributions. To test these predictions, we performed double-digest restriction-site associated DNA (ddRAD) sequencing and inferred patterns of population genetic structure, genetic diversity and inbreeding. Using these foundational data, we reconstruct the history of population size changes, gene flow and the timing of diversification (testing Predictions a and b) and examine genome-wide signatures of local ecological adaptation (testing Prediction c). Our analyses demonstrate a complex demographic and ecological history underlying Astragalus diversification, which we discuss in the context of the evolution and conservation of rare and endemic plants inhabiting heterogeneous landscapes.

| Study system
We studied two previously described Astragalus species in south- Welsh (stage milkvetch; hereafter vehiculus). These taxa were originally described primarily on the basis of relatively minor morphological differences in flower size, flower colour and pod size (Welsh, 1974(Welsh, , 1998. For instance, pods are generally smaller in iselyi (25-32 mm long; Welsh, 1974) compared to vehiculus and sabulosus (20-48 mm long; Welsh, 1998). Astragalus species are reportedly all monoecious and some species show a high degree of self-compatibility (Karron, 1989), although the degree of self-compatibility is unknown among iselyi, sabulosus and vehiculus. Seed dispersal and pollination vectors in these taxa are also unknown, though research from other Astragalus species suggests seed dispersal is primarily gravity-or water-mediated (Morris et al., 2002;Vicente et al., 2011) and that various insects (particularly bees) serve as pollinators (Green & Bohart, 1975;Karron, 1987;Sugden, 1985). Each taxon is restricted to particular geological formations and altitudinal zones: iselyi occurs on seleniferous and uranium-rich soils of the Morrison, Paradox and Mancos formations between elevations of 1525-2010 m (Goodrich et al., 1999), sabulosus is restricted to the Mancos and Morrison formations between 1,300 and 1,600 m, and vehiculus is restricted to the Morrison formation between 1,370 and 1,465 m (Welsh, 1998).
In addition to their small ranges (Figure 1a), recent censuses revealed small population sizes (N c ; iselyi N c = 3,192, sabulosus N c = 4,692, ve-hiculus N c = 1,831; Rita Reisor pers. comm.). Moreover, in addition to climate change, all taxa appear vulnerable to mining, oil and gas development, grazing, and off-highway vehicle use (Franklin, 2003) and are being petitioned for listing under the United States Endangered Species Act. We know of no published genetic or ecological studies in this group. Cytological studies suggest that polyploidy is extremely rare in New World Astragalus (Wojciechowski et al., 1993); thus, we assume these species are diploid. F I G U R E 1 (a) Sampling localities of iselyi (n = 5; purple), sabulosus (n = 6; teal), and vehiculus (n = 2; maroon) populations in southeastern Utah. (b) The first two principal components (PC1 = 13.6% of variation; PC2 = 6.5% of variation) of a PCA showing three genetic clusters corresponding to iselyi, sabulosus, and vehiculus and their corresponding pairwise
At sites with fewer than 15 individuals, we sampled leaves from every individual (separated by at least 5 m). When more than 15 individuals were present at a site, we sampled a representative subset of as many as 18 individuals spaced up to 100 m apart. Leaf tissues were stored in silica desiccant prior to DNA extraction.
DNA was extracted using DNeasy 96 Plant Kits (Qiagen) following the manufacturer's protocol. Genomic DNAs were individually barcoded and processed into ddRAD libraries (Peterson et al., 2012).
Briefly, DNA was digested with EcoRI and MspI restriction enzymes, followed by the ligation of Illumina adaptor sequences and barcodes.
Ligation products were pooled and amplified using 18 polymerase chain reaction cycles. A Pippin Prep (Sage Science) was used to select amplicons between 400 and 600 base pairs (bp). The resulting library was sequenced on a HiSeq 4000 (Illumina) at the University of Oregon's Genomics and Cell Characterization Core Facility to generate single-end 100 bp reads.

| Genomic data preprocessing
We processed and aligned raw sequence data with Stacks v2.41 (Rochette et al., 2019). Reads were removed if they contained uncalled bases or if the average phred-scaled quality score fell below 22 within a sliding window of 15 bp. The last 5 bp were then removed from high-quality reads. We followed the procedure outlined by Rochette and Catchen (2017) to determine appropriate parameter settings in Stacks for de novo assembly. Briefly, we processed a test data set of 12 high-coverage samples (five sabulosus, two vehiculus, five iselyi) while varying the parameters M (the maximum nucleotide distance allowed between stacks) and n (the number of mismatches allowed between loci during catalog construction) from 1 to 9 and while keeping M = n and setting m = 3 (minimum depth of coverage to create a stack). Based on these tests, we chose to process the full data set with M = 3 and n = 3 because the total number of loci retained in 80% of individuals (r80) and the distribution of the number of single nucleotide polymorphisms (SNPs) per locus stabilized at these values (see Figure S1). However, we also tested several other filtering parameters to determine the robustness of our data set to filtering criteria (see Table S1). Using vcftools v1.14 (Danecek et al., 2011), we filtered out variant sites represented by less than 50% of individuals and with a Hardy-Weinberg departure p-value < .0001. Singletons were also excluded because they were more common in individuals used to generate the RAD loci catalog (cstacks program in Stacks), which is expected (Rochette & Catchen, 2017). Finally, we generated a separate "unlinked" data set, consisting of only one randomly selected SNP per locus, to minimize the number of non-independent SNPs.

| Estimating population structure and diversity
Using our unlinked data set, we inferred genetic structure patterns with a principal components analysis (PCA) in EIGENSOFT v6.1.4 (Patterson et al., 2006) and a population cluster analysis in fast-Structure v1.0 (Raj et al., 2014). For fastStructure, we tested K values (representing the number of population clusters) from 1 to 8 and selected the model with the highest log-marginal likelihood lower bound. Based on results from these analyses, we used vcftools to calculate pairwise genetic differentiation (F ST ) between genetic clusters both genome-wide and per-locus. For each genetic cluster, we calculated mean inbreeding coefficients (F IS ) and nucleotide diversity (π) using the populations program in Stacks. To determine whether any genetic cluster represents an admixed or hybrid population, we performed F 3 tests in treemix v1.13 (Pickrell & Pritchard, 2012) and used a block-jackknife approach with sets of 100 random variants to estimate standard errors.

| Inferring demographic history
Using the full data set, we next examined the history of divergence among taxa by inferring a maximum-likelihood phylogenetic tree with a general time reversible (+GAMMA) model in RAxML v8.2.8 (Stamatakis, 2014). We ran 1,000 replicates to estimate bootstrap support values for each node. With the unlinked data set, we then used the program ∂a∂i (Gutenkunst et al., 2009) to estimate the timing of divergence, effective population sizes and levels of gene flow based on patterns of the folded site frequency spectrum (SFS). We used a folded SFS (i.e., major/minor allele calls rather than derived/ ancestral) because we lacked a closely related reference genome.
To determine the evolutionary relationships among taxa, we tested models of three possible bifurcating tree topologies. For each different topology, we tested a model with no migration and a model with one symmetrical migration rate (m, migrants per generation) between all taxa. Although migration rates may differ between taxa, we only tested a single migration rate in order to avoid over-parameterizing our model. Under each model, we estimated divergence times (t, in generations) between taxa, the ancestral population size (N anc ) and the population sizes for each taxon (N) estimated as the relative change in population size (ν) compared to N anc (e.g., if population size increased 10× compared to N anc , then ν = 10). Here, population size changes were modelled as instantaneous changes following divergence.
For each model, we performed 20 independent runs starting with parameter values sampled randomly across a uniform distribution (0.01 < ν < 100, 0 < 2N anc t < 20, 0 < 2N anc m < 5) and determined parameter values with the highest log-likelihood. We then examined which of the maximum-likelihood models produced the best overall fit to the data using a composite-likelihood ratio test with the Godambe Information Matrix (Coffman et al., 2016).
Confidence intervals (CI) for each parameter were estimated with the Godambe Information Matrix with 100 bootstrap data sets comprising 50% of SNPs randomly selected from the full data set. As parameters t and m are reported in coalescent units, we scaled these parameters by θ, assuming the Arabidopsis genomewide SNP mutation rate (μ = 7e−9; Ossowski et al., 2010) and a 132:1 ratio of callable sites:SNPs (based on output from Stacks).
We validated the final model by comparing the predicted SFS to the observed SFS for each population. We also simulated 10,000 SNPs under the final model with msms (Ewing & Hermisson, 2010) and generated expected per-SNP F ST distributions with msstats (Thornton, 2003) to compare against the observed F ST distributions calculated with vcftools.

| Testing geographic and ecological correlates of genetic differentiation
To determine the relative roles of geography and ecology in shaping genome-wide patterns of differentiation within and among Astragalus taxa, we used BEDASSLE v1.5 (Bradburd et al., 2013) implemented in R (R Core Team, 2021). BEDASSLE is a Bayesian statistical approach that models the covariance in allele frequencies between populations as a function of geographic and ecological distance. Parameter values are estimated using a Markov Chain Monte Carlo (MCMC) algorithm. The crucial parameter is αE/αD, which measures the effect of ecological distance (αE) on genetic differentiation relative to geographic distance (αD). We created a genetic distance matrix based on sampling localities using pairwise Weir-Cockerham F ST values estimated with the unlinked SNP data set in vcftools. We generated a geographic distance matrix between sampling localities using a publicly available script written by Peter Rosenmai (https://eurek astat istics.com/calcu latin g-a-dista nce-matri x-for-geogr aphic -point s-using -r/). To generate ecological distance matrices, we first downloaded 11 physical and chemical soil variables at 250 m resolution from soilgrids.org and 19 climatic variables from Worldclim.org (Hijmans et al., 2005; see Table S2 for variable information). We then tested for multicollinearity between all soil and climatic variables using a custom R script and removed highly colinear variables (p < .01). The final set included nine uncorrelated soil and four uncorrelated climate variables ( Figure S2). To further reduce the dimensionality of the data set, we performed a PCA in R (prcomp in "stats" package; R Core Team, 2021) on the ecological variables scaled to have unit variance. We then generated ecological distance matrices in R based on the climate and soil variables.
Prior to running the final analyses with BEDASSLE, we performed short runs of 10,000 generations and tuned MCMC parameter values until acceptance rates were ~20%-50%. For our final analyses, we ran BEDASSLE with a beta-binomial model for two million generations, sampling every 1,000 generations.
We performed three separate runs to ensure convergence on a single posterior distribution. We verified that parameter and log-likelihood values stabilized over the course of each run to evaluate model performance. Based on these evaluations, we discarded the first 1.5 million generations as burn-in and used the remaining samples to generate posterior distributions of parameter values.

| Uncovering putative locally adaptive loci
To identify signatures of selection on SNPs within each population while accounting for the unique history of the clade, we identified  (Yi et al., 2010). Here, we used msms (Ewing & Hermisson, 2010) to simulate 10,000 SNPs under the inferred population history. Output from msms was analysed using the msstats software (Thornton, 2003)

| Genetic structure and diversity
We obtained a total of 10,241 high-quality SNPs distributed across 5,269 loci. Our unlinked data set (one SNP per locus) therefore consisted of 5,269 SNPs. Across our final data set, we obtained a mean coverage of 23.1× ± 7.0× standard deviation (SD) per individual with an average missing data of 24.9% ± 9.0% SD per individual (see Table S3 for additional information). Data generated from this study are available at the USGS ScienceBase-Catalog (Jones et al., 2021).
We uncovered three distinct genetic clusters corresponding to previously described taxonomic groups A. iselyi, A. sabulosus var. sabulosus and A. sabulosus var. vehiculus with a PCA on unlinked SNPs ( Figure 1b). We found strong evidence for the same genetic groups and no evidence of gene flow or continuous structure between them with fastStructure (K = 3 best model; Figure 1c; see Figure S3 (Table 2), indicating shared ancestral polymorphism or gene flow. We also found relatively high genetic diversity (π ≈ 0.002-0.003) and low inbreeding coefficients (mean F IS ≈ 0) within each taxon (Table 1).

| Evolutionary history of divergence, population size and gene flow
Phylogenetic analyses revealed iselyi, sabulosus and vehiculus individuals each formed their own monophyletic clades with 100% bootstrap support values (Figure 2a). However, demographic modelling tests of alternative evolutionary relationships among taxa were ambiguous. Maximum-likelihood values were not significantly different between models with either iselyi (log-likelihood = −2648), sabulosus (log-likelihood = −2651), or vehiculus (log-likelihood = −2647) as the earliest diverging group. Moreover, parameter estimates were essentially identical between alternative models (Table S4).
Therefore, for the remainder of the paper, we present only results from the (iselyi, (sabulosus,vehiculus)) topology model. Although we cannot be certain of the precise tree topology for these taxa, a topology with sabulosus and vehiculus as reciprocally monophyletic is consistent with current taxonomy and the number of private and shared SNPs between taxa (see Table 2).
Under our best split scenario, we found support for a model with

| Landscape genetics of population differentiation
We next explored differences in the environments occupied by Astragalus taxa. A PCA of climatic and soil variation showed that the first two principal components explained 33.2% and 21.7% of TA B L E 1 Population genetic summary statistics for each Astragalus taxon, including the three-population admixture statistic (F3), mean individual inbreeding coefficient (F IS ) and nucleotide diversity (π) Note: Diagonal numbers indicate the number of private SNPs found in the respective taxon. Non-diagonal numbers show the number of SNPs uniquely shared between two taxa. Here, SNPs were called separately for each taxon in Stacks, which is why the total number of SNPs (shown in parentheses) differs from the numbers reported in the full and unlinked SNP data sets. We found a strong correlation between genetic differentiation and geographic distance (adjusted R 2 = .71; p < 2.2e−16; Figure 3b), temperature seasonality distance (adjusted R 2 = .57; p < 2.2e−16; Figure 3c) and silt content distance between sites (adjusted R 2 = .15; p = 1.2e−7; Figure 3d). However, geographic distance was strongly correlated with temperature seasonality (adjusted R 2 = .70; p < 2.2e−16) F I G U R E 3 (a) The first two principal components from a PCA of climatic and soil variables (PC1 explains 33% and PC2 explains 22% of variation) at each of the 13 sampling localities for iselyi (n = 5 sites; purple), sabulosus (n = 6 sites; teal) and vehiculus (n = 2 sites; maroon  Figures S7b and S9), indicating that an approximately 0.3%-0.4% difference in soil silt content between sites has the same effect on genetic differentiation as 1 km. Therefore, the total effect on genetic differentiation is equivalent to an additional ~22-28 km of distance between vehiculus and sabulosus sites (9% difference in silt content) and ~17-21 km of geographic distance between vehiculus and iselyi sites (6.7% difference in silt content).

| Putative adaptive loci
We next used our "linked" data set to search for loci and SNPs with significant levels of differentiation between populations relative to our null demographic model.

| The demographic history of Astragalus diversification
Diversification of angiosperms is hypothesized to be strongly shaped by population fragmentation and limited gene flow between local populations (Niklas et al., 1986), factors that may be strongly influenced by large-scale environmental change. We therefore tested whether the timing of divergence among Astragalus taxa coincides with known periods of environmental change, which would provide evidence for a role of temporal environmental heterogeneity in shaping diversification.
Our study indicates diversification of Astragalus taxa occurred ~130 thousand years ago (kya; assuming one generation every 2 years, as is common in many graminoids and forbs occupying extreme environments that do not reproduce in the year that seeds germinate and establish, e.g., Körner, 2003), roughly corresponding to the beginning of the Last Glacial Period (~110 kya) during the Pleistocene. The Pleistocene was a biogeographically dynamic epoch defined by climate change and cyclical glaciations, which facilitated large-scale shifts of species' ranges, population sizes, and gene flow or hybridization between different populations or species (Galbreath & Cook, 2004;Graham, 1986;Hofreiter & Stewart, 2009;Swenson & Howard, 2005). Changes in both demographic and environmental conditions during this period provided opportunities for adaptive and non-adaptive diversification across diverse clades (Egan & Crandall, 2008;Galbreath et al., 2009;Graham, 1986;Knowles, 2001;Rickart, 2001;Spellman & Klicka, 2007). Thus, in agreement with our prediction, we find evidence that large-scale climate change may have played an important role in generating the demographic or ecological conditions for diversification of Astragalus taxa.
Plant diversification during environmental change can be triggered by numerous biogeographical processes including range fragmentation (i.e., vicariance; Kropf et al., 2006), range shifts (Knowles & Massatti, 2017) and rare long-distance dispersal (Särkinen et al., 2007), leading to the geographic and genetic isolation of populations. To test whether such processes may have played a role in the diversification of Astragalus taxa (Prediction b, in agreement with Sanderson & Wojciechowski, 1996), we used demographic modelling to examine historic population sizes and rates of gene flow between iselyi, sabulosus and vehiculus. We found that divergence between the lineages leading to the present taxa occurred nearly simultaneously from a common ancestor with a small population size (N anc = 1,644; Figure 2b). We suggest these patterns may be best explained by founder events caused by fragmentation of the ancestral range into three parts, rather than other processes like independent long-distance dispersal events (e.g., see Schneeweiss & Schönswetter, 2010). Consistent with this hypothesis, we find evidence of extremely low gene flow between Astragalus taxa (approximately 1 migrant every 167 generations; Figure 2b) and limited dispersal across the landscape (i.e., isolation-by-distance, Figure 3b).
The substantial increases in population size following divergence may be the result of the opening of new habitat or local adaptation following environmental change (see below for more discussion of local adaptation). Although additional genomic research of this genus is needed, our findings are consistent with the hypothesis that vicariance events and genetic isolation of local demes contributes to rapid diversification of Astragalus following climate change (Gavrilets & Vose, 2005;Simões et al., 2016), similar to diversification patterns in other rare and narrowly distributed plant species (e.g., Backes et al., 2019;Godt et al., 1996;Kropf et al., 2006). Under a vicariance model, we would generally expect strong genetic drift within small, isolated populations to decrease intrapopulation genetic diversity. Yet, in contrast to this prediction, we uncovered relatively high nucleotide diversity (Table 1; comparable to wild and cultivated soybeans; Lam et al., 2010), high shared ancestral polymorphism (Table 2), and large long-term effective population sizes (~7-20× larger than N anc ; Figure 2b). In fact, studies have revealed relatively high levels of genetic diversity across many rare Astragalus taxa (Alexander et al., 2004;Baskauf & Burke, 2009;Harrison et al., 2019;Karron et al., 1988;Neel, 2008;Walker & Metcalf, 2008) and some rare plant populations more generally (Barrett & Kohn, 1991;Ellstrand & Elam, 1993). We suggest that aspects of Astragalus seed biology may lend important insights into these somewhat counterintuitive results.
Astragalus seeds are large, nutrient dense, and tend to exhibit high rates of predation and low rates of germination (Combs et al., 2013;Green & Palmbald, 1975;Kaye, 1999). Collectively, these characteristics are thought to lead to inefficient dispersal, low population growth rates and small population sizes (Combs et al., 2013;Kaye, 1999). This is consistent with the low census estimates of these populations (iselyi N c = 3,192, sabulosus N c = 4,649, vehiculus N c = 1,831) and low rates of gene flow across the landscape. Astragalus is also known to have genetically diverse and strongly age-structured seed banks (Molnár et al., 2015;Morris et al., 2002;Van Buren & Harper, 2003), which may maintain genetic variation through time despite small census population sizes and limited gene flow. Seed banks are hypothesized to contribute to the maintenance of variation in rare plant species more generally (Ellstrand & Elam, 1993); however, this hypothesis remains largely untested. Alternatively or in addition, the high N e :N c ratio among iselyi, sabulosus and vehiculus could be explained by recent population declines linked to habitat destruction (e.g., oil and gas development, livestock grazing, camping, and off-highway vehicle use; Ouren et al., 2007). Recent changes in population sizes are incredibly difficult to detect with the SFS (Beichman et al., 2018;Robinson et al., 2014), and our data set likely lacks the power to test this hypothesis. However, we found no signatures of inbreeding (Table 1), which may be expected with recent severe population declines (e.g., Jones et al., 2020). Although somewhat speculative, we suggest that limited seed dispersal and the maintenance of healthy seed banks are parsimonious explanations for widely observed patterns of rarity and high genetic diversity among Astragalus species.

| The ecological drivers of diversification
Given the apparent absence of adaptive innovations defining Astragalus (and the more-inclusive Astragalean clade), diversification of this genus has been largely attributed to demographic factors (Sanderson & Wojciechowski, 1996;Stebbins, 1981). Yet, many Astragalus species are restricted to narrow edaphic or climatic regions, suggesting they may be locally adapted ecological specialists (Baskin & Baskin, 1974Rundel et al., 2015). In fact, rare plant species occupying spatially heterogenous regions are often habitat specialists with narrow ranges (Cowling et al., 1996), suggesting the persistence of rare species may be facilitated by local adaptation.
Patterns of limited gene flow and high genetic diversity that are relatively common in Astragalus and other rare species likely facilitate ecological adaptation (alternatively, ecological adaptation may contribute to limited gene flow); however, empirical tests of local adaptation are challenging and rare. Here, we leveraged landscape and population genomic approaches to test for signatures of local adaptation within Astragalus taxa.
Environmental variables are often the primary ecological factors shaping adaptive variation among a wide variety of plant species (Fischer et al., 2013;Manel et al., 2012), and Pleistocene glaciations are known to cause populations to experience novel environmental conditions (Jackson & Overpeck, 2000;Massatti & Knowles, 2020;Williams & Jackson, 2007). Soil and climate data show that iselyi, sabulosus and vehiculus occupy distinct ecological niches relative to one another (Figure 3a). We found evidence that iselyi (which occurs at higher elevations) occupies colder, wetter and less seasonal environments compared to sabulosus and vehiculus. Moreover, iselyi habitat appears to be differentiated from the other taxa based on certain soil characteristics, including higher soil organic carbon stock (OCSTHA; Figure 3a). Habitat for vehiculus appears to be comprised of soils with relatively high sand content and low silt content.
We explored whether these environmental features are associated with patterns of genome-wide differentiation, which would be indicative of local adaptation. As predicted, we found that genome-wide differentiation across the landscape appears to be strongly shaped by both climatic and soil variation (i.e., isolationby-environment; Wang & Bradburd, 2014;Figures 3c,d and S7).
We should note that, although genome-wide differentiation is associated with contemporary environmental data, it remains possible that differentiation could be more strongly associated with historic environments. Given our ignorance about the historic distributions of these taxa, this is difficult to test and we therefore assume that the associations we uncovered reflect recent or contemporary adaptation. In support of local adaptation, we also identified several highly differentiated, putative adaptive loci in the focal taxa, even after controlling for the demographic and diversification history of the clade (Figure 4). Additional research is needed to validate these selective signatures and uncover the possible functional genetic elements and phenotypes associated with these genetic changes; however, given the nature of the reduced representation data set used herein, there are likely many more putative adaptive loci to be uncovered using targeted approaches (Jones & Good, 2016). The evidence of local adaptation we uncovered does not necessarily support a role of ecological opportunity in initiating divergence between Astragalus, as adaptation may have occurred well after taxa became isolated from one another. Regardless, spatial climatic and soil variation appear to be important factors shaping local adaptation and possibly endemism within Astragalus species.
Our study demonstrates that a combination (and potentially an interaction) of demographic and ecological forces have contributed to diversification of Astragalus taxa. The timing of divergence coincides with a major historic climate change event, which potentially contributed to the geographic or genetic isolation of populations in novel environments. Moreover, the demographic factors shaping low gene flow and high genetic diversity in Astragalus (e.g., seed biology) may have also been conducive to allowing local adaptation to these novel environments, which could further reduce gene flow.
However, our insights into the diversification process of this clade are relatively coarse, and further research is needed to more rigorously test these hypotheses.

| Conservation implications
Our investigation of the iselyi-sabulosus-vehiculus clade in southern Utah revealed evidence for three highly differentiated genetic clusters that represent the currently defined taxonomic units (Figure 1; Barneby, 1989;Welsh, 1974Welsh, , 1998. Demographic modelling suggests that these taxa diverged almost simultaneously from a common ancestor (Figure 2b). More detailed ecological and phylogenetic studies (including closely related congeners) are needed to determine whether these taxa should each be considered their own species or as varieties of a single species; other genetic studies of Astragalus clades suggest that either scenario may be acceptable based on observed values of genetic differentiation (e.g., Alexander et al., 2004;Baskauf & Burke, 2009;Massatti et al., 2018;Travis et al., 1996).
Regardless of taxonomic classification, conservation and management efforts may best be served by maintaining the unique genetic and adaptive histories of each taxon. While they do not appear to currently face high risk of population decline due to genetic factors (e.g., low genetic diversity or inbreeding depression), all taxa face potential risks of decline due to land use (Franklin, 2003). A crucial goal of future research will be to discern whether small census sizes are typical for these taxa or are the result of recent population declines. In addition, modelling habitat availability in anticipation of climate and land use change may help to establish conservation priorities by identifying whether certain Astragalus populations are disproportionately vulnerable.
Although future research is required, Astragalus seed banks may provide a crucial natural resource for conservation and restoration of vulnerable species. Collectively, our study provides critical insights into the concerted influence of demography and adaptation during the diversification of an endemic clade, the details of which likely have important implications for the conservation of many rare species across the region.

ACK N OWLED G EM ENTS
We would like to thank Shannon Lencioni, Jennifer Shostrand, Mindy Wheeler and Robert Fitts for assistance with field collections.
We thank Gery Allan, Lela Andrews and Erica Sukovich at Northern

Arizona University's Environmental Genetics and Genomics
Laboratory for expert support with laboratory procedures. We would also like to thank Jennifer Lewinsohn, Rita Reisor, and the U.S.
Geological Survey Ecosystem Mission Area for providing the opportunity for this research. Finally, this manuscript has been greatly improved by comments from Martin Wojciechowski and three anonymous reviewers. Any use of trade, product or firm names is for descriptive purposes only and does not imply endorsement by the U.S.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13288.

DATA AVA I L A B I L I T Y S TAT E M E N T
Soil and climate variables (found in Table S2) can be downloaded from the Soil Grids database (http://www.soilg rids.org) and Worldclim database (http://www.world clim.org). Files used in analyses are archived in ScienceBase (Jones, Winkler, et al., 2021).