Skip to main content
Download PDF

Genomic insights into local adaptation and vulnerability of Quercus longinux to climate change

BMC Plant Biology volume 24, Article number: 279 (2024) Cite this article

Abstract

Background

Climate change is expected to alter the factors that drive changes in adaptive variation. This is especially true for species with long life spans and limited dispersal capabilities. Rapid climate changes may disrupt the migration of beneficial genetic variations, making it challenging for them to keep up with changing environments. Understanding adaptive genetic variations in tree species is crucial for conservation and effective forest management. Our study used landscape genomic analyses and phenotypic traits from a thorough sampling across the entire range of Quercus longinux, an oak species native to Taiwan, to investigate the signals of adaptation within this species.

Results

Using ecological data, phenotypic traits, and 1,933 single-nucleotide polymorphisms (SNPs) from 205 individuals, we classified three genetic groups, which were also phenotypically and ecologically divergent. Thirty-five genes related to drought and freeze resistance displayed signatures of natural selection. The adaptive variation was driven by diverse environmental pressures such as low spring precipitation, low annual temperature, and soil grid sizes. Using linear-regression-based methods, we identified isolation by environment (IBE) as the optimal model for adaptive SNPs. Redundancy analysis (RDA) further revealed a substantial joint influence of demography, geology, and environments, suggesting a covariation between environmental gradients and colonization history. Lastly, we utilized adaptive signals to estimate the genetic offset for each individual under diverse climate change scenarios. The required genetic changes and migration distance are larger in severe climates. Our prediction also reveals potential threats to edge populations in northern and southeastern Taiwan due to escalating temperatures and precipitation reallocation.

Conclusions

We demonstrate the intricate influence of ecological heterogeneity on genetic and phenotypic adaptation of an oak species. The adaptation is also driven by some rarely studied environmental factors, including wind speed and soil features. Furthermore, the genetic offset analysis predicted that the edge populations of Q. longinux in lower elevations might face higher risks of local extinctions under climate change.

Peer Review reports

Introduction

Trees profoundly influence global carbon cycles and ecosystem stabilization [1,2,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,5,6,7]. However, the impacts of climate change on focal species vary [8,9,10,11,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,15,16,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,23,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,42,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,45,46,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,50,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.

The ecological function of oaks in subtropical forests is essential, and it will be further important to evaluate the adaptation and their vulnerability when climate change is believed to alter their distribution and deteriorate their survival. The first step to conservation is revealing their population structure and adaptation pattern and identifying potential genetic sources and vulnerable populations. Therefore, 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. Second, we aimed to identify environmental features that affect spatial genetic variation and adaptive genetic divergence. Third, we aimed to use adaptive genetic variants to evaluate the vulnerability of populations to climate change.

Materials and methods

Sampling, sequencing, read mapping, and variant calling

All samples analyzed in this study were collected and identified by the authors. We collected 26 populations spanning all known distributional regions of Q. longinux in Taiwan (Fig. S1a; Table S1). Three populations in Southern Taiwan (SK, LZ, and GS) were morphologically identified as Q. longinux var. kuoi; three populations (TP, NS, JSY) were collected near the locations with documented Q. longinux var. lativiolaciifolia [48]; and the rest populations were identified as Q. longinux var. longinux. 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. Voucher specimens of this study have been deposited at National Taiwan Normal University Herbarium (TNU) under deposition number TNU057201–TNU057214.

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, at the Technology Commons in College of Life Science (National Taiwan University, Taiwan).

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].

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].

To provide an in-depth investigation of the evolution history between Q. longinux var. kuoi and the rest of the populations, we simulated and compared four scenarios (Fig. S1b) with different patterns of migration using fastsimcoal2 [68]. All models were simulated and optimized using the composite likelihoods calculated with SFSs. Because the computation of fastsimcoal2 is demanding, observed SFSs were projected to five individuals per simulated population. A mutation rate of 1.01 × 10− 8 per bp per year was assumed for estimating all parameters with a generation time of 50 years [67]. To select the best model, we first used the Akaike information criterion (AIC), which represents the difference between the observed and expected SFSs, to rank the optimized models. Second, the parametric-bootstrap approach was used to calculate 100 likelihoods and the 95% CI for the parameters estimated by the best models and to select the best-fitting models according to the lowest bootstrapping AIC.

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 [69] and two aridity-associated layers [70]. Eight physicochemical soil characteristics downloaded from SoilGrids 2.0 were used as the soil variables [71]. 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 [72]. 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 [73]. The boundaries of Taiwan and nearby islands were downloaded from GADM database (http://www.gadm.org/). All maps in this study were depicted by the authors using R 4.2.3 [73].

Assessment of ecological and morphogenetic divergence

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 [74]. 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 [75].

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) [76,77,78]. 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 [79]. 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 [78]. 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 [80]. 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 [81].

