Contrasted patterns of local adaptation to climate change across the range of an evergreen oak, Quercus aquifolioides

Abstract Long‐lived tree species are genetically differentiated and locally adapted with respect to fitness‐related traits, but the genetic basis of local adaptation remains largely unresolved. Recent advances in population genetics and landscape genomic analyses enable identification of putative adaptive loci and specific selective pressures acting on local adaptation. Here, we sampled 60 evergreen oak (Quercus aquifolioides) populations throughout the species' range and pool‐sequenced 587 individuals at drought‐stress candidate genes. We analyzed patterns of genetic diversity and differentiation for 381 single nucleotide polymorphisms (SNPs) from 65 candidate genes and eight microsatellites. Outlier loci were identified by genetic differentiation analysis and genome–environment associations. The response pattern of genetic variation to environmental gradient was assessed by linear isolation‐by‐distance/environment tests, redundancy analysis, and nonlinear methods. SNPs and microsatellites revealed two genetic lineages: Tibet and Hengduan Mountains–Western Sichuan Plateau (HDM‐WSP), with reduced genetic diversity in Tibet lineage. More outlier loci were detected in HDM‐WSP lineage than Tibet lineage. Among these, three SNPs in two genes responded to dry season precipitation in the HDM‐WSP lineage but not in Tibet. By contrast, genetic variation in the Tibet lineage was related to geographic distance instead of the environment. Furthermore, risk of nonadaptedness (RONA) analyses suggested HDM‐WSP lineage will have a better capacity to adapt in the predicted future climate compared with the Tibet lineage. We detected genetic imprints consistent with natural selection and molecular adaptation to drought on the Qinghai–Tibet Plateau (QTP) over a range of long‐lived and widely distributed oak species in a changing environment. Our results suggest that different within‐species adaptation processes occur in species occurring in heterogeneous environments.


