1. Introduction

Global environmental change has profoundly influenced biodiversity on earth (Shivanna, 2022). Under global climate changes, unraveling the genetic basis of adaptation to changing environments has attracted widespread attentions from evolutionary geneticists for decades as the capacity to adapt to environmental change affects the long-term survival of species (Lasky et al., 2022). Local adaptation occurs when individuals have higher fitness in their native habitats than elsewhere as various environmental factors exert different selection pressures across a species range, and is commonly found within plants and animals (Kawecki and Ebert, 2004). Despite its prevalence, the genetic underpinnings of local adaptation are still poorly understood (Lascoux et al., 2016). Studying local adaptation is of great importance for understanding how species adapt to environmental change, designing and predicting the effects of environmental change on species distributions (Aguirre-Liguori et al., 2021; Sang et al., 2022).

Traditional approaches, such as reciprocal-transplant experiment and common garden experiment have been used to uncover local adaptation by measuring genetic signals associated with adaptive phenotypes (Capblancq et al., 2023; Wadgymar et al., 2022). However, these methods pose challenges and are usually unfeasible for non-model plants due to long life history and experimental intractability (Sang et al., 2022). Landscape genomics, which integrates landscape ecology and population genomics, offers a relatively new approach to determine the effects of environmental factors shaping adaptive genetic variations and find variants driving local adaptation through genotype-environment associations (GEAs) (Zhang et al., 2024). With the rapid advancement of next-generation sequencing technologies, landscape genomics has greatly enhanced our understanding of local adaptation by focusing on identifying outliers or candidate genes associated with environmental factors, and evaluating current and prospective genotype-environment mismatched scenarios to ascertain the extent of genomic plasticity with regards to climate change (Gugger et al., 2021; Heraghty et al., 2023; Wang et al., 2023; Zhang et al., 2023a).

Previous studies about adaptation to climate are heavily biased towards annual crops and tree species with social and economic benefits (Sang et al., 2022; Yuan et al., 2023). However, as an important component of subtropical and tropical forest ecosystem, perennial lianas (woody vines) were largely ignored for a long time, although they comprise up to 10% of woody individuals and species in subtropical forests and usually have a large impact on restoration and succession in forests (Cai and Song, 2000; Gentry, 1992; Schnitzer and Bongers, 2002). Investigation of adaptive evolution in such taxon will provide additional perspectives into how lianas species response to climate change. Here, we addressed Actinidia eriantha Benth (2n = 2x = 58, Actinidia, Actinidiaceae), a perennial deciduous liana with ecological and economic values (Sun et al., 2015). As an endemic species, A. eriantha is widely distributed in subtropical China with an elevational distribution ranging from 100 to 1500 m (Testolin et al., 2016), and wild individuals show abundant variations in fruit shapes and flower colors (Liao et al., 2021). However, it faces threats from predatory exploitation and climate change now, with its natural habitats rapidly decreasing (Zhang et al., 2023b). Previous phylogeographical study showed that A. eriantha originated from the eastern China with 'oceanic' adaptation and expanded to western China during the glacial periods (Guo et al., 2022). Prior GEA analysis based on restriction site-associated DNA sequencing (RAD-seq) and limited populations demonstrated that local climate factors such as temperature and precipitation imposed strong selection pressures and leaded to adaptation to climatic change for A. eriantha (Zhang et al., 2023b). However, it should be noticed that RAD-seq typically represented only a fraction of the genome and yielded low-density SNPs, which might lose crucial genetic footprints involved in local adaptation (Feng and Du, 2022; Hoffmann et al., 2021). In addition, given that genetic offsets are sensitive to the number of populations especially for species with high genetic differentiation, fewer than 15 populations used to predict maladaptation may compromise the accuracy of genomic offset projections (Aguirre-Liguori et al., 2023). Moreover, other environmental variables such as solar radiations, soil components and metal ion, which could also play a critical role in local adaptation, have not been analyzed in previous study (Zhang et al., 2023b), which hinder a comprehensively understanding of genetic basis of local adaptation in A. eriantha.

Therefore, this study aims to illustrate the genetic basis of adaptation to diverse environments in Actinidia eriantha through whole-genome resequencing and landscape genomics. To this end, we sampled more independent natural populations along environment gradients of A. eriantha and performed GEA analysis to identify genetic variations (single nucleotide polymorphisms and insertions/deletions) associated with environmental variations. Our results highlighted that environment factors driven the evolutionary processes of local adaptation and plants adapted to heterogeneous environments through changes in allele frequencies along with environments in native habitats, which provided an additional insight on how liana species adapted to spatially heterogeneous environments.

