Association of genetic and climatic variability in giant sequoia, Sequoiadendron giganteum, reveals signatures of local adaptation along moisture‐related gradients

Abstract Uncovering the genetic basis of local adaptation is a major goal of evolutionary biology and conservation science alike. In an era of climate change, an understanding of how environmental factors shape adaptive diversity is crucial to predicting species response and directing management. Here, we investigate patterns of genomic variation in giant sequoia, an iconic and ecologically important tree species, using 1,364 bi‐allelic single nucleotide polymorphisms (SNPs). We use an F ST outlier test and two genotype–environment association methods, latent factor mixed models (LFMMs) and redundancy analysis (RDA), to detect complex signatures of local adaptation. Results indicate 79 genomic regions of potential adaptive importance, with limited overlap between the detection methods. Of the 58 loci detected by LFMM, 51 showed strong correlations to a precipitation‐driven composite variable and seven to a temperature‐related variable. RDA revealed 24 outlier loci with association to climate variables, all of which showed strongest relationship to summer precipitation. Nine candidate loci were indicated by two methods. After correcting for geographic distance, RDA models using climate predictors accounted for 49% of the explained variance and showed significant correlations between SNPs and climatic factors. Here, we present evidence of local adaptation in giant sequoia along gradients of precipitation and provide a first step toward identifying genomic regions of adaptive significance. The results of this study will provide information to guide management strategies that seek to maximize adaptive potential in the face of climate change.


