1. Introduction

As ecosystems undergo unprecedented alterations in the face of changing climates, understanding the genetic basis of species' adaptability becomes increasingly important (Hosius et al., 2001; Cote and Darling, 2010; Nolan et al., 2018). Introgression describes the transfer of alleles between species through hybridization and repeated backcrossing, and can play a key role in to how species adapt in the face of global change (Harrison and Larson, 2014; Sun et al., 2014; Suarez-Gonzalez et al., 2018; Aguillon et al., 2022; Xiao et al., 2023). Although introgressive hybridization is often stochastic, there is increasing evidence that adaptive introgression, where there is selective retention of beneficial introgressed alleles, can play a key role in enhancing species' survival and ecological resilience (Rieseberg, 2001; Genner and Turner, 2012; Twyford and Ennos, 2012; Savolainen et al., 2013; De et al., 2015; Slatkin and Racimo, 2016; Chen et al., 2018a; Oziolor et al., 2019; Leroy et al., 2020; Fu et al., 2022). Such beneficial alleles may influence traits of recipient species like stress tolerance or phenological shifts. Once integrated into the recipient species' genomic background, these alleles can interact epistatically with native loci, accelerating a population's ability to adapt to novel environments (Barton, 2001; Whitney et al., 2015; Suarez-Gonzalez et al., 2018; Gagalova et al., 2022; Wang et al., 2023a). Thus, gene flow among species not only increases within-species genetic diversity, but can also promote evolutionary innovation, offering a crucial mechanism of adaptation.

Despite widespread recognition of the importance of introgression in facilitating local adaptation, particularly in the context of long-lived organisms like trees (Petit and Hampe, 2006; Savolainen and Pyhajarvi, 2007), our current understanding of the ramifications, genetic framework, and the specific traits influenced by adaptive introgression remains relatively constrained. However, recent developments in bioinformatics have substantially enhanced our ability to detect and characterize introgression (Than et al., 2008; Durand et al., 2011; Pickrell and Pritchard, 2012; Martin et al., 2015; Pease and Hahn, 2015). These advancements have revealed a number of valuable insights, such as the ability to identified conserved alleles shared between crops and wild relatives (Stewart et al., 2003), adaptive introgression between archaic and modern humans (Racimo et al., 2015) and pervasive hybridization during diversification of Heliconius butterflies (Edelman et al., 2019). Indeed, there is considerable evidence from a number of tree species suggesting that between-species introgression can facilitate adaptation to local environments (Suarez-Gonzalez et al., 2018; Ma et al., 2019; Numaguchi et al., 2020; Fu et al., 2022). The use of genotype–environment association analyses can help identify relationships between genotypes and environmental factors (Forester et al., 2018) and reveal genes important for adaptation (Zhou et al., 2014a; Feng et al., 2023a, 2023b; Liu et al., 2024; Lu et al., 2024). When coupled with introgression analysis, this allows precise identification of introgressed alleles that facilitate adaptation, empowering strategies for biodiversity conservation and climate change mitigation (Brancourt-Hulmel, 2000).

Among widespread forests dominated by coniferous trees, members of the genus Picea are both ecologically dominant and genetically diverse, making them a compelling model for studying introgression and adaptation (Lockwood et al., 2013; Ru et al., 2018; Wang et al., 2023b). In North American spruces, for example, there is evidence for both hybridization-associated adaptive clines (e.g., in P. glauca × P. sitchensis hybrids) (Hamilton et al., 2013) and molecular signatures of local adaptation across climatic gradients in P. mariana, including evidence of convergent evolution in genes associated with phenology and stress tolerance (Prunier et al., 2012). However, less is known about the role of adaptive introgression in facilitating environmental adaptation among Picea species from other parts of the world, presenting an opportunity for further investigation in this climate sensitive biodiversity hotspot.

Here, our primary focus rests on the dynamic interplay of adaptive introgression among Picea meyeri Rehder & E.H. Wilson and its closely related counterparts, P. crassifolia Kom. and the more recently diverged P. asperata Mast. (Bi et al., 2016; Feng et al., 2019, 2023a, 2023b). P. meyeri thrives in northern Hebei and Shanxi Provinces, as well as the southeastern part of the Inner Mongolia Autonomous Region (Han et al., 2023). P. meyeri is listed as Near Threatened on the IUCN Red List and dominates the cold evergreen coniferous forest belt, ranging from 1850 to 2700 m above sea level. This species is well-adapted to moderately moist and cool habitats, typically at relatively lower elevations where conditions are neither excessive dry nor overly humid. It prefers well-drained, slightly acidic to neutral soils and exhibits considerable ecological flexibility (Liu et al., 2024). Nearby, P. crassifolia flourishes in drier, colder environments, while P. asperata is adapted to wetter, warmer climates (Feng et al., 2023a). These distinct ecological preferences and their overlapping distributions make this trio of species an ideal system for investigation how adaptive introgression may have contributed to species divergence and environmental resilience. Moreover, all three species hold significant economic and ecological value, contributing to forest ecosystem stability (Xue et al., 2007; Wang et al., 2020a, 2021).

Our research aims to illuminate the processes of adaptive introgression within this complex ecological framework, focusing on the interactions among Picea meyeri and its close relatives, P. crassifolia and P. asperata. We seek to quantify adaptive introgression alleles and identify the specific traits they influence, with particular attention to how introgression aids local adaptation and divergence in these species. Our results will enhance our understanding of the genetic mechanisms underlying adaptive introgression and its role in the evolution of these ecologically important spruce species. Ultimately, our work contributes to a broader understanding of the mechanisms driving ecological adaptation, especially through introgression among closely related taxa, in topographically complex regions.

2. Materials and methods 2.1. Material and RNA sequencing