2. Materials and methods 2.1. Plant samples and environment factors collection

In 2021, we collected young fresh leaves of 311 individuals originated from 25 populations covering the distribution range of Actinidia eriantha (Fig. 1a). The fresh leaves were dried with silica gel for DNA extraction and the detailed sampling information was shown in Table S1 and Table S2.

A total of 117 environment factors spanning 40 years (1970–2010) of 25 provenance sites (resolution 2.5 arcmin) were extracted with R package raster (https://CRAN.R-project.org/package=raster). Detailed explanations of environmental factors and raw environment data were provided in Table S3 and Table S4. Among these raw data, six UV-B variables were extracted from glUV database (Beckmann et al., 2014) and eight soil factors were obtained from Soil Sub Center, National Earth System Science Data Center, National Science & Technology Infrastructure (Liu et al., 2022). The remained climate data were downloaded from the World Climate database (WorldClim2) (Fick and Hijmans, 2017).

2.2. DNA extraction, library construction and resequencing

Genomic DNA were extracted via a modified cetyltrimethylammonium bromide (CTAB) method (Doyle JJT and Doyle JL, 1990). After DNA extraction, 1 μg genomic DNA was randomly fragmented to an average size of 200–400 bp to construct the libraries. The qualified libraries were sequenced by DNBSEQ platform with a target depth of 10X. The raw sequencing data was processed using SOAPnuke (Chen et al., 2018) with the parameters "-n 0.01 -l 20 –q 0.3 –adaMR 0.25 –ada_trim –polyX 50". Finally, about 7 Gb clean data per sample was obtained.

2.3. Variation calling and annotation

All clean paired-end data were aligned to the reference genome of Actinidia eriantha (Yao et al., 2022) via the algorithm "bwa-mem" implemented in Burrows-Wheeler Aligner (BWA, v.0.7.17-r1198-dirty) (Li and Durbin, 2009) with the additional parameter '−M'. Aligned reads were sorted and converted into BAM format from SAM format using SAMtools v.1.13 (Danecek et al., 2021). Picard v.2.24.0 (https://github.com/broadinstitute/picard) was used to remove the duplicated reads due to PCR amplification with default parameters. The sequencing depths and coverage rates were calculated by SAMtools. The translated BAM files were used for variant calling with Genome Analysis Toolkit v.4.2.0.0 (McKenna et al., 2010).

The SNP detection was performed based on the following pipeline. Firstly, a genomic variant call format (GVCF) file in ERC mode for each sample was obtained using "HaplotypeCaller" program of GATK (DePristo et al., 2011). Then all GVCF files merged into a single file with the tool "CombineGVCFs". Finally, primary SNPs and Indels were detected at the population level by "GenotypeVCFs". To obtain high-quality variants, a hard filtered standard was formulated using "VariantFiltration". For SNPs, the parameters were set as follows: QD < 2.0, MQ < 40.0, FS > 60.0, SOR > 3.0, MQRankSum < −12.5, ReadPosRankSum < −8.0. Indels (shorter than 50 bp) were selected with the recommended parameters: QD < 2.0, FS > 200.0, SOR > 10.0, MQRankSum < −12.5, ReadPosRankSum < −8.0. To improve the accuracy of variations, all variations were further filtered using VCFtools v.0.1.16 (https://vcftools.github.io/man_latest.html) with the following parameters '–maf 0.05 –max-missing 0.9 –min-alleles 2 –max-alleles 2 –minQ 30'. The remained sites were considered as the core subset for subsequent population analysis and annotated by SnpEff v.5.0 (Cingolani et al., 2014).

2.4. Population structure and genetic diversity analysis

iTools implement in Reseqtools (He et al., 2013) was employed to delineate the genomic positions of 4DTv (fourfold synonymous transversion) sites. 4DTvs were isolated with BCFtools (https://github.com/samtools/bcftools) based on their positional coordinates and used for building a phylogenetic tree and population structure analysis. For phylogenetic relationship inference, we firstly searched for the best substitution model by ModelFinder algorithm (Kalyaanamoorthy et al., 2017) among 283 DNA models and the best-fit model "GTR + F + R6 + ASC" was chosen according to Bayesian Information Criterion (BIC). The ultimate phylogenetic tree was built with maximum-likelihood (ML) method using IQ-TREE v.1.6.12 (Nguyen et al., 2015) with 1000 bootstrap replicates. Principal component analysis (PCA) was performed via PLINK v.1.90b6.24 (Purcell et al., 2007). Additionally, ADMIXTURE v.1.3.0 (Alexander et al., 2009) was used to conduct the population structure analysis. We performed ten independent runs for each K value from 2 to 10.

We used VCFtools to calculate the nucleotide diversities (π), Tajima's D value and pairwise F_ST with a 100-Kb window size and a step size of 10 Kb for groups and populations, respectively. To elucidate the relationships between spatial distances and genetic distances (F_ST/(1-F_ST)) (isolation by distance, IBD), and between environment distances and genetic divergences (isolation by environment, IBE), a simple Mantel test and a partial Mantel test were performed by vegan v.2.6–4 packages with 9999 permutations, respectively (Dixon, 2003).

2.5. Inference of population demographic history

The effective population size (Ne) was inferred using SMC++ (Terhorst et al., 2017). To improve reliability, genome regions for each chromosome where reads were not uniquely mapped were masked according to SNPable pipeline (http://lh3lh3.users.sourceforge.net/snpable.shtml). Ten representative and deeply sequenced individuals were selected from each cluster (Table S2). The estimate mutation rate was set to 1 × 10⁻⁸ per bp per generation and the generation time was defined as 7 years (Liu et al., 2017).

2.6. Estimation of the importance about environments driving genetic variations

Gradient Forest (GF) was employed to assess the most important environmental factor influencing the genetic variations, utilizing nonparametric and machine learning regression trees (Ellis et al., 2012). We ran GF using 20, 000 random SNPs selecting from the whole dataset with five replications to reduce the random error. The parameter listed as following: ntree = 500, maxLevel = log2(0.368n)/2 (n is the number of candidate SNP), and correlation threshold was set to 0.5.

To avoid the multicollinearity of predictors, environmental variables with Pearson correlation coefficients |r| < 0.7 were retained for subsequent analyses. Partial redundancy analysis (pRDA) was performed to estimate the influences of independent environmental factors, geography and population genetic structure on genetic variations. In this study, the first three PCs were sought as natural population structures, while latitudes and longitudes of individuals represented the geography. Each of three categories, named geography, genetic structure and environmental factors, was separately used as explanatory variable, with the remaining two variables analyzed as conditioning variables in partial model. Additionally, all variable regard as explanatory variables in full model. All models were carried out with genotypes as response variables.

2.7. Identification of candidate genetic variations associated with environments

We used two approaches, latent factor mixed model (LFMM) and redundancy analysis (RDA), to detect the candidate loci associated with environments using variants without missing data. LFMM introduced population structure as the latent factor to control false positive rate (Frichot et al., 2013). In this study, we adopted the K value inferred by ADMIXTURE. The association value between the frequency of each variant site and each environmental variable was calibrated using genomic inflation factors (GIF). P-value for all sites were transformed into q-values using the R package qvalue (https://github.com/StoreyLab/qvalue), with sites having q-value < 0.01 considered significantly associated with environmental variables. Identified variants were retained if they were correlated with at least one environmental variable.

RDA was considered as a knife with high true-positive and low false-positive rates and widely used in identifying adaptive loci (Capblancq and Forester, 2021). We performed the pRDA with first 3 PC as conditional variants and first three significant RDA axes accounting for most variants were selected to detected candidate loci. Loci with a standard deviation cutoff of 3.5 were considered as outliers associated with environment predictors. Finally, SNPs or Indels identified by two methods were considered as the core adaptive outliers for downstream analysis.

2.8. GO enrichment analysis

GO enrichment analysis were performed with R package Clusterprofiler (Wu et al., 2021). The candidate protein sequences were aligned to the local eggnog-MAPPER database (Cantalapiedra et al., 2021) and annotated with default parameters based on the mmseq2 (Steinegger and Soding, 2017). GO database was constructed using AnnotationForge (https://bioconductor.org/packages/AnnotationForge).

2.9. Genetic mismatches between current and future climates

To predict the potential adaptability and change of Actinidia eriantha populations under future climate accurately, we considered four emission scenarios of the shared socioeconomic pathway (SSP126, SSP245, SSP370 and SSP585) under four different climate models (BCC-CSM2-MR model, ACCESS-CM2 model, CMCC-ESM2 model and GISS-E2-1-G model), respectively. We used nine environmental variables including bio1, bio2, bio3, bio7, bio15, bio18, prec03, prec09, and tmin07 under both current and future conditions and adaptive loci to estimate genomic vulnerability. Besides, all loci were also employed to assess the risk of non-adaptedness (RONA) with Euclidean distances between current and future allele frequencies for each natural population under changing climate scenarios, following the methods outlined by Rellstab et al. (2016).

3. Results 3.1. SNP genotyping, population structure analysis and demographic history

Resequencing of 311 individuals of Actinidia eriantha yielded a total of 2.2 Tb clean data with an average depth of 11.74X and a mean coverage of 92.22% (Table S2). A core subset of 11, 278, 868 high-confidence SNPs and 1, 420, 355 high-quality Indels corresponding to an average density of 18 SNPs and 2.5 Indels per kilobase were identified, respectively (Fig. 2). A large portion of the SNPs and indels were located in non-coding regions, with 50% of SNPs in intergenic regions and 36% of indels in intron regions. Only a small fraction of SNPs (5%) and indels (3%) were found within the exon regions, respectively (Table S5).

As inferring the genetic basis of adaptation needs to consider the potential influence of neutral processes that can shape genetic variations, we first elucidated the population structure and demographic history of Actinidia eriantha. Population structure analyses using the fourfold degenerate synonymous sites (4DTvs, = 188, 329) showed that all individuals were divided into four main subgroups denoted as cluster East, Middle, West and Southeast according to the ML tree corresponding to Wuyi Mountains, Luoxiao Mountains and Nanling Mountains, Xuefeng Mountains, and southeastern coastal areas, which is further supported by PCA analysis, with 21.2%, 12.6% and 5.6% of the variation explained by three principal components axes, respectively (Fig. 1a–c). Cluster East predominantly contained natural populations from Zhejiang, Jiangxi, Fujian and Guangdong provinces, while middle group consisted of five populations including two from western of Jiangxi province, one from north of Guangdong province and two populations in Hunan province. West cluster contained seven populations mainly collected from the western China (Guizhou, Guangxi, and Hunan provinces). Notably, seven marginal populations (WH, XW, YC, JGS, MS, JH and ZY populations) within genetic clusters exhibited significant mixtures. In addition, the Southeast cluster consisted of a single population located in Hua'an country in Fujian province that exhibited phenotypic divergence in hair color and strong genetic differentiation from the other clusters (Figs. 1a, e–g and S1). Multiple genetic bottlenecks were observed in these genetic clusters except the East characterized by a gradual reduction in effective population size, followed by a brief recovery and a steady decline (Fig. S2). Specifically, the East, West, and Middle groups experienced a sharp decrease in effective population size around 0.5–1 million years ago. While between 20, 000 and 70, 000 years ago, the population sizes of three groups decreased except for the East group.

3.2. Genetic diversity and genetic differentiations

The East group had highest genetic diversity, with a coefficient of 3.44 × 10⁻³, followed by the Middle group and the West group, with coefficients of 3.40 × 10⁻³ and 2.49 × 10⁻³, respectively. The Southeast group displayed the lowest genetic diversity, with a coefficient of only 1.72 × 10⁻³. Tajima's D results were consistent with the patterns of genetic diversity. Besides, Southeast group exhibited significant differentiation from other three groups (F_ST > 0.3) (Fig. 1d). It was further supported by the clustering analysis of wild individuals, which revealed a greater genetic distance between the Southeast group and the other three groups (Fig. 1b). At the population level, nucleotide diversity ranged from 1.51 × 10⁻³ (GD) to 3.60 × 10⁻³ (GX) (Fig. S3). Notably, the genetic diversity was significantly correlated with longitude (Fig. S4).

In addition, the Mantel's test revealed a significant correlation between geographic distance and genetic distance (r = 0.62 and P < 0.001). This was further confirmed by a partial Mantel's test (r = 0.58 and P < 0.001), which accounted for environmental distance. Although a weaker correlation between genetic distance and environmental distance was observed (r = 0.32 and P < 0.001), it was not significant when controlling for geographic distance (r = 0.13 and P = 0.11) (Table S6).

3.3. Primary environmental factors and drivers of genetic variations

To identify the primary environmental factors influencing genetic variations, we employed GF to rank the importance of environmental variables. Prec_03 emerged as the most crucial environmental driver of allele frequency variations within Actinidia eriantha populations, followed by solar radiation (Figs. 3a and S5). We retained 23 independent environmental factors including five temperature factors (bio_1, bio_2, bio_3, bio_7, tmin_07), four precipitation variables (bio_15, bio_18, prec_03, prec_09), seven soil indexes (cec1530, cf1530, ph1530, soc1530, tk1530, tn1530, tp1530), three UVB factors (UVB2, UVB3, UVB5), three wind elements (wind_01, wind_04, wind_11) and one solar radiation index (srad_03) to characterize the climatic features of A. eriantha for subsequent analysis (Fig. 3b). We further assessed the proportion of variation explained by environmental factors, geography and genetic structure using pRDA to infer the drivers of patterns of genetic diversity. The results demonstrated that three explanatory variables totally explained 65% of the overall genetic variations across A. eriantha (Table 1). Specifically, environment variables explained 19% of variations when controlling for geography and genetic structure, implying that climate was the primary factor driving the patterns of genetic diversity. However, when controlling for another two factors, genetic structure accounted for only 0.4% of total variations. Similarly, the effect of geography explained 0.14% of genotypic variations controlling for climate and genetic structure. In addition, approximately 34.78% variations remained unexplained by climate, geography and genetic structure.

3.4. Genetic variations associated with environmental adaptation

To reduce the influence of population structure on association analysis, we performed GEA analysis using two high-quality datasets (SNPs and Indels), both with and without HA population. Remarkably, a total of 458, 446 and 82, 470 SNPs associated with at least one environment were identified by LFMM and RDA across all populations, along with 57, 467 and 8708 indels highlighted, respectively (Fig. S6). However, excluding the HA population, we identified a total of 196, 626 SNP sites and 15, 661 Indel sites in the remained populations using the LFMM method. Conversely, the RDA method identified a greater number of SNPs (122, 847) and Indels (11, 159) compared to the analysis included HA population (Fig. S6). In total, we obtained 17, 022 SNPs and 1829 Indels associated with environments identified simultaneously by aforementioned four methods. These loci were associated with 1566 genes including the members of NAC transcription factors (NAC010, NAC102), ERF gene family (ERF3, ERF14 and ERF096, etc.), and the homologs of calcium-dependent protein kinases (CPK5, CPK11, and CPK15) (Table S7). GO analysis suggested these genes were predominantly enriched in processes including positive regulation of long-day photoperiodism and flowering, histone modification, DNA replication (Fig. S7), indicating they might be involved in the adaptation to heterogeneous environments.

Importantly, we pinpointed two missense mutations (Chr26:17099817 and Chr26:17100362) exhibiting low genetic diversity yet significant genetic differentiation between East and West cluster, which caused the changes of amino acids (Fig. 4a and b). These loci were situated within the coding regions of an environment-associated gene MHch26g015340, which shared homology with the ERF110 gene in Arabidopsis thaliana (Table S7). The allele frequency of locus Chr26:17100362 increased with longitude (Fig. 4d) and showed a similar trend with the precipitation in March (Prec_03), while the other (Chr26:17099817) was significantly associated only with solar radiation (Srad_03) (Fig. 4e–g). In addition, we also noticed that the reference genotype of locus Chr26:17100362 was primarily distributed in low longitude areas (West of China), whereas the alternative allele was prevalent at high longitude areas (East of China) (Fig. 4c). Regarding locus Chr26:17099817, the allele T was mainly distributed in the regions with weak solar radiation, while allele G was almost fixed in areas experiencing high solar radiation (Fig. 4f). Both showed that adaptive loci frequently exhibited environmental and geographic signals.

3.5. Genetic maladaptation in the future climate models and application on protections

To evaluate whether the adaptive loci identified through GEA analysis were associated with the environmental variables using random variants, we re-ran GF again and precipitation showed the high R² weight importance, of which prec_03 was still the most important factor (Fig. S8). Taking prec_03 (precipitation) and bio_7 (temperature) as examples to assess the cumulative importance of allelic frequency, the consequences revealed that a faster turnover rate appeared in prec_03 than bio_7, which indicated that precipitation had a greater impact on changes in allele frequencies than temperature (Fig. S9). Meanwhile, we estimated the risk of maladaptation to future climate changes for natural populations and inferred the potential change of geographic distribution using 18, 851 adaptive loci (Fig. S6). We firstly predicted the genetic offsets reflecting the overall influences on allele frequency of environment using GF with five temperature factors and four precipitation variables mentioned in section 3.3 under BCC-CSM2-MR climate model. The results demonstrated that higher genetic offset values appeared in populations from Middle and East clusters, implying they are more vulnerable to future climate changes (Figs. 5 and S10). To improve the accuracy of prediction, we further used other three climate models and reached the similar conclusions (Figs. S11-S13). It suggested an increasing trend over time and with rising emissions among the four future models. Besides, we also calculated the RONA values of nine environment factors under four climate models, and found that the trends of average values of RONA changed under different scenarios were similar with genetic offset analysis (Table S8 and Fig. 6). It was worth mentioning that populations in middle and east clusters showed higher genomic vulnerability especially from 2081 to 2100 under SSP245 and SSP585 scenarios revealed by both methods, which indicated these populations would be not suitable to future environments and need more attentions for a long time. Besides, for the most important environment factor Prec_03, population MS showed the highest RONA value and low adaptive potential under all SSP models in the future (Table S8).

4. Discussion

Local adaptation promotes genetic divergence of populations, which can enable populations to sustain productivity during climate changes (Wadgymar et al., 2022). Unraveling the adaptive mechanisms is crucial for survivals and reproductions of organisms. Here, we applied landscape genomics to reveal how environments shaped the patterns of genetic variation of Actinidia eriantha and genetic signatures of adaptation to heterogeneous environments. Our results provided novel insights into the adaptation on how plants adapted to their local habitats and offered the guides on applications in conservation efforts.

4.1. Evolutionary lineages shaped by geography and genetic diversity pattern alongside longitudinal gradients

Genetic structure is shaped by multiple factors including gene flow, selection, population demographic history and landscape ecology (Chen et al., 2023). With PCA and clustering analysis, we observed that natural populations of Actinidia eriantha were divided into four main clusters from east to west mirroring their geographic origins (Fig. 1). This pattern of population differentiation had also been observed in other plants such as ginkgo (Zhao et al., 2019), Emmenopterys henryi (Xu et al., 2021), and Phyllostachys edulis (Zhao et al., 2021). In addition, we noted a longitudinal correlation with genetic diversities gradually increasing from west to east, with the highest genetic diversity observed in Cluster East, particularly among populations proximate to Wuyi Mountains, and the lowest in the Southeast group (Figs. S3 and S4). The higher genetic diversity in Cluster East also reflected the demographic history of A. eriantha revealed by chloroplast and microsatellite markers, as the Wuyi Mountains act as a refugium of A. eriantha during the glacial periods (Guo et al., 2022). Intriguingly, we found evidence suggesting a decline in effective population size during the Mid-Pleistocene Transition and the late glacial maximum (10, 000–20, 000 years ago) in all groups except the east cluster (Fig. S2), further supporting the role of the areas surrounding Wuyi Mountains as a historical refuge for A. eriantha (Guo et al., 2022; Huang et al., 2023; Ma et al., 2022).

4.2. Precipitation and solar radiation driven genetic variation and contributed to local adaptation

To our knowledge, natural populations spanning wide geographic ranges are subjected to selection pressures imposed by heterogeneous environments, which often shape the patterns of genetic diversity. For instance, precipitation, one crucial climate factor, exerts contrasting impacts on plants, with appropriate amounts and timing playing pivotal roles in growth, reproduction, and distribution (Li et al., 2019; Xiao et al., 2024). Previous studies have suggested that Actinidia eriantha is sensitive to precipitation, and optimizing water absorption while maintaining an appropriate water balance is crucial for preventing both water deficits and waterlogging due to the species' root system structure and spatial distribution (Xiloyannis et al., 2023). Additionally, Guo et al. (2022) found that the radical expansion of habitats for A. eriantha occurred during the Last Glacial Maximum, while the species faced its most restricted distribution during the Last Interglacial, similar to Taxus cuspidata (Luo et al., 2024). It identified precipitation played a significant role in the distribution and adaptation of kiwifruit to some extent. In our study, the outcome of GF analysis also provided additional support for precipitation as the primary driver of genetic variations (Figs. 3 and S5). Besides, we further identified several homologous genes associated with precipitation or drought stresses, including GBF3, PUBs and ARF1 (Table S7) (Ramegowda et al., 2017; Song et al., 2017; Verma et al., 2022). Among these genes, U-box proteins have been demonstrated to play a significant role in adaptation and are correlated with drought stress in A. eriantha (Zhang et al., 2023b).

Besides precipitation, sunlight is also one of the most important sources of energy, tightly associated with plant growth and development (Dotto et al., 2018). It has been proven that solar radiation serves as a crucial factor influencing genetic variations and adaptation (Bosco et al., 2020; Wang et al., 2008). For example, Triticum dicoccoides between sun and shade habitats exhibited a significant difference of overall variations and showed a non-ignorable genetic differentiation, implying that solar radiation driven natural selections in wheat (Ren et al., 2017). In our study, we not only observed that solar radiation significantly impacts genetic variations as indicated by GF (Figs. 3 and S5), but also found associated homologous genes such as RAD54, Lhcb1.1, and SPL (A.-H.-Mackerness et al., 1998; Hirakawa et al., 2017; Xie et al., 2020). Hence, it is reasonable to infer that precipitation and solar radiation variables are likely to be important factors driving genetic divergence and promoting adaptation within Actinidia eriantha.

4.3. Genomic signatures of adaptation to heterogeneous environments

Natural selection consistently leads to the genetic and phenotype differentiation due to heterogeneous environments especially for trees with wide distributions, as well as resulted in local adaptation of populations, and numerous genes were involved in this process (Dorman et al., 2009; Morente-Lopez et al., 2022; Nocchi et al., 2023; Yuan et al., 2023). An effective tool to deeply detected adaptation-associated loci along environmental gradients is GEA, which could illustrate the genetic basis of local adaptation. However, we obtained different numbers of variant sites using datasets that include and exclude the HA population (Fig. S6). This significant difference in the number of outlier loci suggested populations at distribution limits influence their adaptive capacity to environmental changes, and high genetic differentiation between them and other populations had great influences on detection effectiveness in GEAs as well (Mendez-Cea et al., 2023).

Totally, we identified 18, 851 outliers linked to the environmental adaptation and uncovered 1, 566 genes (Table S7), some of which played vital roles in abiotic and biotic stresses, phenotypic traits regulation. For instance, AP2/ERF transcript factors were considered as participating in growth, development of plants, and they were also involved in flowering development in Actinidia eriantha (Jiang et al., 2022). We identified multiple members of AP2/ERF family including ERF3, ERF14 and ERF096. Thereinto, ERF096 has been implicated in regulation of Abscisic Acid (ABA) response and selenium tolerance in Arabidopsis thaliana (Jiang et al., 2020; Wang et al., 2015), indicating several ERF genes have similar functions to promote adaptation to stresses in A. eriantha. Besides, we also demonstrated that calcium-dependent protein kinases, including CPK5, CPK11, and CPK15 (Table S7). Their homologous genes, which have been identified in A. thaliana, play crucial roles in plant signal transduction, thereby mediating plant development and environmental responses (Yip Delormel and Boudsocq, 2019).

4.4. Spatial pattern of adaptation loci shaped by environment factors

More importantly, among these adaptive loci, we discovered clear genetic footprints of divergent selection at two significant loci of AeERF110 which exhibited divergent genotypes between populations of the East and the West. The frequencies of alleles at these loci in the exon were tightly correlated with environmental factors, particularly extremes in solar radiation and precipitation (Fig. 4). This unique distribution pattern of allele frequencies has been identified not only in kiwifruit but also in other plant species, such as Populus koreana (Sang et al., 2022). A similar phenomenon was also reported in Populus cathayana, which showed one adaptive SNP of UVR8 exhibited different alleles in regions with high and weak solar radiation, which played an important role in high-altitude adaptation (Xiang et al., 2024). They also reaffirmed the significance of precipitation and solar radiation in shaping patterns of genetic variation and provided reliable evidence of local adaptation in Actinidia eriantha. Besides, all the cases discussed above indicated that environmental gradients could exert different and strong selective pressures on plant populations, and selections due to heterogeneous environments likely shaped the spatial and temporal distributions of genetic variations across natural populations of plants, potentially driving local adaptation.

It also needs to be pointed out that ERF110 as a crucial gene in flowering time control has been documented in various species, such as chrysanthemum (Huang et al., 2022) and Morusa-tropurpurea (Dai et al., 2023). In angiosperms, flowering time as one of important traits is tightly associated with potentiality of adaptation to local environment and adjusting optimal flowering time based on ambient environments is an important strategy for plants to adapt to climate changes (Price et al., 2020; Shim et al., 2017). Flowering time was associated with many factors, of which rainfall and solar radiation were also sought as important factors (Moore and Lauenroth, 2017). For example, UV-B as a main component of solar radiation can delay flowering time in Arabidopsisthaliana through regulating expressions of MSI1 and CLF (Dotto et al., 2018). Therefore, we speculated precipitation and solar radiation might modulate flowering time through genotype variations, which promoted the adaptation of A. eriantha in local habitats. However, further works including point mutations, transgenic experiment, stress experiments and RNA-seq, etc. are needed to be conducted in the future for illuminating detailed mechanisms about functional differentiation of ERF110 alleles in mediating adaptation to local environments. Besides, in this study, allelic variation in ERF110 may modulate flowering time in A. eriantha, which presents an opportunity to take genome editing technologies such as prime editing in kiwifruit germplasm conservation (Yin et al., 2024; Zhou et al., 2024).

4.5. Genetic mismatches to the changing climates and managements to protect the Actinidia eriantha

Changes of local environments would bring challenges for organisms, especially for sessile plants. As climate warming leads to the loss of original habitats, plants are forced to migrate to new areas, such as higher altitudes or latitudes, in search of suitable environments (Aguirre-Liguori et al., 2019; Wu et al., 2022). In the face of potential decrease of individual numbers and distribution areas under environmental change, efforts should be made to explore populations with pre-adaptive genotypes and enhance the gene flow between populations, which may help these risk populations mitigate the impact of future climate change (Kremer et al., 2012). However, more researches now predict the changes based on species distribution models (SDM), which ignore evolution history and genetic change of population and limited in estimated the adaptive potentials of populations (Shen et al., 2022; Thuiller et al., 2019; Wang et al., 2021). To overcome the gaps of knowledge, we used Gradient Forest to value the genomic vulnerability, which was widely used to other species (Bay et al., 2018; Jiang et al., 2021). Results showed that populations in middle and east regions were more vulnerable and sensitive to future environment changes in our study, which was urgent to propose suitable rescue and protection measures. It is important to emphasize that effective biodiversity conservation should firstly build upon traditional methods such as land conservation, species protection, and policies aimed at preventing deforestation. Additionally, it should purposefully integrate local adaptation strategies to manage gene flow, as well as focus on replanting, reintroduction, and restoring resilient populations in areas where local extinctions have occurred. In regions with high risks, particularly in the middle and east populations, utilizing a composite seed source that includes both native (seeds from plants adapted to local environments) and non-native seeds (seeds introduced from foreign areas) can better boost population diversity and resilience. This multi-source approach may effectively enhance diversity in central populations and provide a buffer against climate change, thereby increasing genetic diversity and reducing the risk of inbreeding depression (Aitken and Bemmels, 2016; Meek et al., 2023). Furthermore, through assessment of adaptive capacity, locally-adapted individuals may also be translocated to areas where they can survive in suitable habitat climate conditions (Thurman et al., 2022). These strategies are also applicable to other species in similar circumstances and offer valuable guidance for conservation efforts.

Although genetic offset is an increasingly popular approach to assess genetic maladaptation under changing environments and used to provide implications for biodiversity conservation, it should be also emphasized the limitations and caveats of genetic offset. First, genetic offset analysis does not provide the evidence that outlier loci actually influence the fitness in relation to climate. Besides, it assesses the risk of maladaptation based on the current climate and genetic data, but could not predict how populations respond to climate change in future. In addition, genetic offset analysis assumes that the relationships between environments and genetics remains unchanged from the present to the future (Aguirre-Liguori et al., 2023; Capblancq et al., 2020; Rellstab, 2021; Rellstab et al., 2021). Given these limitations, the results of genetic offsets should be treated cautiously, and further research is needed to validate its assumptions with experiments.

5. Conclusion

In the present study, large-scale genome resequencing was used to shed light on population demographic history, genetic diversity and genetic underpinnings of local adaptation. We constructed a detailed genetic variation map encompassing both SNPs and indels in kiwifruit. In addition, we pointed out that the environmental factors driven the process of local adaptation and genotype-environment associations revealed 1566 genes contributed to the adaptation to various environments. Among them, AeERF110 showed strong geographic and climatic signatures implicated in local adaptation. Moreover, we predicted that populations from east and middle cluster showed higher genetic offsets under future climate changes, which could contribute to designing effective management and conservation strategies. Overall, our study offered a novel understanding of how species adapt to heterogeneous environments and provided the fundamental conservation strategies.

Acknowledgements

This study was funded by the National Natural Science Foundation of China (grants number 32070377 and 31770374) and Science Fund for Creative Research Groups of the Natural Science Foundation of Hubei Province (2024AFA035).

CRediT authorship contribution statement

Quan Jiang: Writing – review & editing, Writing – original draft. Yufang Shen: Software. Lianhai Wu: Resources. Zhengwang Jiang: Resources. Xiaohong Yao: Writing – review & editing, Conceptualization.

Data availability

Raw sequence reads are deposited in the SRA under BioProject PRJNA916086.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.02.003.