Landscape genomics of the streamside salamander: Implications for species management in the face of environmental change

Abstract Understanding spatial patterns of genetic differentiation and local adaptation is critical in a period of rapid environmental change. Climate change and anthropogenic development have led to population declines and shifting geographic distributions in numerous species. The streamside salamander, Ambystoma barbouri, is an endemic amphibian with a small geographic range that predominantly inhabits small, ephemeral streams. As A. barbouri is listed as near‐threatened by the IUCN, we describe range‐wide patterns of genetic differentiation and adaptation to assess the species’ potential to respond to environmental change. We use outlier scans and genetic‐environment association analyses to identify genomic variation putatively underlying local adaptation across the species’ geographic range. We find evidence for adaptation with a polygenic architecture and a set of candidate SNPs that identify genes putatively contributing to local adaptation. Our results build on earlier work that suggests that some A. barbouri populations are locally adapted despite evidence for asymmetric gene flow between the range core and periphery. Taken together, the body of work describing the evolutionary genetics of range limits in A. barbouri suggests that the species may be unlikely to respond naturally to environmental challenges through a range shift or in situ adaptation. We suggest that management efforts such as assisted migration may be necessary in future.


| INTRODUC TI ON
A central goal of molecular ecology and evolutionary biology is to understand the genetic basis of local adaptation (Hoban et al., 2016;Storz, 2005). Landscape genomics has emerged as an analytical framework to achieve this goal (e.g., Joost et al., 2007;Manel & Holderegger, 2013;Storfer et al., 2018), which is particularly pertinent in light of the potential influence of climate change on species' geographic distributions (Razgour et al., 2019). Indeed, many species are already shifting their geographic distributions toward the poles (Chen et al., 2011;Hickling et al., 2006;Parmesan et al., 1999;Parmesan & Yohe, 2003;Pecl et al., 2017), and approximately a quarter of species assessed by the IUCN face the threat of extinction (IUCN, 2020 impacts on endemic species with small geographic distributions because such species tend to have narrow tolerance ranges to climatic factors such as temperature and precipitation (Sheth & Angert, 2014;Urban, 2015;Yu et al., 2017).
Species range limits often exist, in theory, because populations at the range edge are small and lack genetic diversity to adapt beyond the edge (Eckert et al., 2008;Sexton et al., 2009). Edge populations are thus often the first to disappear when a species undergoes a range-wide decline (Doherty et al., 2003). Landscape genomics studies elucidate patterns of adaptive genomic diversity that management initiatives can use to conserve the capacity of taxa to respond to environmental change (Funk et al., 2019). For example, assisted migration has emerged as one strategy to relocate populations of species threatened by climate change to areas of suitable habitat (Hällfors et al., 2014;Vitt et al., 2010). Another strategy is genetic rescue, or the process of introducing individuals from more genetically variable populations into those that are genetically depauperate to enhance their adaptive potential and fitness (Fitzpatrick et al., 2020;Pimm et al., 2006;Whiteley et al., 2015).
When specific threats are known (e.g., increasing average temperatures), conservation efforts may benefit from genotyping individuals and translocating those identified as harboring putatively beneficial alleles (Hohenlohe et al., 2019).
Amphibians appear particularly sensitive to threats imposed by climate change (Blaustein et al., 2010;Lourenço-de-Moraes et al., 2019), owing to their limited dispersal abilities to track climate shifts and frequent reliance on ephemeral habitats for breeding and larval survival (Blaustein et al., 2010). The streamside salamander, Ambystoma barbouri (Kraus & Petranka, 1989), is an endemic species with a small geographic range primarily restricted to central Kentucky, southeastern Indiana, and southwestern Ohio (Kraus & Petranka, 1989;Petranka, 1998). A. barbouri is listed as nearthreatened by the IUCN, which cites urban and agricultural development, forest harvesting, invasive species, and climate change as drivers (Hammerson, 2004). A. barbouri adults lay eggs with consequent larval development primarily in small, ephemeral streams that are prone to early drying in hot years (Petranka, 1998). Larval activity varies with stream permanence, and experiments suggest a strong genetic component that limits the extent of phenotypic plasticity of developmental timing, which may experience selection (Micheletti & Storfer, 2020;Storfer & Sih, 1998).
On a range-wide scale, A. barbouri conforms to the expectations of the central-marginal hypothesis (Brown, 1984;Eckert et al., 2008). That is, there is a negative correlation of genetic diversity, effective population size, and population connectivity moving from the core of the species' distribution to the edge (Micheletti & Storfer, 2015). Additionally, despite asymmetric gene flow along transects from the range core to range periphery, a reciprocal transplant field experiment suggests that edge populations show evidence of local adaptation (Micheletti & Storfer, 2020). The geographic range of A. barbouri appears to be restricted owing to a combination of relatively cryptic environmental variables elucidated by a landscape genetics study, including limited limestone availability and increases in growing season precipitation that create high resistance to gene flow (Micheletti & Storfer, 2017). To date, however, genetics studies on A. barbouri have been conducted on small numbers of presumably neutral loci, including allozymes (Storfer, 1999) and microsatellites (Micheletti & Storfer, 2015. Collectively, the relatively strong genetic isolation among study sites (Storfer, 1999) and evidence for local adaptation at a range-wide scale warrants further investigation of the genomic basis of adaptation and the capacity to respond to environmental change in this species.
The major goal of this study was to identify genetic markers showing large signatures of divergent selection among study sites spanning the geographic range of A. barbouri, including a site at the core, as well as sites at the northern, southern, and western range edges ( Figure 1). Specifically, we used double-digest restriction siteassociated DNA (ddRAD) sequencing (Peterson et al., 2012) to: (1) describe genetic differentiation among range core and peripheral study sites; (2) evaluate genomic evidence for local adaptation; and (3) gain insight into the potential for adaptive responses to environmental change. The collection localities included fishless, ephemeral streams in the geographic core of the distribution and at the extreme northern, southern, and western edges. We selected these specific edge sites because previous research suggests that they belong to genetically differentiated clusters (Micheletti & Storfer, 2017) and would allow the characterization of range-wide genetic variation. The tissue samples from each site consisted of tail tips from larval salamanders, which were collected among different stream pools within single stream basins to avoid sampling full siblings. Tissues were stored in 95% ethanol.

