1. Introduction

Elucidating the mechanisms by which natural selection maintains genetic divergence in the face of gene flow remains a key question in evolutionary biology. Given that variable environments promote ecological, physiological, and genetic differentiation across gradients (Johannesson et al., 2024), long-term divergent selection can subsequently shape distinct lineages by favoring genotypes that are locally adapted to their specific habitats (Qiu et al., 2011; Jiang et al., 2025). Under such scenarios, local adaptation reduces migrant fitness and establishes strong extrinsic barriers to gene flow (i.e., isolation-by-environment, IBE), which may ultimately lead to reproductive isolation and subsequent ecological speciation (White and Butlin, 2021; Borthakur et al., 2022; Li et al., 2024a). Prior to complete reproductive isolation, however, locally adapted lineages may hybridize, forming hybrid zones either by primary intergradation or following climate-driven secondary contact (Rolán-Alvarez et al., 1997; Abbott, 2017). These hybrid zones produce recombinant genotypes that facilitate the identification of loci exhibiting excess ancestry and potential selection signatures (Sung et al., 2018), thus serving as crucial natural laboratories for studying how evolutionary forces and ecological conditions maintain lineages distinctiveness despite ongoing hybridization (Abbott et al., 2013; Stankowski et al., 2021).

Subtropical China is acknowledged as one of the world’s biodiversity hotspots, harboring approximately one-third of China’s vascular plant diversity (Mi et al., 2021). Despite ongoing debate, this region is increasingly recognized as an evolutionary cradle for numerous plant species (López-Pujol et al., 2011; Zhu and Tan, 2024), a status attributed to its extensive forest habitats (i.e., 22–34°N, 98–123°E) with heterogeneous topographic and climatic conditions that promote local adaptation (Qian and Ricklefs, 2000; Qiu et al., 2011). Meanwhile, this distinct geological history (e.g., orogeny and/or climatic fluctuations) has also created well-known hybrid zones. For instance, recent admixture triggered by interglacial intervals gave rise to a hybrid zone in Quercus acutissima Carruth., but the western and eastern lineages of this species maintain genetic and phenotypic differentiation due to ecological adaptation (Zhang et al., 2018; Gao et al., 2021; Yuan et al., 2023). Similar patterns have been reported in multiple plant species, indicating that hybridization during species or lineage divergence (speciation with gene flow) in subtropical China may be more prevalent than previously recognized (Bock et al., 2023; Hu et al., 2025). These observations suggest that this region represents an ideal natural laboratory to investigate how genetic divergence is maintained despite ongoing gene flow. However, systematic genomic investigations of such hybrid zones remain scarce. Elucidating the processes of local adaptation and hybridization is thus critical for understanding the mechanisms underlying the origin of biodiversity in subtropical China (Meek et al., 2023; Hu et al., 2025).

Chinese walnut (Juglans cathayensis Dode) is an economically and ecologically significant forest tree widely distributed across the subtropical evergreen broad-leaved forest (EBLF) of China (Zhou et al., 2017) (Fig. 1A). The species comprises two morphologically and geographically distinct varieties: J. cathayensis var. formosana and J. cathayensis var. cathayensis. The former is found in eastern China, bearing smoother, bi-ridged nuts with only faint wrinkles and lacking spines or deep pits; whereas the latter distributes in central–western China, producing nuts with 6–8 prominent longitudinal ridges and sharp spines (Kuang and Lu, 1979) (Fig. 1B). Interestingly, previous phylogeographic analyses using both nuclear and chloroplast DNA datasets have consistently revealed an east–west differentiation (Bai et al., 2014, 2016; Geng et al., 2024). This divergence aligns with the biogeographic boundary between two major East Asian temperate forest subkingdoms: the Sino-Japanese Forest in eastern warm-to-cold temperate lowlands under Pacific monsoon influence, and the Sino-Himalayan Forest in western cold–dry regions dominated by the Indian monsoon (Qiu et al., 2011; Bai et al., 2014; Geng et al., 2024). These patterns hint that the two varieties of J. cathayensis might have evolved under distinct selective pressures associated with different environmental conditions. In addition, although pollen-mediated gene flow is prevalent in J. cathayensis, topographic fragmentation and climatic gradients across subtropical China likely prevent the homogenization of its two varieties, thus maintaining a narrow genetic mixing zone between them (Bai et al., 2014; Geng et al., 2024). However, due to the limited resolution and neutral evolution of molecular markers used in previous studies (Bai et al., 2014), how local adaptation maintains genetic differentiation in J. cathayensis despite ongoing hybridization has yet to be systematically investigated at the genomic level.

When hybrids occur naturally within wild populations, identifying outlier loci subject to spatially divergent selection is a key point (Barraclough, 2024; Wu et al., 2024). Recent improvements in whole-genome resequencing and statistical methods have rendered this exploration feasible (Filipe et al., 2022; Zhang et al., 2023). For example, genome scans for genetic differentiation (e.g., F_ST and D_XY; termed as “genomic islands”) (Bock et al., 2023) and for selection signatures (e.g., XP-CLR) are employed to detect loci linked to population-specific adaptation (Bourgeois and Warren, 2021; Feng et al., 2024). Moreover, genotype–environment associations (GEAs) analyses provide a robust approach for detecting adaptive loci linked to heterogeneous environmental landscapes (Forester et al., 2018). Although such analyses enhance our understanding of how ecological factors shape genetic differentiation (Yuan et al., 2025), they lack power to determine whether the identified loci exhibit greater resistance to gene flow in nature (Bock et al., 2023). In contrast, cline models (e.g., Bayesian genomic cline methods, BGC) can quantify gene flow patterns across hybrid zones and infer selection strength acting on specific loci (reviewed in Gompert et al., 2017). Two parameters, α and β outliers, are employed to characterize the deviation of introgression for individual loci from the genomic background and selection patterns within hybrid zones (Gompert and Buerkle, 2012). Integrating BGC with genome scans and GEA outliers provides a feasible approach to test whether loci involved in hybridization and/or introgression are associated with local adaptation (Liu et al., 2024). To our knowledge, however, such an approach has been rarely applied to tree species characterized by high outcrossing rates, extensive gene flow, and complex evolutionary histories (but see Menon et al., 2018; Guo et al., 2023a).