We collected fresh mature needles from first-year branches from natural populations of Picea asperata, P. crassifolia, and P. meyeri from China. In total, we sampled a total of 114 individuals from 21 populations (Fig. 1 and Table S1). Sampling locations were recorded using an eTrex GIS (Garmin, Germany). Immediately after collection, we froze needles in liquid nitrogen and stored at −80 ℃ until RNA extraction. To reduce sampling bias in estimates of genetic diversity, we sampled individuals within each population at least 100 m apart.

Total RNA was extracted from the needles of each individual using a modified CTAB method (Ma et al., 2015; Amraee and Rahmani, 2020). For RNA sequencing, we adopted the Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA) to build libraries and perform paired-end sequencing of the extracted RNA. A minimum of 8 Gb of raw sequence data was generated for each individual.

2.2. Reads mapping and variant calling

Raw reads were preprocessed using FastQC v.0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for initial quality assessment and fastp v.0.20.0 (Chen et al., 2018b) for quality filtering, using the following parameters: -c -q 15 -u 40 -g -n 5. High-quality paired-end reads that remained after filtering were mapped to a revised transcriptome of Picea abies, excluding previously reported fungal transcripts (10.5281/zenodo.15009537) (Nystedt et al., 2013; Delhomme et al., 2015; Ru et al., 2016; Feng et al., 2023a, 2023b) using the MEM algorithm of BWA v.0.7.17 (Li and Durbin, 2009) with default settings. The resulting binary alignment map (BAM) files were sorted using SAMtools v.0.1.19 (Li et al., 2009). After sorting, duplicate reads were marked using Picard Tools v.2.21.5 (https://broadinstitute.github.io/picard) and excluded from further analyses. Indels were detected and realigned using the "RealignerTargetCreator" and "IndelRealigner" tools from The Genome Analysis Toolkit (GATK) v.3.8 (DePristo et al., 2011).

Single nucleotide polymorphisms (SNP) calling was conducted using the "mpileup" command in SAMtools v.0.1.19 (Li et al., 2009). The resulting SNPs dataset was filtered based on the following criteria: (1) Only genotypes with a Phred-scaled likelihood score of 20 or higher were retained to ensure a genotyping accuracy rate of at least 99%; Genotypes with lower quality scores were treated as missing values; (2) Only biallelic loci were retained for subsequent population-level analyses; (3) SNPs located within indels (including a 5 bp buffer around the indel) were excluded to minimize false positive rates; (4) SNPs with a non-reference (ALT) allele frequency of less than 1% were removed; SNPs with a minor allele frequency (MAF) of less than 5% were removed; (5) Sites with depth (DP) < 10 per individual were considered missing, and SNPs were included only if the missing rate within each species was < 50%. After applying these filters, a subset of high-quality, eligible SNPs was retained for further analysis. The filtered SNP dataset was annotated using SnpEff v.4.3t (Cingolani et al., 2012) and the reference genome annotation file was obtained via Transdecoder v.5.7.1 (https://github.com/TransDecoder/TransDecoder) (Kim et al., 2016).

2.3. Analysis of population structure, phylogenetic and genetic diversity

To investigate the population structure of each of the three Picea species, we performed a principal component analysis (PCA) on the high-quality SNPs. To do so, the variant call format (VCF) file was first converted to binary PED format using VCFtools v.0.1.14 (Danecek et al., 2011) and PLINK v.1.07 (Purcell et al., 2007) software with the parameters -indep-pairwise 50 5 0.2 to prune SNPs based on linkage disequilibrium (Purcell et al., 2007). We conducted the PCA using Eigensoft-smartpca v.6.1.4 (Price et al., 2006), based on the SNPs dataset with default parameters and visualized it using R v.4.4.1 software. We determined the significance levels of the principal components using the Tracy–Widom test. To further infer the population structure, we used ADMIXTURE v.1.3.0 (Alexander et al., 2009) with the number of genetic clusters (K) ranging from 1 to 10. The optimal K value was identified by minimizing the cross-validation error through a 10-fold cross-validation procedure.

In addition, with the unpruned high-quality SNPs data, we constructed a maximum likelihood (ML) phylogenetic tree for the three Picea species using IQ-TREE v.2.0.3 (Nguyen et al., 2015) with 1000 bootstrap replicates to ensure robust support for the inferred relationships. We identified the TVM + F + R10 substitution model as the best-fit model by using model testing and comparison with the built-in ModelFinder tool (Kalyaanamoorthy et al., 2017), which we then used for tree construction. P. breweriana was used as an outgroup to root the tree to help interpret evolutionary relationships among the species. To characterize genetic diversity within each species and genetic differentiation among the species, we calculated nucleotide diversity (π) and pairwise genetic differentiation (F_ST) per locus using VCFtools v.0.1.14 (Danecek et al., 2011).

2.4. Analysis of gene flow 2.4.1. ABBA-BABA test

To detect gene flow between species, we calculated Patterson's D statistic using the Dtrios module in Dsuite v.0.5 (Malinsky et al., 2021). This analysis assumes a phylogenetic relationship among the three taxa (P1, P2, P3) and an outgroup (O), constructed as {[(P1, P2), P3], O}, where P1 is closely related to P2 but more distantly related to P3. For biallelic SNPs, the ancestral allele carried by the outgroup (O) is denoted as A, and the derived allele is denoted as B. In this context, an ABBA site represents a pattern where P1 carries the outgroup allele (A), while P2 and P3 share the derived allele (B). Conversely, a BABA site corresponds to a pattern where P1 and P3 share the derived allele (B) and P2 carries the outgroup allele (A). Significant results from the ABBA-BABA test can indicate both gene flow and population subdivision in the ancestral species (Slatkin and Pollack, 2008; Durand et al., 2011).

Patterson's D specifically measures the extent of gene flow: if the D statistic is significantly greater or less than 0, this suggests that there is significant gene flow between P3 and either P2 or P1 (Green et al., 2010; Patterson et al., 2012; Freitas et al., 2021). Patterson's D is calculated as D = (nABBA-nBABA)/(nABBA + nBABA), where nABBA is the total count of the ABBA patterns and nBABA is the total count of the BABA patterns. To determine the significance of the D statistic, we calculated Z-scores for each contig, with a Z-score absolute value > 3 considered significant (Durand et al., 2011; Novikova et al., 2016). We applied the ABBA-BABA test to the SNPs dataset in order to evaluate potential gene flow between Picea asperata (P1), P. crassifolia (P2) and P. meyeri (P3), with P. breweriana serving as the outgroup (O) (Martin et al., 2015).

2.4.2. TreeMix

TreeMix is a method used to estimate population splits and migration events by building tree models incorporating division and gene flow among populations. This method builds a maximum likelihood graph based on the covariance matrix of allele frequencies between populations, inferring relationships among populations and their common ancestors. It also accounts for admixture events (referred to as "migration") to improve the accuracy of the inferred tree. Using the allele frequency data, we used TreeMix v.1.13 (Pickrell and Pritchard, 2012) to generate a ML tree for distinct populations of the three Picea. By mapping gene flow events onto the ML tree, we could visualize and analyze historical gene exchange between populations. By examining the distribution and direction of these gene flow events, we could identify potential cases of historical admixture among target species.

2.5. Screening for introgressed alleles

To identify introgressed alleles between pairs of species, we calculated the modified f-statistic (fdM) (Malinsky et al., 2015; Martin et al., 2015) across the transcriptome. The fdM is a statistical indicator of introgression that measures changes in allele sharing between P3 and P2 or between P3 and P1. We used ABBABABAwindows.py in genomics_general v.0.4 (https://github.com/simonhmartin/genomics_general) to calculate fdM. This approach uses sliding windows across a dataset of SNPs (-w 10000 -m 10 -s 1000) to assess allele introgression at the genetic level. For the analysis, we designated P1 as Picea asperata, P2 as P. crassifolia, P3 as P. meyeri, and P. breweriana as the outgroup (O).

To explicitly define introgressed alleles, we applied the following criteria: (1) The absolute value of fdM exceeded 0.2 (|fdM| > 0.2), reflecting excess allele sharing between P3 and either P1 or P2; (2) At least three consecutive 10000 bp windows showing fdM values above this threshold, minimizing the influence of random variation; (3) Interspecific F_ST values for loci with introgressed alleles (calculated using VCFtools v.0.1.14) (Danecek et al., 2011) were significantly lower than those without, suggesting reduced genetic differentiation between species at these loci (Feng et al., 2023b). For simplicity, we defined allele sharing between P1 and P3 as AM_introgression, and between P2 and P3 as CM_introgression. We annotated AM_introgression and CM_introgression using SnpEff v.4.3t (Cingolani et al., 2012) and counted the number of synonymous and non-synonymous mutations they contained.

Finally, to determine whether any functional classes of loci were overrepresented among the candidate introgression loci, we performed a functional annotation analysis using EggNOG-mapper v.2.1.12 (Cantalapiedra et al., 2021) with an E-value < 1.0 × 10⁻³. Gene ontology (GO) enrichment analysis was performed using TBtools v.2.096 (Chen et al., 2020). The Benjamini-Hochberg FDR (false discovery rate) was used to correct p-values, and GO terms with corrected p-values ≤ 0.05 were considered enriched.

2.6. Detection of selective sweeps

To quantify polymorphism (θπ, pairwise nucleotide variation) within each species, as well as genetic differentiation (F_ST) between each species pair, we used VCFtools v.0.1.14 (Danecek et al., 2011) with 10 kb sliding windows with a 1 kb step size. All analyses were conducted using a pruned SNPs dataset to minimize the effects of linkage disequilibrium. We transformed the ratios of θπ using log₂ conversions for further analysis. Using the Picea crassifolia-P. meyeri pair as an example, selective signals for P. meyeri and P. crassifolia were identified based on the top 5% of F_ST values and the top 5% of log₂(θπcra/θπmey) and log₂(θπmey/θπcra), respectively (Li et al., 2013; Freitas et al., 2021; Zhang et al., 2022b).

2.7. Identification and analysis of adaptive introgressed alleles

To identify and characterize adaptive introgressed alleles, we analyzed the following loci: (1) Introgressed alleles overlapping with positively selected genes (PSGs): These include introgressed loci between Picea crassifolia and P. meyeri (PCM_introgression) and between P. asperata and P. meyeri (PAM_introgression) that overlap with genes under positive selection. (2) Environment-associated introgressed PSGs: These are positively selected introgressed loci associated with environmental factors, identified through redundancy analysis (RDA). This category includes environment-associated PCM_introgression alleles (PRCM_introgression) and PAM_introgression alleles (PRAM_introgression). The RDA was performed using the vegan v.2.6.8 (Dixon, 2003) package in R. To investigate the functional roles of these loci, we aligned their sequences with the Arabidopsis thaliana protein database using BLASTX, implemented by DIAMOND v.0.9.14 (Buchfink et al., 2015), with an E-value < 1.0 × 10⁻⁵.

2.8. Isolation by distance and isolation by environment

To examine how geographic and environmental factors influence genetic composition, we evaluated isolation by distance (IBD) and isolation by environment (IBE) effects. We calculated genetic distances among the three species using Arlequin v.3.5 (Excoffier and Lischer, 2010) and environmental distances as Euclidean distances using the vegan v.2.6.8 (Dixon, 2003) R package. We derived pairwise geographic distances from sampling site coordinates using the Geosphere v.1.5–20 package (https://github.com/rspatial/geosphere) in R. To disentangle the effects of IBD and IBE, we used a partial Mantel test to quantify IBD while controlling for environmental effects using the vegan v.2.6.8 (Dixon, 2003) package in R. The significance of the partial Mantel tests was assessed through 9999 permutations.

2.9. Analysis of species predicted distributions

We extracted a total of 92 environmental variables under current conditions from the World Bioclim Database (Worldclim, http://www.worldclim.org) with a resolution of 2.5 arc-min. These variables included 19 bioclimatic variables, as well as average monthly climate data for precipitation, solar radiation, wind speed, water vapor pressure, elevation, and minimum, mean, and maximum temperatures. We used presence data from 21 sites for Picea asperata, 35 sites for P. crassifolia, and 19 sites for P. meyeri for species distribution modeling. To avoid multicollinearity among environmental variables, we calculated Pearson's correlation coefficients between each pair of variables. For any pair with |r| > 0.7, we eliminated the variable with the lower relative contribution to the model (Dormann et al., 2013; Wang et al., 2023c). We estimated paleoclimate data for the mid-Holocene and the Last Glacial Maximum based on three different GCMs (CCSM4, MIROC, and MPI). We used species distribution modelling (SDM) to simulate the distribution of each species using MaxEnt v.3.3.4 (Phillips and Dudik, 2008). Suitable habitat maps were created using ArcGIS v.10.6.1 (ESRI, Redlands, CA, USA) to visualize model results and assess potential distribution shifts across climate contexts.

3. Results 3.1. Transcriptome sequencing and genetic variants

In all, we achieved a mean sequencing coverage of 72% and a mean mapping ratio of 73% (Table S2). In addition, using strict quality filters, we identified 170, 406 high quality SNPs from 114 individuals from the three Picea species. These SNPs were distributed across 14, 765 genes. Among them, 60, 427 were synonymous mutations, and 66, 125 were non-synonymous mutations. Specifically, we identified 124, 054 SNPs from 25 individuals across 6 populations of P. asperata, 136, 204 SNPs from 59 individuals across 12 populations of P. crassifolia and 122660 from 30 individuals across 3 populations of P. meyeri. We annotated SNPs for each population (Table S3) and found that the number of total variants (synonymous and non-synonymous) was highest in P. crassifolia and lowest in P. meyeri. In addition, the number of shared SNPs was greater between P. crassifolia paired with either P. asperata (105, 770) or P. meyeri (99, 730) than between P. asperata and P. meyeri (96, 743). This indicates closer genetic relationships between P. crassifolia and the other two species.

3.2. Population genetic structure and genetic diversity

Using the transcriptome of Picea breweriana as an outgroup, we constructed genealogies for all 114 individuals using ML methods. The phylogenetic analysis based on SNP variation showed that P. meyeri samples grouped into a single clade, while P. asperata and P. crassifolia were further subdivided, so that there were three distinct clades (Fig. 2A). The results from ADMIXTURE analysis were consistent with this observation. When we performed the analysis with K = 2, P. meyeri diverged first from the other two species and formed a separate cluster (Fig. 2A). When K = 3, the lowest cross-validation error (0.3859) resulted in the division of individuals from the three species into three distinct genetic clusters (Fig. 2A and Fig. S1).

PCA clustering results corroborated the findings from both the ADMIXTURE analysis and the phylogeny tree. The first principal component clearly distinguished Picea meyeri from two other species (Fig. 2B). The second principal component further separated P. asperata from P. crassifolia (Fig. 2B). The three species had similar genetic diversity: P. asperata (0.003532 ± 0.003276), P. crassifolia (0.003043 ± 0.003044) and P. meyeri (0.003474 ± 0.003277) (Table 1). We found moderate genetic differentiation (pairwise F_ST) among the three species. 0.0820 ± 0.1090 between P. asperata and P. crassifolia; 0.1040 ± 0.1320 between P. crassifolia and P. meyeri; 0.1290 ± 0.1540 between P. asperata and P. meyeri (Table 1).

3.3. Gene flow among three species and analysis of gene introgression

Our ABBA-BABA tests yielded a D value of 0.0430 and a Z-score of 10.9176, indicating that Picea meyeri shares more allelic loci with P. crassifolia than with P. asperata (Fig. S2). The results from the TreeMix analysis suggest the ML tree closely aligns with that constructed using IQ-TREE. This topology revealed that P. asperata and P. crassifolia were positioned on the same major branch, with P. meyeri diverging earliest (Fig. 3B). We also detected gene flow among the three species, with particularly strong evidence between P. crassifolia and P. meyeri.

The fdM analysis detected 2246 candidate loci of introgressed alleles between Picea asperata and P. meyeri (AM_introgression), which is fewer than the 2467 loci of introgressed alleles identified between P. crassifolia and P. meyeri (CM_introgression) (Fig. 3A). Subsequently, we calculated F_ST for both introgressed and non-introgressed allelic loci in the species pairs P. asperata and P. meyeri, P. crassifolia and P. meyeri, finding that the F_ST values of both species pairs the introgressed alleles were consistently lower than those of the non-introgressed allelic loci (Table S4). Among the introgressed allelic loci, P. asperata and P. meyeri exhibited higher F_ST values than P. crassifolia and P. meyeri, which was similar to that in the non-introgressed allelic loci (Table S4). Further detection of the number of mutations in the introgressed alleles, we found that there were 2121 synonymous and 2221 non-synonymous mutations in the AM_introgression, whereas there were 2328 synonymous and 2271 non-synonymous mutations in the CM_introgression.

GO enrichment analyses were conducted on AM_introgression and CM_introgression datasets. In the AM_introgression dataset, 53 allelic loci exhibited significant enrichment. Among these, 10 loci were assigned to the "Nucleotide phosphorylation" biological process pathway. Furthermore, the "Photosynthesis" and "Photosynthesis, light reaction" biological process pathways were enriched (GO: 0015979, GO: 0019684) (Fig. S3). In the CM_introgression dataset, nine loci were assigned to the biological process category. Notably, some of these loci are related to root morphology (GO: 0010015) (Fig. S3), suggesting that Picea crassifolia likely has more extensive root development than the other species, potentially due to its greater adaptation to arid environments. Moreover, biological processes involved with "metallo-sulfur cluster assembly" and "iron-sulfur cluster assembly" were significantly enriched, indicating long-term adaptation to metal-contaminated environments or drought conditions (Abdirad et al., 2022).

3.4. Positive selection and adaptive introgression in the three Picea species

In selective scanning analysis, we detected a total of 77 PSGs in Picea asperata (log₂(θπmey/θπasp) > 2.38) and 171 PSGs in P. crassifolia (log₂(θπmey/θπcra) > 2.42) when using P. meyeri as the reference species. In contrast, we identified 79 PSGs in P. meyeri when using P. asperata as the reference (log₂(θπasp/θπmey) < −2.35), and 43 PSGs in P. meyeri when using P. crassifolia as the reference (log₂(θπcra/θπmey) < −2.41) (Fig. 3C and D; Table 2). In combination with these findings, we detected nine candidates adaptive introgressed alleles in P. asperata, 36 alleles in P. crassifolia, four alleles in P. meyeri with P. asperata as the reference, and three alleles in P. meyeri with P. crassifolia as the reference (Table 2).

Alignment of the candidate adaptive introgressed alleles with the Arabidopsis thaliana proteome revealed several environmentally relevant genes across the three Picea (Fig. 3C and D; Table S5). In P. asperata, we identified three genes (PRAM_introgression) related to drought and salt resistance, consistent with the species' well-documented adaptation to arid conditions (Ingram and Bartels, 1996; Duan et al., 2005). Additionally, we detected a gene associated with flowering time and light response (comp86109_c0_seq1) (Table S5). This suggests that P. asperata may enhance reproductive success by modulating flowering times in response to environmental cues, which could influence resilience under changing environmental conditions.

In Picea crassifolia, we identified a total of 17 environmentally responsive genes (PRCM_introgression), most of which are involved in responses to drought, salt, and cold damage (Yu et al., 2014; Zhang et al., 2022a). Notably, we also detected genes associated with DNA damage repair (comp48457_c0_seq1, comp90811_c0_seq1, and comp95493_c0_seq1), which are likely related to mitigating ROS-DNA damage resulting from environmental stress (Zhang et al., 2017). Additionally, we detected a gene associated with flowering time (comp75039_c0_seq1) (Table S5).

In Picea meyeri, when using P. asperata as the reference, we identified two genes —RBOHD (comp72876_c0_seq1) and RHD2 (comp90619_c0_seq1) — involved in salt resistance and immune response. On the other hand, we identified three environmentally responsive genes with P. crassifolia as the reference, one of these genes — ERH1 (comp50475_c0_seq1) is involved in regulating the plant's response to pathogen attacks. This may be because P. meyeri typically inhabits moist conditions where the ability to respond to pathogens is particularly crucial, as excessive moisture is often associated with increased pathogen pressure (Oberhuber et al., 1999; Siitonen et al., 2005). Furthermore, we detected a late embryogenesis abundant hydroxyproline-rich glycolprotein (comp76784_c0_seq1), which is known for its role in drought, salt, and cold stress tolerance. This gene is likely to have originated from P. crassifolia and subsequently introgressed into P. meyeri (Goyal et al., 2005; Tolleter et al., 2010), thereby potentially enhancing P. meyeri's tolerance to environmental stresses.

3.5. Genotype-environmental association analysis

To investigate the effects of IBE on genetic structure, we selected 11 environmental variables (BIO4: temperature seasonality, BIO5: mean temperature of the warmest month in the coldest season, BIO13: precipitation in the wettest month, BIO14: precipitation in the driest month, BIO15: precipitation seasonality, SRAD10: solar radiation in October, TMAX3: maximum temperature in March, WIND1: wind speed in January, WIND6: wind speed in June, WIND12: wind speed in December, ALT: altitude) for further analyses. Subsequently, we used the partial Mantel test to examine the relationship between geographic and environmental distance and genetic distance. We found a significant relationship between genetic distance and both geographic distance and environmental distance. For the combined populations of Picea asperata, P. crassifolia, and P. meyeri, genetic distance was significantly associated with environmental distance (R² = 0.1019, p < 0.05), suggesting the impact of environmental gradients on local adaptation (Fig. S4). These findings confirm that both geographic and environmental distances significantly contribute to the genetic divergence observed in these species (Fig. S4).

To further explore the genetic basis of local adaptation, we used redundancy analysis (RDA) to identify candidate SNPs associated with the selected environmental variables. In Picea asperata and P. meyeri, we identified several environmental factors (BIO1, BIO2, BIO4, BIO13, and BIO15) as the primary drivers of genetic differentiation (Fig. S5). Likewise, in the RDA for P. crassifolia and P. meyeri, genetic differentiation was strongly influenced by several of the same variables (BIO1, BIO2, BIO4, and BIO15), and some different (BIO12, BIO14, and BIO15) variables (Fig. S5). Across all three species (P. asperata, P. crassifolia, and P. meyeri), BIO4, BIO6, BIO13, BIO14, and BIO15 were the most influential variables shaping genetic differentiation.

In Picea asperata and P. meyeri, we identified a total of 2143 potential candidate genes using RDA, we obtained four genes (PRAM_introgression) by intersecting with the previously identified PAM_introgression alleles, one of which is responsive to flowering time (comp86109_c0_seq1). In P. crassifolia and P. meyeri, we identified 1761 potential candidate genes, 17 of which were identified as PRCM_introgression alleles, by intersecting with the PCM_introgression alleles. Among them, we detected a gene associated with flowering time (comp75039_c0_seq1) (Table S5). Notably, three of these PRCM_introgression alleles were related to drought stress and DNA repair (comp48457_c0_seq1, comp90811_c0_seq1, and comp95493_c0_seq1).

3.6. Distribution models for the three Picea species

We identified differences in the key environmental factors influencing the spatial distribution of the three Picea species. We found that precipitation plays a strong role in the distribution of P. crassifolia, temperature influences the distribution of P. asperata, and both temperature and precipitation affect the distribution of P. meyeri (Table S6).

The results from three global climate models (GCMs: CCSM4, MIROC, and MPI) for the mid-Holocene and the Last Glacial Maximum indicate that the potential distribution range of Picea meyeri was larger than that of P. crassifolia and P. asperata during both periods (Fig. 4). The distribution ranges of all three species were smaller in the mid-Holocene than during the Last Glacial Maximum. However, since the mid-Holocene, P. meyeri, P. crassifolia, and P. asperata have each experienced significant habitat expansion. Moreover, the predicted contemporary distributions of these species align well with their actual distributions. The suitable habitat of P. asperata tends to be concentrated towards the southwest of its current range, while P. crassifolia is towards the northwest, and P. meyeri is towards the northeast (Fig. 4).

4. Discussion

In this study, we used high-quality transcriptome sequencing data from three closely related spruce species — Picea asperata, P. crassifolia, and P. meyeri — to explore genetic divergence among the species and characterize selection-driven introgression. By examining SNPs variation, genetic diversity, population structure, gene flow, and adaptive introgression, we provide critical insights into the evolutionary dynamics of these species. Our results help us to understand genetic mechanisms of adaptation and can also inform management strategies for spruce populations in the face of environmental changes.

4.1. Genetic differentiation and gene flow

Incomplete reproductive isolation among spruce species leads to gene flow, a well-documented phenomenon (Bouille and Bousquet, 2005; Zou et al., 2013; Ru et al., 2018; Feng et al., 2019; Zhu et al., 2024). In particular, population genomic studies of Picea asperata and P. crassifolia reveal substantial genetic exchange. This suggests that recurrent climatic oscillations that alternatively isolate and reconnect populations may have played a crucial role in shaping their evolutionary trajectories (Feng et al., 2023a, 2023b). This pattern is further supported by evidence of cryptic speciation in P. brachytyla (Lu et al., 2024). In these studies, genes exhibiting interspecific divergence were strongly associated with local environmental adaptations, suggesting that adaptation to distinct ecological niches can occur despite ongoing gene flow.

Our phylogenetic, ADMIXTURE, and PCA analyses revealed distinct genetic clusters among the three spruce species, reflecting their genetic distinctiveness. Moreover, our ABBA-BABA tests and TreeMix analyses provided robust evidence of gene flow among these species, particularly between Picea meyeri and P. crassifolia (Fig. 3A and B). In addition, our analyses of pairwise F_ST values revealed that lowest differentiation was between P. asperata and P. crassifolia, while the highest was between P. asperata and P. meyeri that between P. crassifolia and P. meyeri was intermediate (Table 1). In all, these findings suggest that these species maintain obvious genetic distinctiveness despite ongoing gene flow (Fig. 2) (Hendry and Taylor, 2004; Hendry et al., 2007; Feng et al., 2019). Numerous studies have demonstrated that adaptive divergence between closely related species can proceed despite gene flow, and that significant geographic isolation is not always required (Nosil, 2008; Thibert-Plante and Herdry, 2011; Richards et al., 2019). Furthermore, divergence between closely related species is often accompanied by ongoing gene flow and the presence of hybrid individuals in regions where their distributions overlap (Hoskin et al., 2005; Martin et al., 2013; Ma et al., 2018; Gao et al., 2020; Zhu et al., 2024).

4.2. Environmental influences on species distribution

Our analysis of key environmental factors that affect the spatial distribution of the three spruce species reveals clear ecological differentiation among the species. The distribution of Picea asperata is primarily temperature-dependent, with this species thriving in stable, moderate climates typical of high altitudes (Fig. 4 and Table S6). In contrast, P. crassifolia is strongly influenced by precipitation, preferring relatively environments that are essential for its growth in otherwise arid regions (Feng et al., 2023a). P. meyeri, however, seems less constrained by either temperature or precipitation, exhibiting greater ecological plasticity and occupying a broader range of habitats (Fig. 4 and Table S6) (Liu et al., 2024).

Global climate models (GCMs: CCSM4, MIROC, and MPI) for the mid-Holocene and the Last Glacial Maximum support these patterns, indicating that Picea meyeri has long maintained a broader ecological niche than the other two species, likely due to its adaptive flexibility. Cooler and drier conditions during the Last Glacial Maximum likely restricted the ranges of P. asperata and P. crassifolia, while P. meyeri persisted across a wider area. Interestingly, the distributions of all three species were smaller during the mid-Holocene than during the Last Glacial Maximum (Feng et al., 2023a, 2023b; Kimmitt et al., 2023; Korznikov et al., 2023). Since the mid-Holocene, all three species have undergone substantial habitat expansion because post-glacial climatic amelioration provided more favorable conditions for growth and dispersal (Fig. 4) (Gerardi et al., 2010; Zani et al., 2023).

This study highlights the complex interplay between climate and species distribution. Understanding these dynamics is crucial for predicting species' responses to future climate change (Zhu et al., 2021). Accordingly, conservation strategies must consider these environmental dependencies to preserve suitable habitats for these species amid changing climates.

4.3. Signals of adaptive introgression in three closely related spruce species

As climate change accelerates, tree species face challenges adapting through migration or phenotypic plasticity — mechanisms often constrained by long lifespans and limited dispersal abilities (Kijowska-Oberc et al., 2020). Instead, genetic adaptation, particularly via interspecific gene flow, represents a critical avenue for acquiring adaptive variation (Zhou et al., 2014b; Hamilton and Miller, 2016). This process enhances evolutionary potential and is essential for maintaining resilience in long-lived species subjected to ongoing environmental stress (Mallet, 2005; Edmands, 2007).

We find asymmetric introgression, with more extensive gene flow occurring between Picea crassifolia and P. meyeri (CM_introgression), with more synonymous and non-synonymous mutations than between P. asperata and P. meyeri (AM_introgression). At the same time, the number of non-synonymous mutations is similar to the number of synonymous mutations in the two species pairs (Table S3). The pattern of adaptive introgressed alleles is consistent with this observation, with a greater number between P. crassifolia and P. meyeri (PCM_introgression) than between P. asperata and P. meyeri (PAM_introgression) (Table 2). This biased introgression —both in overall gene flow and the subset of adaptive introgressed alleles — may have facilitated the adaptive divergence between P. crassifolia and P. asperata. This could reflect stronger or more localized environmental selection acting on these two species due to their more specialized habitat requirements. In parallel, introgression from P. asperata and P. crassifolia into P. meyeri may have contributed to P. meyeri's enhanced tolerance to a wider range of environmental conditions, thereby supporting its ecological plasticity.

GO enrichment analyses of introgressed alleles further support this conclusion, revealing functions related to environmental adaptation, such as response to drought, temperature and salt stress (Fig. S3). Key adaptive introgressed alleles, including those encoding calponin-like domain protein, Rho GTPase activating protein, and ABC transporter proteins, illustrate the functional significance of positively selected introgressed allele in ecological adaptation (Fig. S3).

Our identification of environment-associated introgressed positively selected genes underscores the role of environmental pressures in shaping genetic diversity and adaptive responses. Among these, genes associated with flowering time (e.g., WNK1) (Wang et al., 2008) and stress tolerance reflect targeted genetic responses to specific environmental stimuli. These findings not only advance our comprehension of the genetic mechanisms underlying adaptive introgression, but also offer valuable insights for conservation strategies, emphasizing the importance of preserving genetic variation in the face of environmental change.

Importantly, our results align with broader patterns observed in other conifer species. For instance, genomic analyses of 13 closely related Eurasian pines revealed complex gene introgression across multiple species, with clear niche differentiation despite extensive gene flow (Zhao et al., 2024). Similarly, genome assemblies of pine haplotypes have identified genes linked to flowering regulation and abiotic stress tolerance (Jang et al., 2024), mirroring the adaptive patterns we observed here. Altitude-associated SNPs near flowering-time genes (e.g., ENSA and FTL2) in Pseudotsuga menziesii echo our findings on flowering-time genes and confirm a widespread photoperiod adaptation mechanism in conifers (Gugger et al., 2010).

Evidence from other studies on spruce further supports the idea that shared environmental pressures can lead to convergent adaptive strategies. For instance, adaption to drought in several coniferous species involves expansion or positive selection of gene families associated with water-use efficiency (Xue et al., 2022; Fu et al., 2024; Yin et al., 2024). Thus, our findings on selection-driven introgression are consistent with broader evolutionary trends in forest trees.

Adaptive introgression likely also contributed to morphological and physiological differentiation among species. The three spruce species differ substantially in a number of traits, such as needle structure, cone morphology, and bark texture, that can strongly influence their survival and reproduction (Wang et al., 2020b). For example, in Picea asperata, exogenous abscisic acid (ABA) application enhances leaf mass ratio, root/stem ratio and fine/total root ratio under drought (Duan et al., 2007), while P. meyeri exhibits altitudinal conduit variation suited for moist, temperate climates (Wang et al., 2019; Li et al., 2022). P. crassifolia, in contrast, shows pronounced stomatal plasticity, reflecting its specialized adaptation to drought (Wang et al., 2019).

Physiological traits further emphasize ecological divergence. Picea asperata has lower stomatal conductance and higher photosynthetic efficiency, allowing it to maintain water and energy balance in environments with extreme moisture and temperature variation (Duan et al., 2005, 2007). P. meyeri, however, possess high water use efficiency, enabling sustained growth in variable climates (Wang et al., 2019). P. crassifolia can survive under high evaporative demand and low precipitation by reducing water loss and increasing water use efficiency, demonstrating higher drought stress tolerance as evidenced by stomatal conductance and chlorophyll content under water-limited conditions (Duan et al., 2005; Wang et al., 2019). Coincidentally, the mean annual precipitation in the native ranges of P. meyeri) (491 mm) and P. crassifolia (417 mm) is considerably lower than that in the native range of P. asperata (714 mm).

Together, these morphological and physiological distinctions underscore how environmental pressures, coupled with introgression of adaptive alleles, have contributed to species divergence. Gene flow plays a dual role by enhancing genetic diversity to facilitate adaptation, while promoting niche differentiation, which reduces competition. For example, these species can coexistent regionally because Picea asperata is suited to warm, humid zones, while P. crassifolia and P. meyeri thrive in colder or moderately humid regions.

To better understand the spatial dimension of introgression, we also examined population-level patterns of introgression across species ranges. Populations in closer proximity to the donor species show broader F_ST distributions and lower median values, which is consistent with ongoing gene flow. In contrast, more geographically isolated populations exhibit higher genetic differentiation, likely due to limited connectivity (Fig. S6). Furthermore, certain introgressed alleles (e.g., WNK1, RBOHD) (Table S5) appear to be under positive selection in some populations, indicating their role in fine-scale environmental adaptation. These findings underscore the importance of spatial context in shaping both the frequency and functional relevance of introgression.

In addition to enhancing adaptability in some instances, gene flow also carries risks. Excessive introgression can erode species boundaries and reduce habitat adaptations. In regions with steep environmental gradients, maladaptive gene flow could compromise local fitness. Therefore, it is essential that conservation and management strategies weigh both the benefits and drawbacks of gene flow. Doing so will help preserve, genomic integrity, adaptive potential and long-term ecological function in these coniferous ecosystems.

4.4. Conservation and management

Our findings underscore the important implications of genetic diversity and adaptive potential for the conservation of long-lived tree species, particularly species under rapidly changing environmental conditions. In particular, our results underscore the importance of preserving genetic connectivity for maintaining the resilience of spruce populations against environmental fluctuations. Facilitating gene flow can significantly enhance this connectivity, promoting a more flexible and evolutionarily responsive gene pool capable of withstanding future climatic and ecological challenges (Flores et al., 2023).

To support this, conservation efforts should prioritize the preservation and restoration of habitat corridors — strips of land that connect fragmented habitats and allow for the movement of individuals and their genes between populations (Gregory et al., 2014; Zhang et al., 2024). These corridors are vital both for the everyday survival and reproduction of species, but also for maintaining long term genetic diversity. In parallel, identifying and mitigating barriers to gene flow, such as physical obstructions, inbreeding, or anthropogenic disturbances, is necessary to sustain population connectivity (Statham et al., 2022).

Alongside habitat-focused strategies, long term genomic monitoring is essential to assess the impacts of gene flow and introgression over time and understand their role in shaping evolutionary trajectories and ecosystem dynamics (Ma et al., 2019). However, to fully interpret the adaptive significance of these patterns, it is equally important to validate the roles of candidate adaptive alleles. Determining how these genes contribute to critical traits — such as stress tolerance and flowering time regulation — will provide the mechanistic insight necessary for designing targeted and effective conservation strategies (Flanagan et al., 2018; Shi et al., 2024).

Comparative studies across closely related species can reveal both conserved and lineage-specific mechanisms of adaptation, offering broader insights into the genomic basis of environmental responses (De Lafontaine et al., 2015; Brunton et al., 2024). Genes associated with climate adaptation, such as those identified in our study and others (e.g., Feng et al., 2023a), serve as valuable molecular markers in applied breeding programs aimed at enhancing environmental resilience. To ensure the effectiveness and sustainability of such programs, it is essential to preserve genetic diversity by selecting diverse parental lines, avoiding inbreeding (Jones et al., 2006; Kiwuka et al., 2021), and maintaining population connectivity through habitat corridors (Gilbert-Norton et al., 2010). At the same time, careful monitoring of gene flow is needed to avoid disrupting local adaptation — particularly when managing introgression across steep environmental gradients (Sexton et al., 2014). Regular genetic assessments and trait performance evaluations will help guide outcomes toward both ecological compatibility and climate resilience.

In addition to in situ approaches, ex situ conservation strategies (e.g., gene and seed banks) are indispensable for preserving adaptive genetic material and supporting restoration under future climate scenarios (Sochacki et al., 2024). These resources can serve as reservoirs for targeted breeding to increase stress resistance and ensure adaptive potential across populations. Finally, successful conservation and breeding efforts must consider ecological interactions, including those with pollinators and competing species, which influence reproductive success and long-term population viability (Potts et al., 2010). Integrating these ecological dimensions is crucial for maximizing the effectiveness and sustainability of adaptive conservation frameworks.

4.5. Study limitations and further perspectives

It is important to note several limitations that could influence our interpretations. First, our transcriptome-based alignments, while effective for detecting coding variation, have limited sensitivity to low-expression genes and may miss genes expressed only in specific tissues or developmental stages. In addition, post-transcriptional RNA modifications and alternative splicing events can introduce alignment ambiguities, complicating the interpretation of sequence variants. Moreover, aligning reads to transcripts rather than the full genome can obscure gene family relationships and hinder the assignment of orthologous or paralogous sequences, especially in conifer genomes with large, repetitive gene families (Gagalova et al., 2022).

Second, because our analysis focused primarily on coding regions, it may have excluded important adaptive variation in non-coding regulatory elements (e.g., enhancers, promoters, and untranslated regions). Convergent or parallel adaptations frequently occur in non-coding sequences and can play important roles in evolutionary responses to environmental pressures (Novella et al., 2004; Nakayama and Makino, 2024). To address these limitations, we suggest that next steps would be to employ whole-genome sequencing approaches that include coding and non-coding regions and to assess the potential generality of adaptive mechanisms in spruce and other coniferous species.

A third limitation involves the geographic scope of our sampling. The complex terrain of the study region constrained accessibility and limited our ability to sample populations in certain areas. This was particularly for Picea meyeri populations. Nonetheless, our population structure analyses suggest that low sample size did not substantially bias our core findings (Fig. 2A). Even so, next steps would certainly be to sample from broader geographic coverage to improve the representativeness and robustness of population-level inferences.

Finally, while our analysis focused on species-specific patterns of adaptive introgression, we acknowledge that some adaptive alleles may represent reused ancestral polymorphisms rather than newly derived variants. Previous research in Picea (e.g., Feng et al., 2019) has documented widespread trans-specific polymorphisms, suggesting that standing genetic variation — persisting through speciation events — can serve as a reservoir for adaptation. Although we did not directly test for the presence of such trans-specific or reused alleles here, future analyses incorporating temporal allele frequency shifts or broader lineage comparisons may help elucidate the evolutionary origins of these adaptive loci.

Acknowledgements

This work was supported by the Project of Qinghai provincial central government guides local funds for science and technology development (2024ZY005).

CRediT authorship contribution statement

Shuo Feng: Writing – review & editing, Writing – original draft, Funding acquisition, Formal analysis. Haixia Ma: Writing – original draft, Software, Methodology, Formal analysis, Data curation. Yu Yin: Methodology, Formal analysis. Wei Wan: Software, Methodology, Formal analysis. Kangshan Mao: Writing – review & editing, Validation, Investigation. Dafu Ru: Writing – review & editing, Validation, Supervision, Investigation, Conceptualization.

Data availability

The sequencing data have been deposited to SRA database at NCBI under BioProject PRJNA846694.

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