1 Introduction
Trees profoundly influence global carbon cycles and ecosystem stabilization [1-3]. The effects of anthropogenic climate change on precipitation and temperature patterns have altered the distribution, community composition, and phenology of forest trees worldwide [4-7]. However, the impacts of climate change on focal species vary [8-12]. For example, warming temperatures may have a positive effect on the growth rate of trees from colder environments but a negative effect on trees from tropical regions [10, 13]. Higher temperatures and increased aridity are predicted to increase stress reactions in forest trees under drought conditions, especially tropical and subtropical tree species [14-17].
Studies of the potential response to climate change based on field experiments have emphasized the importance of phenotypic variation for adaptation to local climate and acclimation to drought and cold stresses [18, 19]. However, field experiments are unrealistic for species with long lifespans, such as trees, or endangered species with limited populations. An alternative approach is genotype–environment association (GEA) analysis, which identifies environmental factors that select for genetic characteristics [20, 21]. By predicting the vulnerability of forest trees to climate change and environmental factors that may restrict future distribution, GEA analyses can provide a foundation for conservation projects and forest management [22-24].
Although GEA studies have implied that standing genetic variation may enable some tree species to cope with climate change, the long lifespans and low germination rates of tree species may limit the pace of adaptation to acute and drastic environmental changes [25]. Moreover, the applicability of the results of GEA studies to conservation projects may be confounded by several factors. For example, the geographical distance between source and sink populations and the composition of natural barriers to gene flow may influence pollen dispersal direction and germination rates for introgression between populations [26]. These factors could reduce the spread of beneficial alleles (i.e., genetic rescue effects) [27]. Incorporating landscape analyses into GEA studies could improve the evaluation of genetic vulnerability, assessment of vulnerable populations, and inference of potential dispersal routes under climate change [28, 29].
The tree family Fagaceae is ecologically and economically important and has been widely studied in landscape genetics research [30, 31]. The varied natural habitats of Fagaceae species offer ideal study systems for exploring the effects of environment on genetic diversity, local adaptation, and the response to climate change [32, 33]. Studies have shown that environmental factors, such as precipitation and temperature, significantly affect the adaptive divergence of Fagaceae species [34], but the impacts of other abiotic factors, such as wind, topology, and soil, have largely been ignored [35, 36]. Field experiments have shown that wind and soil factors are prominent drivers of local adaptation of plant morphology and physiology [37, 38], but the potential impacts of these factors on genomic architecture and phenotypic variation across the heterogeneous landscape of Fagaceae species have not been comprehensively established [19, 37, 39, 40]. This is particularly true for Fagaceae species on the subtropical island of Taiwan, where the mountainous terrain and diverse climate have created a range of habitats [41-43].
The rapid development of novel analytic methods in landscape genomics provides unprecedented opportunities to examine new hypotheses and assess vulnerable populations using different statistical assumptions in the context of climate change [2]. If climate change disrupts the allele frequencies that underlie current genetics–environment relationships, vulnerable populations may become less resilient or even extinct locally [44-47]. In this study, we adopted Quercus longinux Hayata, an endemic Fagaceae species in Taiwan, as the study species to evaluate the complex effects of environment on local adaptation and the response to climate change. Q. longinux inhabits mountainous ranges across Taiwan, from low to middle elevations (500 m to 2,200 m above sea level) [48]. Based on its wide distribution along latitudinal and altitudinal gradients, Q. longinux is classified into three varieties (i.e., var. longinux, var. kuoi, and var. lativiolaciifolia), and large morphological variations in fruits and leaves between habitats have been observed [48]. Q. longinux var. longinux and var. lativiolaciifolia grow sympatrically at low to middle elevations in Taiwan. By contrast, Q. longinux var. kuoi is allopatric with the other varieties and is limited to southeastern Taiwan. Compared with the other two varieties, Q. longinux var. kuoi has longer, broader leaves that are green on both sides when fresh and have a non-violet abaxial surface when dried [48]. Environmental variations may influence several leaf traits associated with photosynthetic efficiency and acclimation to abiotic stresses [49-51]. Identifying variations in not only genetic data but also leaf morphology that are associated with environment could shed light on the processes that gave rise to the unique Q. longinux var. kuoi in southern Taiwan.
This work was guided by three objectives. First, we aimed to infer genetic structure and assess the consistency of relationships between genetic data, ecological niches, and phenotypic traits for Q. longinux var. kuoi, a special allopatric variety limited to southern Taiwan. We hypothesized that the distinct habitats in southern Taiwan influence population structures and result in differentiated genetic characteristics that correspond with morphological traits. Second, we aimed to identify environmental features that affect spatial genetic variation and adaptive genetic divergence. We hypothesized that classic temperature and precipitation factors influence genetic divergence and that wind- and soil-related factors participate in adaptation to heterogeneous habitats in Taiwan. Third, we aimed to use adaptive genetic variants to evaluate the vulnerability of populations to climate change. We hypothesized that edge populations and populations subject to the most severe environmental changes will experience the highest genetic offsets and face the highest risk of local extinction.
2 Materials and methods
2.1 Sampling, sequencing, read mapping, and variant calling
We collected 26 populations spanning all known distributional regions of Q. longinux in Taiwan (Fig. S1; Table S1). From each selected tree, a branch with more than ten mature leaves was collected for morphological measurements. Fresh leaves were stored in silica gel at 4°C until DNA extraction.
DNA was extracted using the modified CTAB method based on Doyle [52]. DNA quality and quantity were evaluated with a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States), respectively. After quality control, DNA from 205 individuals was used for dd-RAD library construction. The library was double digested by Sbf1 and Msp1 and ligated to Illumina sequencing adapters with individual barcodes and library indices. Fragments of 250-500 bp were selected and amplified by polymerase chain reaction (PCR). Finally, the fragments were sequenced using Illumina HiSeq X (Illumina, USA) to generate 150-bp paired-end reads.
We trimmed the raw reads to remove adapters, reads < 50 bp, and reads with quality < 30 using fastp [53]. The clean reads were mapped to the reference genome of Q. robur [54] using the mem algorithm of BWA [55]. Aligned reads with a mapping quality score < 20 were discarded. PCR duplicates were marked with the Picard Toolkit (http://broadinstitute.github.io/picard/). We further performed realignment around indels using ABRA2 [56]. Summary statistics, including read number, mapping rate, and coverage, were calculated from SAMtools [57]. Variants were called with BCFtools [58]. We first generated genotype likelihoods at each genomic position with coverage using the BCFtools mpileup command and called the variants with default parameters. After removing indels, we retained biallelic SNPs with missing rates < 0.4, minor allele frequencies > 0.01, genotypes with quality scores >30, and mean minimum depths > 3 as filtering thresholds in VCFtools [59].
2.2 Genetic diversity, genetic structure, and demographic analysis
To evaluate and visualize the genetic clusters of Q. longinux, we first performed principal component analysis (PCA) using a genetic matrix in which missing data were replaced with mean values of each site in the R package adegenet [60]. Next, pairwise FST, population-specific FST, and permutational multivariate analysis of variance (PERMANOVA) were performed with the R packages hierfstat [61] and vegan [62] to evaluate differentiation among defined genetic clusters. A neighbor-joining (NJ) tree was also constructed based on the Nei distances from the R package ape [63]. We then employed StrAuto [64] analysis to assess the pattern of admixture among populations using 100,000 steps of MCMC chains with 25,000 steps burn-in and 10 iterations from K = 1 to 8. The results were uploaded to the CLUMPAK [65] server to generate consensus plots and evaluate the best K according to the ΔK value.
A stairway plot was used to infer historical changes in the effective population size (Ne) of
Q. longinux [66]. We first generated the site frequency spectra (SFSs) from easySFS (GitHub repo:
https://github.com/isaacovercast/easySFS) using all SNPs and selected the projection values with the highest segregating sites for each group. We then ran bootstrapping 1,000 times to estimate median sizes and 95% CIs, assuming a mutation rate of 1.01 × 10
-8 per bp per generation and a generation time of 50 years for oak species [67].
2.3 Environmental variables and GIS processing
Three categories of environmental variables were collected: 105 climate variables, eight soil variables, and four topological variables (Table S2). The climate variables included BIO01-19, seven monthly precipitation factors, temperature, solar radiation, and wind speed downloaded from WorldClim2 [68] and two aridity-associated layers [69]. Eight physicochemical soil characteristics downloaded from SoilGrids 2.0 were used as the soil variables [70]. For topography, the altitude layer was downloaded from GOV.DATA.TW (https://data.gov.tw/), and three altitude-derived layers (i.e., roughness, aspect, and slope) were computed with the R package raster [71]. All variables were resampled to a resolution of 1 km2 for the downstream analyses. To address the issue of multicollinearity, variables with VIF > 10 and Pearson correlation coefficient |r| > 0.8 were discarded, leaving 20 environmental variables among the three categories (Table S3; Supplementary data) for use in the subsequent analyses. We compared differences in the environmental variables between groups using analysis of variance (ANOVA) and Tukey's honestly significant difference test applied in the R package stats [72].
2.4 Assessment of ecological and morphogenetic divergence
2.4.1 Niche differentiation assessment
To harness the power of GEA analysis to evaluate a species’ ability to cope with future climate change, it is crucial to include non-genetic factors that may be associated with fitness to local environments, such as niche characteristics and phenotypic traits. To investigate and quantify the niche overlap and diversification between the three (eastern, western, and southern) genetic groups revealed in our results, we first assessed the overlap of the realized niche using the first two axes from the environmental PCA performed with the 20 retained environmental variables in the R package ecospat [73]. Five hundred pseudo-absence points were added to the analyses as background. In addition, we calculated two niche overlap indices: Schoener’s D and the standardized Hellinger-transformed Warren's I based on occurrence density grids. Next, we conducted the equivalency test with 100 permutations and the background test with 100 runs to evaluate the significance of niche differentiation at p < 0.05. The equivalency test and background test were implemented with the R package ENMtools [74].
2.4.2 Phenotypic trait analysis
To quantify the contributions of genetic and environmental factors to phenotypic variations, we measured 12 morphological leaf traits that are commonly used in morphological studies of Fagaceae and other plants (Table S3; Fig. S2; Supplementary data) [75-77]. A total of 826 dried leaves of 191 individuals from 26 populations (1–5 leaves per individual) were recorded with a digital camera and measured using ImageJ [78]. Because the measured traits may be highly correlated with leaf size, an additional six computational traits that are not typically correlated with leaf size were added [77]. We first performed PCA and PERMANOVA with the six computational traits (i.e., shape-related traits and specific leaf area (SLA)) and abaxial surface color to quantify morphological differences among groups using the R package vegan [62]. Biplots of the attributes and contributions were visualized with the R package factoextra [79]. Differences in traits between groups were compared using ANOVA and Tukey's post hoc test. Next, we performed partial redundancy analysis (RDA) to evaluate the contributions of geography, demography, and environment to phenotypes using all scaled phenotypic traits as the response.
To further investigate the associations between morphological variables and environmental predictors, we constructed univariate generalized linear models (GLMs) with one of the leaf traits as the response (Table S4) and one of the environmental variables as the predictor in the R package MASS [80].
2.5 Detection of local adaptation and underlying gene functions
2.5.1 Identification of environment-associated genetic variants
We integrated two FST-based outlier detection methods and one genotype–environment association (GEA) approach to detect environment-associated SNPs. First, we implemented BayeScan [81] and PCAdapt [82] to identify SNPs with a significant departure of allele frequency from neutrality while controlling for background genetic structure. BayeScan was run with posterior odds (PO) = 100, and SNPs with q < 0.01 were retained as putative outliers. PCAdapt was performed with the number of K principal components (PCs) selected by the decreasing order of the percentage of variation explained by each PC. We further implemented Benjamini–Hochberg correction [83] on the results and retained the SNPs with a false-discovery rate (FDR) < 10% as candidate outliers.
Second, we utilized a univariate approach, the latent factor mixed models (LFMM) algorithm [84, 85], to discover SNPs that were significantly correlated with environmental variables while accounting for population genetic structure. Imputation of missing data and determination of latent factors (K) were performed by the snmf function in the R package LEA [84]. LFMM was performed in the R package LEA [84] using the retained environmental variables with ten runs, and Z-scores from each run were combined. To further adjust the p-values, the Benjamini-Hochberg procedure [83] was used, and FDR < 10% was used to identify putative adaptive outliers. SNPs that overlapped between significant environment-associated outliers and the BayeScan or PCAdapt outliers were regarded as putative adaptive SNPs for use in downstream analyses.
To investigate and compare the roles of geography, demography, and environment in shaping the genetic variation of adaptive and all SNPs, we decomposed the relative contributions of each group of predictors using three methods. First, Mantel tests were used to test for associations between FST/(1-FST) and geographic distance (isolation by distance, IBD) and environmental distance (isolation by environment, IBE) using the R package vegan [62]. Environmental distances were represented by the Euclidean distances of scaled climatic, soil, and topographical variables. Values of pairwise FST between populations were calculated with the R package hierfstat [61]. Differences in FST between adaptive and all SNPs were compared by ANOVA in the R package stats [72].
Second, we performed partial redundancy analysis (RDA) using three predictor datasets: (1) proxies of geographic structure obtained by converting the geographic coordinates to uncorrelated axes using principal coordinates of neighborhood matrices (PCNM) in the R package vegan [62] and retaining half of the positive axes; (2) proxies of demographic history obtained by performing PCA using the genetic matrix of all SNPs and retaining the PC scores from PC1 to PC2; and (3) the 20 retained environmental variables. The genetic metrics calculated from the adaptive or all SNPs were used as the response variables. To further quantify the variation explained by each group of environmental variables, we performed another partial RDA using the 20 environmental variables, which were classified into (1) climatic, (2) soil, and (3) topographical variables, as predictors of the two genetic metrics as the response. The significance of the full models and pure fractions was assessed using 999 permutations with the function anova.cca() of the R package vegan [62].
Third, we used GradientForest (GF), a non-linear and machine-learning approach derived from the random forest algorithm [86]. GF was performed with the genetic metrics of the adaptive or all SNPs. The number of regression trees was set to 500, the maximum number of splits was set using the formula log2 (0.368 × number of individuals/2), and the correlation threshold r was set to 0.5, as suggested in a previous study [86]. The importance of each group of predictors was evaluated by the percentage of weighted R2.
2.5.2 Gene annotation and enrichment analysis
To further investigate the potential gene functions of putative environment-associated SNPs, we retrieved the closest genes with a physical distance < 10 Kbps as potential underlying genes using BEDOPS [87]. The closest genes were further annotated to the Arabidopsis thaliana genome using KOBAS-I [88] with the criterion of e-value < 10-5. After retrieving the homologous genes, gene enrichment analysis was performed with the GO and KEGG databases to detect enriched pathways using KOBAS-I. Significance was assessed with Fisher's exact test [89], and FDR was corrected with the Benjamini and Hochberg procedure [83]. FDR < 5% indicated significantly enriched GO terms or pathways.
2.6 Investigating the influences of current landscape variables on genetic differentiation
To compare landscape characteristics that shaped the differentiation of adaptive or all SNPs, we constructed IBD, IBE, and isolation by resistance (IBR) models using FST/(1-FST) calculated from adaptive or all SNPs as a response. Six predictor matrices were generated, including a geographical Euclidean distance metric (IBD), three environmental Euclidean distances calculated from the climate, soil, and topology factors (IBEclimate, IBEsoil, and IBEtopology), and two circuit-theory-based resistance layers calculated with CIRCUITSCAPE [90] using distribution maps of ecological niche modeling from the factors of climate and topography (e.g., IBRclimate, IBRtopography).
We used two complementary approaches to compare the models. First, reciprocal causal modeling (RCM) [91] was performed to evaluate the relative importance of each model based on the relative values of partial Mantel tests estimated from each focal and alternative model. All partial Mantel tests were computed using the R package vegan [62] with 999 permutations. Second, we implemented the maximum likelihood population mixed-effects model (MLPE) [92] to compare each predictor based on the AIC values. Models with DAIC > 2 were considered to have different model fits. MLPE was implemented with mlpe_rga in the R package ResistanceGA [93]. In addition, we performed an estimated effective migration surfaces (EEMS) analysis [94] to visualize geographic regions deviating from the assumption of IBD. Missing data were imputed and converted to genetic distance using the function bed2diffs. Next, EEMS was executed with RunEEMS_SNPS with a Deme number of 500. A total of 500,000 Markov chain Monte Carlo (MCMC) iterations were run after a burn-in of 150,000 and a thinning interval of 9,999.
2.7 Investigating the impacts of future climate change on genetic adaptedness
2.7.1 Ecological niche modeling
To investigate the current potential distribution range of Q. longinux, we performed ecological niche modeling (ENM) using an ensemble approach based on the true skill statistic (TSS)-weighted combination of six methods in the R package sdm: maxent, glm, svm, gam, mda, and mlp [95]. Five hundred pseudo-absence data were generated and added to the analyses. A total of 60 runs with 10-fold cross-validation analyses were performed. We then used the area under the receiver-operator curve (AUC) to evaluate model performance. Finally, the habitat prediction was transformed into a binary map that classified suitable and unsuitable regions using the threshold of maximum test sensitivity plus specificity.
2.7.2 Genetic offset assessment
We applied five temperature and precipitation variables (i.e., BIO1, BIO3, BIO12, PREC4, and PREC10) selected from the obtained factors with accessible prediction layers to assess the influence of climate change on Q. longinux. We downloaded and averaged the predictions from three future models (CCSM4, MIROC-ESM, and MIROC5) in 2070 to account for model variability. We also considered two contrasting representative concentration pathways (RCPs): a low-emission model (RCP2.6) and a high-emission model (RCP8.5) from CMIP5 for 2070. Three complementary approaches were used to evaluate genetic offset. We first calculated the risk of non-adaptedness (RONA), which represents the theoretical changes needed in allele frequency to maintain the linear environment-genotype relationships with correction on weighted R2, implemented in pyRONA [96]. We retained the top three representative climatic variables with the highest number of significant outliers.
Second, as a complementary method to RONA, we used a random forest algorithm to model the non-linear relationships between adaptive SNPs and the same climatic variables as RONA in the R package GradientForest [86]. PCA was performed on the GF model to visualize the prediction of genetic variation in spatial regions, and the first three principal components (PCs) were assigned to a red–green–blue palette. The genetic offset in the GF model was calculated as the Euclidian distance between the current and future genetic composition at each grid. The binary map generated by ENM was used as a mask to limit prediction within suitable habitats.
Third, following Gougherty, Keller [97], we integrated migration to predict maladaptation to future climate using a generalized dissimilarity modeling (GDM) algorithm [2] with local, forward, and reverse genetic offsets. Local offset represents the predicted in situ change in allele frequencies with no migration. Forward offset was calculated by selecting the minimum offset between each grid in the current range and all grids in the future climate, and reverse offset was calculated by identifying the minimum offset between each grid within the current range in the future and all grids in the current climate. For forward offset, the distance required to migrate to the place with minimum forward genetic offset and the direction (bearing) of migration were also recorded. No dispersal limitation was assumed for forward and reverse offsets to any locations in the suitable habitats. Finally, to simultaneously visualize local, forward, and reverse offsets, we mapped the three metrics as red, green, and blue in an RGB image.
3 Results
3.1 Population structure and demographic history
On average, 70.9% of reads per sample were successfully mapped to the reference genome of Q. robur. We recovered 624,451 SNPs and retained 1,933 high-quality SNPs with an average individual missing rate of 11%. StrAuto revealed three genetic clusters (Fig. 1a): the eastern and western clusters, which are mainly separated by the central mountain range (CMR), and a cluster in the Henchung Peninsula that was limited to southern Taiwan, denoted as HC. Multiple individuals showed admixture with other clusters (Fig. 1b). PCA and the neighbor-joining (NJ) tree produced consistent patterns (Fig. 1c; Fig. S3a-b). FST was higher between genetic clusters than within groups, and divergence was highest between HC and either the eastern or western cluster (Fig. S4c). Similar results were obtained by PERMANOVA (Table S5). According to the stairway plot (Fig. S4d), Ne has differed between HC and the other two clusters since 1 million years ago (Fig. S4d). The apparent population decline of HC 1 million years ago was preceded by a period of larger population size. By contrast, the eastern and western clusters underwent a relatively recent population decline in the last 100,000 years.
3.2 Influence of landscape factors on population genetic divergence
We identified 182 FST outlier SNPs (PCAdapt + BayeScan) and 540 SNPs that were significantly associated with one or more environmental variables by LFMM (Fig. 2). A total of 105 putative adaptive SNPs (identified by the overlap between LFMM and FST outliers) were found in 35 genes (Table S6) that were widely distributed across the genome and did not cluster in specific regions.
FST was stronger and more significant when the adaptive SNPs were used instead of all SNPs (mean FST all loci = 0.09; mean FST outlier = 0.59; p < 0.05), indicating that selection may drive the spatial genetic differentiation between populations. Compared to all loci, the adaptive loci exhibited stronger IBD and IBE patterns (Fig. 3a). Using reciprocal causal modeling (RCM) and a maximum likelihood population mixed-effects model (MLPE), IBE was consistently identified as superior to other competitive models using adaptive SNPs when controlling population structure, indicating that genetic differentiation of the adaptive SNPs was mainly influenced by environmental variation (Fig. 3b-c; Table 1). Furthermore, the divergence of the genetic structure of HC from those of the other clusters was greater when divergence was assessed by adaptive SNPs than by neutral SNPs (Fig. S4a-b). By contrast, IBR based on topographical resistance was selected as the best model by MLPE and RCM using all sites, suggesting that the overall genetic differentiation of the species was mainly influenced by topographical dispersal barriers (Fig. 4a-b; Table 1). The conductance layers generated by the CIRCUITSCAPE algorithm and the estimated effective migration surface (EEMS) analysis also supported topographical resistance across regions, with greater barriers to dispersal at higher elevations (Fig. 4c-d).
The results of partial RDA indicated that pure environment was the most significant variable affecting the genetic variation of adaptive SNPs among three pure-effect portions; whereas pure demography was the most important pure-effect portions of all SNPs when controlling other factors (Fig. S5; Table 2). Forward selection also identified population structure as the most significant variable affecting the genetic variation of adaptive SNPs (Table S7). Environmental variables contributed to 80% of the explained variation in adaptive loci but 60% of the explained variation in all loci (Fig. S5). For adaptive loci, the three groups of variables explained a large proportion of the joint effects (45% of explained variation) (Fig. S5; Table 2). Another partial RDA that decomposed the contributions of climate, soil, and topography to the genetic variation in either the adaptive or all SNPs revealed that pure climate and soil variables explained significantly more of the genetic variation in adaptive SNPs than topography variables (Fig. S5). In agreement with the results of partial RDA, GF analysis demonstrated that demography was the top variable with the highest weighted R2 (Table S8; Fig S6), although geography explained a higher proportion of variation after summing the total weighted R2 (demography: 32%; geography: 47%).
3.3 Functional annotation of adapted loci
Several Q. longinux loci involved in climate and soil adaptation were associated with environmental variables (Table S6). Two functional pathways were significantly enriched in adaptive loci: oxidative phosphorylation and photosynthesis.
Furthermore, the allele frequencies of the annotated loci were associated with environmental gradients (Fig. 5a), indicating adaptation signals. Wind and soil gradients were significantly associated with the allele frequencies of ATPD, NF-Y19, cob, rpoB, and ABCG34 (Fig. 5b; Table S6). Some outliers, such as ADH1 orthologs, were associated with annual mean temperature and precipitation in spring (Fig. 5a). Other genes, including AT2G40435, the NET1D ortholog, and the HSP70 ortholog, were also involved in precipitation- and temperature-associated adaptation (Fig. 5b; Table S6).
3.4 Analyses of ecological niche and morphometric distinctiveness
Consistent with the genetic data, niche analysis of the first two axes of the environmental PCA and both overlap indices showed that the eastern and western clusters had a higher degree of niche overlap with each other than with HC (Fig. S7). HC was characterized by niches with higher mean annual temperature and lower precipitation in spring compared to those of the other clusters (Fig. S7). The equivalency test, background test, and ANOVA demonstrated that the three clusters occupied significantly different ecological niches (Fig. S7; Fig. S8).
The PCA and RDA of phenotypic traits were also congruent with the genetic analyses and showed that HC was morphologically and ecologically differentiated from the other clusters (Fig. 6a-b; Fig. S9; Fig. S10), which was classified as Q. longinux var. kuoi. Although the eastern and southern clusters had a high degree of overlap according to the first two PC axes, PERMANOVA indicated significant differentiation between the two clusters (Table S9). Partial RDA using all phenotypic traits as responses demonstrated that pure geography contributed more variation than pure environment, and a large intersection between their interactions was found (Fig. 6c-d; Table S10). The GLMs revealed that leaf traits were significantly associated with environment, but the directions of the relationships (i.e., positive or negative) differed depending on the trait and environmental variable (Fig. S11). For example, leaf thickness was positively correlated with annual temperature (r = 0.22, p < 0.01), whereas leaf length was negatively correlated with annual precipitation (r = −0.22, p < 0.01).
3.5 Genetic offset and prediction of the response to future climate change
The high AUC value (average AUC = 0.81) suggested a good model fit for the predicted distribution of Q. longinux (Table S11; Fig. S12). As the predictions under different RCPs were highly correlated (Table S12), we inferred RONA based on the average values from both predictions. Substantial variations in the RONA estimates between populations and the top three representative environmental variables were observed (Fig. 7a-c; Table S10). The RONA estimates were larger in regions with greater differences between current and future climates (Fig. 7a-c). Whereas the eastern and western populations had relatively low RONA values (< 0.2) for the three variables, the northern populations were predicted to suffer from severe winter rainfall (precipitation in October, Fig. 7a) in the future and had much higher RONA values (> 0.6).
The GF model constructed with five climatic variables suggested that precipitation in winter was the most influential variable with the highest weighted R2 (Fig. S13). The focal species exhibited strong spatial patterns, indicating adaptation to local climate conditions (Fig. 7d). Consistent with the results of the RONA analysis, the GF model estimated highly correlated results between RCPs with similar genetic offset patterns (Fig. 7e-f). Both the RONA analysis and the GF model indicated that the northern populations were the most vulnerable to future climate change (Fig. 7e-f).
Although the predicted patterns of local, forward, and reverse offsets varied throughout the range of Q. longinux, these offsets were consistently predicted to be highest in the northernmost and southeastern populations (Fig. 8a-b). More migration events were predicted for the northern populations, with longer distances to minimize forward offsets in the RCP8.5 model than in the RCP2.6 model (Fig. 8c-f).
4 Discussion
We analyzed genomic data and phenotypic traits to explore the genetic architecture of several environment-associated adaptations of Q. longinux, a dominant evergreen forest tree species on the subtropical island of Taiwan in East Asia. We identified several SNPs with strong effects on adaptation to environmental factors, including some factors that have rarely been discussed in GEA studies (e.g., soil- and wind-related factors). We found that leaf traits were influenced by the interaction of demographic and environmental factors. Moreover, we determined that the populations in northern and southeastern Taiwan are the most vulnerable to future climate change. Finally, we identified populations with unique genetic and phenotypic characteristics in southern Taiwan. These populations are potential targets for conservation efforts in forest management to preserve unique and adaptive genetic resources.
4.1 Distinct genetic separation of southern populations from eastern and western populations in Taiwan
Three putative genetic clusters were classified using PCA and StrAuto. The eastern and western clusters were mainly separated from each other by the CMR and were mixed in northern and southern Taiwan. Similar patterns of east–west divergence have been observed for other plants in Taiwan [98, 99], implying that the mountain ridges and rugged topography act as profound barriers to gene flow and contribute to the divergence of species occupying low-to-middle elevations in Taiwan. The third cluster, HC, was limited to the Henchung Peninsula in southern Taiwan, which has been identified as the main glacial refuge in Taiwan for other Fagaceae species, such as Q. glauca and Castanopsis carlesii [100, 101]. The EEMS analysis revealed higher genetic diversity in HC and genetically admixed populations in southern and northern Taiwan. The populations of HC were genetically differentiated from other populations, as suggested by higher pairwise FST and distinct trajectory of changes in Ne in the past 1,000 years. Taken together, our results reveal significant genetic differentiation between the eastern cluster, western cluster, and HC and suggest that HC, morphologically classified as Q. longinux var. kuoi, in the southernmost part of Taiwan may have diverged earlier than the other two clusters.
A distinct group of Q. longinux populations, group HC, with unique leaf and fruit traits in tropical marine climates in southern Taiwan was previously classified as Q. longinux var. kuoi. This variety has no whitish epicuticular wax coating on abaxial surfaces [48]. The classification of this variety is consistent with the genetic and morphological separation of the southernmost population from other populations of Q. longinux that we observed in the present study. Niche analyses also demonstrated that the habitats harboring Q. longinux var. kuoi were significantly different from those of the other populations, with no ecological overlaps, indicating that Q. longinux var. kuoi may face different environmental pressures. Moreover, we found evidence of adaptive divergence between Q. longinux var. kuoi and other populations. For example, the low spring precipitation and high clay content in habitats in southeastern Taiwan may act as strong environmental stresses that initiate genetic and phenotypic adaptation (e.g., drought resistance) in response to local conditions. In plant species, hostile environmental conditions in edge populations prompt local adaptation processes [102]. The substantial divergence, relatively high genetic diversity, and high offsets of the populations in southern Taiwan indicate that Q. longinux var. kuoi is a conservation unit that should be prioritized for protection as a source of adaptive genetic variations related to high temperature and drought resistance under climate change.
4.2 Topography as a key driver of genetic differentiation in Quercus longinux
The IBR model caused by the topographical barrier was the most influential factor restricting gene flow and driving the overall genetic differentiation of Q. longinux. The conductance layers of topographical barriers and the EEMS analysis also revealed significantly higher resistance in regions at high elevations. The diverse topography and mountainous areas of Taiwan have led to genetic diversification of various plants and animals [103, 104]. Elevation has been recognized as the main driver of genetic differentiation in other Quercus species distributed in temperate regions [105]. We found that topographical resistance in alignment with mountain ridges was the most profound factor contributing to genetic differentiation on Taiwan Island. Fagaceae typically disperse their seeds through rodents or birds, while their pollens are transported by wind [106-109]. Rodents are the primary carriers and predators of several Quercus species in temperate ecosystems [110, 111] and subtropical forests [112]. Studies suggest that genetic clusters of rodents found at mid-elevations in Taiwan, where most Fagaceae species grow, are separated by mountain ridges and rugged topography [113-115]. Therefore, the dispersal distance of Quercus seeds in Taiwan may be limited by the dispersal abilities of sympatric rodents. Monteiro, Veiga [116] also suggested that elevational barriers are a common obstacle to gene flow for animals in tropical and subtropical regions. Animals in these areas face difficulties in crossing high elevations due to their inability to acclimate to sudden changes in temperature [117].
4.3 Environmental heterogeneity drives adaptive genetic divergence and phenotypic variation
IBE was the most strongly supported model based on putative adaptive loci, whereas IBR mainly drove genetic differentiation based on all SNPs. PCA also revealed less genetic admixture when genetic differentiation was assessed by GEA outliers compared to neutral SNPs, indicating that the three genetic clusters (i.e., eastern, western, and HC) were exposed to different environmental pressures and had undergone adaptive divergence. However, the partial RDA revealed a large intersection of explained variation (45%) shared by environment, geography, and colonization history, suggesting that environmental variation highly covaried with other confounding factors. Similar covariations have been observed for other evergreen species. For example, the south-to-north postglacial expansion of the red spruce along the Appalachian Mountains created high collinearity between genetic structure, climate gradients, and geographic distributions [118]. Similarly, evergreen subtropical trees in Taiwan underwent south–north expansions after the Last Glacial Maximum and may have developed adaptations to latitudinal gradients of temperature and precipitation, resulting in confounding relationships between geography, genetic structure, and adaptation. Consequently, it was challenging to attribute and disentangle the genetic variation explained by each category of predictors, leading to a non-significant impact of pure climate variables.
Leaf shape is affected by various environmental factors, such as precipitation and temperature, which maximizes photosynthetic efficiency and adaptation to harsh environments [75, 119-121]. We found that several leaf traits of Q. longinux were influenced by environmental factors (Fig. S11; Table S10). The negative correlation between leaf length and annual precipitation contradicts previous findings [51, 122]. However, we observed a positive correlation between annual precipitation and wind speed in winter and a negative correlation between annual precipitation and solar radiation in summer (Fig. S14). Strong wind and weaker solar radiation may counteract the effects of increased annual precipitation on leaf length. Similar confounding effects of climatic factors on leaf growth and elongation were observed in Fagus sylvatica [76]. We demonstrated that the interaction of environment and geographic relationships mainly contributed to the explained variation in leaf traits. Demographic history provided only a limited contribution to leaf variation, suggesting that phenotypic plasticity or local adaptation contributed by local climate surpasses the impact of demographic history on leaf traits.
Temperature and precipitation are known drivers of genetic adaptation in Fagaceae [34]. This study expands the environmental factors considered compared to previous GEA studies and demonstrates that significant selection is initiated by multiple environmental factors in an endemic Quercus species. First, we found that wind speed in cold seasons influenced leaf traits and adaptive genetic variation (Table S7; Fig S11). Wind intensity has been shown to reduce plant growth and height and increase stem thickness [123, 124]. Wind also influences transpiration rates and stomatal conductance, indirectly affecting photosynthetic efficiency and water requirements [125, 126]. Some genes correlated with wind speed were also significantly associated with precipitation-related factors (e.g., ATPD, NF-YCP), implying that elevated evaporation rates caused by strong wind may result in drought-like stress, which plants may respond to through similar genetic pathways. Second, we determined that soil-related variables contributed 60% of the variation in adaptive divergence, indicating critical roles of these variables in local adaptation of Q. longinux. The characteristics of soil particle sizes represent the potential water content and salt stress in local soils. In general, soil particle size is negatively correlated with water availability, implying potential abiotic pressure from water deprivation during dry seasons [36, 127]. Consistent with this conclusion, we observed a correlation between genes involved in the response to drought and grid size. In summary, our findings provide a new perspective for future GEA studies by indicating that some environmental variables with important but rarely tested physiological impacts can be used to unravel the intricate mechanisms by which plants respond and adapt to heterogeneous environments.
4.4 Adaptive genes underlying local adaptation in Quercus longinux
We identified several genes involved in adaptation to the local environment (Table S6). We also identified two significantly enriched functional pathways, oxidative phosphorylation and photosynthesis, which have been implicated in plant adaptation and physiological responses to abiotic stress [128, 129]. Oxidative phosphorylation helps regulate reactive oxygen species (ROS) generation by plant mitochondria under abiotic stresses. The efficiency and regulation of photosynthesis strengthen the plant and sustain its growth and development under stressful or unfavorable conditions [130]. These results suggest that the identified adaptive SNPs underlie the response to different abiotic stresses. We also observed associations between the allele frequencies of the annotated genes and environmental gradients. For example, loci of ADH1 are strongly associated with annual mean temperature and precipitation in spring. ADH1 is responsive to multiple abiotic stimuli, including low temperature, hypoxia, flooding, salt, and dehydration [131-133]. In legumes, ADH1 is a target of miRNA regulation under water‐deficit to coordinate ROS levels [134]. Under stressful cold situations, ADH1 plays a crucial role in maintaining the stability of membrane structure to enhance cold resistance in plants [128]. Other genes involved in precipitation- and temperature-associated adaptation include orthologs of AT2G40435, which is involved in responses to biotic and abiotic stresses [135, 136]; NET1D, which is expressed in root tissues and mediates uptake in response to stress [137]; and HSP70, which stabilizes and refolds heat-inactivated proteins to protect cells from heat damage [138].
Although we used several procedures to evaluate the putative adaptive SNPs identified in this study, the interpretations of the constructed GEA relationships and identified outliers are subject to limitations. First, we scanned the outliers using relatively few loci (approximately 2,000 SNPs), which may bias the identified loci and elevate the false-positive rate in GEA approaches [139]. Second, a large proportion of the variation was explained by several confounding factors simultaneously, making it challenging to detect true adaptive signals from spatial autocorrelation and colonization history [118]. Controlling for genetic structures may reduce the false-positive rate but cause many or most adaptive loci to be missed [118, 139]. To minimize the false-positive rate, we used two other FST-based outlier analyses in addition to LFMM. This approach strongly reduces the false-positive rate but may bias the scanning toward strong selective sweeps by excluding several small-effect loci that may be equally important in shaping adaptive variation.
4.5 Assessing genetic vulnerability and climate adaptation in Quercus longinux
The projections of genetic offset from GF and RONA revealed that the populations in northern Taiwan might experience the most considerable turnover of genetic composition to cope with future climate change (Fig. 7). Winter precipitation in northern Taiwan is expected to more than double according to both emission models (current: 302 mm, RCP2.6: 682 mm, RCP8.5: 635 mm). The drastic increase in winter precipitation will negatively impact forest productivity [140] because the complex relationship between precipitation and water availability affects plant growth and phenology [140, 141]. Winter precipitation significantly impacts the phenology of oaks, including the onset and duration of flowering, bud bursting, and leaf flushing, and thus may greatly affect the likelihood and extent of gene flow between populations [142]. Considering the long generation time of oaks and the difficulty of juvenile growth in occupied forests, the expected changes in allele frequency to adapt to increased winter precipitation in the northern populations (RONA > 0.6) may not be achievable through standing genetic variation alone.
Under the emission model of intensified global warming (RCP8.5), southward migration over long distances (> 200 km) will increase to minimize forward genetic offset. Although we accounted for migration in estimating genetic offset, the northernmost and southeastern populations consistently showed relatively high local, forward, and reverse offsets. Our results indicate that no current populations across the distribution of Q. longinux are preadapted to future climates in these regions. The northern and southeastern populations are crucial genetic sources for climatic adaptation in other regions and should be prioritized in conservation strategies and protection efforts. Moreover, the rugged topography in mountainous regions may further hamper the movement of populations to higher elevations. Assisted gene flow from other populations preadapted to future climates may help marginal populations mitigate the effects of climate change [143]. For example, southern populations (e.g., HC) may act as potential sources of adaptation to high temperatures and seasonal arid climates for populations at higher altitudes and latitudes where higher future temperatures are predicted. However, the source populations must be selected carefully to ensure genetic compatibility with the sink populations and new environments [143].
Acknowledgments
The authors express their gratitude to Yun-Hsin Lai for her invaluable contributions to sample collection and leaf morphology measurement. Special thanks to Dr. Chih-Kai Yang and Huan-Ching Lin for their assistance in sample collection. The authors also thank Dawn Schmidt for her assistance in English editing.
References
1. Aitken, S.N., et al., Adaptation, migration or extirpation: climate change outcomes for tree populations. Evolutionary applications, 2008. 1(1): p. 95-111.
2. Fitzpatrick, M.C. and S.R. Keller, Ecological genomics meets community‐level modelling of biodiversity: Mapping the genomic landscape of current and future environmental adaptation. Ecology letters, 2015. 18(1): p. 1-16.
3. Rellstab, C., B. Dauphin, and M. Exposito‐Alonso, Prospects and limitations of genomic offset in conservation management. Evolutionary Applications, 2021. 14(5): p. 1202-1212.
4. Seidl, R., et al., Forest disturbances under climate change. Nature climate change, 2017. 7(6): p. 395-402.
5. Kramer, K., I. Leinonen, and D. Loustau, The importance of phenology for the evaluation of impact of climate change on growth of boreal, temperate and Mediterranean forests ecosystems: an overview. International Journal of Biometeorology, 2000. 44(2): p. 67-75.
6. Piao, S., et al., Plant phenology and global climate change: Current progresses and challenges. Global change biology, 2019. 25(6): p. 1922-1940.
7. Flannigan, M.D., B.J. Stocks, and B.M. Wotton, Climate change and forest fires. Science of the total environment, 2000. 262(3): p. 221-229.
8. Kirschbaum, M.U. and A. Fischlin, Climate change impacts on forests. 1996.
9. Lindner, M., et al., Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. Forest ecology and management, 2010. 259(4): p. 698-709.
10. Way, D.A. and R. Oren, Differential responses to changes in growth temperature between trees from different functional groups and biomes: a review and synthesis of data. Tree physiology, 2010. 30(6): p. 669-688.
11. Wang, T., G.A. O'Neill, and S.N. Aitken, Integrating environmental and genetic effects to predict responses of tree populations to climate. Ecological applications, 2010. 20(1): p. 153-163.
12. Pedlar, J.H. and D.W. McKenney, Assessing the anticipated growth response of northern conifer populations to a warming climate. Scientific Reports, 2017. 7(1): p. 1-10.
13. Saxe, H., et al., Tree and forest functioning in response to global warming. New phytologist, 2001. 149(3): p. 369-399.
14. Waring, R.H. and W. Schlesinger, Forest ecosystems. Analysis at multiples scales, 1985. 55.
15. Camarero, J.J., et al., Growth response to climate and drought change along an aridity gradient in the southernmost Pinus nigra relict forests. Annals of forest science, 2013. 70(8): p. 769-780.
16. Lines, E.R., et al., Predictable changes in aboveground allometry of trees along gradients of temperature, aridity and competition. Global Ecology and Biogeography, 2012. 21(10): p. 1017-1028.
17. Waring, R.H., Characteristics of trees predisposed to die. Bioscience, 1987. 37(8): p. 569-574.
18. Ramírez-Valiente, J.A. and J. Cavender-Bares, Evolutionary trade-offs between drought resistance mechanisms across a precipitation gradient in a seasonally dry tropical oak (Quercus oleoides). Tree physiology, 2017. 37(7): p. 889-901.
19. Sork, V.L., Genomic studies of local adaptation in natural plant populations. Journal of Heredity, 2018. 109(1): p. 3-15.
20. Waldvogel, A.-M., et al., Climate change genomics calls for standardized data reporting. Frontiers in Ecology and Evolution, 2020. 8: p. 242.
21. Rellstab, C., et al., A practical guide to environmental association analysis in landscape genomics. Molecular Ecology, 2015. 24(17): p. 4348-4370.
22. Gugger, P.F., et al., Landscape genomics of Quercus lobata reveals genes involved in local climate adaptation at multiple spatial scales. Molecular Ecology, 2021. 30(2): p. 406-423.
23. Cavender‐Bares, J., Diversification, adaptation, and community assembly of the American oaks (Quercus), a model clade for integrating ecology and evolution. New Phytologist, 2019. 221(2): p. 669-692.
24. Gao, J., et al., Combined genotype and phenotype analyses reveal patterns of genomic adaptation to local environments in the subtropical oak Quercus acutissima. Journal of Systematics and Evolution, 2021. 59(3): p. 541-556.
25. Sork, V., et al., Putting the landscape into the genomics of trees: approaches for understanding local adaptation and population responses to changing climate. Tree Genetics & Genomes, 2013. 9(4): p. 901-911.
26. Holderegger, R., et al., Landscape genetics of plants. Trends in plant science, 2010. 15(12): p. 675-683.
27. Richardson, J.L., et al., Navigating the pitfalls and promise of landscape genetics. 2016, Wiley Online Library.
28. Razgour, O., Beyond species distribution modeling: a landscape genetics approach to investigating range shifts under future climate change. Ecological Informatics, 2015. 30: p. 250-256.
29. Razgour, O., et al., Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proceedings of the National Academy of Sciences, 2019. 116(21): p. 10418-10423.
30. Pollegioni, P., et al., Landscape genetics of Persian walnut (Juglans regia L.) across its Asian range. Tree Genetics & Genomes, 2014. 10(4): p. 1027-1043.
31. Mattioni, C., et al., Landscape genetics structure of European sweet chestnut (Castanea sativa Mill): indications for conservation priorities. Tree Genetics & Genomes, 2017. 13(2): p. 39.
32. Petit, R.J., et al., Fagaceae trees as models to integrate ecology, evolution and genomics. New Phytologist, 2013. 197(2): p. 369-371.
33. Kim, B.Y., et al., RADseq data reveal ancient, but not pervasive, introgression between Californian tree and scrub oak species (Quercus sect. Quercus: Fagaceae). Molecular Ecology, 2018. 27(22): p. 4556-4571.
34. Müller, M. and O. Gailing, Abiotic genetic adaptation in the Fagaceae. Plant Biology, 2019. 21(5): p. 783-795.
35. Aguirre‐Liguori, J.A., et al., Divergence with gene flow is driven by local adaptation to temperature and soil phosphorus concentration in teosinte subspecies (Zea mays parviglumis and Zea mays mexicana). Molecular ecology, 2019. 28(11): p. 2814-2830.
36. Rellstab, C., et al., Signatures of local adaptation in candidate genes of oaks (Quercus spp.) with respect to present and future climatic conditions. Molecular Ecology, 2016. 25(23): p. 5907-5924.
37. Macel, M., et al., Climate vs. soil factors in local adaptation of two common plant species. Ecology, 2007. 88(2): p. 424-433.
38. Cavender-Bares, J. and J.A. Ramírez-Valiente, Physiological evidence from common garden experiments for local adaptation and adaptive plasticity to climate in American live oaks (Quercus Section Virentes): implications for conservation under global change, in Oaks physiological ecology. Exploring the functional diversity of genus Quercus L. 2017, Springer. p. 107-135.
39. Smith, D.S., et al., Soil-mediated local adaptation alters seedling survival and performance. Plant and Soil, 2012. 352(1): p. 243-251.
40. Byars, S.G., W. Papst, and A.A. Hoffmann, Local adaptation and cogradient selection in the alpine plant, Poa hiemata, along a narrow altitudinal gradient. Evolution: International Journal of Organic Evolution, 2007. 61(12): p. 2925-2941.
41. Fang, J.-Y., et al., Divergent selection and local adaptation in disjunct populations of an endangered conifer, Keteleeria davidiana var. formosana (Pinaceae). PLoS One, 2013. 8(7): p. e70162.
42. Hsieh, Y., et al., Historical connectivity, contemporary isolation and local adaptation in a widespread but discontinuously distributed species endemic to Taiwan, Rhododendron oldhamii (Ericaceae). Heredity, 2013. 111(2): p. 147-156.
43. Huang, C.-L., et al., Genetic relationships and ecological divergence in Salix species and populations in Taiwan. Tree Genetics & Genomes, 2015. 11(3): p. 1-17.
44. Foster, P., The potential negative impacts of global climate change on tropical montane cloud forests. Earth-Science Reviews, 2001. 55(1-2): p. 73-106.
45. Dirnböck, T., F. Essl, and W. Rabitsch, Disproportional risk for habitat loss of high‐altitude endemic species under climate change. Global Change Biology, 2011. 17(2): p. 990-996.
46. Still, C.J., P.N. Foster, and S.H. Schneider, Simulating the effects of climate change on tropical montane cloud forests. Nature, 1999. 398(6728): p. 608-610.
47. Cazzolla Gatti, R., et al., Accelerating upward treeline shift in the Altai Mountains under last-century climate change. Scientific reports, 2019. 9(1): p. 1-13.
48. Huang, T., Flora of Taiwan, 2nd edn, Vols 1–5. Editorial Committee of the Flora of Taiwan, Taipei, 1994.
49. Zhong, M., et al., Leaf morphology shift of three dominant species along altitudinal gradient in an alpine meadow of the Qinghai-Tibetan Plateau. Polish Journal of Ecology, 2014. 62(4): p. 639-648.
50. Vaca-Sánchez, M.S., et al., Genetic and functional leaf traits variability of Quercus laurina along an oak diversity gradient in Mexico. European Journal of Forest Research, 2021. 140(5): p. 1211-1225.
51. Peppe, D.J., et al., Sensitivity of leaf size and shape to climate: global patterns and paleoclimatic applications. New phytologist, 2011. 190(3): p. 724-739.
52. Doyle, J., DNA protocols for plants, in Molecular Techniques in Taxonomy. 1991, Springer. p. 283-293.
53. Chen, S., et al., fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 2018. 34(17): p. i884-i890.
54. Schmid-Siegert, E., et al., Low number of fixed somatic mutations in a long-lived oak tree. Nature Plants, 2017. 3(12): p. 926-929.
55. Li, H., Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:1303.3997, 2013.
56. Mose, L.E., C.M. Perou, and J.S. Parker, Improved indel detection in DNA and RNA via realignment with ABRA2. Bioinformatics, 2019. 35(17): p. 2966-2973.
57. Li, H., et al., The sequence alignment/map format and SAMtools. bioinformatics, 2009. 25(16): p. 2078-2079.
58. Li, H., A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics, 2011. 27(21): p. 2987-2993.
59. Danecek, P., et al., The variant call format and VCFtools. Bioinformatics, 2011. 27(15): p. 2156-2158.
60. Jombart, T., adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics, 2008. 24(11): p. 1403-1405.
61. Goudet, J., Hierfstat, a package for R to compute and test hierarchical F‐statistics. Molecular ecology notes, 2005. 5(1): p. 184-186.
62. Oksanen, J., et al., Package ‘vegan’. Community ecology package, version, 2013. 2(9): p. 1-295.
63. Paradis, E., J. Claude, and K. Strimmer, APE: analyses of phylogenetics and evolution in R language. Bioinformatics, 2004. 20(2): p. 289-290.
64. Chhatre, V.E. and K.J. Emerson, StrAuto: automation and parallelization of STRUCTURE analysis. BMC bioinformatics, 2017. 18(1): p. 1-5.
65. Kopelman, N.M., et al., Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Molecular ecology resources, 2015. 15(5): p. 1179-1191.
66. Liu, X. and Y.-X. Fu, Stairway Plot 2: demographic history inference with folded SNP frequency spectra. Genome biology, 2020. 21(1): p. 1-9.
67. Kleinschmit, J. Intraspecific variation of growth and adaptive traits in European oak species. in Annales des sciences forestières. 1993. EDP Sciences.
68. Fick, S.E. and R.J. Hijmans, WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International journal of climatology, 2017. 37(12): p. 4302-4315.
69. Trabucco, A. and R.J. Zomer, Global aridity index and potential evapotranspiration (ET0) climate database v2. CGIAR Consort Spat Inf, 2018. 10.
70. Poggio, L., et al., SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. Soil, 2021. 7(1): p. 217-240.
71. Hijmans, R.J., et al., Raster package in R. 2013, Version.
72. R Core Team, R., R: A language and environment for statistical computing. 2013.
73. Di Cola, V., et al., ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography, 2017. 40(6): p. 774-787.
74. Warren, D.L., R.E. Glor, and M. Turelli, ENMTools: a toolbox for comparative studies of environmental niche models. Ecography, 2010. 33(3): p. 607-611.
75. Sun, M., et al., Variations in leaf morphological traits of Quercus guyavifolia (Fagaceae) were mainly influenced by water and ultraviolet irradiation at high elevations on the Qinghai-Tibet Plateau, China. Int. J. Agric. Biol, 2016. 18: p. 266-273.
76. Meier, I.C. and C. Leuschner, Leaf size and leaf area index in Fagus sylvatica forests: competing effects of precipitation, temperature, and nitrogen availability. Ecosystems, 2008. 11(5): p. 655-669.
77. Kremer, A., et al., Leaf morphological differentiation between Quercus robur and Quercus petraea is stable across western European mixed oak stands. Annals of Forest Science, 2002. 59(7): p. 777-787.
78. Abràmoff, M.D., P.J. Magalhães, and S.J. Ram, Image processing with ImageJ. Biophotonics international, 2004. 11(7): p. 36-42.
79. Kassambara, A. and F. Mundt, Package ‘factoextra’. Extract and visualize the results of multivariate data analyses, 2017. 76(2).
80. Ripley, B., et al., Package ‘mass’. Cran r, 2013. 538: p. 113-120.
81. Foll, M., BayeScan v2. 1 user manual. Ecology, 2012. 20(10).
82. Luu, K., E. Bazin, and M.G. Blum, pcadapt: an R package to perform genome scans for selection based on principal component analysis. Molecular ecology resources, 2017. 17(1): p. 67-77.
83. Thissen, D., L. Steinberg, and D. Kuang, Quick and easy implementation of the Benjamini-Hochberg procedure for controlling the false positive rate in multiple comparisons. Journal of educational and behavioral statistics, 2002. 27(1): p. 77-83.
84. Frichot, E. and O. François, LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 2015. 6(8): p. 925-929.
85. Frichot, E., et al., Testing for associations between loci and environmental gradients using latent factor mixed models. Molecular biology and evolution, 2013. 30(7): p. 1687-1699.
86. Pitcher, C., et al.,
Example analysis of biodiversity survey data with R package gradientForest. R vignette. Available at
http://gradientforest. r-forge. r-project. org/biodiversity-survey. pdf [Verified 27 March 2017], 2011.
87. Neph, S., et al., BEDOPS: high-performance genomic feature operations. Bioinformatics, 2012. 28(14): p. 1919-1920.
88. Bu, D., et al., KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic acids research, 2021. 49(W1): p. W317-W325.
89. Upton, G.J., Fisher's exact test. Journal of the Royal Statistical Society: Series A (Statistics in Society), 1992. 155(3): p. 395-402.
90. McRae, B.H. and V.B. Shah, Circuitscape user’s guide. The University of California, Santa Barbara, 2009.
91. Cushman, S.A., et al., Re-evaluating causal modeling with Mantel tests in landscape genetics. Diversity, 2013. 5(1): p. 51-72.
92. Peterman, W.E., et al., A comparison of popular approaches to optimize landscape resistance surfaces. Landscape Ecology, 2019. 34(9): p. 2197-2208.
93. Peterman, W.E., ResistanceGA: An R package for the optimization of resistance surfaces using genetic algorithms. Methods in Ecology and Evolution, 2018. 9(6): p. 1638-1647.
94. Petkova, D., J. Novembre, and M. Stephens, Visualizing spatial population structure with estimated effective migration surfaces. Nature genetics, 2016. 48(1): p. 94-100.
95. Naimi, B. and M.B. Araújo, sdm: a reproducible and extensible R platform for species distribution modelling. Ecography, 2016. 39(4): p. 368-375.
96. Pina‐Martins, F., et al., New insights into adaptation and population structure of cork oak using genotyping by sequencing. Global Change Biology, 2019. 25(1): p. 337-350.
97. Gougherty, A.V., S.R. Keller, and M.C. Fitzpatrick, Maladaptation, migration and extirpation fuel climate change risk in a forest tree species. Nature Climate Change, 2021. 11(2): p. 166-171.
98. Yu, T.-L., H.-D. Lin, and C.-F. Weng, A new phylogeographic pattern of endemic Bufo bankorensis in Taiwan Island is attributed to the genetic variation of populations. PloS one, 2014. 9(5): p. e98029.
99. Huang, S.F., et al., Phylogeography of Trochodendron aralioides (Trochodendraceae) in Taiwan and its adjacent areas. Journal of Biogeography, 2004. 31(8): p. 1251-1259.
100. Shih, F., et al., Partial concordance between nuclear and organelle DNA in revealing the genetic divergence among Quercus glauca (Fagaceae) populations in Taiwan. International Journal of Plant Sciences, 2006. 167(4): p. 863-872.
101. Cheng, Y.P., S.Y. Hwang, and T.P. Lin, Potential refugia in Taiwan revealed by the phylogeographical study of Castanopsis carlesii Hayata (Fagaceae). Molecular Ecology, 2005. 14(7): p. 2075-2085.
102. Yang, Y.-Z., et al., Parallel adaptation prompted core-periphery divergence of Ammopiptanthus mongolicus. Frontiers in Plant Science, 2022. 13: p. 956374.
103. Lin, H.-D., Y.-R. Chen, and S.-M. Lin, Strict consistency between genetic and topographic landscapes of the brown tree frog (Buergeria robusta) in Taiwan. Molecular Phylogenetics and Evolution, 2012. 62(1): p. 251-262.
104. Li, Y.-S., et al., Testing the effect of mountain ranges as a physical barrier to current gene flow and environmentally dependent adaptive divergence in Cunninghamia konishii (Cupressaceae). Frontiers in genetics, 2019. 10: p. 742.
105. Gharehaghaji, M., et al., Effects of landscape features on gene flow of valley oaks (Quercus lobata). Plant Ecology, 2017. 218(4): p. 487-499.
106. Liepelt, S., R. Bialozyt, and B. Ziegenhagen, Wind-dispersed pollen mediates postglacial gene flow among refugia. Proceedings of the National Academy of Sciences, 2002. 99(22): p. 14590-14594.
107. Craft, K.J. and M.V. Ashley, Landscape genetic structure of bur oak (Quercus macrocarpa) savannas in Illinois. Forest Ecology and Management, 2007. 239(1-3): p. 13-20.
108. Roth, J.K. and S.B. Vander Wall, Primary and secondary seed dispersal of bush chinquapin (Fagaceae) by scatterhoarding rodents. Ecology, 2005. 86(9): p. 2428-2439.
109. Wang, J., et al., Effects of mast seeding and rodent abundance on seed predation and dispersal of Quercus aliena (Fagaceae) in Qinling Mountains, Central China. Plant Ecology, 2017. 218(7): p. 855-865.
110. Yu, F., et al., Seed predation patterns favor the regeneration of dominant species in forest gaps compared with the understory in an oak-pine mixed forest. Acta theriologica, 2014. 59(4): p. 495-502.
111. Gómez, J.M., D. Garcıa, and R. Zamora, Impact of vertebrate acorn-and seedling-predators on a Mediterranean Quercus pyrenaica forest. Forest ecology and management, 2003. 180(1-3): p. 125-134.
112. Xiao, Z., Z. Zhang, and C.J. Krebs, Long‐term seed survival and dispersal dynamics in a rodent‐dispersed tree: testing the predator satiation hypothesis and the predator dispersal hypothesis. Journal of Ecology, 2013. 101(5): p. 1256-1264.
113. Yu, H.-T., Patterns of diversification and genetic population structure of small mammals in Taiwan. Biological Journal of the Linnean Society, 1995. 55(1): p. 69-89.
114. Yuan, S.L., L.K. Lin, and T. Oshida, Phylogeography of the mole‐shrew (Anourosorex yamashinai) in Taiwan: implications of interglacial refugia in a high‐elevation small mammal. Molecular Ecology, 2006. 15(8): p. 2119-2130.
115. Hsu, F.-H., F.-J. Lin, and Y.-S. Lin, Phylogeographic structure of the Formosan wood mouse, Apodemus semotus Thomas. ZOOLOGICAL STUDIES-TAIPEI-, 2001. 40(2): p. 91-102.
116. Monteiro, W.P., et al., Everything you always wanted to know about gene flow in tropical landscapes (but were afraid to ask). PeerJ, 2019. 7: p. e6446.
117. Janzen, D.H., Why mountain passes are higher in the tropics. The American Naturalist, 1967. 101(919): p. 233-249.
118. Capblancq, T., et al., From common gardens to candidate genes: exploring local adaptation to climate in red spruce. New Phytologist, 2023. 237(5): p. 1590-1605.
119. Joel, G., G. Aplet, and P.M. Vitousek, Leaf morphology along environmental gradients in Hawaiian Metrosideros polymorpha. Biotropica, 1994: p. 17-22.
120. Liu, W., L. Zheng, and D. Qi, Variation in leaf traits at different altitudes reflects the adaptive strategy of plants to environmental changes. Ecology and Evolution, 2020. 10(15): p. 8166-8175.
121. Hovenden, M.J. and J.K. Vander Schoor, Nature vs nurture in the leaf morphology of Southern beech, Nothofagus cunninghamii (Nothofagaceae). New Phytologist, 2004. 161(2): p. 585-594.
122. Wiemann, M.C., et al., Estimation of temperature and precipitation from morphological characters of dicotyledonous leaves. American Journal of Botany, 1998. 85(12): p. 1796-1802.
123. Biddington, N.L., The effects of mechanically-induced stress in plants—a review. Plant growth regulation, 1986. 4(2): p. 103-123.
124. Smith, V. and A. Ennos, The effects of air flow and stem flexure on the mechanical and hydraulic properties of the stems of sunflowers Helianthus annuus L. Journal of Experimental Botany, 2003. 54(383): p. 845-849.
125. Anten, N.P., et al., Wind and mechanical stimuli differentially affect leaf traits in Plantago major. New Phytologist, 2010. 188(2): p. 554-564.
126. Onoda, Y. and N.P. Anten, Challenges to understand plant responses to wind. Plant Signaling & Behavior, 2011. 6(7): p. 1057-1059.
127. Zhang, X., et al., Relationship between soil water content and soil particle size on typical slopes of the Loess Plateau during a drought year. Science of the total environment, 2019. 648: p. 943-954.
128. Zsigmond, L., et al., Arabidopsis PPR40 connects abiotic stress responses to mitochondrial electron transport. Plant Physiology, 2008. 146(4): p. 1721-1737.
129. Vashisth, D., et al., Transcriptome changes induced by abiotic stresses in Artemisia annua. Scientific Reports, 2018. 8(1): p. 3423.
130. Muhammad, I., et al., Mechanisms regulating the dynamics of photosynthesis under abiotic stresses. Frontiers in plant science, 2021. 11: p. 615942.
131. Christie, P.J., M. Hahn, and V. Walbot, Low-temperature accumulation of alcohol dehydrogenase-1 mRNA and protein activity in maize and rice seedlings. Plant Physiology, 1991. 95(3): p. 699-706.
132. de Bruxelles, G.L., et al., Abscisic acid induces the alcohol dehydrogenase gene in Arabidopsis. Plant Physiology, 1996. 111(2): p. 381-391.
133. Komatsu, S., et al., Characterization of a novel flooding stress-responsive alcohol dehydrogenase expressed in soybean roots. Plant Molecular Biology, 2011. 77: p. 309-322.
134. De la Rosa, C., A.A. Covarrubias, and J.L. Reyes, A dicistronic precursor encoding miR398 and the legume‐specific miR2119 coregulates CSD1 and ADH1 mRNAs in response to water deficit. Plant, cell & environment, 2019. 42(1): p. 133-144.
135. Liu, J.X., et al., Salt stress responses in Arabidopsis utilize a signal transduction pathway related to endoplasmic reticulum stress signaling. The Plant Journal, 2007. 51(5): p. 897-909.
136. Lin, F., et al., Molecular response to the pathogen Phytophthora sojae among ten soybean near isogenic lines revealed by comparative transcriptomics. BMC genomics, 2014. 15(1): p. 1-13.
137. Hawkins, T.J., et al., The evolution of the actin binding NET superfamily. Frontiers in Plant Science, 2014. 5: p. 254.
138. De Maio, A., Heat shock proteins: facts, thoughts, and dreams. Shock, 1999. 11(1): p. 1-12.
139. Ahrens, C.W., et al., The search for loci under selection: trends, biases and progress. Molecular ecology, 2018. 27(6): p. 1342-1356.
140. Zeppel, M., J.V. Wilks, and J.D. Lewis, Impacts of extreme precipitation and seasonal changes in precipitation on plants. Biogeosciences, 2014. 11(11): p. 3083-3093.
141. Lipton, D., et al., Ecosystems, ecosystem services, and biodiversity. 2018, US Global Change Research Program.
142. Armstrong-Herniman, W. and S. Greenwood, The role of winter precipitation as a climatic driver of the spring phenology of five California Quercus species (Fagaceae). Madroño, 2021. 68(4): p. 450-460.
143. Aitken, S.N. and M.C. Whitlock, Assisted gene flow to facilitate local adaptation to climate change. Annual review of ecology, evolution, and systematics, 2013. 44: p. 367-388.
Data Accessibility
The following primary datasets associated with this paper have been deposited on Figshare under the DOI: 10.6084/m9.figshare.24276946:
1. VCF files encompassing all individuals subjected to analysis in this study.
2. Leaf morphological traits of each sample measured in this study.
3. Binary habitat layer estimated with ENM.
4. Environmental data corresponding to each sampling site.
5. Topographical and climate resistance layers simulated in this study.
Funding
This work was supported by the National Science and Technology Council under project ID NSTC 112-2621-B-003-001-MY3.
Competing Interests
The authors have no relevant financial or non-financial interests to disclose.
Author contribution
P.C.L. designed and supervised the project. P.W.S. collected samples and performed laboratory experiments and statistical analyses. P.W.S. and P.C.L. drafted the manuscript. J.T.C. and M.X.L. assisted with data interpretation, provided valuable feedback, and revised the manuscript. All authors have read and approved the final version of the manuscript.