In the present study, we aimed to determine the mechanisms that maintain intraspecific divergence in J. cathayensis under persistent hybridization in subtropical China. For this purpose, we used population genetics, landscape genomics, and ecological niche modeling to evaluate how evolutionary history shaped genomic variation patterns between the two varieties of J. cathayensis and to identify the genomic basis of local adaptation in these varieties. Specifically, we examined if loci involved in hybridization are related to local adaptation. This study not only elucidates intraspecific diversification and hybridization in J. cathayensis, but also offers novel perspectives on the formation and maintenance of plant diversity in subtropical China.

2. Materials and methods 2.1. Sampling, sequencing and variant calling

We newly sequenced 30 individuals from 2021 to 2023 and downloaded 21 sequences from the SRA database (Zhang et al., 2019a, 2022). These individuals, from 28 populations, essentially cover the geographical range of Juglans cathayensis in subtropical China (Fig. 1A). Two J. regia L. individuals were included as outgroups (Table S1). All fresh leaf samples were dried with silica gel and stored at −20 ℃ in the Lushan Botanical Garden (Jiujiang, Jiangxi). Total genomic DNA was extracted from leaf tissue using a modified CTAB method (Doyle and Doyle, 1987). Paired-end sequencing libraries were constructed with an insert size of 400 bp and whole-genome resequencing (> 15× coverage) was performed on the MGI DNBSEQ-T7 platform.