| INTRODUC TI ON
In an era of unprecedented climate change, the adaptive potential of populations has become an increasingly important topic to conservation biologists, raising questions of landscape partitioning of adaptive variation and management strategies to maintain population viability. Given the rapid rate of climate change, new beneficial mutations are expected to play a limited role for species with low mutation rates and long generation times. Therefore, adaptive evolution under climate change for many species will depend on standing genetic variation (Aitken, Yeaman, Holliday, Wang, & Curtis-McLane, 2008; Barrett & Schluter, 2008) that may vary across the landscape and include alleles that gain adaptive value as selection pressures change (Olson-Manning, Wagner, & Mitchell-Olds, 2012). Reliance on standing genetic variation is likely to be particularly true for long-lived sedentary species, such as forest trees that are characterized by adaptive constraints that can limit their evolutionary response to rapid environmental change: extended generation times that result in local persistence, increased rates of genetic drift associated with overlapping generations (Rogers & Prügel-Bennett, 2000), and limited rates of migration due to long generation times. For these species, understanding the distribution of adaptive diversity in relation to recent or past climatic gradients is a critical first step in promoting the adaptive potential of populations in hopes of maintaining future viability (Aitken & Whitlock, 2013;Holderegger, Kamm, & Gugerli, 2006).
The rich history of field research on phenotypic traits in plants (common gardens and reciprocal transplant studies) provides evidence for abundant heritable variation for quantitative traits that are organized along environmental clines (Morgenstern, 1996;Savolainen, Pyhäjärvi, & Knurr, 2007). Until recently, determining the molecular basis of this variation has been less tractable. However, the rapid advancement of genome sequencing, including methods that use reduced genomic complexity (e.g., genotyping by sequencing (GBS), restriction-site associated DNA sequencing (RADseq)), has opened the door to more comprehensive assessments of population-level diversity and allowed for the detection of regions under selection. Although some instances of strong selection on single or few gene loci have been noted (Akey, 2009;Linnen, Kingsley, Jensen, & Hoekstra, 2009;Sella, Petrov, Przeworski, & Andolfatto, 2009), many traits of adaptive importance in plants are believed to be polygenic in nature (Holland, 2007;Le Corre & Kremer, 2012;Pritchard & Di Rienzo, 2010;Yeaman et al., 2016). Under selection, these traits can exhibit subtle changes in frequency across many loci of small effect. Further, demographic processes can shape genetic diversity in ways that mimic selective gradients, as geographic distance and climatic gradients are often autocorrelated. As a result, imprints of selection within the genome can be difficult to detect (Yeaman, 2015), and it is necessary to parcel out the contribution of geographic space in order to successfully identify regions of functional importance (Excoffier, Hofer, & Foll, 2009;Rellstab, Gugerli, Eckert, Hancock, & Holderegger, 2015).
Determining the presence of adaptively important genetic variation and its distribution across a species range is crucial to predicting species' responses to global climate change and directing biodiversity conservation and management efforts (Aitken & Whitlock, 2013;Alberto et al., 2013;Funk, McKay, Hohenlohe, & Allendorf, 2012;Sgrò, Lowe, & Hoffmann, 2011;Sork et al., 2013). This has become an urgent challenge in California, where a protracted drought has resulted in massive tree mortality (USDA, 2016). The Sierra Nevada of California is a high mountain range that collects precipitation from the Pacific Ocean mostly in the form of winter rain and snowfall.
The slow release of water from snowmelt in the spring is an important source of moisture for seedling growth and establishment.
Sierra snowpack has declined in recent years (Fyfe et al., 2017) and high-resolution regional climate models suggest that spring snow water equivalent will decline by 73% by the end of the century, with midelevations (1,500-2,500 m) experiencing the greatest declines (Sun, Berg, Hall, Schwartz, & Walton, 2018).
This elevational range includes the extant groves of the iconic, long-lived conifer, giant sequoia (Sequoiadendon giganteum [Lindl.] Buchholz) that occur in a highly disjunct range consisting of ~70 groves spanning approximately 400 km north to south (Figures 1 and   2). Currently, most giant sequoia populations are in protected areas as this species is valued both culturally and for ecotourism. However, despite this protected status, a key question is whether populations of giant sequoia will remain viable under changing climate.
Our previous work has shown very restricted gene flow (DeSilva & Dodd, 2020), suggesting that natural dispersal outside of existing groves will be unlikely. Long generation times (~305 years; Dodd & DeSilva, 2016) will slow the expansion of new variants that may arise through mutation, which underscores the role of standing genetic variation in determining the future viability of giant sequoia populations. In our landscape genetics study of microsatellite variation, we found some evidence for isolation by environment (IBE), linking genetic divergence at putatively neutral loci to dissimilarity in precipitation, and temperature-related variables (DeSilva & Dodd, 2020).
Although IBE is consistent with local adaptation, it is dependent on a reduction of gene flow from divergent habitats due to selection against nonadapted immigrants, and therefore, patterns of IBE may not be reflective of local adaptation when gene flow is low or absent (Nosil, Vines, & Funk, 2005;Wang & Bradburd, 2014), as is likely the case in sections of giant sequoia range (DeSilva & Dodd, 2020). A recent common garden study reported provenance variation in growth performance, providing support for the existence of adaptive genetic variation across the species' range (Valness, 2016). Yet, to date, no studies have investigated local adaptation in giant sequoia using genomic data. Our ultimate goal was to detect populations that may be genetically responsive to anticipated climate change, so here, we build upon this earlier work by reporting on genomic signatures of selection using a range-wide genotyping-by-sequencing dataset.
Specifically, we utilize an F ST outlier test and gene-environment association methods (LFMM and RDA), to find signatures of local adaptation among giant sequoia populations and locate potential genomic regions under selection.