| Next-generation sequencing
We performed double-digest restriction site-associated DNA sequencing (ddRADseq; Peterson et al., 2012), using the restriction enzymes EcoRI and PstI, both of which use six-base-pair recognition sites. We multiplexed individuals with 30 unique adapter barcodes. We size-selected 500-600 bp fragments using a Pippin

| Variant calling and filtering
We carried out an integrated alignment procedure (Paris et al., 2017) by first assembling RAD loci de novo in Stacks Version 2.52 (Rochette et al., 2019) and then using the mem option in BWA (Li & Durbin, 2009) to align catalogued de novo loci to the June 2018 chromosome-level genome assembly of A. mexicanum (Smith et al., 2019), the most closely related salamander genome available. The de novo RAD locus catalog was created using the Stacks script denovo_ map.pl, with options set to require a minimum of three reads to form a stack (m = 3), allow three mismatches between read stacks within individuals (M = 3), and allow four mismatches between read stacks across individuals (n = 4). These parameter values were optimized by following Paris et al. (2017). For alignment, the clipping penalty for BWA-mem was specified as 10 and all other settings were left as default. Mapping information for RAD loci with MAPQ > 10 was integrated with the Stacks RAD locus catalog using the Stacks script stacks-integrate-alignments.pl.
We filtered the resulting single nucleotide polymorphisms (SNPs) using an iterative filtering procedure modified from O'Leary et al.
(2018). Briefly, iterative filtering involves alternating and progressively increasing the stringency of filters, a process that may generate a final dataset containing a larger number of retained individuals and higher-quality markers than a filtering process that applies each filter at a strict level only once (O'Leary et al., 2018). The iterative filtering stage of our pipeline involved cycling between removing SNPs with high missingness in at least one of the four collection localities and removing individuals with high missingness. We also highlight that we employed a number of filters intended to reduce the abundance of RAD loci representing collapsed paralogs, which may otherwise be prevalent given the large, repetitive genomes of salamanders. The full filtering scheme is detailed in Table S1.
We evaluated the genomic distribution of SNPs in the final dataset by correlating numbers of SNPs that mapped to each reference chromosome with reference chromosome size; a strong correlation would suggest proportional representation of chromosomes in our dataset. We also characterized whether SNPs were within or between gene models of the A. mexicanum genome. We first downloaded BED files containing gene model annotations for each chromosome of the June 2018 A. mexicanum genome assembly from the ucsc.edu) and then used the BEDTools closest command (Quinlan & Hall, 2010) to identify the gene models nearest each SNP.