| INTRODUC TI ON
Climatic oscillations can have a profound impact on biodiversity and sustainability at every level of the biota, from genes to ecosystems (de Lafontaine, Napier, Petit, & Hu, 2018;Hewitt, 2000;IPCC, 2013). In response to rapid climate change, natural populations can migrate to new favorable locations, adapt locally to novel environments, or become extinct (Aitken, Yeaman, Holliday, Wang, & Curtis-McLane, 2008;Savolainen, Pyhäjärvi, & Knürr, 2007;Sork et al., 2010). Among these responses, local adaptation to altered environments seems especially crucial for sedentary and long-lived organisms such as forest tree species, because moving to favorable locations through propagule dispersal might not be fast enough to cope with the rate of ongoing climate change (Corlett & Westcott, 2013;Hughes, Inouye, Johnson, Underwood, & Vellend, 2008;Kremer et al., 2012;Sork et al., 2013).
Long-lived tree species with large population size are typically locally adapted to different conditions within highly heterogeneous environments (Savolainen, 2011). Recent advances in ecological genomics of nonmodel species have started to unravel the molecular genetic basis for local adaptation in tree species (Sork et al., 2013). Numerous methods are available, and all of them ultimately rely on the rejection of the neutral model of evolution to detect genomic imprints of natural selection (Hoban et al., 2016).
One approach infers natural selection by comparing the strength of population structure among loci (F ST outlier analysis) (Lotterhos & Whitlock, 2014). Another approach estimates the strength of associations between environmental variables in structuring genetic variation by relying on genotype-environment associations (GEAs), tests of isolation-by-environment, constrained ordination techniques, or generalized dissimilarity modeling (Sork, 2018). Jointly using these various methods to study the molecular imprint of local adaptation is likely to provide the best inferences because these approaches are conceptually distinct, with different set of assumptions and computational frameworks.
The Qinghai-Tibet Plateau (QTP) in southwest China is the largest and highest plateau (average elevation > c. 4,000 m) in the world (Zhang, Ye, & Sun, 2016), with the Himalaya-Hengduan Mountains considered as one of the World's biodiversity hotspots (Marchese, 2015;Myers, Mittermeier, Mittermeier, Da Fonseca, & Kent, 2000). In recent years, numerous phylogeographical surveys have examined more than 80 plant species in the QTP and its surroundings (reviewed by Yu et al., 2018). These studies focused on neutral genetic variation to investigate population dynamics and biogeographic history, while overlooking plant adaptive potential in relation to the local environment. Yet, organisms within natural populations inhabiting the QTP are likely to be locally adapted to the high altitude, an extreme environment of the plateau. Rapid desertification, directly caused by an increasingly drier and warmer climate (Xue, Guo, Han, Sun, & Liu, 2009) in the QTP, should trigger detectable adaptive responses in plants. Specifically, drought (water deficit) represents a critical condition, which imposes a strong selective pressure fostering rapid local adaptation in plants (e.g., Eveno et al., 2008;Petit et al., 1999). Investigating drought adaptive response of long-lived trees is timely because most climate-change scenarios suggest a general increase in aridity worldwide (e.g., Park, Sur, Kim, & Lee, 2018).
Quercus L. (oak) is one of the most diverse and ecologically important tree genera in the Northern Hemisphere, with high species diversity in Central America and South-East Asia (Denk, Grimm, Manos, Deng, & Hipp, 2018). Because of its tolerance to a wide range of environments and contrasting habitats (Xu, Dimitrov, Shrestha, Rahbek, & Wang, 2019), the oak genus has provided useful insights on evolutionary mechanisms underlying local adaptation (Kremer, 2016;Petit et al., 2013). Landscape genomic studies have revealed patterns of local adaptation in European white oaks (Quercus petraea, Q. pubescens, Q. robur;Rellstab et al., 2016) and North America valley oak (Q. lobata; Sork et al., 2010Sork et al., , 2016. The recent release of a high-quality assembled genomic sequence of Q. robur (Plomion et al., 2018) (Lyu et al., 2018). Studies specifically addressing local adaptation of oaks in Asia are thus urgently needed.
The Himalaya-Hengduan Mountains biodiversity hotspot of the QTP hosts no less than five evergreen highland oak species (Meng et al., 2017). Among these, the endemic Quercus aquifolioides is the most widely distributed and occupies the highest elevation, reaching the tree line in some areas (Du, Hou, Wang, Mao, & Hampe, 2017).
According to an updated classification of oaks, Q. aquifolioides belongs to subgenus Cerris in section Ilex (Denk et al., 2018). The species displays exceptional environmental adaptations, including broad tolerance from subtropical humid to cold dry environments in steep high solar radiation slopes of rugged mountains, at an elevation ranging from 1,900 to 4,600 m a.s.l. (Huang, Zhang, & Bartholomew, 1999;Tang, 2006). Because of its wide geographical distribution and ecological amplitude, Q. aquifolioides might provide useful insights on the genetic mechanism of adaptive variation. Here, we relied on a joint analysis from state-of-the-art approaches to identify genetic imprints consistent with natural selection and molecular adaptation across the range of Q. aquifolioides. To this end, we targeted candidate gene sequences known to be involved in adaptation to drought stress in other oak species from the EvolTree database (http://www.evolt ree.eu/index.php/candi date-genes -db). We hypothesized that the combination of F ST outlier detection methods, SNP annotations, and various procedures testing associations between genetic and environmental variables could identify molecular signatures consistent with local adaptation to drought, a major climatic driver limiting plant growth in the Himalaya-Hengduan Mountains biodiversity hotspot. Furthermore, in order to determine how adaptation would contribute to Q. aquifolioides response to ongoing climate change, we predicted the "adaptedness" (i.e., adaptive capacity; Foden et al., 2013) of Q. aquifolioides to future local climate F I G U R E 1 Geographical distribution, location of study sites, and mapping of the Bayesian genetic clusters in Quercus aquifolioides. (a) Species range (gray shadow) and sampling locations (black dots). (b) Two main clusters indicated by red (Tibet lineage) and green (HDM-WSP lineage) based on Bayesian cluster analysis of eight nuclear microsatellite markers from 587 individuals in 60 populations based on climate data modeled under a scenario of future climate warming. Our study lays the groundwork for further investigations unraveling ecological adaptation of oaks in Asia and provides new insights into local adaptive response to climatic selective pressure in tree species.

| Field sampling, DNA isolation, and microsatellite genotyping
We sampled leaf material from 60 study sites spaced >30 km apart and covering the range of Q. aquifolioides in southwest China as described in the Chinese Virtual Herbarium. At each site, foliage of 7-11 individuals spaced >100 m apart was sampled and stored in silica gel, reaching a total of 587 individuals ( Figure 1; Table S1).
Total genomic DNA was isolated from dry leaf tissue using a Plant Genomic DNA Extraction Kit (Tiangen). All 587 individuals were genotyped at eight neutral nuclear microsatellite (SSR) markers described by Du et al. (2017). This dataset served as a control, indicative of neutral genetic variation in our sample.

| Climatic variables
For each study site, a total of 103 climate variables for current conditions (~1970-2000) and future predictions (2050 and 2070) were extracted from WorldClim Version2 raster layers at 30 s (~1 km 2 ) of resolution (Fick & Hijmans, 2017; http://www.world clim.org/ver-sion2). These included the full suite of 19 mean annual bioclimatic variables along with average monthly climate data for precipitation, solar radiation, wind speed, water vapor pressure, and minimum, mean, and maximum temperature. To avoid biased estimates of model coefficients and spurious significance levels resulting from multicollinearity, we excluded highly correlated climate variables with the threshold values of 0.7 using a variance inflation factor (VIF) test in "usdm" R package (Naimi, Hamm, Groen, Skidmore, & Toxopeus, 2014;R Core Team, 2018). Along with the two geographic variables (latitude and longitude), four climatic variables-isothermality (bio03, a measure of monthly temperature oscillation relative to annual temperature oscillation), mean temperature of the driest quarter (bio09), precipitation during the dry season (prec01), and precipitation during the wet season (prec06)-were finally retained for subsequent analyses (Table S1).

| Candidate gene fragment selection and annotation
A set of 180 drought-related stress candidate genes in Quercus petraea and Quercus robur were downloaded from the EvolTree database (Table S2; http://www.evolt ree.eu/index.php/candi date-genes -db). We used Primer3web (Untergasser et al., 2012;http://prime r3.ut.ee) to design PCR primer pairs for each candidate gene. We tested each pair of primers on eight individuals under various PCR conditions. The PCR products were visually checked on agarose gel and validated by the Sanger sequencing (ABI X3730XL DNA analyzer). We retained 65 polymorphic gene fragments for further genotyping. Sequences were submitted to GenBank and used as a reference in further analyses.
We designed pair-end barcode sequences to distinguish each individual after pool sequencing ( Figure S1). Specifically, 20-bp barcode sequences (10 bp was connected to the 5′ end of forward and reverse primers, respectively) with a 3-bp difference were designed to distinguish among 100 individuals through pairwise coupling. We amplified the 65 candidate gene fragments and pooled PCR products in equimolar ratios into six libraries, each including PCR products from ca. 100 individuals. Libraries were sequenced using an Illumina MiSeq System with a paired-end sequence of 250-nt reads, at the Center of Biomedical Analysis (Tsinghua University).
Illumina raw reads were split according to barcode information.
Low-quality sequences and Illumina-specific adapters were removed by Trimmomatic 0.36 (Bolger, Lohse, & Usadel, 2014), and the clean reads were mapped onto the original reference sequences by DNA mapping algorithm Burrow-Wheeler Aligner, BWA-MEM (Li & Durbin, 2009). Duplicates were marked in each aligned file by Picard Tools (http://broad insti tute.github.io/picard), and the reads around any insertion/deletion (indel) were realigned by the Genome Analysis Toolkit (GATK; Mckenna et al., 2010). The initial SNPs were called, filtered, a new round of SNP calling was performed based on the initial SNPs, and the desired phased haplotypes were produced from our variants. Using VCFtools (Danecek et al., 2011), we retained high-quality (quality score >30) biallelic SNPs with a minor allele frequency (MAF) ≥2.5% (Rellstab et al., 2016).
We downloaded amino acid sequences for European white oak, Q. robur (Plomion et al., 2018) from the OAK GENOME SEQ UENCING website (http://www.oakge nome.fr/), and performed a blastx search of the candidate genomic sequences in Q. aquifolioides against Q. robur amino acid sequences. If the high scoring segment pairs (HSPs) of the blast top hit contained no stop codons, we assumed that the segments were likely to be coding sequences (CDSs). We used SnpEff (Cingolani et al., 2012) to annotate likely gene functions of the variants.
We estimated mean frequency of the most frequent allele at each locus (P), mean observed heterozygosity (H O ), mean expected heterozygosity (H E ), and mean nucleotide diversity (π) for each lineage based on all SNPs and non-F ST outlier SNPs using STACKS 1.47 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013).
For SSRs, the mean effective population size (N E ), mean observed heterozygosity, mean expected heterozygosity, and mean Shannon index (I) for each lineage were calculated by GenAlEx 6.5 (Peakall & Smouse, 2012). For each summary statistic, 2-group Mann-Whitney U test was used to evaluate the significance of the difference between the two lineages.
Contrary to F ST outlier and GEAs, multivariate approaches do not search for molecular imprints of locus-specific adaptive effects.
Instead, they can provide insights on the role of adaptation by testing the multivariate relationships between environmental gradients and genetic structure across populations, while accounting for spatial genetic structure caused by neutral evolutionary processes.

| F ST outlier detection
We combined two methods to identify outlier loci with extremely high F ST values, assumed to reflect a locus-specific imprint of diversifying selection. First, BayeScan relies on a Bayesian approach that directly estimates the posterior probability that a given locus is under selection, taking the effective population size and migration rate into account, thus reducing false positives (Foll & Gaggiotti, 2008). We computed 20 pilot runs with a run length of 5 × 10 3 and a burn-in of 5 × 10 4 followed by 5 × 10 3 iterations with a thinning interval of 10.
Second, FDIST2 implemented in Arlequin 3.5 was used to establish, based on 10 5 simulations, a null distribution of F ST across loci as a function of heterozygosity, reflecting migration-mutation-drift balance with no selection. Outlier loci correspond to those SNPs with observed F ST values falling outside the 95% confidence interval (CI) neutral envelope (Beaumont & Nichols, 1996;Excoffier et al., 2009).

| Genome-environment association (GEA) outlier detection
We integrated two approaches to identify GEA outlier loci displaying significant statistical associations between allele frequency and climate variables. Bayesian generalized linear mixed models (BayEnv, Coop et al., 2010) relies on a Bayesian approach to test whether environmental factors improve fit over a null model (Günther & Coop, 2013). First, the non-F ST outlier SNP subset was used to estimate the null model of allele frequencies covariance across populations. Second, a Bayes factor (BF) is calculated for each SNP, representing the ratio of probabilities of a correlation between the allele frequency variation in the SNP and an environmental factor versus the null model given by the covariance matrix alone. Ten independent runs of 10 5 MCMC cycles were produced for each SNP, and the means of the results were used to estimate correlations between allele frequencies and environmental variables. SNPs with log 10 (BF) >0.5 have strong support for associations between allele frequencies and environmental gradients (Jeffreys, 1998).
Latent factor mixed models (LFMMs) can improve detection accuracy using latent variables and taking the population structure into account (Frichot et al., 2013). Ten independent runs with 10 5 iterations after a 5 × 10 4 burn-in step were produced to compute correlations between allele frequencies and climate variables in "LEA" R package (Frichot & François, 2015). The latent factors were set to two for analyses including all populations and one for analyses within each lineage (Tibet and HDM-WSP). SNPs with p-values < .05 are significantly associated with climate.

Linear relationships
Mantel tests of isolation by distance (IBD) and isolation by environment (IBE) were performed to assess linear relationships between geographic distance (IBD; pairwise Euclidean distance) or environmental distance (IBE; Bray-Curtis distance) and genetic distance (pairwise F ST ) using R package "ecodist" (Goslee & Urban, 2007). To disentangle the effect of IBD and IBE, a partial Mantel test was used to measure IBD by controlling the impact of environment and IBE by controlling the impact of geography. In addition, multiple regression on distance matrix (MRM) analysis was performed to test a multiple regression of genetic distance on geographic and environmental distances. Statistical significance of Mantel tests and MRM was evaluated from 10 4 permutations.
Redundancy analyses (RDAs) and partial redundancy analyses (pRDAs) were used to detect linear relationships between genetic variations (from SNPs and SSRs datasets) and multivariate climatic gradients, using "vegan" R package (Oksanen et al., 2017). SSRs and SNPs with the four climate variables were considered as constrained factors and geographic variables (longitude and latitude) as conditioned factors. Statistical significance was evaluated from 999 permutations.

Nonlinear relationships
Generalized dissimilarity modeling (GDM, Ferrier, Manion, Elith, & Richardson, 2007) was performed using the "gdm" package (Ferrier et al., 2007) to identify nonlinear relationships between genetic distance and environmental and geographic distances by fitting splines (Fitzpatrick & Keller, 2015). Geographic distance was based on Euclidean distance between sites, and the matrices of genetic distance were based on the mean allele frequency of each SNP. Spline shape and height describe the allelic compositional change along the environmental gradient and the importance of the environmental variable, respectively.

| Risk of nonadaptedness under future climatic scenarios
We performed a risk of nonadaptedness (RONA, Rellstab et al., 2016) to predict the adaptive potential of Q. aquifolioides to future local climate using default settings in PYRONA v0.

| Genetic structure and diversity
The 65 (Table S2). Of the 65 successfully amplified polymorphic loci, 64 had at least one blast hit against Q. robur amino acid sequences (Table S3). We identified 71 coding sequences (CDSs) from 60 loci (Table S4). From the 51,623,432 paired-end Illumina reads obtained from six sequencing libraries, 831 high-quality biallelic SNPs were called by the GATK pipeline.
Principal component analysis using all SNPs clearly assigned the 587 sampled individuals to two distinct genetic lineages (Tibet and HDM-WSP) that could not be distinguished on the basis of neutral markers (non-F ST outlier SNP subset or SSRs) alone (Figure 2).
Using all SNPs, AMOVA indicated 13% of the molecular variance is explained by grouping populations in the two lineages (Table 1).
Overall, based on non-F ST outlier SNPs and especially SSRs, diversity indices indicate lower level of genetic diversity in Tibet than in HDM-WSP (most p-values ≤ .05; Table 2, Table S6).  Table S7).

| GEA outlier detection
Details of GEA outlier SNPs are shown in Figure 3b, Figures S3 and   S4, and Tables S7, S8 and S10. In the Tibet lineage, 19 candidate TA B L E 1 Hierarchical analyses of molecular variance (AMOVA) of Q. aquifolioides populations based on all SNPs, non-F ST outlier SNPs, and SSRs long-chain acyl-CoA synthetase 4, LACS4; Figure S5). In total, 12 candidate genes containing GEA outlier SNPs were detected by both GEA approaches in the two lineages ( Figure 3b; Table S7). Among these, 28 SNPs were significantly associated with the same climatic variable (prec01), and the remainder were associated with the other three climate variables (bio03, bio09 and prec06) (Table S7).

| Linear relationships
Mantel and partial Mantel tests revealed spatial genetic structure consistent with IBD in Tibet and HDM-WSP as well as for all populations combined. Correlation between genetic and environmental distances, consistent with IBE, was detected in HDM-WSP and for all populations combined, but not in Tibet (Figure 4; Table S11). More specifically, genetic distance was significantly associated with precipitation during the dry season (prec01) in HDM-WSP and for all populations (Table S12). These findings were corroborated by optimal MRM models that yielded similar results (Table S13).
The percentages of variance explained by RDA and pRDA were similar, and we thus report results for pRDA. In Tibet, SNP variation was not associated with geography or climate variables, whereas in HDM-WSP and for all populations, temperature during the driest quarter (bio09) and precipitation during the dry season (prec01) contributed the most to SNP variation (Table 3; Figure S6). In HDM-WSP, 11% of the SNP variance was explained by climate, while geography only accounted for 6%. For all populations combined, 8% of the explained SNP variance was due to climate and 9% to geography (Table 3). For the SSRs, pRDA revealed that precipitation during the dry season (prec01) explained most genetic variance in Tibet.
In HDM-WSP and all populations, prec06 and prec01 were the two environmental variables most explanatory, respectively (Table 3; Figure S7). Partitioning of the total SSR variance revealed that 5%, 3%, and 2% of the explained variance were due to climate, while 3%, 2%, and 4% were due to geography in Tibet, HDM-WSP, and all populations, respectively (Table 3).

| Nonlinear relationships
Geographic and environmental distances for the variables considered with GDM explained 82% and 6% of genetic variation for all populations, 55% and 14% of variation in Tibet, and 13% and 24% in HDM-WSP (Table S14). Hence, while geography was the most important predictor in Tibet and for all populations combined, prec01 was the most important environmental factor related to genetic differentiation in HDM-WSP lineage ( Figure 5). These results were corroborated by I-spline analysis ( Figures S8-S10).

| Risk of nonadaptedness under future climatic scenarios
The most represented environmental factor identified by RONA analyses was the precipitation during the dry season (prec01) in Tibet and HDM-WSP lineages. The RONA value of Tibet lineage was on average higher than that of HDM-WSP lineage under both RCP26 and RCP85 predictions in 2050 and 2070 ( Figure 6; Figures S11 and S12 and Tables S15 and S16), indicating that the HDM-WSP lineage has higher potential to adapt to future local climate than Tibet

| D ISCUSS I ON
Our study integrated numerous population genetics and landscape genomic methods to detect putative adaptive genetic variation in Q. aquifolioides populations throughout the species' range in the Himalaya-Hengduan Mountains. We first relied on F ST outlier and GEA methods to reveal SNP markers consistent with locus-specific imprints of natural selection in each lineage individually and throughout the species range. Based on annotations of these potentially adaptive SNP variants and their predicted coding effects, we can compare the functional differences between genes identified in the distinct lineages and further analyze the effects of these differences on the ecological adaptation of Q. aquifolioides in different lineages. We also performed linear and nonlinear multivariate analyses by integrating genetic and environmental factors to identify complex relationships between SNP frequency gradients with climate variables and determine key climate variables most likely to trigger adaptive response in Q. aquifolioides. We finally predicted the capacity of Q. aquifolioides to adapt to future local climate, based on two modeled scenarios of climate warming.
We found that genetic variation in Q. aquifolioides showed contrasted patterns of local adaptation in the two lineages. In Tibet lineage, the genetic variation was mainly driven by geographical distance, whereas in HDM-WSP lineage, the climatic variables, especially precipitation, were instrumental in shaping a genetic imprint consistent with adaptive genetic variation.

| Patterns of genetic differentiation and diversity in Q. aquifolioides
Results from AMOVA and PCA revealed higher level of genetic differentiation when all SNPs were considered compared with reduced marker datasets including only non-F ST outlier SNPs or SSRs.
Diversifying selection can reduce effective gene flow and increase divergence across populations (Guichoux et al., 2013;Nosil, Funk, & Ortiz-Barrientos, 2009). In addition, gene flow across oak populations is limited due to the complex topography in the Himalaya-Hengduan

F I G U R E 4 Mantel tests of genetic distance (F ST /(1 − F ST )) against geographic and environmental distances in each lineage and for all populations combined
TA B L E 3 Summary and partitioning of the variance associated with climate and geographic variables based on redundancy analysis (RDA) and partial RDA (pRDA) in the Tibet lineage, the HDM-WSP lineage, and all populations combined for SNPs and SSRs   Table S14 Mountains (Meng et al., 2017), and hence, local adaptation in López-Pujol, Zhang, Sun, Ying, & Ge, 2011). The reduced genetic diversity with fewer private alleles in the Tibet lineage likely reflects extinction-recolonization dynamics in this area, whereby only a subset of the Hengduan Mountain refugium gene pool recolonized Tibet (Du et al., 2017;Meng et al., 2017;Qiu, Fu, & Comes, 2011).

| Functional interpretation of outlier candidate gene SNPs
Accurately identifying loci most likely to be under selection is a crucial step in candidate gene studies. All methods present their own limitations and advantages (de Lafontaine et al., 2018;Hoban et al., 2016;Sork, 2018). In the present study, a combination of F ST -based and genotype-environment association (GEA) methods was used, as this should increase the robustness of outlier detection. Another challenge is to find the most likely causes of outlier syndrome. In this respect, the recent release of the oak genome (Plomion et al., 2018) provides useful information to help interpreting gene functions. In addition, thanks to a large number of genomic studies focusing on F I G U R E 6 Risk of nonadaptedness plot (RONA) for three environmental factors under RCP26 (top panels) and RCP 85 (lower panels) prediction models in 2070. Bars represent weighted means (by R 2 value), and lines represent standard error for each population in the Tibet and HDM-WSP lineages. The exact R 2 values are shown in Table S16 drought stress in plants, our results can be contextualized in light of available knowledge about these stress-related genes.
Genotype-environment associations identified 12 genes associated with climate variables (mostly precipitation during the dry season-prec01) in both Tibet and HDM-WSP lineages, suggesting some level of parallel adaptation to drought in the species. Moreover, GEAs identified two genes displaying lineage-specific association with climatic variables in HDM-WSP (CL6004CT6724_02 and CL9715CT14526_03).
Three SNPs within these two genes were associated with temperature during the driest quarter (bio09) and precipitation during the dry season (prec01). These three markers were also flagged as positive F ST  (Shockey & Browse, 2011). In Arabidopsis, LACS4 inactivation resulted in a significant overaccumulation of tryphine lipids and displayed morphological anomalies of the pollen grains (Jessen et al., 2011).

| Contrasted patterns of genetic variation across the range
Castanopsis (Sun, Surget-Groba, & Gao, 2016). Here, we investigated these patterns in two lineages and revealed strikingly contrasted signatures of evolutionary dynamics across the range.

| CON CLUS ION
Combining the results from population genetics and landscape genomic methods, we detected contrasted patterns of genetic variation in response to environmental gradients in Quercus aquifolioides, a naturally distributed forest species. Selectively neutral evolutionary processes (e.g., isolation by distance) are key drivers of the genetic variation in the Tibet lineage, whereas adaptive processes (e.g., isolation by environment) prevail for shaping genetic diversity of the HDM-WSP lineage. More specifically, we identified two genes that were correlated with climate gradients only in the HDM-WSP lineage, supporting lineage-specific signals of local adaptation. Our study thus provides valuable insights on local adaptation in trees and their main environmental drivers. In the changing climatic environment associated with global warming, the HDM-WSP lineage may provide valuable genetic resources for the species. In the near future, more studies integrating genomic, phenotypic, and environmental data are required to gain further insights into the mechanisms of oak adaptation.

ACK N OWLED G EM ENTS
The authors would like to thank Prof. Victoria L. Sork

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

AUTH O R CO NTR I B UTI O N S
FD designed the research; TRW and YYW performed the experiments; FD, TRW, YYW, and SU performed the analysis; FD, SU, and GdL wrote the manuscript; and all authors revised the manuscript.