Raw sequence reads were cleaned using fastp v.0.12.4 (Chen et al., 2018) to remove adapters and low-quality bases. The resultant high-quality reads were aligned to the J. mandshurica Maxim. reference genome (a synonym of J. cathayensis native to northern China) (Li et al., 2022; Geng et al., 2024) with BWA v.0.7.17 (Li and Durbin, 2009). The aligned reads were sorted using SAMtools v.1.10 (Li et al., 2009) and PCR duplicates were removed using Picard v.3.0 (https://github.com/broadinstitute/picard). Raw SNPs were identified with HaplotypeCaller, CombineGVCFs and GenotypeGVCFs of GATK v.4.2.3.0 (McKenna et al., 2010), producing a merged VCF file. Subsequently, SNPs were filtered based on the following criteria: QD < 2.0; MQ < 40.0; FS > 60.0; SOR >3.0; MQRankSum < −12.5; ReadPosRankSum < −8.0. Finally, all variations were further filtered using VCFtools v.0.1.16 (Danecek et al., 2011) with the following parameters: –max-missing 0.8 –maf 0.05 –min-alleles 2 –max-alleles 2 –maxDP 40 –minQ 3.

2.2. Population structure and genetic diversity

LD-pruning was performed using PLINK v.1.90b6.24 (Purcell et al., 2007) with the settings ‘indep-pairwise 50 10 0.2’, resulting in 578,735 independent SNPs for subsequent population structure analysis. Population structure was inferred with ADMIXTURE v.1.3.0 (Alexander and Lange, 2011) by testing a range of K-values from 2 to 10. PCA was performed using the pca function implemented in PLINK v.1.90b6.24. In addition, a maximum likelihood (ML) phylogenetic tree was constructed with FastTree v.2.1.10 (Price et al., 2010) and displayed using iTOL v.6.0 (Letunic and Bork, 2024). Finally, we carried out identity-by-descent (IBD) analysis in BEAGLE v.5.2 (Browning and Browning, 2007) with the following settings: window = 100,000; overlap = 10,000; ibdtrim = 100; ibdlod = 10; ibd = true.

Nucleotide diversity (π) and pairwise genetic differentiation (F_ST) were estimated for each lineage pair with VCFtools v.0.1.16, applying 20 kb non-overlapping sliding windows.

2.3. Hybridization detection

We employed multiple methods to detect hybridization within the Admixed lineage (see below) at the population level, as ADMIXTURE and/or PCA alone cannot distinguish admixture from isolation-by-distance (IBD) (Wiens and Colella, 2025). To detect gene flow between parents and admixed lineages, we computed Patterson’s D-statistics (ABBA-BABA) using the Dtrios function in Dsuite v.0.5.r44 with default parameters, relying on its automated inference of population relationships (Malinsky et al., 2021). Strong gene flow was detected between two populations (P3 and P2) based on a predetermined four-taxon tree topology (((P1, P2), P3), O) (Z > 3 rejects null hypothesis). Notably, D-statistics exhibited higher false positives under lineage-specific rate variation and lower power with an increasing number of hybridization events, as multiple events within a small subset of taxa can mask each other (Frankel and Ané, 2023). We therefore employed TreeMix v.1.13 (Pickrell and Pritchard, 2012) to cross-validate the D-statistics results, running 10 iterations for each of the 0–10 migration events (m). All the above analyses used two J. regia accessions as outgroups. Finally, the R package INTROGRESS (Gompert and Alex Buerkle, 2010) was used to examine hybrid index (Hi) for each individual and visualized through a triangle plot. The Hi ranges from 0 (the East lineage) to 1 (the West lineage), with intermediate values indicating varying degrees of hybridization.

2.4. Demographic history inference

PopLDdecay v.3.4.2 (Zhang et al., 2019b) was used to measure and compare linkage disequilibrium (LD) patterns for each lineage. The demographic history of divergent lineages was estimated by three complementary methods. For PSMC analysis (Li and Durbin, 2011), one representative individual per population was selected and analyzed using the following parameters: -N30 -t15 -r4 -p ‘4 + 25 × 2 + 4+6’. The results of 100 bootstraps were combined to plot. We then inferred recent effective population size (N_e) history of all individuals using SMC++ v.1.11.1 (Terhorst et al., 2017) with the following parameters: –spline cubic –knots 20 –timepoints 1 200,000.

To infer historical split times and migrations rates, we generated a two-dimensional joint site frequency spectrum using easySFS.py (Gutenkunst et al., 2009) based on four-fold degenerate sites (4DTVs) (He et al., 2013) and then analyzed it in fastsimcoal2 (Excoffier et al., 2013). The best-fitting demographic scenario was identified among the four classical speciation models (Nielsen and Wakeley, 2001) as well as twelve derived variants exhibiting distinct population dynamics (Fig. S1). Parameter estimation and the best-fitting model selection were carried out as described in Li et al. (2024a). All demographic history analyses used a mutation rate of 2.06 × 10⁻⁹ per site per year and a generation time of 30 years following Zhang et al. (2022). Because recent admixture events can inflate estimates of ancestry proportions and gene flow (Meier et al., 2017), admixed individuals identified by ADMIXTURE (0.01 < Q < 0.99) and INTROGRESS (0.01 < Hi < 0.99) were excluded before demographic analyses.

2.5. Niche modeling and ecological divergence

The potential distributions of J. cathayensis were projected under the current, the Last Glacial Maximum (LGM), the Mid-Holocene (MH), and the Last Interglacial (LIG) climatic periods. After removing non-natural populations and spatial thinning to a 5-km distance (to mitigate autocorrelation), 44 occurrence records sourced from the Chinese Virtual Herbarium (https://www.cvh.ac.cn/) and this study (Table S2) were retained. Paleoclimate simulations for the MH and LGM periods were based on the MIROC and CCSM models. Nineteen bioclimatic variables were downloaded from the WorldClim website (http://www.worldclim.org/) (Table S3), with a 2.5 arc minute resolution for each environmental layer. Only climatic variables with a pairwise Pearson correlation coefficient of |r| < 0.7 were retained for projections of species distribution (Bio1, annual mean temperature; Bio2, mean diurnal range; Bio3, Isothermality; Bio7, temperature annual range; Bio12, annual precipitation; Bio15, precipitation seasonality). The climate niches of J. cathayensis were modeled using MAXENT v.3.4.1 (Phillips et al., 2006) with default parameters. An AUC value exceeding 0.9 was set as the threshold for acceptable model discrimination.

For ecological divergence analysis, we used the geographical locations of the 28 populations from our resequencing data. Based on the lineage assignments from the ADMIXTURE analysis (see Fig. 1A), we calculated Schoener’s D and standardized Hellinger distance (I) between lineages using niche overlap and identity tests in ENMTools v.1.4.4 (Warren et al., 2008, 2010).

2.6. Identifying genomic regions of differentiation and selection

To avoid confounding effects of ancestral admixture on inter-lineages divergence and selection analyses, admixed individuals were excluded prior to conducting these tests. We calculated the genome-wide distribution of F_ST, π, and Tajima’s D with VCFtools v.0.1.16, while absolute genetic divergence (D_XY) was computed with the Python script popgenWindows.py (https://github.com/simonhmartin/genomics_general). The population-scaled recombination rate (Rho = 4N_er) was calculated for each lineage using FastEPRR (Gao et al., 2016). All parameters were estimated in 20-kb non-overlapping sliding windows (excluding those with < 10 SNPs).

We followed Ma et al. (2018) to obtain a null F_ST distribution via 5,000,000 times genome-wide SNP permutations for each of the 20-kb non-overlapping windows. Subsequently, window-based F_ST estimates were compared to null distributions to obtain P-values and adjusted using the false discovery rate (FDR). Windows residing in the top 5% of the empirical F_ST distribution and exhibiting an FDR < 0.01 were identified as outlier genomic windows. Adjacent outlier windows were further merged into a larger divergent regions (i.e., genomic islands). In addition, we calculated the cross-population composite likelihood ratio (XP-CLR) scores (Chen et al., 2010) to detect potential selective signals in each lineage with 20-kb non-overlapping windows. The top 5% windows identified by both F_ST and XP-CLR methods were considered as putative positively selected genes (PSGs) for each lineage. Finally, the differences of D_XY, π, Tajima’s D, Rho, and XP-CLR scores between genomic islands and the rest of genome (genomic background) were examined using the Mann–Whitney U test. Pairwise correlations of these parameters with F_ST were then evaluated through Spearman’s correlation analysis.

To further validate sample size adequacy, we performed a saturation analysis via random subsampling, retaining only one individual per population (East: 8; West: 6) to recalculate genetic divergence (F_ST and D_XY), selection signatures (Tajima’s D and XP-CLR), and genetic diversity (π) across non-overlapping 20-kb windows. The correlation between these genetic parameters derived from the full and subsampled datasets was then evaluated per window using Spearman correlation analysis.

2.7. Genotype-environment association analyses

To identify candidate loci associated with multivariate environmental axes (Capblancq and Forester, 2021), we performed redundancy analysis (RDA) in the R package vegan (Dixon, 2003) using the six uncorrelated environmental variables selected for niche modeling. Outlier loci were identified using a standard deviation cut-off of 3 along one or more RDA axes. As an alternative approach, we conducted GEAs using a univariate latent-factor linear mixed model (LFMM) implemented in the R package LEA v.3.3.2 (Frichot and François, 2015), specifying two latent factors to control for population structure. Significant associations between allele frequencies and the 19 bioclimatic variables were identified with a 5% FDR threshold. Genes concurrently identified by RDA and LFMM were designated as “core adaptive genes” underlying local adaptation.

To assess the effects of geography and environment on genetic variation for both neutral (LD-pruned SNPs) and adaptive loci, we first examined the correlation between environmental (Euclidean distance) and geographical distance. Then, Partial Mantel tests were applied to investigate the correlation between genetic distance and geographic (isolation-by-distance, IBD; after controlling environment influence) and environmental (isolation-by-environment, IBE; after controlling geography influence) distance with significance determined using 999 permutations in the R package vegan. Population genetic distances were represented by a pairwise F_ST/(1-F_ST) matrix constructed using the R package hierfstat v.0.5.11 (Goudet, 2005).

2.8. Gene ontology (GO) analyses

For GO enrichment analysis, positively selected and core adaptive genes were first extracted from the J. mandshurica genome and functionally annotated using eggNOG-mapper (Cantalapiedra et al., 2021). These annotations were then analyzed with the R package ClusterProfiler (Wu et al., 2021), using an FDR cutoff of 5% for significance.

2.9. Bayesian genomic cline analysis

To quantify genome-wide variation in introgression among admixed individuals from the inter-lineage hybrid zone, we performed the Bayesian Genomic Cline (BGC) analysis using the program bgc (Gompert and Buerkle, 2012). Briefly, BGC quantifies the hybrid index (Φ) and infers the locus-specific ancestry origin from a parental lineage (here, West or East lineage). The genomic cline parameters (α and β) are then used to estimate the posterior probability of ancestry at each locus within the admixed lineage (Gompert and Buerkle, 2011). The parameter α denotes the position of the cline center and infers the direction of introgression as excess ancestry from one reference lineage (positive α) or the other reference lineage (negative α). The parameter β denotes the slope of the cline, with positive and negtaive β outliers indicating reproductive isolation and adaptive introgression, respectively (Gompert and Buerkle, 2012). The analysis consisted of two independent MCMC runs of 50,000 steps each. After a burn-in of 5,000 steps, samples were recorded every 10th iteration. Outputs of the two runs combined after being visually inspected for convergence, and the results were plotted using the R package ClineHelpR (Martin et al., 2021).

Subsequently, we quantified the overlap by mapping both positive selection genes and the core adaptive genes to the candidate introgressed regions detected through BGC analysis. These observed overlap proportions were compared to those expected from the random expectation, generated by 1000 permutations of randomly sampled gene sets using R v.4.1.2 software (R CoreTeam, 2021). These analyses provide a fruitful way to test whether loci associated with local adaptation (potential reproductive isolation barriers) also exhibit patterns of excess ancestry (Sung et al., 2018; Liu et al., 2024).

3. Results 3.1. Population structure, genetic diversity and hybridization

Whole-genome resequencing of 49 samples from 28 populations yielded 513 Gb of data, with an average sequencing depth of 30.95×. Using the chromosome-level Juglans mandshurica reference genome (Li et al., 2022), the average mapping rate of raw reads was 97.22% (Table S1). After filtering, 5,910,597 high-quality SNPs were retained for subsequent analyses. Although sample sizes per site were limited, individuals were grouped into genetic clusters (see below) for population genomic analyses. The resulting sample sizes per analytical unit, combined with genome-wide high-depth sequencing and millions of SNPs, provided sufficient statistical power to estimate allele frequencies and detect selection signatures (Willing et al., 2012; Nazareno et al., 2017).

ADMIXTURE analysis of the LD-pruned SNPs divided all individuals into two genetic clusters when K = 2: the eastern cluster corresponded to J. cathayensis var. formosana (East lineage, n = 17), while the western cluster comprised both genetically pure J. cathayensis var. cathayensis (West lineage, n = 13) and admixed individuals (Admixed lineage, n = 19) (Fig. 1 and Fig. S2; Table S4). The relationships of these lineages were further supported by principal component analysis (PCA) (Fig. 1C) and the ML phylogenetic tree (except for four admixed individuals nested within the West lineage) (Fig. S3), where the Admixed lineage exhibited an intermediate position between the West and East lineages. Moreover, the much more extensively shared IBD (identity-by-descent) blocks between the Admixed lineage and the East/West lineage indicated recent ancestral admixture of the two varieties (Fig. 1D). Nucleotide diversity (π) was slightly lower in the West lineage (2.52 × 10⁻³) than in the East (2.67 × 10⁻³) and Admixed lineages (2.65 × 10⁻³). Pairwise genetic differentiation (F_ST) were notably lower between the Admixed lineage and East (0.032) or West lineage (0.013) compared to that between the East and West lineages (0.054) (Fig. 1E).

D-statistics analysis revealed only one predominant topology and exhibited a high level of positive Z-score (8.70), indicating close relatedness between the West and Admixed lineages and significant introgression from the East lineage into the Admixed lineage (Fig. 1F). The result of TreeMix confirmed the same phylogenetic relationships and introgression patterns (Fig. 1G). Hybrid index (Hi) estimates for the Admixed lineage ranged from 0.3 to 0.7 (Fig. 1H).

3.2. Demographic history

The decay of linkage disequilibrium showed a similar pattern in the East and West lineages (Fig. S4). PSMC analysis indicated both lineages reached a high effective size approximately 4 million years ago (Ma), followed by a prolonged decline in most individuals (Fig. 2A). Notably, subsets of individuals from each lineage showed population expansions around 0.025 Ma, which temporally coincided with the LGM (Fig. 2A). SMC++ analysis showed that the East and West lineages experienced population expansions at about 0.01 Ma (Fig. S5).

The best-fitting demographic scenario based on Akaike Information Criterion (lowest AIC = 28904.544, w_i = 1; Table S5) indicated that the ancestral population (N_e = 670,174) diverged into the East and West lineages around 4 Ma, with continuous bidirectional gene flow followed by isolation at about 0.6 Ma (Fig. 2B and Table S6). Migration rate from the West lineage to the East (0.0035) was higher than vice versa (0.0012). Additionally, both lineages underwent substantial shifts in population size around 2.4 Ma, with the East lineage expanding dramatically (N_e from 2686 to 366,087), whereas the West lineage experienced a moderate contraction (N_e from 648,783 to 313,505) (Fig. 2B).

3.3. Species distribution models and niche divergence

In the MAXENT analysis, six uncorrelated bioclimatic variables were retained (Fig. S6) and all niche models exhibited high discriminative power, with AUC values exceeding 0.9. Distribution modeling of J. cathayensis showed repeated range fluctuations since the LIG (Fig. 2C–F and Fig. S7), suggesting periodic isolation–reconnection cycles between lineages. The ENMTools analysis showed that the observed niche similarity (D and I) differed significantly from null distributions between the West and East lineages, but not between the Admixed lineage and either pure lineage (Fig. 2G–I).

3.4. Heterogeneous genomic divergence and positive selection analysis

Of 1289 outlier genomic windows, 689 genomic islands were identified that showed elevated differentiation (mean F_ST = 0.2928 ± 0.0957) relative to the genomic background (mean F_ST = 0.0563 ± 0.0467) (Fig. 3A and Fig. S8; Table S7). These islands displayed strongly elevated D_XY and XP-CLR scores, but significantly reduced π, Rho and Tajima’s D (Fig. 3B and Fig. S8). In addition, F_ST showed strong positive correlations with D_XY and XP-CLR, and strong negative correlation with π, Rho, and Tajima’s D (Fig. 3C). Subsampling saturation analysis demonstrated significant positive correlation between population genetic parameters estimated from the full and subsampled datasets across all 20-kb non-overlapping windows (Fig. S9), indicating the robustness and stability of these parameter patterns even with reduced sampling.

F_ST and XP-CLR tests simultaneously identified a total of 1146 candidate genes under positive selection. Of these, 783 genes were detected in the East lineage and 766 genes in the West lineage (Fig. S10; Tables S8 and S9), with 403 PSGs shared between lineages. These PSGs are primarily involved in adaptive responses to diverse abiotic and biotic stressors, including cold or heat tolerance (e.g., SIZ1 and GA2ox3) (Miura and Nozawa, 2014; Li et al., 2019), salt resistance (e.g., SOS1 and UPL1) (Yang and Guo, 2018; Marczak et al., 2025), and disease resistance (NAC4) (Yuan et al., 2019). Some genes play diverse roles in plant growth and development (e.g., MOR1 and RAPTOR1) (McCready et al., 2020; Liu and Yu, 2023). Additionally, we identified two genes associated with flowering time regulation (EBF2 and ATR2) (Qureshi et al., 2011; Li et al., 2015), which may play a crucial role in the prezygotic isolation between the two lineages. GO enrichment analysis revealed that the PSGs in each lineage were enriched for terms related to various functions, including organ development, environmental adaptation, reproductive isolation, and others (Figs. S11 and S12; Tables S10 and S11).

3.5. Genomic variants associated with environmental adaptation

RDA analysis with six uncorrelated variables identified 214,681 candidate adaptive SNPs (Fig. 4A). RDA1 axes distinguished three lineages (explaining 25.48% of the variance), with three significant contributors to intraspecific divergence: Isothermality, annual precipitation, and precipitation seasonality (Fig. 4B). We then identified a total of 73,065 SNPs significantly associated with at least one environmental variable by LFMM, of which 3652 SNPs (in 307 genes) overlapped with those identified by RDA (Fig. 4A). These shared genes were considered as “core adaptive genes” for environmental adaptation. GO enrichment analysis showed that these core adaptive genes are primarily associated with activity and metabolic processes, suggesting their role in local adaptation in J. cathayensis (Fig. S13 and Table S12). Due to the strong autocorrelation between environmental and geographical distances (Fig. S14), Partial Mantel tests were used to assess patterns of isolation-by-distance (IBD) and isolation-by-environment (IBE) for the neutral and adaptive SNPs, respectively (Fig. 4C). The results showed robust IBE signals but an absence of significant IBD in both neutral and adaptive variants.

We selected Bio7 (temperature-related) and Bio15 (precipitation-related) to illustrate how significant genes and their Arabidopsis thaliana homologs are associated with representative climate variables (Fig. 4D; Table S13). For example, UGT85A1, associated with Bio7, is preferentially expressed in the leaf epidermis and guard cells in A. thaliana and plays a role in drought tolerance (Rehman et al., 2018). CKX3 regulates the activity of reproductive meristems, flower organ size, and seed yield (Bartrina et al., 2011). AGO2 was found to be involved in salinity stress resistance (Wang et al., 2019a). G-TMT is a key gene in the vitamin E biosynthesis pathway and regulates the spectrum of accumulated tocopherols in plants (Koch et al., 2003). Additionally, EIL1, associated with Bio15, is involved in the ethylene signaling pathway and plays a crucial role in the anther development (Zhu et al., 2022). AGD4 is involved in vesicle trafficking and plays a crucial role in the development of pollen (Yang, 2018). ATG9 is known to be important for autophagy (Zhuang et al., 2017).

3.6. Introgression loci associated with local adaptation

The BGC-estimated hybrid index ranged from 0.3 to 0.7 (Fig. 5A), consistent with results from INTROGRESS. Genomic admixture patterns were heterogeneous across the genome, with the cline center parameter α being more variable (−1.171–1.049) than the cline slope parameter β (−0.006–0.006) (Table S14). In addition, we detected 6677 outlier loci exhibiting significant locus-specific admixture that deviated from the genome-wide average (i.e., excess ancestry), including 3004 positive α outliers and 3673 negative α outliers. The positive and negative α outliers represented excess ancestry in the West and East lineages, respectively. We did not identify any loci that were β outliers or had excess ancestry indicated by β (Table S14).

We subsequently examined the association between these 6,677 α outliers (distributed across 853 genes) and both PSGs and core adaptive genes. Of these, only 6.02% of genes (69) overlapped with the PSGs (Fig. 5B; Tables S8 and S9), and 5.54% of genes (17) overlapped with the core adaptive genes (Fig. 5C; Table S13). The observed overlaps were significantly lower than random expectations (permutation test, P < 0.001; Fig. S15). However, none of these shared genes were identified as β outliers, suggesting an absence of adaptive introgression (Guo et al., 2023a).

4. Discussion 4.1. Demographic history and hybridization between the two varieties of Juglans cathayensis

In subtropical China, historical climatic changes have profoundly shaped the population genetic structure and demographic history of extant species (Hewitt, 2000; Qiu et al., 2011). In this study, range-wide whole-genome resequencing of Juglans cathayensis revealed two genetic clusters corresponding to its two varieties, which were further subdivided into three lineages: an East lineage of J. cathayensis var. formosana individuals, and West and Admixed lineages of J. cathayensis var. cathayensis individuals (Fig. 1). This result is consistent with previous studies using microsatellite and nuclear genomic data (Bai et al., 2014; Geng et al., 2024). We further estimated the lineage divergence time and demographic history (Fig. 2). PSMC analysis indicated a persistent decline in effective population size from approximately 4 Ma, likely promoting lineage divergence. This timing was corroborated by fastsimcoal2 analysis and published chloroplast haplotype data (Geng et al., 2024). Such a pattern may have resulted from the intensification of East Asia monsoon during the Pliocene, which amplified climatic differences between eastern and western China and potentially drove genetic divergence through ecological selection (Wang et al., 2019b; Yuan et al., 2023). Alternatively, the sharp decline in effective population size may have reduced genetic diversity while enhancing the effects of genetic drift and inbreeding, further facilitating lineage divergence in J. cathayensis (Jing et al., 2025). Notably, similar east–west divergences have been documented in multiple plant species across subtropical China, such as Gaultheria crenulata Kurz (Li et al., 2024b), Pseudotaxus chienii (W.C. Cheng) W.C. Cheng (Kou et al., 2020), Castanopsis carlesii (Hemsl.) Hayata (Sun et al., 2016), and C. fargesii Franch (Sun et al., 2014), suggesting a common biogeographic response to regional historical climatic forces.

Based on fastsimcoal2 simulations, bidirectional gene flow persisted between the East and West lineages from their initial divergence until about 0.6 Ma (Fig. 2B), a period that coincides with increased amplitude of glacial-interglacial variations and intensified monsoon cyclicity during the mid-Pleistocene transition (MPT, 0.6–1.25 Ma) (Ao et al., 2023). Climatic divergence between eastern and western China likely reinforced genetic differentiation through long-term local adaptation, ultimately leading to complete cessation of gene flow (Qiu et al., 2011). These possibilities are supported by niche identity tests and significant IBE signals (Fig. 3C). In addition, 19 individuals from intermediate regions and niche characteristics exhibited admixed ancestry, indicating recent hybridization and the formation of a hybrid zone in the ecotone between parental habitats (Fig. 1, Fig. 2). Demographic and niche modeling have suggested that these hybridization events may have been triggered by range fluctuations (e.g., vertical habitat shifts in response to Quaternary climate oscillations) between the East and West lineages (An et al., 2001; Dang et al., 2025) (Fig. 2), and/or long-distance pollen-mediated gene flow (Bai et al., 2014). It should be noted that the admixed individuals were excluded from our fastsimcoal2 analysis to focus on the East-West lineage relationship. The potential impact of this exclusion on inferred gene flow strength and divergence time warrants future investigation. Despite this caveat, our findings indicate that pronounced Pliocene–Holocene climatic fluctuations likely contributed to the emergence of the three extant lineages in J. cathayensis (Ren et al., 2024), with both historical gene flow and contemporary hybridization occurring throughout its divergence history.

4.2. Genomic islands of divergence and local adaptation

Comparisons of genome-wide variation in recently diverged species or lineages often reveal a heterogeneous landscape of genetic divergence marked by genomic islands (Han et al., 2017). Although the overall pairwise genetic differentiation (F_ST) between the East and West lineages of J. cathayensis was relatively low (0.054), we identified 689 genomic islands exhibiting elevated differentiation (Table S7). The reductions in nucleotide diversity (π), recombination rate (Rho), and Tajima’s D within these regions, coupled with their negative correlation to genetic divergence (Fig. 3), implicate linked selection as a key process in generating genomic islands (Noor and Bennett, 2009; Ke et al., 2022; Cao et al., 2023; Shi et al., 2024). More specifically, the elevated XP-CLR scores within islands and their strong positive association with F_ST further indicates significant effects of positive selection (i.e., divergence hitchhiking) in shaping the genomic divergence landscape of J. cathayensis during its ecological adaptation (Xu et al., 2024).

Consistent with previous studies on interspecific and intraspecific divergence (Ma et al., 2018; Hu et al., 2022; Guo et al., 2023b; Xu et al., 2024), we found that absolute genetic divergence (D_XY) increased markedly in the genomic islands and there was a positive correlation between D_XY and F_ST (Fig. 3). Generally, the elevated values of both F_ST and D_XY in these highly divergent regions likely result from divergent sorting of ancient polymorphisms and/or divergent selection with gene flow (Ravinet et al., 2017; Bock et al., 2023). In our study, the sharing of 403 PSGs scattered across the genome between lineages of J. cathayensis (Tables S8 and S9) indicates that these genes likely originated from divergent selection of the ancient haplotypes in the ancestral populations (see also Hu et al., 2022). Furthremore, demographic history and hybridization analyses suggest that both historical gene flow and a hybrid zone existed between the East and West lineages (Fig. 1, Fig. 2). The observed reduction in recombination rate within these divergent regions further indicates that genetic barriers to gene flow arise from both recombination rate variation and the presence of barrier loci in low-recombination genomic regions (Xu et al., 2024) (Fig. 3B). Therefore, gene flow may also have contributed to the formation of these genomic islands (Guo et al., 2023b).

The evergreen broad-leaved forest (EBLF) of subtropical China comprises two major subregions: the eastern moist EBLF, under the Pacific monsoon, occupies warm-to-cold temperate hills and lowlands (< 500 m) with 1000–2000 mm annual rainfall. In contrast, the western semi-moist EBLF is dominated by the Indian monsoon, featuring colder and drier conditions on plateaus and basins (>1000 m) with 900–1200 mm of annual precipitation (Wang et al., 2007; Qiu et al., 2011; Li et al., 2024b). In response to these contrasting environments, the eastern and western varieties of J. cathayensis likely experienced divergent selection, leading to adaptive divergence and consequent allele frequency shifts at multiple loci (Zou et al., 2024). Using F_ST and XP-CLR methods, we identified a total of 1146 protein-coding genes that have undergone positive selection. These genes are implicated in ecological adaptation, organ development, and flowering time, indicating that they play critical roles in reproductive isolation between lineages (Hu et al., 2022) (Tables S8–S11). Moreover, the 307 core adaptive genes identified through GEAs showed functional enrichment for environmental adaptation (Tables S12 and S13). Partial Mantel tests further indicate that genetic divergence among populations is shaped more strongly by environmental factors than by geographic distance. Taken together, our results demonstrate that these environmental adaptation genes have played pivotal roles in the genetic differentiation of J. cathayensis across subtropical China, which may underlie the observed ecological and phenotypical divergence among lineages (Gao et al., 2021; Liu et al., 2023). Notably, most core adaptive genes were not identified as PSGs, consistent with the polygenic barriers hypothesis via the accumulative addition of small individual effects (Hu et al., 2022; Sang et al., 2022; Xu et al., 2024).

Climate change is reducing crop yields and shifting species distributions, making it critical to identify adaptive genes that can enhance crop resilience, ensure species survival, and safeguard global food security (Campbell et al., 2025). Beyond identifying genomic signatures of adaptation, our study pinpointed specific candidate genes with direct relevance to environmental stress resistance, including SIZ1, GA2ox3, SOS1, UPL1, NAC4, UGT85A1, and AGO2 (Miura and Nozawa, 2014; Rehman et al., 2018; Yang and Guo, 2018; Li et al., 2019; Wang et al., 2019a; Yuan et al., 2019; Marczak et al., 2025). More importantly, two temperature-associated genes identified as adaptive to J. cathayensis are homologous to Arabidopsis genes (CKX3 and G-TMT) that play roles in seed yield and tocopherol accumulation (Koch et al., 2003; Bartrina et al., 2011). As a wild relative of cultivated walnut, J. cathayensis harbors these candidate adaptive genes that represent promising targets for walnut marker-assisted breeding. Its varieties thus provide an invaluable reservoir of adaptive traits, making them valuable genetic assets for enhancing cultivated walnut’s resilience to adverse climates (Zhou et al., 2017; Bernard et al., 2018; Campbell et al., 2025). Further functional validation is required to confirm the roles of these candidate genes.

4.3. Limited introgressed loci associated with local adaptation in the hybrid zone

When hybrid fitness involves multiple loci with different impacts, those under strong selection are likely to be significant β outliers, whereas weakly selected loci will show excess ancestry based on α outliers (Gompert and Buerkle, 2012; Guo et al., 2023a). In this study, Bayesian genomic cline (BGC) analysis identified 6,677 α outlier loci showing pervasive introgression in the Admixed lineage of J. cathayensis (Table S14), implying moderate divergent natural selection within the hybrid zone (Menon et al., 2018). Previous studies have shown that adaptive alleles transferred across species boundaries through hybridization and introgression occur in various taxa, including herbs (Sung et al., 2018), trees (Guo et al., 2023a), rodents (Horníková et al., 2024), fishes (Kakioka et al., 2021) and zooparasites (Fairfax et al., 2022), potentially enhancing environmental adaptation in the recipient lineages. However, our study detected thousands of α outliers but no β outliers (Fig. 5A), indicating the absence of strong adaptive introgression between the two J. cathayensis varieties (Guo et al., 2023a). This pattern raises concerns about the adaptive evolutionary potential of the Chinese walnut. Future studies should integrate genetic offset and fitness-related traits to evaluate population maladaptation in J. cathayensis varieties, which could help predict future habitat suitability and inform conservation strategies. The lack of β outliers further suggests that this viable hybrid zone is likely maintained primarily by extrinsic ecological factors, rather than pervasive intrinsic incompatibilities (Fouet et al., 2017; Ryan et al., 2017), a phenomenon also observed in other tree species such as white pine (Menon et al., 2018).

Hybrid zones may become isolated from parental lineages, collapse, expand, or stabilize, with their dynamics depending on the interplay between dispersal and selection pressures (Gompert et al., 2017; Stankowski et al., 2021). Thus, studying hybrid zones between intraspecific lineages can illuminate how natural selection maintains divergence by identifying which loci or phenotypic traits exhibit restricted gene flow relative to genome-wide expectations during the process of speciation (Wait et al., 2025). Given that substantial pollen-mediated gene flow may prevent genetic divergence in J. cathayensis populations, Bai et al. (2014) concluded that the lineages of Chinese walnut are not at an incipient stage of allopatric speciation. Despite this, only a small portion of α outliers overlapped with genes related to local adaptation in the present study (Fig. 5B and C), indicating that the genomes of these speciating lineages (at the so-called “semi-isolated species” stage) are semipermeable (Monnet et al., 2025), with gene flow between lineages varying among genomic regions, i.e., reduced in genomic regions containing genes associated with local adaptation, but increased in regions less constrained by selective forces (e.g., neutral loci unlinked to those selected loci) (White and Butlin, 2021; Freedman et al., 2023). As speciation proceeds, genomic regions with reduced gene flow are expected to expand through a process called divergence hitchhiking, leading to greater overall reproductive barriers (Bock et al., 2023). Therefore, we can expect that long-term local environmental adaptation through habitat preference would further accumulate genetic divergence and greatly limit hybridization between the varieties of J. cathayensis, ultimately culminating in complete reproductive isolation and speciation (Dai et al., 2024; Hou et al., 2025). Considering the absence of obvious geographic separation by physical barriers between populations (Fig. 3C), we propose that these varieties represent an incipient stage of ecological speciation or parapatric speciation involving extensive hybridization (Coyne and Orr, 2004; Menon et al., 2018), a process that may contribute to the formation and maintenance of biodiversity in subtropical China as a whole.

5. Conclusion

Based on the integration of whole-genome resequencing data and ecological niche modeling, our findings reveal that Chinese walnut comprises two major genetic lineages (East and West), with an admixed lineage forming a hybrid zone in the ecotone. Following initial divergence in the Middle Pliocene, the East and West lineages maintained sustained gene flow until the mid-Pleistocene, a process likely mediated by long-term adaptation to distinct regional environments. Landscape genomics analyses further confirmed this adaptive divergence, identifying positively selected and environment-associated loci that have facilitated intraspecific differentiation. Notably, Bayesian genomic cline analysis revealed that these adaptive loci exhibit restricted introgression in the hybrid zone, demonstrating that natural selection serves to maintain divergence in locally adapted genomic regions. Although this study is limited by a relatively small sample size from each population, our findings provide important insights into how local adaptation sustains genetic divergence in the presence of gene flow. Future research should expand individual sampling and incorporate additional genomic variants (e.g., chromosomal inversions) to improve the reliability and generality of these inferences, ultimately advancing our understanding of the evolutionary mechanisms shaping plant diversity in subtropical China.

Acknowledgements

We gratefully thank Dr. Rengang Zhang from the Kunming Institute of Botany, Chinese Academy of Sciences, for his valuable assistance with data analysis. We are grateful to the handling editor and the two anonymous reviewers for their constructive comments and suggestions on our manuscript. This work was financially supported by the National Natural Science Foundation of China (grant numbers 32460063 and 41961009), the Innovation Leading Talent Program in Jiangxi Province (JXSQ2023101107), the National Key R&D Program of China (2024YFF1307405) and the Key R&D Program of Jiangxi Province (20252BCF320039).

CRediT authorship contribution statement

Wan Hu: Methodology, Formal analysis, Visualization, Writing–original draft, Funding acquisition. Caixia Han: Methodology, Formal analysis, Investigation. Li Li: Methodology, Investigation, Data curation. Chunhua Zang: Methodology, Investigation, Formal analysis. Runzi Li: Methodology, Formal analysis. Hong Nie: Investigation, Data curation. Jianbo Nie: Investigation, Data curation. Yong Shi: Conceptualization. Chen Feng: Conceptualization, Funding acquisition. Jie Zhang: Conceptualization. Yixuan Kou: Writing–review & editing. Zhiyong Zhang: Writing–review & editing. Dengmei Fan: Conceptualization, Writing–review & editing. Danqi Li: Conceptualization, Investigation, Data curation, Supervision, Writing–review & editing, Funding acquisition.

Data availability

All re-sequencing data of population have been deposited in Genome Sequence Archive (GSA) of National Genomics Data Center (NGDC) under project PRJCA048367.

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