| Summary statistics and population genetic structure
We used fastStructure (Raj et al., 2014) and the find.clusters function in the R package Adegenet 2.1.3 (Jombart, 2008) to infer the number of genetic clusters represented in the dataset. For fast-Structure, we assessed support for K = 1-10 using 20 replicates for each value of K and subsequently used the python script chooseK included with fastStructure to identify the value of K with the F I G U R E 1 Map indicating the four collection localities. The continuous geographic range of Ambystoma barbouri is indicated as the white outline in the main plot and blue shaded region in the inset. Tissue samples of 28 larval A. barbouri were collected from each of four locations, marked with white points and location labels, corresponding to the range center (C) and the extreme northern, western, and southern range edges (N, W, and S, respectively). Map features are plotted over mean annual temperature greatest statistical support. We visualized admixture proportions of individuals using the software Distruct v1.1 (Rosenberg, 2004). The Adegenet find.clusters function transforms data using principal components analysis (PCA) before using a successive K-means clustering algorithm to evaluate support for different values of K. We assessed results for K = 1-10 and selected the value of K with the lowest value of the Bayesian Information Criterion. We also employed PCA and discriminant analysis of principal components (DAPC) in Adegenet (Jombart, 2008). Prior to performing PCA and DAPC, we used the impute function from the R package LEA  to impute missing data based on ancestry and genotype frequency estimates obtained from the best run of 500 replicates of the snmf function. We used xval optimization to determine the number of principal components to be used in DAPC and retained K discriminant axes. We evaluated both PC and DAPC scores of individuals because PCA characterizes variation occurring both within and between genetic groups, whereas DAPC focuses on variation between groups (Jombart et al., 2010). We additionally evaluated the correlation between PCs, latitude, and longitude using Pearson's correlation test (e.g., Novembre et al., 2008).
We used the R package diveRsity (Keenan et al., 2013) to calculate Weir and Cockerham's estimators of the inbreeding coefficient F IS and F ST (Weir & Cockerham, 1984). We estimated F IS for each population using the basicStats function and estimated 95% confidence intervals using 1000 bootstraps. We estimated 95% confidence intervals (CI) for pairwise F ST based on 1000 bootstraps using the bias-corrected bootstrapping method implemented in the diffCalc function; we carried out bootstrapping across SNPs, as recommended by Weir and Cockerham (1984). We also recalculated pairwise F ST after removing SNPs determined to have significant signatures of selection in univariate tests (see below). We further characterized genetic diversity for each locality by calculating individual observed heterozygosity and per-SNP nucleotide diversity (π) in VCFtools (Danecek et al., 2011).
We estimated effective population size (N e ) for each locality using the linkage disequilibrium-based method in NeEstimator v2.1 (Do et al., 2014); N e estimates were obtained after removing SNPs with significant signatures of selection. We repeated effective population size estimates using allele frequency thresholds of 0.05, 0.02, and 0.01.

| Environmental variables
We initially downloaded GIS layers for the 19 bioclimatic variables available from the WorldClim database (Fick & Hijmans, 2017) and six additional variables, briefly outlined here (Table S2). The bioclimatic variables are all related to temperature and precipitation, which are related to known stressors in amphibians (e.g., desiccation; Daszak et al., 2005). Temperature seasonality may select for the width of thermal tolerance (i.e., the difference between upper and lower thermal limits) in amphibians (Snyder & Weathers, 1975).
The bioclimatic variables also include quarter-specific environmental conditions, which appear to impact habitat suitability for A. barbouri (Micheletti & Storfer, 2015).
We also included solar radiation averaged across months and elevation in our analyses (Fick & Hijmans, 2017) because ultraviolet radiation exposure negatively correlates with aspects of fitness in various amphibian species (Bancroft et al., 2008). A. barbouri also favors forested habitat (Kraus & Petranka, 1989), so we included percent forest canopy cover (Coulston et al., 2012) and land cover (Homer et al., 2020)  Soil pH and percent soil organic carbon both correlate with many biophysicochemical properties of soil (Brady & Weil, 1999;Smith & Doran, 1997), so we included both factors, predicted for a depth of 5 cm, as additional variables (Ramcharan et al., 2018). We converted pH to the corresponding concentration of hydronium ions prior to analyses.
Using QGIS Desktop 3.12.1 (QGIS development team, 2009), we defined a 2.5 km radius buffer zone around each collection locality; the larvae sampled are unlikely to be offspring of individuals that dispersed from locations outside of this buffer (e.g., Orloff, 2011;Rittenhouse & Semlitsch, 2007). We then extracted the mean value of each environmental variable within the buffer to determine broad differences in environmental conditions among sites. For the binary variable of forest habitat vs. nonforest habitat, we calculated the proportion of the buffer zone corresponding to forest habitat. The full set of variables and values for each collection locality can be found in Table S2. Note that the geographic range of A. barbouri is small, so the species' environmental tolerance ranges are likely narrow (Micheletti & Storfer, 2015;Sheth & Angert, 2014;Yu et al., 2017). Although environmental factors may appear to have little variation among collection localities, this variation may be important from a biological perspective. For example, the maximum and minimum elevations (differing by 110 m) observed among our sites are similar to those identified from species presence data used in Micheletti and Storfer (2015).
For genetic-environment association (GEA) tests, we reduced the number of highly correlated environmental variables using a twostage approach. In the first stage, we grouped logically related variables (e.g., temperature-related variables) and retained only those variables correlated at r ≤ 0.90 within each group. We then pooled the remaining variables and again retained only those correlated at r ≤ 0.90. A matrix of Pearson's correlation coefficients for every pair of variables can be found in Table S3. Although this process filtered out annual precipitation, we retained this variable because it is relevant for explaining the distribution of A. barbouri (Micheletti & Storfer, 2017) and generally relevant to various aspects of amphibian ecology (e.g., aquatic breeding, aquatic larval stage, and vulnerability to desiccation); imperfectly correlated variables also may reveal nonredundant signatures of selection. This yielded a final set of six variables: mean annual temperature, temperature seasonality, mean temperature of the wettest quarter, annual precipitation, elevation, and percent soil organic carbon (Table 1). Using these six variables, we additionally quantified environmental distance among the sampling locations. Environmental variables were demeaned and divided by their standard deviations; we conducted PCA on the resulting centered and scaled environmental data and then calculated pairwise Euclidean distances among the four sites using all principal components (e.g., Wang, 2013). We also assessed correlations among pairwise environmental distance, geographic distance, and Rousset, 1997). We additionally assessed correlations between PC1 and PC2 from the genetic PCA (see above) and each of the six environmental factors independently.