Detection of local adaptation and underlying gene functions

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 [82] and pcadapt [83] 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 [84] 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 [85, 86], 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 [85]. LFMM was performed in the R package LEA [85] 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 [84] 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 [73].

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 [87]. 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 [87]. The importance of each group of predictors was evaluated by the percentage of weighted R2.

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 [88]. The closest genes were further annotated to the Arabidopsis thaliana genome using KOBAS-I [89] 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 [90], and FDR was corrected with the Benjamini and Hochberg procedure [84]. FDR < 5% indicated significantly enriched GO terms or pathways.

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 [91] 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) [92] 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) [93] to compare each predictor based on the AIC values. Models with ΔAIC > 2 were considered to have different model fits. MLPE was implemented with mlpe_rga in the R package ResistanceGA [94]. In addition, we performed an estimated effective migration surfaces (EEMS) analysis [95] 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.

Investigating the impacts of future climate change on genetic adaptedness

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 [96]. 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.

Genetic offset assessment

Because the accessible prediction layers are limited, we only applied five temperature and precipitation variables (i.e., BIO1, BIO3, BIO12, PREC4, and PREC10) selected from the obtained factors 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 [97]. We retained and discussed 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 [87]. 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 [98], 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 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.

Results

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. From the result of StrAuto, individuals in each cluster were genetically admixed with other clusters (Fig. 1b). PCA and the neighbor-joining (NJ) tree also suggested that the boundaries of the three clusters were not clearly defined (Fig. 1c; Fig. S3a-b). FST was higher between genetic clusters than within groups, and divergence was highest between HC and the eastern cluster (Fig. S4c). Similar results were obtained by PERMANOVA (Table S5). The results of the estimated effective migration surface (EEMS) analysis also suggested the connectivity and geographic distribution of genetic diversity were mainly separated by CMR (Fig. 1d-e). 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 Ky. The fastsimcoal2 analysis suggested no apparent gene flow between group HC and the eastern and western clusters since their divergence of 19.6 Kya (Table S6). Recent bottleneck events were also revealed at 11.5 Kya from the best models for all simulated populations.

Fig. 1
figure 1

Spatial genetic structure and diversity of Q. longinux in Taiwan. (a) Sampling populations annotated with ancestral components inferred by StrAuto according to K = 3. (b) Results of StrAuto from K = 2 to 3. The height of each colored segment represents the possible ancestry of each individual derived from inferred ancestors. (c) Results of PCA with different colors reflecting different groups. (d) Results of the EEMS analysis of regions significantly deviating from isolation by distance (IBD). The blue, white, and orange colors illustrate regions with high dispersal rates (dispersal corridor), IBD, and low dispersal rates (dispersal barrier), respectively. (e) Results of the effective diversity (i.e., the modeled dissimilarity between pairs of individuals from the same location) estimated by the EEMS analysis. Blue and orange represent regions of high and low genetic diversity, respectively. The black dots indicate the sampling locations

Full size image

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 S7) that were widely distributed across the genome and did not cluster in specific regions. FST estimate from adaptive SNPs was higher and more significant than 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 S8). 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. 1d-e; Table S8). The conductance layers generated by the CIRCUITSCAPE algorithm also supported topographical resistance across regions, with greater barriers to dispersal at higher elevations (Fig. S4e-f).

Fig. 2
figure 2

Results of the genetic scans. The Venn diagram shows the overlap of the loci identified by each method

Full size image
Fig. 3
figure 3

Landscape genetics of Q. longinux and factors influencing genetic divergence. (a) Isolation by distance (IBD) analysis and three isolation by environment (IBE) analyses based on adaptive sites (green dots and upper lines) or all SNP sites (red dots and lower lines). (b) Results of reciprocal causal modeling using all SNP sites or (c) adaptive loci. The y-axis represents focal models, and the x-axis represents alternative models. The colors of the heatmap indicate differences in Mantel’s R between the focal and alternative models. Fully supported focal models are marked with asterisks

Full size image

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 portion of all SNPs when controlling other factors (Fig. S5; Table S9). Forward selection also identified population structure as the most significant variable affecting the genetic variation of adaptive SNPs (Table S10). Environmental variables contributed to 80% of the explained variation in adaptive loci but 60% 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 S9). 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 S11; Fig. S6), although geography explained a higher proportion of variation after summing the total weighted R2 (demography: 32%; geography: 47%).

Functional annotation of adapted loci

Several Q. longinux loci involved in climate and soil adaptation were associated with environmental variables (Table S7). 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. 4a), indicating adaptation signals. Wind and soil gradients were significantly associated with the allele frequencies of ATPD, NF-Y19, cob, rpoB, and ABCG34 (Fig. 4b; Table S7). Some outliers, such as ADH1 orthologs, were associated with annual mean temperature and precipitation in spring (Fig. 4a). Other genes, including AT2G40435, the NET1D ortholog, and the HSP70 ortholog, were also involved in precipitation- and temperature-associated adaptation (Fig. 4b; Table S7).