| DNA extraction, GBS library preparation, and data processing
Foliage was collected from 6 to 9 trees within each of 18 populations of giant sequoia distributed throughout the range (Figure 2).
To reduce the potential of sampling related individuals, we aimed to sample individual trees >40 m apart. However, this was not possible in some small and highly clustered populations. In this latter case, we attempted to maximize the distance between sampled individuals, with the exception of the PLAC population, where all individuals were sampled. Our goal was to maximize the capture of variation across the range of our study species. Thus, we prioritize increased sampling of populations across the S. giganteum range, with the trade-off of limited sampling within each population. Appropriate permits were obtained for all sampling.
High-purity genomic DNA from 143 individuals was isolated from leaf tissue using Plant/Fungi DNA Isolation kits (Norgen Biotek). We constructed three sequencing libraries using a double-digest restriction enzyme-associated genotyping-by-sequencing (GBS) protocol outlined in Peterson, Dong, Horback, and Yong-Bi (2014). Genomic DNA was digested using SbfI and EcoRI restriction enzymes (New England Biolabs). The resulting product was ligated to barcoded adapters and purified, and 46-48 individuals per library were then pooled and subjected to PCR amplification using Phusion High-Fidelity PCR Kit (New England Biolabs) and an automated size selection for fragments between 430 and 570 bp using Pippen Prep. The resulting three libraries were sequenced on an Illumina HiSeq 4000 platform using 150 bp pair-end reads. Sequence data were then demultiplexed using the process_radtags module within the STACKS pipeline (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013), during which reads with a phred quality score < 10 were removed.
Sequences were then aligned to the giant sequoia reference genome v1.0 (Redwood Genome Project, 2019), using the software Bowtie 2 and SAMtools (Langmead & Salzberg, 2012;Li et al., 2009). Variable sites were called using FreeBayes (Garrison & Marth, 2012) and filtered to remove low-quality reads, potential sequencing errors, and paralogs. Data filtering steps included removing loci with uneven mapping quality and those with average read depth >200, requiring a minimum read depth of 5× and a minor allele count >3, removal of loci with more than 80% missing data, and a thinning step that retains one SNP per DNA fragment to remove potentially linked loci (Appendix S1). This filtering protocol resulted in a final dataset of 1,364 bi-allelic SNPs used for outlier tests and environmental association analyses and to obtain genetic diversity statistics.

| Environmental data
To characterize the climatic conditions for each population, we used the spatial centroid of each population to extract and compile twentyone environmental variables at a spatial resolution of approximately 1 km 2 . Nineteen climate variables were obtained from the WorldClim database (Fick & Hijmans, 2017), and elevation and climate water deficit (CWD) were obtained from the California Basin Characterization Model (Flint, Flint, Thorne, & Boynton, 2013). CWD provides an indication of aridity that is important for Mediterranean climate systems, such as in California (Stephenson, 1998). We conducted a principal component analysis (PCA) on the full environmental dataset (21

F I G U R E 1 A giant sequoia tree in Giant Forest, Sequoia and
Kings Canyon National Park, CA, USA variables) after standardization, to reduce dimensionality in the climate data. We retained the first two axes (hereafter PC1 and PC2), which together explained 82% of the climate variation (Appendices S2 and S3). PC1 was driven predominantly by temperature and elevation variables, with a small contribution from annual and winter precipitation, whereas PC2 was determined mostly by precipitation-related variables and CWD with a minor contribution from variables related to temperature seasonality (Appendices S2 and S3).