| Tests for selection
We carried out several tests for selection that broadly belong to two categories: differentiation outlier scans and GEA analyses. Outlier scans identify loci showing unusual differentiation among populations relative to background genetic differentiation, while GEA analyses identify loci with allele frequencies that correlate significantly with environmental variables while accounting for neutral population structure (Hoban et al., 2016;Rellstab et al., 2015). We used stringent significance thresholds recommended by the authors of each test. When authors of tests did not clearly designate significance thresholds, we consulted the literature for common cutoffs.
We performed outlier detection analyses using the software hapflk v1.4 (Bonhomme et al., 2010), pcadapt (Luu et al., 2017), and the Bayenv2 X T X statistic (Günther & Coop, 2013). We used hapflk v1.4 to compute the FLK statistic for each SNP and assessed outlier significance at a false discovery rate threshold of q<0.05 (Benjamini & Hochberg, 1995). Pcadapt calculates z-scores for each SNP by regressing each SNP against principal components; SNPs that emerge as outliers on the basis of Mahalanobis distances calculated from z-scores are interpreted as putatively experiencing divergent selection among populations. We used the conservative Bonferroni correction of 0.01/10,527 (the number of SNPs) to assess the significance of each SNP (Luu et al., 2017). We additionally used the X T X statistic implemented in Bayenv2, which is a measure of subpopulation differentiation analogous to F ST . Bayenv2 accounts for neutral genetic structure using a population covariance matrix, which we obtained by calculating the median matrix across five replicate runs of Bayenv2 for 500,000 iterations each. Using the empirical covariance matrix, we then ran five independent replicates of Bayenv2 for 750,000 iterations and calculated the median value of X T X across replicates for each SNP. To calibrate the significance threshold of the X T X statistic, we used the BayPass accessory R function simulate.baypass to generate a pseudo-observed dataset (POD) of neutral SNPs based on the empirical covariance matrix estimated in Bayenv2 (Gautier, 2015). The POD contained the same number of SNPs as the real dataset, and sample sizes for each SNP were drawn with replacement from the sample sizes in the real dataset. We calculated the X T X statistic for each SNP in the POD using Bayenv2, TA B L E 1 Final set of environmental factors used in genetic-environment association (GEA) analyses with options identical to the analysis carried out on the real dataset.
We used the X T X value corresponding to the top 0.01% quantile of the POD as the threshold for significance in the real dataset (Gautier, 2015).
LFMM uses a number of latent factors (equal to the K inferred genetic clusters in the data) to account for neutral genetic structure when testing for associations between alleles and environmental factors. We accounted for multiple comparisons using a Bonferronicorrected significance threshold of 0.01/10,527 (Frichot et al., 2013;Yang et al., 2020). Bayenv2 calculates a Bayes factor (BF) for each locus-environmental factor combination, which describes the support for a model which includes the environmental factor relative to a null model which excludes it. Bayenv2 also estimates Spearman's Rho, a rank-based correlation coefficient describing the association between each SNP and environmental factor. Spearman's Rho is robust to outlier populations (Günther & Coop, 2013); Günther and Coop (2013) recommend the simultaneous use of a relatively permissive BF threshold and relatively stringent Spearman's Rho threshold when outlier populations are a concern. We ran five independent replicates of Bayenv2 for 750,000 iterations using the same covariance matrix used in the estimation of X T X and calculated the median