Fig. 4
figure 4

Results of environment-dependent outlier detection. (a) Associations between allele frequencies and environmental variables. The black dots represent the average allele frequencies of all adaptive loci correlated with WIND10, PREC04, CLYPPT, and BIO01 for each population. (b) Manhattan plots of SNP sites associated with WIND10, PREC04, CLYPPT, and BIO01. The orange dots represent significant loci at p < 0.05. Selected candidate genes are labeled in the plots at their respective genomic positions

Full size image

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. 5a-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 S12). In contrast, Q. longinux var. lativiolaciifolia was genetically and morphologically admixed with Q. longinux var. longinux, and no apparent structure could be assigned between these two varieties (Fig. 1a-b). 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. 5c-d; Table S13). 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).

Fig. 5
figure 5

Variations in leaf traits and results of partial RDA. (a) Results of PCA colored by different genetic groups. (b) Partial RDA partitioning sources of phenotypic variation into different environmental factors. Environmental factors are depicted with red arrows. (c) Partial RDA of environmental factors, geography, and demography. (d) Partial RDA partitioning explained variation into climatic, soil, and topographical variables. The values are the explained variation. Values < 0 are not shown

Full size image

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 S14; Fig. S12). As the predictions under different RCPs were highly correlated (Table S15), 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. 6a-c; Table S15). The RONA estimates were larger in regions with greater differences between current and future climates (Fig. 6a-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. 6a) in the future and had much higher RONA values (> 0.6).

Fig. 6
figure 6

Prediction of genetic offset under future climate change based on (a, b, c) RONA and (d, e, f) GF methods. (a), (b), and (c) reflect the RONA values estimated for PREC10, BIO12, and BIO03. The color gradients represent the average differences in climatic variables between current conditions and 2070 climate change scenarios. (d) RGB map of the first three PC axes based on the GF prediction, which depicts the genetic turnover of adaptive loci. (e) and (f) illustrate the genetic offset throughout the range of Q. longinux under the RCP2.6 (c) and RCP8.5 (d) scenarios in 2070. The colors of the cells represent the values of genetic offset estimated based on the GF procedure

Full size image

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. 6d). Consistent with the results of the RONA analysis, the GF model estimated highly correlated results between RCPs with similar genetic offset patterns (Fig. 6e-f). Both the RONA analysis and the GF model indicated that the northern populations were the most vulnerable to future climate change (Fig. 6e-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. 7a-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. 7c-f).

Fig. 7
figure 7

Map of the GDM-predicted genetic offset across the distribution of Q. longinux under the RCP2.6 (a, c, e) and RCP8.5 (b, d, f) climate change scenarios. (a) and (b) show the RGB maps of local (red), forward (green), and reverse (blue) offsets under RCP2.6 (a) and RCP8.5 (b) in 2070. Brighter cells represent relatively high values along each of the three axes, whereas darker cells indicate relatively low values. (c) and (d) depict the direction of forward offsets among the distribution range of Q. longinux under RCP2.6 (c) and RCP8.5 (d) in 2070. (e) and (f) represent the estimated distances between source cells with the lowest forward offset to sink cells under RCP2.6 (e) and RCP8.5 (f) in 2070

Full size image

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.

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 [99, 100], 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 [101, 102]. 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. The optimal demographic scenario revealed HC diverged from the rest of the populations dating back to the Last Glacial Maximum (LGM; ∼20Kya) when most of the low-elevation oaks were believed to retreat to the refuge in southern Taiwan [101, 102]. Also, implying colonization resulting from Pleistocene climate oscillations could prompted the diversification and local adaptation of plants in Taiwan Island. 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 have diverged from the two clusters without apparent gene flow since their divergence.

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]. 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 [103]. Distinct environmental pressure, along with a lack of migration between HC and the rest of the Q. longinux population, could speed up the fixation of alleles in a relatively smaller group, further facilitating the genetic and phenotypic adaptation [25, 104]. 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.

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 is highly covaried with other confounding factors. Similar covariations have been observed in the south-to-north postglacial expansion of the red spruce along the Appalachian Mountains which created high collinearity between genetic structure, climate gradients, and geographic distributions [105]. Similarly, evergreen subtropical trees in Taiwan underwent south–north expansions after the LGM 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 maximize photosynthetic efficiency and adaptation to harsh environments [76, 106,107,108]. We found that several leaf traits of Q. longinux were influenced by environmental factors (Fig. S11; Table S13) and a negative correlation between leaf length and annual precipitation contradicts previous findings [51, 109]. 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 [77]. 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 S10; Fig. S11). Wind intensity has been shown to reduce plant growth and height and increase stem thickness [110, 111]. Wind also influences transpiration rates and stomatal conductance, indirectly affecting photosynthetic efficiency and water requirements [112, 113]. 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, 114]. 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 related to plant adaptation.