| Genomic signatures of selection: F ST outliers and gene-environment association tests
To detect F ST outliers that are candidates for selection, we utilized the Bayesian likelihood approach implemented in BayeScan v.2.1 (Foll & Gaggiotti, 2008). This method scans the genome for highly differentiated SNPs that potentially have been subjected to divergent selection while accounting for neutral genetic structure (Narum & Hess, 2011).
BayeScan was run using the false discovery rate (FDR) set to 0.05 under the following parameters: 20 pilot runs of 5,000 with an additional burn in of 50,000 iterations and a subsequent run with 5,000 iterations and a thinning interval of 10. The prior odds for the neutral F I G U R E 2 Range map of giant sequoia (black) showing the gradient of precipitation of driest quarter (mm) across a section of California. Sampled populations indicated by a population code model were increased to 100 (default is 10) as raising this value has been shown to reduce false positives with little effect on false negatives (Lotterhos & Whitlock, 2014). Loci with log 10 values of the posterior odds >1.0 were retained, as the program documentation suggests these loci show "strong" evidence for selection (Foll, 2010).
To test for associations between genomic variation and environmental factors, we utilized latent factor mixed models (LFMMs; Frichot et al., 2013), as implemented in the LEA package in R (Frichot & François, 2015a). LFMM is a univariate approach that treats each individual locus as a response variable with climate data (PC1 and PC2 separately) as the explanatory variable, while incorporating neutral structure using latent factors (Frichot et al., 2013). In simulation studies, LFMM has demonstrated a good balance between high power and low false-positive rate (Frichot et al., 2013;de Villemereuil et al., 2014). As suggested by Frichot et al. (2013) and Frichot and François (2015a), we used two methods to determine the optimal number of latent factors (K value) that correct for the neutral genetic structure of our data. First, we ran a principal component analysis on the individual allele frequencies. We then determined the number of components that explain the genetic variance, based on the Tracy-Widom test on the eigenvalues, as an estimate of K (Frichot & François, 2015a). Second, we utilized the Bayesian clustering algorithm STRUCTURE that estimates the number of genetic clusters (K) without prior information about geographic origin (Pritchard, Stephens, & Donnelly, 2000). The best K value was determined using the ∆K statistic as suggested by Evanno, Regnaut, and Goudet (2005). We used four replicates and a burnin of 300,000 and 1,000,000 MCMC repeats after burnin for K = 2-12.
We ran LFMM to test for associations between SNP′s and two composite climate variables (PC1 and PC2) using ten independent replications at 50,000 iterations after a burnin period of 25,000 with the number of latent factors (K) ranging from 8 to 12, as the methods outlined above suggested K equal to 10 and 9, respectively. We chose high run length parameters because of the relatively small number of individuals and loci. LFMM uses the z-scores to indicate the strength of the gene-environment association (Frichot & François, 2015b). As suggested by the authors, we calculated the median z-score from ten replicate runs, re-adjusted the p-values, controlled for FDR using the q-value of 0.05, and determined candidate SNPs based on the Benjamini-Hochberg procedure (Frichot & François, 2015b).
We also utilized RDA, a multivariate GEA method, to test for more subtle polygenic signatures of adaptation and detect outlier loci as candidates of functional importance. Redundancy analysis (RDA) is an extension of multiple regression to multivariate response variables (Legendre & Legendre, 2012). In finding the ideal combination of predictor and explanatory variables, RDA has shown high power to detect potential signals of polygenic adaptation (Forester et al., 2018;Harrisson et al., 2017). For these analyses, Hellingertransformed allele frequencies (Legendre & Gallagher, 2001) were treated as response variables. Because RDA models do not allow missing data, we imputed allele frequency data using probabilistic principal component analysis (ppca) as implemented in the "pcaMethods" Ppca uses a decomposition of SNP frequencies to create principal components; the components with the largest eigenvalues are then used to impute the missing data. We evaluated space and climate as explanatory variables. Space was defined by distance-based Moran's eigenvector maps (dbMEMs; Borcard & Legendre, 2002;Dray, Legendre, & Peres-Neto, 2006) based on Euclidean distances between all giant sequoia groves (sampled and unsampled) and extracting the values that correspond to our sample sites. Then, we conducted backward model selection, using the "ordistep" function within the vegan package for R (Oksanen et al., 2013), to reduce the number of dbMEM vectors. For climate, we reduced the twenty-one untransformed environmental variables described above, first by removing highly correlated environmental variables, (|r| < 0.7), and subsequently by using the "ordistep" function for backwards model selection to remove variables lacking explanatory power. The above process resulted in climate being represented by "isothermality" (ISO), a measure diurnal and annual temperature fluctuation, "precipitation of driest quarter" (PDQ), a measure of summer precipitation in Mediterranean climates, and "climate water deficit" (CWD), a measure of aridity, in all RDA models, and space represented by two dbMEM vectors, MEM3 and MEM5. All variables were centered and standardized before use in each model.
We set up multiple RDA models to determine the relative amount of variation in allele frequency explained by climate after correcting for geographic space as a signature of local adaptation (Harrisson et al., 2017;Lasky et al., 2012;Sork et al., 2016). First, to elucidate the major factors shaping genetic variation and to detect potential signals of local adaptation, we set up three models for comparison: a full RDA model where allele frequencies were associated with both climate and spatial explanatory variables, a partial RDA in which the effects of climate were conditioned on geography (dbMEMs), and a second partial RDA, where the effects of geography were conditioned on climate. Next, to detect outlier loci, allele frequencies were associated with climate predictors after removing the effects of spatial predictors (Forester et al., 2018;Harrisson et al., 2017;Lasky et al., 2012). Using the first constrained axis, we identified candidate SNPs of potential adaptive importance as those with loadings in the tails of a 95% confidence interval from the mean or 2.0SD from the mean loadings. One risk of using such a low cutoff is an elevated rate of false positives.
However, we chose this to maximize the number of SNPs detected, as we did not expect to find single loci that would be under very strong selection for climate variation in the range of giant sequoia. Moreover, we also identified the climate predictor with the highest correlation to each indicated SNP. In all RDA models, we assessed model and constrained-axis significance using 999 permutations.