| Additive polygenic scores
Polygenic scores have emerged in the landscape genomics literature as a method to summarize genotype-environment associations across SNPs at the individual level (e.g., Babin et al., 2017;Xuereb et al., 2018), and this implementation is distinct from polygenic scores used in relation to genome-wide association and quantitative genetics studies that establish links between genetic data and phenotypes and/or fitness (e.g., Wray et al., 2014); accordingly, the use of polygenic scores in the present study concerns their implementation in landscape genomics research and does not involve direct measures of phenotypes or fitness. Additive polygenic scores provide evidence for multilocus or polygenic local adaptation through the co-occurrence of alleles across multiple significant SNPs within individuals (Gagnaire & Gaggioti, 2016). We created six sets of significant SNPs, one set per environmental variable. The sets contained significant SNPs identified using a GEA test and at least one other GEA or outlier test. Restricting the polygenic scores analysis to SNPs identified by multiple tests for selection may reduce the influence of false positives.
Following Babin et al. (2017) and Xuereb et al. (2018), for each significant SNP, we identified the allele exhibiting a positive association between its frequency and an environmental factor. We then summed, within a given individual, the dosage (0, 1, or 2) of the positively associated allele across all significant SNPs, taking this value as a raw additive polygenic score. We calculated raw additive polygenic scores in R using a modified version of the script provided by We also calculated polygenic scores in relation to each environmental factor using random subsets of the entire SNP dataset. Each random subset contained an identical number of SNPs as the set of SNPs significantly associated with the environmental factor of interest. For each environmental factor, we repeated this process 1000 times and calculated individual percent polygenic scores. Using each replicate, we created a GLMM as described previously. We then compared the beta coefficient for the environmental factor using the observed polygenic scores with the distribution of beta coefficients obtained from the random SNP subsets to evaluate whether the observed coefficient is greater than expected by chance. This provides a basis for comparing the pattern of polygenic scores observed from the significant SNPs to the background genome.

| Candidate gene identification
We used the BEDTools closest command (Quinlan & Hall, 2010)

| Summary statistics and population genetic structure
FastStructure and the find.clusters algorithm of Adegenet found the greatest support for four genetic clusters corresponding to the four sampling localities ( Figure S2); we therefore grouped individual samples based on collection locality for all subsequent analyses.
FastStructure found little evidence for admixture among the four genetic clusters (Figure 2A)  (Table S4). DAPC highlighted strong isolation between the northern locality and remaining three localities ( Figure 2C).  (Table S5).