We identified two significantly enriched functional pathways, oxidative phosphorylation and photosynthesis, which have been implicated in plant adaptation and physiological responses to abiotic stress [115, 116]. Oxidative phosphorylation helps regulate reactive oxygen species (ROS) generation by plant mitochondria under abiotic stresses [115, 116]. The efficiency and regulation of photosynthesis strengthen the plant and sustain its growth and development under stressful or unfavorable conditions [117]. These results suggest that the identified adaptive SNPs underlie the response to different abiotic stresses. Several genes demonstrating significant associations with environments also have potential functions for adaptation. 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 [118,119,120]. In legumes, ADH1 is a target of miRNA regulation under water-deficit to coordinate ROS levels [121]. Under stressful cold situations, ADH1 plays a crucial role in maintaining the stability of membrane structure to enhance cold resistance in plants [115]. Other genes involved in precipitation- and temperature-associated adaptation include orthologs of AT2G40435, which is involved in responses to biotic and abiotic stresses [122, 123]; NET1D, which is expressed in root tissues and mediates uptake in response to stress [124]; and HSP70, which stabilizes and refolds heat-inactivated proteins to protect cells from heat damage [125].

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. 6). 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 [126] because the complex relationship between precipitation and water availability affects plant growth and phenology [126, 127]. 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 [128]. 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 [129]. 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 [129].

Data availability

Sequencing data have been deposited on NCBI BioProject PRJNA1049506 under GenBank accession number SAMN38699773–SAMN38700060. The Variant Call Format (VCF) file along with additional files used for the analyses have been deposited on Figshare under DOI 10.6084/m9.figshare.24276946.