| Genomic context of outlier loci
To gain insights into the potential adaptive significance of outlier loci, we obtained the flanking sequence of each outlier SNP locus from the giant sequoia reference sequence (Redwood Genome Project, 2019). Since the giant sequoia reference genome is not annotated, functional annotation was performed using the online BLAST (Basic Local Alignment Search Tool) database. Using a 601bp sequence (300 bp upstream and downstream of the SNP site), we searched the NCBI database using BLASTn with an e-value cutoff set to 1 × 10 5 and the requirement of >70% sequence similarity.

| Genetic diversity and differentiation
Genetic diversity and differentiation differed substantially across the eighteen sampled populations (   Hierarchical AMOVA found a small, but significant variance due to regions (Fct = 0.02, p = .000) and a larger portion of genomic variation distributed among populations (F ST = 0.15, p = .000; Appendix S4).

| Candidate genomic regions associated with climate variables
Univariate environmental association analyses (LFMM with K = 9) indicated a total of 58 loci with strong correlations to composite TA B L E 1 Population information and genetic diversity summary statistics calculated for each population. Diversity statistics calculated without minor allele filtering are noted within parentheses.  (Figure 4; Appendix S6). Of these, 51 were correlated with PC2 that was predominantly driven by precipitation, and seven were correlated with the temperature-driven PC1 ( Figure 4; Table 2; Appendix S6). An examination of the adjusted P-values from all runs (K set from 8-12) provided additional support for K = 9 (Appendix S7; Frichot & François, 2015b). Our BLAST analysis was successful in finding functional annotation for three loci that were strongly associated with PC1 and 11 loci that were associated with PC2 (Table 2).
We used the partial RDA model to detect outlier loci as candidates of importance in selection in a multivariate context, where allele frequencies were associated with climate after removing the effect of spatial predictors. Using the SNP loadings on the first RDA axis, we identified 24 outlier loci beyond the 95% confidence that demonstrated strong correlations to environmental variation, all of which were most correlated with precipitation of the driest quarter (PDQ) (Figure 5; Appendices S6 and S8). Our annotation procedure supported functional importance for seven of the 24 loci (Table 2).

| Concordance among tests for signatures of selection
Overall, eight loci were detected as outliers by both RDA and LFMM.
All of these loci were most associated with precipitation-related variables (PC2 in LFMM and PDQ in RDA). Annotation through BLAST identified two of these loci as having a putative function (Table 2).
Overlap was found between a BayeScan (F ST outlier) and LFMM at one locus (1,123; Table 2).

| Partitioning variation between climate and geographic space
The full RDA model explained 45% of the total variation in allele frequency and supported an influence of climate and/or space in shaping allelic variation (p = .001; adjusted R 2 = .22). The first two canonical axes from the full RDA model were significant (p = .002, and 0.013 respectively) and together accounted for 77% of the explained variation ( Figure 5). The partial RDA model, with climate conditioned on space, was significant (p = .019; adjusted R 2 = .09) and constrained 49% of the variance explained by the full model. The first partial RDA axis was significant at the 0.1 level (p = .098) and accounted for 45% of the variation. The partial RDA model with space conditioned on climate accounted for 24% of the explained variation and was nonsignificant (p = .207). The remaining 27% of the explained variation was confounded between climate and geography.

| D ISCUSS I ON
Giant sequoia is a paleoendemic of California that has likely suffered from a long-term demographic decline (Dodd & DeSilva, 2016).

Today, it is limited to a number of restricted groves in the Sierra
Nevada mountain chain. Small grove sizes and limited gene exchange among populations (DeSilva & Dodd, 2020) might be expected to limit its adaptive potential through inbreeding effects and genetic drift. However, through different approaches we have found evidence for a signal of spatially varying local adaptation associated with climate variables and, in particular, along gradients dominated by precipitation. We report here that population genetic structure in giant sequoia has been shaped by local adaptation overlain on historical population processes. From our study of genomic variation, we detected 79 loci as either F ST outliers or loci with strong associations to climate as candidate regions of adaptive importance. Of these, we highlight 26 SNPs, found from multiple methods, or that correspond with functional annotation, as prime candidates for additional research. We emphasize that these outlier loci may include false positives and that experimental studies are needed to demonstrate functional significance of putative adaptive genomic regions (Barrett & Hoekstra, 2011;Kawecki & Ebert, 2004). Here, we present a first step toward understanding local adaptation in an iconic forest tree.

| Evidence for local adaptation
Despite the strong structure among populations, analysis of our genomic data revealed a signal of divergent selection associated Local adaptation is further supported by nine loci detected by multiple methods and 19 candidate loci with functional annotation (Table 2). Overlap in the detection of outlier loci has been reported in numerous field studies (Harrisson et al., 2017;Hess, Zendt, Matala, & Narum, 2016;Sork et al., 2016). A carefully designed simulation study demonstrated that overlap between GEA methods was found more often for actual targets of selection rather than false positives (Forester et al., 2018). In addition, two of the eight loci detected by both RDA and LFMM have relevant functional annotation, (loci 827 and 1,229; Table 2). A BLAST search suggests that Locus 827 is a kinesin-like protein, KIN-13A, which has been found to be involved in trichome morphogenesis (Lu, Lee, Pan, Maloof, & Liu, 2005 Sletvold & Ågren, 2012). Locus 1,229 represents a potential pollen allergen gene, which is thought to be involved in plant responses to stress (Chen et al., 2016). Although our BLAST search suggested functional importance for 20 loci, many of the annotated regions are characterized only as mRNA with further functional roles yet to be determined ( Table 2). The non-annotated outlier loci are promising candidates for future research as they may be of unknown importance, linked to adaptive genes, in regulatory regions, or represent false positives. Any future annotation of the giant sequoia genome will provide valuable clarity as to the specific role of all outlier loci.
Yet, we emphasize the candidate loci noted in this study demonstrate only strong associations with climate and identifying the exact targets of selection involves rigorous experimental research. Taken together, outlier loci with functional annotation or those detected by multiple methods provide strong support for adaptive variation across the range of giant sequoia.

| Outlier SNPs driven by precipitation
Interestingly, variables associated with precipitation appeared to be the major drivers of local adaptation, which perhaps reflects the strong gradients of water relations on the western slope of the Sierra Nevada. Using gene-environment association methods (LFMM and RDA), we found evidence for adaptive differentiation across giant sequoia populations in response to gradients in precipitation and a more limited signal of local adaptation to temperature.
For Mediterranean type climates, PDQ is of particular relevance for desiccation sensitive species as it equates to summer precipitation, which may represent a vital water source for giant sequoia seedlings during a vulnerable establishment phase.

| Limitations, opportunities, and future implications
It is important to note that GEA methods can suffer from low power or high false-positive rate under some demographic scenarios. Although LFMM has been shown to be robust to various demographic scenarios, including those that create high levels of population structure, this method can have elevated false discovery rates (FDR) under scenarios that create IBD (Forester, Jones, Joost, Landguth, & Lasky, 2016;Lotterhos & Whitlock, 2015;de Villemereuil et al., 2014). Here, we do not find a significant signal of IBD in our data. Yet, the role of IBD or other neutral factors affecting population structure in giant sequoia has not been fully elucidated.
Previous research has indicated isolation by distance (IBD) and/or ancient divergence separating the northern and southern popula- Here, a large portion of the explained variance in our data (27%) was confounded between climatic variation and geographic space, which is perhaps due, to the strictly north-south range of giant sequoia, making it inherently difficult to decouple distance from environmental gradients that vary latitudinally. Therefore, LFMM results should be treated with some caution due to the potential contribution of IBD to population structure. In contrast to LFMM, RDA models show high power and low false-positive rates under IBD (Forester et al., 2016). Simulation studies indicate that the performance of RDA also remains high when population structure and selective gradients are explicitly correlated (Capblancq et al., 2018). Yet, RDA is not without limitations, as it can have low power under island demographic models (Forester et al., 2018) or when selective pressures are highly clustered (Capblancq et al., 2018). Given that each GEA method has particular limitations, we believe the outlier loci detected by both LFMM and RDA, as well as outliers with functional annotation remain strong candidate regions of adaptive importance.
There is an ongoing need for future studies to provide additional clarity on the distribution of adaptive variation and genetic architecture of local adaptation in giant sequoia. To our knowledge, this study represents the first investigation of adaptive variation using genome-wide data. Yet, the results presented here are based on a small subset of the genomic variation within the species, as the giant sequoia genome is very large (8.5 Gb, Redwood Genome Project). More comprehensive sampling of the genome as well as an incorporation of phenotypic information will greatly improve our understanding of local adaptation in this species. In addition, future annotation of the giant sequoia genome will provide opportunities to better understand the genetic underpinnings of many phenotypic traits.
A trend toward increased aridity along midelevation Sierra Nevada forests could undermine the long-term persistence of giant sequoia. With end-of-century predictions for this region that include decreasing snowfall and earlier snowmelt, forests of the Sierra Nevada mountains will likely experience an accentuation of the summer drought that is typical of Mediterranean climates (Fyfe et al., 2017;Stewart, Cayan, & Dettinger, 2004;Sun et al., 2018).
Considering the evidence presented here, we highlight the potential that increased water stress may create maladaptation of giant sequoia populations to their environment. It has been suggested that a long-term decline (over the last ~2My) of giant sequoia is tied to increasing aridity during the development of current climate regimes (Dodd & DeSilva, 2016). Moreover, some giant sequoia populations suffered extensive foliage die back during the drought period from 2012 to 2016 (Stephenson et al., 2018), providing further indication of sensitivity of giant sequoia to arid conditions.

| CON CLUS IONS
We provide evidence of local adaptation along gradients of pre-

ACK N OWLED G M ENTS
We are grateful to Matthew Hughes for assistance with leaf collections and to Prahlada Papper for vital help with data processing. This research was partially funded by Mcintire Stennis grant CA-B-ECO-0165-MS awarded to R.S.D. We also thank reviewers whose feedback helped improve the manuscript.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The full SNP data are available through DRYAD https://doi.