| Tests for selection
LFMM 2 identified 512 SNPs significantly associated with at least one of the six environmental factors, and at least 77 SNPs were identified as having a significant association with each environmental factor. Bayenv2 identified 145 SNPs associated with at least one environmental factor. There was some overlap in the SNPs identified by the two GEA analysis programs for a given environmental factor, as well as overlap in the SNPs identified by each program individually across the six environmental factors ( identified by at least two tests for selection ( Figure S4). As evaluated using a permutation test, this represented significantly greater overlap among tests than expected by random chance (p < 0.0001).
Overlap among pairs of tests was also generally significantly greater than expected by random chance (p < 0.007), bar the overlap for the Bayenv2-LFMM and Bayenv2-PCAdapt comparisons (Table S8).  Figure S5. Significant SNPs were found across all reference chromosomes, even when stratified by the environmental factor they were associated with. The number of significant SNPs found on a given reference chromosome was correlated with chromosome size (Table S9). Numbers of significant SNPs used subsequently in additive polygenic scores analyses and the identification of candidate genes are provided in Table S10.

| Additive polygenic scores
There was generally limited variation in the additive polygenic scores of individuals within populations ( Figure 3A-C Figure 3A-C), and permutation tests showed that beta coefficients corresponding to polygenic scores derived from the significant SNP sets were significantly greater than expected by random chance identification of SNPs (p < 0.001 for annual precipitation, mean annual temperature, and temperature seasonality). Models containing either latitude or longitude were also fit for each set of polygenic scores. R 2 for all models (those including the environmental factor and those including latitude or longitude, along with the random effect of site) were >0.99. For each SNP set, polygenic scores appeared well explained by either the environmental factor of interest or latitude ( Figure S6). The relationship between population means of additive polygenic scores and the environment for percent soil organic carbon, elevation, and mean temperature of the wettest quarter were nonmonotonic ( Figure S7D-F), and we did not fit models for these variables.

| Candidate gene identification
We used BEDtools to identify the A. mexicanum gene models nearest each of the 732 significant SNPs. 49 of the 732 significant SNPs occurred in RAD loci that aligned to contigs of the A. mexicanum genome assembly that have not been assigned chromosomal positions; we therefore could not identify gene models near these SNPs.  (Table S11).

F I G U R E 3 Select polygenic scores.
Additive polygenic scores are shown in relation to (a) annual precipitation, (b) mean annual temperature, and (c) temperature seasonality. Polygenic scores were calculated based on SNPs identified as significant in a GEA analysis and at least one other test for selection. Polygenic scores in different panels are based on different but sometimes overlapping sets of SNPs. Scores are reported as a fraction of the maximum possible summed dosage of positively associated alleles an individual could have, given its missing rate. Light gray crosses and curves represent a single representative replicate of polygenic scores calculated from a random subset of the 10,527 SNPs in the full dataset Many of the candidate genes identified near the significant SNPs are associated with general biological processes (e.g., transcription).
We therefore restricted our discussion to several genes that have associations with more specific ecological drivers of local adaptation; a full table of candidate genes, associated GO terms, and additional data can be found in Table S12.
We note WNT7B and AXIN1, which are both related to multicellular organism development (e.g., Zeng et al., 1997), as well as U2AF1L4, which may be involved in regulating the circadian rhythm (Preußner et al., 2014). Among the genes that mapped near SNPs uniquely associated with mean annual temperature, TEX2 has a signature of positive selection in high altitude sheep (Shi et al., 2019) and both ACER1 and PSMG3 are implicated in hypoxia response (Kharrati-Koopaee et al., 2019;Wang et al., 2013). The genes OTX1 and LVRN were found in relation to both mean annual temperature and temperature seasonality, and they have been implicated in hypoxia response and human altitudinal adaptation, respectively (Yi et al., 2013;Zhang et al., 2020). The MGAT5B gene was near a SNP associated with mean temperature of the wettest quarter, and this gene has previously shown differential expression in association with variation in aquatic habitat permanence (i.e., temporary ponds vs. permanent streams) in larvae of the salamander species Salamandra infraimmaculata (Goedbloed et al., 2017). RIN2, found here in association with mean temperature of the wettest quarter, has been implicated in dehydration response (Kordonowy & MacManes, 2017).
ZCCHC6 is a hypoxia biomarker (Mosqueira et al., 2012) and was identified here in association with both mean annual temperature and elevation. PIK3R5 has previously been implicated in thermal stress response (Kim et al., 2017) and was found in relation to mean annual temperature and temperature seasonality. Allele frequency shifts for SNPs near these candidate genes are provided in Figure   S8. No GO biological process, molecular function, or cellular component terms were significantly over-or under-represented among our list of candidate genes at an FDR of 0.05.

| DISCUSS ION
Landscape genomics has emerged as a widely used framework for understanding the spatial distribution of adaptive variation (Joost et al., 2007;Storfer et al., 2018). In light of accelerating climate change, the identification of putatively adaptive variation provides a starting point for assessing whether species may have the capacity to respond to shifting environmental conditions. Using a suite of suggest that a range contraction may be more likely. We arrive at this conclusion through a combination of: (1) evidence for local adaptation; (2) evidence for strong range-wide genetic differentiation, and (3) a synthesis of prior results suggesting that A. barbouri conforms to the central-marginal hypothesis, which we discuss in tandem.

| Tests for selection and additive polygenic scores
Approximately 7.0% of SNPs were found to have either significant associations with at least one environmental factor or significant genetic differentiation among populations, indicating that environmental heterogeneity likely plays an important role in shaping spatial patterns of genomic variation across the range of A. barbouri.
Additive polygenic scores identify signatures of polygenic adaptation (Gagnaire & Gaggiotti, 2016), and we find a strong association of polygenic scores with annual precipitation ( Figure 3A). This suggests that A. barbouri populations may be locally adapted to precipitation through a polygenic architecture, although this pattern may also be explained by geography. Indeed, precipitation is extremely important for amphibian larval survival in general, particularly for A. barbouri that develop in small ephemeral streams (Petranka, 1998;Sih et al., 2000;Storfer & Sih, 1998).
The northern population has extremely low polygenic scores for mean annual temperature but extremely high polygenic scores for temperature seasonality, relative to the other three populations.
The southern, eastern and western populations, taken together, still show a positive relationship with these two environmental variables. The strong differentiation between the north population and the other three populations may facilitate local adaptation (Garcia-Ramos & Kirkpatrick, 1997;Lenormand, 2002), leading to adaptive fixed differences between the northern population and the other three populations, which are considerably less genetically differentiated from one another. Indeed, the northern population showed evidence of local adaptation, with higher larval survival when reared in northern conditions than core conditions in a reciprocal transplant experiment (Micheletti & Storfer, 2020). Alternatively, it is possible that our GEA analyses may have by chance identified neutral loci exhibiting correlations with environmental factors driven largely by the differentiation of the northern population. Strong population structure and confounding between population structure and the environment likely also contribute to observed patterns of polygenic scores, as polygenic scores can also be explained by geography (i.e., latitude and longitude). Although GEAs account for population structure and our approach to polygenic scores involved the removal of SNPs found by only one test for selection, the influence of false positives may still be substantial (Selmoni et al., 2020), suggesting that future work should focus on validating our results (see Limitations and Future Directions).
In general, the identification of SNPs significantly associated with precipitation and temperature-related variables, as well as evidence suggestive of polygenic adaptation, supports the importance of these factors as selective pressures on amphibian species. Indeed, amphibians are generally vulnerable to desiccation (e.g., Daszak et al., 2005;McMenamin et al., 2008), and drying of breeding habitat negatively impacts recruitment (McMenamin et al., 2008;Semlitsch & Wilbur, 1988). Other ambystomatids become more active during periods of elevated precipitation (e.g., Trenham, 2001), further suggesting that traits related to desiccation resistance or avoidance could improve fitness in areas that are more water-limited.

| Potential to respond to environmental change
Numerous studies point to species' geographic distributions changing in response to climate change. For example, many species have begun to track environmental conditions toward the poles (Chen et al., 2011;Hickling et al., 2006;Parmesan et al., 1999;Parmesan & Yohe, 2003). However, this response is contingent upon species having dispersal capabilities sufficient to match the pace of environmental change (Schloss et al., 2012). Taxa with poor mobility, such as salamanders (Smith & Green, 2005), may have poor capacity to respond to environmental change through range shifts (Angert et al., 2011).
Indeed, we show very strong genetic differentiation across the geographic range of A. barbouri, consistent with prior work (e.g., Storfer, 1999), suggesting that migration rates are low. However, the single individual collected from the western locality that showed evidence for admixture with the central and southern localities suggests that relatively rare, long-distance dispersal may be possible. Additionally, the northern portion of the geographic range is separated from the majority of the range area by the Ohio River, a large river that putatively acts as a barrier to migration (Micheletti & Storfer, 2017).
Resistance to gene flow also becomes more intense toward the range edges (Micheletti & Storfer, 2017), suggesting poor potential for dispersal beyond the current range. As such A. barbouri is not expected to track a climate-induced shift of suitable habitat northward.
Range-wide, A. barbouri largely conforms to the central-marginal hypothesis, with habitat suitability, genetic diversity, and effective population sizes decreasing along transects from the range core to the range periphery (Micheletti & Storfer, 2015). The present results support previous findings of lower genetic diversity at the extreme southern and northern range edges than at the range core, but the western collection locality appears to have highest genetic diversity of the four localities sampled in our study. However, edge populations studied previously do not universally have lower genetic diversity than the core (Micheletti & Storfer, 2015), so the present results should not be construed as invalidating the general trends identified using substantially more sampling localities in previous research.
Reduced genetic diversity may constrain adaptive potential (Gilpin & Soulé, 1986), preventing a species from colonizing new regions with environmental conditions differing from its current distribution (Sexton et al., 2009). Apparent poor dispersal and low genetic diversity at the range periphery together suggest that the range of A. barbouri may be unlikely to shift appreciably from its current expanse.
Although some populations of A. barbouri may be locally adapted to current environmental conditions (Micheletti & Storfer, 2020), the present evidence for strong genetic differentiation and prior evidence in support of the central-marginal hypothesis suggest that edge populations may be unable to adapt in situ as the environment changes. Low genetic diversity at the range edges (Micheletti & Storfer, 2015) may not only constrain a range shift but also hinder in situ adaptation if variation at adaptive loci is limiting. Polygenic scores indicate that the northern population may be strongly locally adapted to relatively low annual precipitation, low mean annual temperature, and high temperature seasonality. Indeed, all individuals sampled from this population have similar polygenic scores (either extremely high or low), indicating the near-absence of alleles putatively adaptive for alternative environmental conditions. To the extent that the polygenic scores may reflect the effects of selection, the northern population may therefore have a poor capacity to respond to changing environmental conditions through in situ allele frequency changes. Migration across the range may also be low enough that the spread of standing adaptive variants may be outpaced by environmental change. Together, potential strong local adaptation, low genetic diversity, and poor migration may contribute to evolutionary mismatch with the environment in future (e.g., Zimova et al., 2016), leading to low population fitness and declining population size in the northern range edge.
Previous studies showed that the southern range edge has the lowest habitat suitability (approximately 15-24%; Micheletti & Storfer, 2015). Larvae from this population also had poor overall survival regardless of rearing conditions in a reciprocal transplant experiment (Micheletti & Storfer, 2020). The elevated F IS estimate for the southern population suggests inbreeding as a possible cause. Mean annual temperature is notably higher in this part of the range, and this variable is negatively associated with gene flow (Micheletti & Storfer, 2017), suggesting that warming associated with climate change may result in increasing isolation among populations at the southern range edge and enhanced genetic drift.
Taken together, these data suggest maladaptation in this apparent sink population, which may presage loss of the southern range edge. Potential strong local adaptation and genetic isolation in the north, as well as maladaptation in the south, suggest that a range contraction may be likely.

| Candidate genes
The recovery of genes involved in development (e.g., AXIN1 and WNT7B) and circadian rhythm (U2AF1L4) may reflect variation in the timing of developmental events, such as metamorphosis. Tests for selection in other species have also recovered genes involved in the circadian rhythm (e.g., Dall'Ara et al., 2016;Geraldes et al., 2014).
Indeed, regulation of developmental timing is essential for amphibian species that undergo metamorphosis Wilbur & Collins, 1973), and for A. barbouri individuals that develop in ephemeral streams.
The repeated identification of hypoxia-related genes associated with mean annual temperature but not elevation suggests that high temperatures toward the southern area of the range may impose selection on A. barbouri through the reduced availability of dissolved oxygen in bodies of water. Numerous genes that mapped near SNPs associated with mean annual temperature have been implicated in adaptation to hypoxia in other species. Research in other ambystomatids has revealed a positive correlation between dissolved oxygen and hatching rate (Sacerdote & King, 2009). Further, timeto-hatching and larval development of two ambystomatids were reduced in hypoxic conditions, with survival following hatching also decreasing for one species (Mills & Barnhart, 1999). While it is possible that another variable strongly correlated with mean annual temperature may be the underlying environmental driver for these results, it appears unlikely to be elevation, which has a correlation with mean annual temperature of only r = −0.52 across our four study sites (Table S3). That multiple genes identified here have also been identified in other studies of adaptation increases our confidence that they are true positives while also providing evidence for common genetic mechanisms underlying evolutionary responses to environmental variation.

| Limitations and future directions
Strong genetic differentiation among subpopulations and confounding between population structure and environmental conditions may increase the false positive rates of various tests for selection (e.g., Hoban et al., 2016;Lotterhos & Whitlock, 2015;Novembre & Di Rienzo, 2009;Rellstab et al., 2015). Small numbers of collection localities also limit power to detect genomic signatures of local adaptation and increase false discovery rate (Selmoni et al., 2020). These circumstances are true for our study, and while our decision to calculate polygenic scores using only SNPs found to have significant signatures of selection in multiple tests may serve to limit the influence of false positives, more data are required to better disentangle selection and genetic drift. Given the limitations of our sampling scheme, future research should extend sampling to a larger number of locations across the distribution of A. barbouri.
It is also important to note that because our procedure for identifying candidate genes near SNPs used a reference genome assembly of A. mexicanum, our set of candidate genes is sensitive to the level of structural conservation between the genomes of A. barbouri and A. mexicanum. Knowledge of the extent of linkage disequilibrium along the A. mexicanum genome is also lacking, although the recombination distance per cM is approximately 5 Mb (Smith et al., 2019).
We suspect that recombination between each candidate SNP and variation in or near its corresponding candidate gene may be low, given that distances between them are all ≤5.82 Mb. Given these limitations, the candidate genes identified here should be acknowledged with caution. The generation of additional salamander reference genomes will enable more confident identification of genes involved in adaptation to environmental heterogeneity. This will in turn facilitate management efforts aimed at conserving adaptive potential in the numerous salamander species that are threatened by climate change, land development, and emerging infectious diseases (Blaustein et al., 2010;Collins & Storfer, 2003;Martel et al., 2013;Stuart et al., 2004).
Given that our results are derived from reduced-representation sequencing, we interrogated only a small fraction of the large genome of A. barbouri. Thus, the putatively adaptive genomic variation identified here provides incomplete insight into the genomic architecture underlying local adaptation in this species (Lowry et al., 2017). GEA and outlier tests are also biased toward detecting largeeffect loci, further compounding this limitation (Hoban et al., 2016).
Future transcriptomic sequencing or targeted DNA sequencing (e.g., whole-exome sequencing) may help overcome the difficulties related to the large, repetitive genomes of salamanders and facilitate more complete characterization of adaptive genomic variation (Weisrock et al., 2018).
Additional research should pursue functional validation of candidate loci, which is generally lacking in population and landscape genomics studies (Li et al., 2017;Storfer et al., 2018). Future work could include reciprocal transplant experiments (sensu Micheletti & Storfer, 2020) using genotyped individuals to link genomic variation with phenotypic variation and fitness in different environments.
Nevertheless, the present study provides an initial characterization of the genomic basis of adaptation across the geographic range of A. barbouri, as well as a number of hypotheses concerning the ecological pressures driving adaptation.

| Conservation management implications
A. barbouri is listed as near-threatened by the IUCN, owing to its limited geographic range, habitat alteration resulting from human land use, and climate change (Hammerson, 2004). Prior genetic work has provided a baseline understanding of neutral evolutionary processes occurring across the range, including evidence for low genetic connectivity among subpopulations (Micheletti & Storfer, 2015, 2020Storfer, 1999). Here, evidence for adaptive differences and strong genetic differentiation among subpopulations, as well as limited genetic diversity within some subpopulations of A. barbouri suggests range shifts in response to climate change are unlikely. Consequently, our results could be used as preliminary data for the development of genotyping panels aimed at monitoring neutral genetic variation or a subset of adaptive variation, a method that has been suggested in relation to conservation of other species (e.g., Hohenlohe et al., 2019). Indeed, future management efforts could genotype individuals at neutral SNPs to identify those best suited for genetic rescue of genetically depauperate subpopulations, as well as candidate SNPs to identify those best suited for assisted migration to areas of varying environmental conditions, especially as anthropogenic land use and climate change intensify (Hällfors et al., 2014;Vitt et al., 2010).
Considering these management actions is especially pertinent given that A. barbouri is an endemic species with a small geographic range; such species are expected to be disproportionately threatened with extinction (Manne & Pimm, 2001;Urban, 2015). Protocol ASAF #04526-001.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.5dv41 ns6w.