References

  1. Aitken SN, et al. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol Appl. 2008;1(1):95–111.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Fitzpatrick MC, Keller SR. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol Lett. 2015;18(1):1–16.

    Article  PubMed  Google Scholar 

  3. Rellstab C, Dauphin B, Exposito-Alonso M. Prospects and limitations of genomic offset in conservation management. Evol Appl. 2021;14(5):1202–12.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Seidl R, et al. Forest disturbances under climate change. Nat Clim Change. 2017;7(6):395–402.

    Article  Google Scholar 

  5. Kramer K, Leinonen I, Loustau D. The importance of phenology for the evaluation of impact of climate change on growth of boreal, temperate and Mediterranean forests ecosystems: an overview. Int J Biometeorol. 2000;44(2):67–75.

    Article  CAS  PubMed  Google Scholar 

  6. Piao S, et al. Plant phenology and global climate change: current progresses and challenges. Glob Change Biol. 2019;25(6):1922–40.

    Article  Google Scholar 

  7. Flannigan MD, Stocks BJ, Wotton BM. Climate change and forest fires. Sci Total Environ. 2000;262(3):221–9.

    Article  CAS  PubMed  Google Scholar 

  8. Kirschbaum MU, Fischlin A. Climate change impacts on forests 1996.

  9. Lindner M, et al. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For Ecol Manag. 2010;259(4):698–709.

    Article  Google Scholar 

  10. Way DA, Oren R. Differential responses to changes in growth temperature between trees from different functional groups and biomes: a review and synthesis of data. Tree Physiol. 2010;30(6):669–88.

    Article  PubMed  Google Scholar 

  11. Wang T, O’Neill GA, Aitken SN. Integrating environmental and genetic effects to predict responses of tree populations to climate. Ecol Appl. 2010;20(1):153–63.

    Article  CAS  PubMed  Google Scholar 

  12. Pedlar JH, McKenney DW. Assessing the anticipated growth response of northern conifer populations to a warming climate. Sci Rep. 2017;7(1):1–10.

    Article  Google Scholar 

  13. Saxe H, et al. Tree and forest functioning in response to global warming. New Phytol. 2001;149(3):369–99.

    Article  CAS  PubMed  Google Scholar 

  14. Waring RH, Schlesinger W. Forest ecosystems Analysis at multiples scales, 1985. 55.

  15. Camarero JJ, et al. Growth response to climate and drought change along an aridity gradient in the southernmost Pinus nigra relict forests. Ann for Sci. 2013;70(8):769–80.

    Article  Google Scholar 

  16. Lines ER, et al. Predictable changes in aboveground allometry of trees along gradients of temperature, aridity and competition. Glob Ecol Biogeogr. 2012;21(10):1017–28.

    Article  Google Scholar 

  17. Waring RH. Characteristics of trees predisposed to die. Bioscience. 1987;37(8):569–74.

    Article  Google Scholar 

  18. Ramírez-Valiente JA, Cavender-Bares J. Evolutionary trade-offs between drought resistance mechanisms across a precipitation gradient in a seasonally dry tropical oak (Quercus oleoides). Tree Physiol. 2017;37(7):889–901.

    Article  PubMed  Google Scholar 

  19. Sork VL. Genomic studies of local adaptation in natural plant populations. J Hered. 2018;109(1):3–15.

    Article  Google Scholar 

  20. Waldvogel A-M, et al. Climate change genomics calls for standardized data reporting. Front Ecol Evol. 2020;8:242.

    Article  Google Scholar 

  21. Rellstab C, et al. A practical guide to environmental association analysis in landscape genomics. Mol Ecol. 2015;24(17):4348–70.

    Article  PubMed  Google Scholar 

  22. Gugger PF, et al. Landscape genomics of Quercus lobata reveals genes involved in local climate adaptation at multiple spatial scales. Mol Ecol. 2021;30(2):406–23.

    Article  CAS  PubMed  Google Scholar 

  23. Cavender-Bares J. Diversification, adaptation, and community assembly of the American oaks (Quercus), a model clade for integrating ecology and evolution. New Phytol. 2019;221(2):669–92.

    Article  PubMed  Google Scholar 

  24. Gao J, et al. Combined genotype and phenotype analyses reveal patterns of genomic adaptation to local environments in the subtropical oak Quercus acutissima. J Syst Evol. 2021;59(3):541–56.

    Article  Google Scholar 

  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 Genet Genomes. 2013;9(4):901–11.

    Article  Google Scholar 

  26. Holderegger R, et al. Landscape genetics of plants. Trends Plant Sci. 2010;15(12):675–83.

    Article  CAS  PubMed  Google Scholar 

  27. Richardson JL, et al. Navigating the pitfalls and promise of landscape genetics. Wiley Online Library; 2016.

  28. Razgour O. Beyond species distribution modeling: a landscape genetics approach to investigating range shifts under future climate change. Ecol Inf. 2015;30:250–6.

    Article  Google Scholar 

  29. Razgour O, et al. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc Natl Acad Sci. 2019;116(21):10418–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Pollegioni P, et al. Landscape genetics of persian walnut (Juglans regia L.) across its Asian range. Tree Genet Genomes. 2014;10(4):1027–43.

    Article  Google Scholar 

  31. Mattioni C, et al. Landscape genetics structure of European sweet chestnut (Castanea sativa Mill): indications for conservation priorities. Tree Genet Genomes. 2017;13(2):39.

    Article  Google Scholar 

  32. Petit RJ, et al. Fagaceae trees as models to integrate ecology, evolution and genomics. New Phytol. 2013;197(2):369–71.

    Article  PubMed  Google Scholar 

  33. Kim BY, et al. RADseq data reveal ancient, but not pervasive, introgression between Californian tree and scrub oak species (Quercus sect. Quercus: Fagaceae). Mol Ecol. 2018;27(22):4556–71.

    Article  CAS  PubMed  Google Scholar 

  34. Müller M, Gailing O. Abiotic genetic adaptation in the Fagaceae. Plant Biol. 2019;21(5):783–95.

    Article  PubMed  Google Scholar 

  35. Aguirre-Liguori JA, 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). Mol Ecol. 2019;28(11):2814–30.

    Article  CAS  PubMed  Google Scholar 

  36. Rellstab C, et al. Signatures of local adaptation in candidate genes of oaks (Quercus spp.) with respect to present and future climatic conditions. Mol Ecol. 2016;25(23):5907–24.

    Article  PubMed  Google Scholar 

  37. Macel M, et al. Climate vs. soil factors in local adaptation of two common plant species. Ecology. 2007;88(2):424–33.

    Article  PubMed  Google Scholar 

  38. Cavender-Bares J, Ramírez-Valiente JA. 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. Oaks physiological ecology. Exploring the functional diversity of genus Quercus L. Springer; 2017. pp. 107–35.

  39. Smith DS, et al. Soil-mediated local adaptation alters seedling survival and performance. Plant Soil. 2012;352(1):243–51.

    Article  CAS  Google Scholar 

  40. Byars SG, Papst W, Hoffmann AA. Local adaptation and cogradient selection in the alpine plant, Poa Hiemata, along a narrow altitudinal gradient. Evolution: Int J Org Evol. 2007;61(12):2925–41.

    Article  Google Scholar 

  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):e70162.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  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):147–56.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Huang C-L, et al. Genetic relationships and ecological divergence in Salix species and populations in Taiwan. Tree Genet Genomes. 2015;11(3):1–17.

    Article  CAS  Google Scholar 

  44. Foster P. The potential negative impacts of global climate change on tropical montane cloud forests. Earth Sci Rev. 2001;55(1–2):73–106.

    Article  Google Scholar 

  45. Dirnböck T, Essl F, Rabitsch W. Disproportional risk for habitat loss of high-altitude endemic species under climate change. Glob Change Biol. 2011;17(2):990–6.

    Article  Google Scholar 

  46. Still CJ, Foster PN, Schneider SH. Simulating the effects of climate change on tropical montane cloud forests. Nature. 1999;398(6728):608–10.

    Article  CAS  Google Scholar 

  47. Cazzolla Gatti R, et al. Accelerating upward treeline shift in the Altai Mountains under last-century climate change. Sci Rep. 2019;9(1):1–13.

    Article  CAS  Google Scholar 

  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. Pol J Ecol. 2014;62(4):639–48.

    Google Scholar 

  50. Vaca-Sánchez MS, et al. Genetic and functional leaf traits variability of Quercus Laurina along an oak diversity gradient in Mexico. Eur J for Res. 2021;140(5):1211–25.

    Article  Google Scholar 

  51. Peppe DJ, et al. Sensitivity of leaf size and shape to climate: global patterns and paleoclimatic applications. New Phytol. 2011;190(3):724–39.

    Article  PubMed  Google Scholar 

  52. Doyle J. DNA protocols for plants, in Molecular techniques in Taxonomy. Springer; 1991. pp. 283–93.

  53. Chen S, et al. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Schmid-Siegert E, et al. Low number of fixed somatic mutations in a long-lived oak tree. Nat Plants. 2017;3(12):926–9.

    Article  PubMed  Google Scholar 

  55. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM arXiv preprint arXiv:1303.3997, 2013.

  56. Mose LE, Perou CM, Parker JS. Improved indel detection in DNA and RNA via realignment with ABRA2. Bioinformatics. 2019;35(17):2966–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Li H, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  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):2987–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Danecek P, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24(11):1403–5.

    Article  CAS  PubMed  Google Scholar 

  61. Goudet J. Hierfstat, a package for R to compute and test hierarchical F-statistics. Mol Ecol Notes. 2005;5(1):184–6.

    Article  Google Scholar 

  62. Oksanen J, et al. Package ‘vegan’ community ecology package. Version. 2013;2(9):1–295.

    Google Scholar 

  63. Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20(2):289–90.

    Article  CAS  PubMed  Google Scholar 

  64. Chhatre VE, Emerson KJ. StrAuto: automation and parallelization of STRUCTURE analysis. BMC Bioinformatics. 2017;18(1):1–5.

    Article  Google Scholar 

  65. Kopelman NM, et al. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015;15(5):1179–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Liu X, Fu Y-X. Stairway plot 2: demographic history inference with folded SNP frequency spectra. Genome Biol. 2020;21(1):1–9.

    Google Scholar 

  67. Kleinschmit J. Intraspecific variation of growth and adaptive traits in European oak species. In Annales des sciences forestières. EDP Sciences; 1993.

  68. Excoffier L, et al. fastsimcoal2: demographic inference under complex evolutionary scenarios. Bioinformatics. 2021;37(24):4882–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37(12):4302–15.

    Article  Google Scholar 

  70. Trabucco A, Zomer RJ. Global aridity index and potential evapotranspiration (ET0) climate database v2. CGIAR Consort Spat Inf; 2018. p. 10.

  71. Poggio L, et al. SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. Soil. 2021;7(1):217–40.

    Article  CAS  Google Scholar 

  72. Hijmans RJ et al. Raster package in R. 2013, Version.

  73. Core Team R. R., R: A language and environment for statistical computing 2013.

  74. Di Cola V, et al. Ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography. 2017;40(6):774–87.

    Article  Google Scholar 

  75. Warren DL, Glor RE, Turelli M. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography. 2010;33(3):607–11.

    Article  Google Scholar 

  76. 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:266–73.

    Article  Google Scholar 

  77. Meier IC, Leuschner C. Leaf size and leaf area index in Fagus sylvatica forests: competing effects of precipitation, temperature, and nitrogen availability. Ecosystems. 2008;11(5):655–69.

    Article  CAS  Google Scholar 

  78. Kremer A, et al. Leaf morphological differentiation between Quercus robur and Quercus petraea is stable across western European mixed oak stands. Ann for Sci. 2002;59(7):777–87.

    Article  Google Scholar 

  79. Abràmoff MD, Magalhães PJ, Ram SJ. Image processing with ImageJ. Biophotonics Int. 2004;11(7):36–42.

    Google Scholar 

  80. Kassambara A, Mundt F. Package ‘factoextra’. Extract Visualize Results Multivar data Analyses, 2017. 76(2).

  81. Ripley B, et al. Package ‘mass’. Cran r. 2013;538:113–20.

    Google Scholar 

  82. Foll M. BayeScan v2. 1 user manual. Ecology, 2012. 20(10).

  83. Luu K, Bazin E, Blum MG. Pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol Ecol Resour. 2017;17(1):67–77.

    Article  CAS  PubMed  Google Scholar 

  84. Thissen D, Steinberg L, Kuang D. Quick and easy implementation of the Benjamini-Hochberg procedure for controlling the false positive rate in multiple comparisons. J Educational Behav Stat. 2002;27(1):77–83.

    Article  Google Scholar 

  85. Frichot E, François O. LEA: an R package for landscape and ecological association studies. Methods Ecol Evol. 2015;6(8):925–9.

    Article  Google Scholar 

  86. Frichot E, et al. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol Biol Evol. 2013;30(7):1687–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. 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.

  88. Neph S, et al. BEDOPS: high-performance genomic feature operations. Bioinformatics. 2012;28(14):1919–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Bu D, et al. KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Res. 2021;49(W1):W317–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Upton GJ. Fisher’s exact test. J Royal Stat Society: Ser (Statistics Society). 1992;155(3):395–402.

    Article  Google Scholar 

  91. McRae BH, Shah VB. Circuitscape user’s guide. Santa Barbara: The University of California; 2009.

    Google Scholar 

  92. Cushman SA, et al. Re-evaluating causal modeling with mantel tests in landscape genetics. Diversity. 2013;5(1):51–72.

    Article  Google Scholar 

  93. Peterman WE, et al. A comparison of popular approaches to optimize landscape resistance surfaces. Landscape Ecol. 2019;34(9):2197–208.

    Article  Google Scholar 

  94. Peterman WE. ResistanceGA: an R package for the optimization of resistance surfaces using genetic algorithms. Methods Ecol Evol. 2018;9(6):1638–47.

    Article  Google Scholar 

  95. Petkova D, Novembre J, Stephens M. Visualizing spatial population structure with estimated effective migration surfaces. Nat Genet. 2016;48(1):94–100.

    Article  CAS  PubMed  Google Scholar 

  96. Naimi B, Araújo MB. Sdm: a reproducible and extensible R platform for species distribution modelling. Ecography. 2016;39(4):368–75.

    Article  Google Scholar 

  97. Pina-Martins F, et al. New insights into adaptation and population structure of cork oak using genotyping by sequencing. Glob Change Biol. 2019;25(1):337–50.

    Article  Google Scholar 

  98. Gougherty AV, Keller SR, Fitzpatrick MC. Maladaptation, migration and extirpation fuel climate change risk in a forest tree species. Nat Clim Change. 2021;11(2):166–71.

    Article  Google Scholar 

  99. Yu T-L, Lin H-D, Weng C-F. A new phylogeographic pattern of endemic Bufo bankorensis in Taiwan Island is attributed to the genetic variation of populations. PLoS ONE. 2014;9(5):e98029.

    Article  PubMed  PubMed Central  Google Scholar 

  100. Huang SF, et al. Phylogeography of Trochodendron aralioides (Trochodendraceae) in Taiwan and its adjacent areas. J Biogeogr. 2004;31(8):1251–9.

    Article  Google Scholar 

  101. Shih F, et al. Partial concordance between nuclear and organelle DNA in revealing the genetic divergence among Quercus glauca (Fagaceae) populations in Taiwan. Int J Plant Sci. 2006;167(4):863–72.

    Article  CAS  Google Scholar 

  102. Cheng YP, Hwang SY, Lin TP. Potential refugia in Taiwan revealed by the phylogeographical study of Castanopsis Carlesii Hayata (Fagaceae). Mol Ecol. 2005;14(7):2075–85.

    Article  CAS  PubMed  Google Scholar 

  103. Yang Y-Z, et al. Parallel adaptation prompted core-periphery divergence of Ammopiptanthus mongolicus. Front Plant Sci. 2022;13:956374.

    Article  PubMed  PubMed Central  Google Scholar 

  104. Manel S, et al. Landscape genetics: combining landscape ecology and population genetics. Trends Ecol Evol. 2003;18(4):189–97.

    Article  Google Scholar 

  105. Capblancq T, et al. From common gardens to candidate genes: exploring local adaptation to climate in red spruce. New Phytol. 2023;237(5):1590–605.

    Article  PubMed  Google Scholar 

  106. Joel G, Aplet G, Vitousek PM. Leaf morphology along environmental gradients in Hawaiian Metrosideros polymorpha Biotropica, 1994: p. 17–22.

  107. Liu W, Zheng L, Qi D. Variation in leaf traits at different altitudes reflects the adaptive strategy of plants to environmental changes. Ecol Evol. 2020;10(15):8166–75.

    Article  PubMed  PubMed Central  Google Scholar 

  108. Hovenden MJ, Vander JK, Schoor. Nature vs nurture in the leaf morphology of Southern Beech, Nothofagus Cunninghamii (Nothofagaceae). New Phytol. 2004;161(2):585–94.

    Article  PubMed  Google Scholar 

  109. Wiemann MC, et al. Estimation of temperature and precipitation from morphological characters of dicotyledonous leaves. Am J Bot. 1998;85(12):1796–802.

    Article  CAS  PubMed  Google Scholar 

  110. Biddington NL. The effects of mechanically-induced stress in plants—a review. Plant Growth Regul. 1986;4(2):103–23.

    Article  CAS  Google Scholar 

  111. Smith V, Ennos A. The effects of air flow and stem flexure on the mechanical and hydraulic properties of the stems of sunflowers Helianthus annuus L. J Exp Bot. 2003;54(383):845–9.

    Article  CAS  PubMed  Google Scholar 

  112. Anten NP, et al. Wind and mechanical stimuli differentially affect leaf traits in Plantago major. New Phytol. 2010;188(2):554–64.

    Article  PubMed  Google Scholar 

  113. Onoda Y, Anten NP. Challenges to understand plant responses to wind. Plant Signal Behav. 2011;6(7):1057–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  114. Zhang X, et al. Relationship between soil water content and soil particle size on typical slopes of the Loess Plateau during a drought year. Sci Total Environ. 2019;648:943–54.

    Article  CAS  PubMed  Google Scholar 

  115. Zsigmond L, et al. Arabidopsis PPR40 connects abiotic stress responses to mitochondrial electron transport. Plant Physiol. 2008;146(4):1721–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  116. Vashisth D, et al. Transcriptome changes induced by abiotic stresses in Artemisia annua. Sci Rep. 2018;8(1):3423.

    Article  PubMed  PubMed Central  Google Scholar 

  117. Muhammad I, et al. Mechanisms regulating the dynamics of photosynthesis under abiotic stresses. Front Plant Sci. 2021;11:615942.

    Article  PubMed  PubMed Central  Google Scholar 

  118. Christie PJ, Hahn M, Walbot V. Low-temperature accumulation of alcohol dehydrogenase-1 mRNA and protein activity in maize and rice seedlings. Plant Physiol. 1991;95(3):699–706.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  119. de Bruxelles GL, et al. Abscisic acid induces the alcohol dehydrogenase gene in Arabidopsis. Plant Physiol. 1996;111(2):381–91.

    Article  PubMed  PubMed Central  Google Scholar 

  120. Komatsu S, et al. Characterization of a novel flooding stress-responsive alcohol dehydrogenase expressed in soybean roots. Plant Mol Biol. 2011;77:309–22.

    Article  CAS  PubMed  Google Scholar 

  121. De la Rosa C, Covarrubias AA, Reyes JL. A dicistronic precursor encoding miR398 and the legume-specific miR2119 coregulates CSD1 and ADH1 mRNAs in response to water deficit. Plant Cell Environ. 2019;42(1):133–44.

    Article  PubMed  Google Scholar 

  122. Liu JX, et al. Salt stress responses in Arabidopsis utilize a signal transduction pathway related to endoplasmic reticulum stress signaling. Plant J. 2007;51(5):897–909.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  123. 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):1–13.

    Article  Google Scholar 

  124. Hawkins TJ, et al. The evolution of the actin binding NET superfamily. Front Plant Sci. 2014;5:254.

    Article  PubMed  PubMed Central  Google Scholar 

  125. De Maio A. Heat shock proteins: facts, thoughts, and dreams. Shock. 1999;11(1):1–12.

    Article  PubMed  Google Scholar 

  126. Zeppel M, Wilks JV, Lewis JD. Impacts of extreme precipitation and seasonal changes in precipitation on plants. Biogeosciences. 2014;11(11):3083–93.

    Article  Google Scholar 

  127. Lipton D, et al. Ecosystems, ecosystem services, and biodiversity. US Global Change Research Program; 2018.

  128. Armstrong-Herniman W, Greenwood S. The role of winter precipitation as a climatic driver of the spring phenology of five California Quercus species (Fagaceae). Madroño. 2021;68(4):450–60.

    Article  Google Scholar 

  129. Aitken SN, Whitlock MC. Assisted gene flow to facilitate local adaptation to climate change. Annu Rev Ecol Evol Syst. 2013;44:367–88.

    Article  Google Scholar 

Download references

Acknowledgements

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.

Funding

This work was supported by the National Science and Technology Council under project ID NSTC 112-2621-B-003-001-MY3. This article was also subsidized by National Taiwan Normal University (NTNU).

Author information

Authors and Affiliations

  1. School of Life Science, National Taiwan Normal University, No. 88 Ting-Chow Rd., Sec. 4, Taipei, 116, Taiwan

    Pei-Wei Sun, Jui-Tse Chang, Min-Xin Luo & Pei-Chun Liao

Authors
  1. Pei-Wei Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jui-Tse Chang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Min-Xin Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Pei-Chun Liao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

P.C.L. designed and supervised the project. P.W.S. collected samples, identified plant materials and specimens, 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.

Corresponding author

Correspondence to Pei-Chun Liao.

Ethics declarations

Ethics approval and consent to participate

All samples analyzed in this study were collected by the authors. P.W.S. identified the plant materials used in this study. The locations of collections are not privately owned or protected areas, and are not involved with endangered or protected species. No permits were required for this study.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, PW., Chang, JT., Luo, MX. et al. Genomic insights into local adaptation and vulnerability of Quercus longinux to climate change. BMC Plant Biol 24, 279 (2024). https://doi.org/10.1186/s12870-024-04942-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-024-04942-8

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us