1. Introduction

High-elevation ecosystems, characterized by steep temperature gradients, intense ultraviolet radiation, seasonal water stress, and low oxygen availability, pose significant adaptive challenges to plants (Liu et al., 2014; Jansen, 2023). These extreme conditions drive distinctive physiological and morphological adaptations (Zandalinas et al., 2018), making alpine taxa ideal models for studying rapid evolution and evolutionary resilience (Hao and Lei, 2022). The Himalaya and Hengduan Mountains (HHM) region exemplifies these dynamics, functioning as a global biodiversity hotspot (Myers et al., 2000; Mittermeier et al., 2011) shaped by sharp ecological gradients and glacial–interglacial cycles that promote ecological divergence and local adaptation (Liu et al., 2014; Favre et al., 2015; Zhang et al., 2023). Simultaneously, adaptive capacity and evolutionary persistence have emerged as a key framework for anticipating how alpine taxa respond to rapid change. The likelihood of persistence depends on standing genetic variation, the strength and direction of selection, dispersal/connectivity, and life history (Bell and Gonzalez, 2011). While high-mountain systems remain understudied in terms of their underlying genomic signals (Carlson et al., 2014; Capblancq and Forester, 2021), the genomic mechanisms driving environmental adaptation in the alpine flora of the highly diverse HHM region also remain largely unclear (Stobdan et al., 2008; Carlson et al., 2014; Zheng et al., 2020; Zhang et al., 2021b).

The Saussurea obvallata complex, holding significant cultural and medicinal value, represents a promising model system for investigating high-elevation adaptation. It comprises eight closely related lineages distributed across a broad elevational and climatic spectrum in the HHM region, occupying diverse habitats such as alpine subnival belts, alpine meadows, montane shrublands, and rocky scree slopes (Raab-Straube, 2017). These heterogeneous environments reflect the fragmented “sky island” landscapes shaped by the uplift of the Qinghai–Tibet Plateau and subsequent climatic fluctuations, which are thought to have driven habitat isolation, ecological divergence, and lineage diversification. The wide ecological breadth of the S. obvallata complex suggests that sky island dynamics have played a key role in promoting local adaptation and evolutionary differentiation. Traditionally, the lineages within the complex have been recognized as distinct species based on morphological traits and geographic ranges (Raab-Straube, 2017). However, recent genomic analyses indicate shallow divergence times (1.5–4.3 Mya), incomplete reproductive isolation, and pervasive gene flow among lineages (Zhang et al., 2021b), supporting the interpretation of the complex as an evolutionary continuum rather than fully differentiated species. While plastome phylogenies occasionally suggest deeper divergences, these patterns are likely the result of incomplete lineage sorting rather than independent evolutionary trajectories. The S. obvallata complex thus exemplifies a recently diverged species continuum, justifying the application of population and landscape genomic approaches (Stapley et al., 2010; Bell and Gonzalez, 2011). This framing allows investigation into how historical divergence and environmental selection shape current genomic variation and adaptive potential. Given that alpine species are particularly vulnerable to climate change, their persistence hinges upon standing genetic variation enabling adaptive capacity and evolutionary persistence (Carlson et al., 2014; Bell and Gonzalez, 2011). Yet, genomic evidence of these processes remains poorly understood in Himalayan taxa.

Addressing this gap, we combine landscape genomic and phylogeographic methods, integrating plastome phylogenies with nuclear SNPs derived from RAD-seq. Plastomes provide insights into maternal lineage history and deep-time divergence (Zhang et al., 2019b), while RAD-seq data provides a genome-wide snapshot of biparental nuclear variation, enabling detection of population structure, gene flow, and associations with environmental variables, even with reduced genomic representation (Catchen et al., 2013; Andrews et al., 2016).

We hypothesize that adaptive evolutionary capacity in Saussurea obvallata complex is shaped by climatic heterogeneity and ecological selection across environmental gradients. Specifically, we (1) characterize population structure, lineage divergence and gene flow, (2) identify key environmental gradients underlying genomic differentiation and candidate adaptive loci, and (3) evaluate how climate change may alter habitat suitability and elevate genomic vulnerability across lineages. We interpret genomic offset as a proxy for vulnerability based on the sample fraction of the genome, not as a direct measure of adaptive potential. Because RAD-seq sparsely samples the genome, we emphasize that polygenic adaptation with many small-effect loci may be under detected.

2. Materials and methods 2.1. Sample collection, library preparation and sequencing

A total of 108 individuals representing seven of the eight species within the Saussurea obvallata complex were collected from 14 distinct geographical locations (Fig. 1A and Table S1). Voucher specimens were deposited at the Herbarium of the Kunming Institute of Botany, Chinese Academy of Sciences (KUN), and the National Herbarium and Plant Laboratories (KATH). Total DNA was extracted from silica gel-dried leaves using a modified CTAB (Cetyltrimethylammonium bromide) protocol (Doyle and Doyle, 1987) for both plastome and RAD-seq sequencing.

For RAD-seq, genomic DNA from all 108 individuals was digested with EcoRI in a 30-μl reaction mixture, followed by ligation of the P1 adapter. Ligated DNA was pooled, sheared, and size-selected (350–550 bp) before ligating a second adapter (P2). The ligation products were purified using a QIAquick Gel Extraction Kit (Qiagen) and amplified via PCR using Illumina TruSeq primers with 12 cycles of amplification. Amplified products were further purified using the QIAquick Gel Extraction Kit and assessed for quality using an Agilent 2100 Bioanalyzer and qPCR. The final libraries were sequenced on two lanes of the Illumina NovaSeq 6000 platform (Novogene, Tianjin, China) with 150 bp paired end reads. Raw reads were processed to remove low-quality bases and adapter sequences using Trimmomatic v.0.39 (Bolger et al., 2014), ensuring high-quality data for downstream analyses. Reads were cleaned by removing adapters, low-quality bases, and reads with ambiguous barcodes, ensuring only reads with an average Phred score above 30 were retained.

For plastome sequencing, 14 representative individuals (one per site) were selected (Table S2). Paired end (2 × 150 bp) sequencing was performed on an Illumina NovaSeq 6000 platform (Novogene, Tianjin, China). Raw reads were filtered using Trimmomatic v.0.39 (Bolger et al., 2014) following the same criteria as applied for RAD-seq data (Table S3). Plastome reads were assembled de novo using GetOrganelle v.1.6.3a (Jin et al., 2020). The plastomes were annotated using Saussurea obvallata (NC_044726) as the reference with Plastid Genome Annotator (Qu et al., 2019; Zhang et al., 2025) and manually verified in Geneious v.7.0.2 (Kearse et al., 2012). See Tables S2 and S3 for summary statistics.

2.2. Paired-end SNP calling

RAD-seq data were do novo assembled using Stacks v.2.55 (Catchen et al., 2013; http://catchenlab.life.illinois.edu/stacks/). Cleaned reads were demultiplexed and filtered using the process_radtags function, retaining only those with correct barcodes and intact EcoRI restriction sites. SNP identification was performed using the denovo_map.pl pipeline, incorporating ustacks, cstacks, and sstacks modules with parameters M = n = 3. Polymorphic loci present in at least 80% of individuals were retained (Table S4), while loci with heterozygosity > 0.5 were excluded. Only biallelic SNPs with a minor allele frequency (MAF) > 0.05 were retained for downstream analysis using the population module in Stacks.

2.3. Phylogenetic analyses and molecular dating

To reconstruct evolutionary relationships and estimate divergence times within the Saussurea obvallata complex, we analyzed 46 complete plastomes, including 14 newly sequenced and 32 GenBank accessions from major clades Ⅰ and Ⅱ (Xu et al., 2019; Zhang et al., 2021a). A concatenated alignment of 79 protein-coding genes was used for both Bayesian Inference (BI) and Maximum Likelihood (ML) phylogenetic reconstruction. The best-fit substitution model (GTR + G + I) was selected using JMODELTEST v.2.0.1 (Darriba et al., 2012) based on the Akaike Information Criterion (AIC). BI was implemented in MrBayes v.3.2.6 (Ronquist et al., 2012), run for 5 million generations with sampling every 1000 generations and a burn-in of 25%. ML analysis was conducted using RAxML v.8.2.12 (Stamatakis, 2014) under the GTRGAMMA model with 1000 bootstrap replicates.

Divergence time was estimated in BEAST2 v.2.6 (Bouckaert et al., 2014) using an uncorrelated lognormal relaxed clock model and Yule speciation prior, suitable for estimating divergence among distinct evolutionary lineages. Two fossil calibration points were applied: the crown age of Saussurea (9.3 ± 1.2 Mya) and the split between clades Ⅰ and Ⅱ (5.4 ± 0.8 Mya) (Zhang et al., 2021a). Analyses run for 100 million generations with sampling every 10,000 steps, and convergence (ESS > 200) was confirmed in TRACER v.1.7.1 (Rambaut et al., 2018). A maximum clade credibility (MCC) tree was summarized in TreeAnnotator v.2.6.3 (Rambaut and Drummond, 2010) and visualized in FigTree v.1.4.4.

To assess historical gene flow and reticulation among lineages, we used TreeMix v.1.13 (Pickrell and Pritchard, 2012) on the nuclear SNP dataset, testing 0 to 10 migration edges. Residual fit and log-likelihood improvements were used to determine the optimal number of migration events. TreeMix inferences were treated as exploratory hypothesis of introgression and cannot alone distinguish hybridization from incomplete lineage sorting (ILS).

2.4. Genomic consequences and ancestry inference

Genetic diversity metrics such as nucleotide diversity (π), observed heterozygosity (H_o), expected heterozygosity (H_e), and inbreeding coefficient (F_IS), were all calculated from the Stacks-derived SNP dataset using the populations module. Pairwise F_ST and absolute divergence (D_XY) were computed at both population and species-complex levels (Table S5). Pairwise F_ST was further visualized using interactive circos plots generated by the chorddiag v.0.1.3 R package (Flor, 2021).

Population genetic structure was assessed using 15,567 SNPs with multiple complementary approaches. We inferred ancestry proportions (Q) with fastSTRUCTURE v.1.0 (Raj et al., 2014) using an LD-pruned, unlinked SNP set (PLINK --indep-pairwise 50 5 0.2, MAF ≥ 0.05). We ran K = 1–8 with multiple seeds, summarized replicate run using CLUMPP v.1.1.2 (Jakobsson and Rosenberg, 2007) and visualized the consensus Q matrices with DISTRUCT v.1.1 (Rosenberg, 2004) to resolve label switching. The optimal K = 7 was chosen with chooseK.py. ADMIXTURE v.1.3.0 (Alexander et al., 2009) was run only to qualitatively cross-check clustering patterns; its Q estimates broadly agreed with fastSTRUCTURE but were not used in figures or downstream analyses. High-resolution co-ancestry pattern were inferred with fineRADstructure v.0.3.2 (Malinsky et al., 2018) and visualized using custom R scripts (fineradstructureplot.r and finestructurelibrary.r, http://cichlid.gurdon.cam.ac.uk/fineRADstructure.html).

To further explore clustering and individual differentiation, we performed Discriminant Analysis of Principal Components (DAPC) and Principal Component Analysis (PCA) using the adegenet v.2.1.8 R package (Jombart, 2008; Jombart and Ahmed, 2011). The optimal number of DAPC cluster was determined by Bayesian Information Criterion (BIC) and retained principal components. Genome-wide linkage disequilibrium (LD) decay was calculated with PopLDdecay v.3.4091 (Zhang et al., 2019a) as the average squared correlation coefficient (r²) between pairwise SNPs with MAF >0.1 over a 100-kbp window. Sample sizes and populations for LD estimation are listed in Table S1. LD curves for small-sample taxa are not interpreted due to upward bias and instability of r² at low n.

2.5. Genotype-environment associations (GEA) analysis

Genotype–environment associations (GEA) were assessed using 15,567 SNPs (MAF > 0.05) and 31 minimally redundant environmental predictors derived from five categories (climatic, topographic, land cover, habitat heterogeneity, and cloud cover metrics) (Hijmans et al., 2015; Tuanmu and Jetz, 2014, 2015; Wilson and Jetz, 2016; Amatulli et al., 2018) after Pearson correlation filtering (|r| < 0.8, Table S7) and PCA. Environmental datasets were obtained at 2.5 arc-minute resolution from WorldClim (www.worldclim.org), SRTM, and EarthEnv sources (www.earthenv.org). These selected variables formed the basis for both GEA and ensemble species distribution modeling.

GEA were conducted using redundancy analyses (RDA) with vegan v.2.5-7 R package (Bourret et al., 2014; Oksanen et al., 2018) and latent factor mixed model (LFMM) with LEA v.0.3.1 R package (Frichot et al., 2013; Frichot and François, 2015). Candidate SNPs in RDA were identified from significant axes (Table S6) using a ±3.5 standard deviation loading threshold. A baseline RDA was also performed to visualize raw patterns before accounting for population structure. LFMM was run with seven latent factors, with significant associations identified after genomic inflation correction and FDR control (threshold: FDR < 0.001).

We further performed partial RDA (pRDA) (vegan v.2.5-7) on Hellinger-transformed SNP genotypes (0/1/2 ÷ 2; Legendre and Gallagher, 2001) with the model Y ~ X + Condition(C). Here, X represented each of five environmental predictor sets, whereas C accounted for neutral genetic structure and spatial effects, captured respectively by PCA axes from neutral SNPs and up to 10 positive dbMEMs derived from geographic coordinates (Dray et al., 2006). Predictors were filtered for collinearity (pairwise |r| < 0.8; retained variables in Table S7). Significant constrained axes (p ≤ 0.05) were assessed with 999 within-populations permutations. Candidate SNPs were identified exclusively from significant axes using a per-axis |z| ≥ 3.5 threshold (counts in Table 1). Variance was partitioned into pure environment, pure structure, shared, and residual components using varpart (Peres-Neto et al., 2006). Result stability was check by permuting environmental labels within populations (500–1000 runs) while holding Y and C fixed.

To distinguish the effects of geographic structure and environmental factors on spatial genetic variation, we performed both Mantel and partial Mantel tests (Mantel, 1967) with 999 permutations using the vegan v.2.5-7 R package. Mantel tests were used to assess correlations between genetic differentiation F_ST/(1-F_ST) and either geographic (isolation by distance, IBD) or environmental distances (isolation by environment, IBE). Partial Mantel tests were then employed to isolate the independent effects of each factor by controlling for the other. Environmental distance was calculated as the Euclidean distance between sites based on standardized environmental variables. To aid interpretation, we also generated residual regression plots, in which genetic distances were regressed against either geographic or environmental distances after accounting for the alternate predictors.

2.6. Spatial prediction of genomic variation

We applied Gradient Forest (GF) modeling (Ellis et al., 2012; Fitzpatrick and Keller, 2015; Rellstab et al., 2021) to predict spatial patterns of adaptive genomic variation and allele frequency turnover across the Saussurea obvallata complex range, using adaptive SNPs identified through genotype–environment association analysis. The same least-correlated variables selected in the GEA steps were used for GF modeling, based on PCA and pairwise correlation thresholds (|r| < 0.8). To provide complementary insights, we constructed two GF models using the same least-correlated predictors as in GEA: i) all least correlated environmental variables, which allowed us to capture present-day adaptive patterns in their full ecological context, and ii) restricted to climatic variables. The latter enables projection to future climates for genomic offset, as only four of the least-correlated climatic predictors have available future climate datasets. Spatial projections were visualized as RGB composite maps (Ellis et al., 2012; Fitzpatrick and Keller, 2015), illustrating allele turnover linked to environment gradients. Additionally, Procrustes analysis was conducted using the procrustes() and protest() functions in the vegan v.2.5-7 R package to assess the concordance between neutral and adaptive genetic variation, helping identify regions of potential local adaptation or mismatch.

Future genomic variation under RCP 4.5 (2061–2080; multi-model median of General Circulation Models [GCMs]; Rana et al., 2021) was projected by applying the GF model to multi-model climate projections, and genomic offset was calculated as the Euclidean distance between present and future allele turnover predictions (Fitzpatrick and Keller, 2015), identifying populations potentially at risk of maladaptation. We interpret genomic offset as a proxy for vulnerability i.e., predicted allele–environment turnover at RAD-seq–sampled loci. Because RAD-seq captures only a fraction of the genome, offset reflects the detectable portion of adaptive variation and should not be treated as a direct measure of total adaptive potential.

Additionally, we used Generalized Dissimilarity Modeling (GDM) implemented in gdm v.1.5-0 R package (Ferrier et al., 2007) to relate genetic dissimilarity (F_ST) against environmental predictors, fitting models with climatic and all environmental variables while accounting for geographic distance (Ingvarsson and Bernhardsson, 2019; Gallegos et al., 2023). Model performance and variable importance were assessed via permutation tests, and spatial maps of genomic turnover were generated to visualize gradients of genetic differentiation (Ferrier et al., 2007; Luqman et al., 2023). Full modeling parameters, SNP filtering criteria, and statistical settings are provided in Supplementary Methods.

2.7. Species distribution modeling

To characterize the shared environmental niche of the Saussurea obvallata complex, we compiled 158 spatially filtered occurrence records for all eight described species from field surveys, herbarium collections (KATH, KUN), and online databases Plants (PE, E, K, BM, TI; acronyms following Thiers 2020, continuously updated) including Global biodiversity information facility (GBIF; www.gbif.org, accessed 25 Sep. 2021) and JSTOR Global. Occurrence data from all eight lineages were pooled to model the collective niche of the complex, consistent with its treatment as a single evolutionary unit. We used 43 least-correlated predictor variables (|r| < 0.8) representing climate, topographic, habitat heterogeneity, land cover, and cloud cover, based on 158 spatially rarefied occurrence points and the same filtering approach as in GEA. Ensemble species distribution models were run with two predictor sets: (i) 11 least-correlated climatic variables, and (ii) all 43 least-correlated variables across all categories (Table S7).

Ensemble species distribution models were constructed using the Biomod2 v.3.3–7.1 R package (Thuiller et al., 2019), integrating eight algorithms spanning climatic envelope, regression, and machine learning approaches. Models were built with climatic and all environmental variables, and ensemble forecasts were generated via weighted mean consensus models. Model performance was evaluated using TSS and AUC scores, and binary suitability maps were derived using a 50% threshold. Suitable habitats were further classified into low 25–50%), moderate (50–75%), and high (> 75%) suitability. Future habitat projections for 2070 under RCP 4.5 were generated using CMIP5 multi-model median climate data. Habitat shifts were classified as reduced, stable, or expanded based on changes in predicted suitability. Detailed modeling parameters, variable selection procedures, and algorithm settings are provided in Supplementary Methods.

3. Results 3.1. Sequence data processing

After quality filtering, a total of 542.93 Gb of clean bases were retained from 108 individuals of the Saussurea obvallata complex, reduced from an initial 613.93 Gb of raw bases (Table S3). De novo assembly yielded 573,414 RAD tags, from which 3782 high-quality SNP loci (15,567 variant sites) were retained after filtering (r = 0.8) for biallelic status, MAF > 0.05, and call rate > 80%. Plastome assemblies from 14 accessions ranged in length from 151,380 to 152,634 bp across (Table S2), and each genome comprised 113 unique genes, including 79 protein-coding genes (CDS), 30 tRNAs, and 4 rRNAs.

3.2. Phylogenetic analysis and molecular dating

Plastome-based phylogenies revealed a polyphyletic Saussurea obvallata complex, with two clades (Fig. S1): Clade 2A (S. hengduanshanensis, S. obvallata, S. glandulosissima, S. sichuanica, and S. septentrionalis) and Clade 2B [S. sikkimensis and S. kawakarpo, along with S. langpoensis, Saussurea sp. (sp9 LX-2019) and S. bracteata]. Divergence of Clade 2A was dated to ~8.45 million years ago (Mya) [95% highest posterior density (HPD): 5.55–9.37 Mya) (node A; Fig. 2A), while most interspecific divergence within the complex occurred during the Pliocene to early Pleistocene (1.57–4.35 Mya; Fig. 2A).

3.3. Landscape genetic dynamics

Population structure analysis using fastSTRUCTURE identified seven genetic clusters (K = 7) (Fig. S2A), corresponding to recognized species (Fig. 1B). All ancestry proportions (Q) reported here are from fastSTRUCTURE (CLUMPP-summarized). ADMIXTURE runs yielded qualitatively similar clustering but were used only for visualization and did not contribute Q values to our analyses. Minimal admixture was observed across clusters, except between Saussurea sichuanica and S. septentrionalis (Fig. S2D–E). PCA, DAPC, and fineRADStructure confirmed this clustering (Figs. 1C–D, S2B, S3). Species such as S. glandulosissima and S. sikkimensis showed high within-population co-ancestry (Fig. S3), while S. septentrionalis exhibited extensive admixture.

Genetic diversity (π) ranged from 0.007 (Saussurea sikkimensis) to 0.051 (S. glandulosissima), with low observed heterozygosity and near-zero inbreeding across all taxa (Table S1). Genetic differentiation was high (pairwise F_ST: 0.123–0.406; Fig. S4A and Table S5), with the greatest divergence between S. septentrionalis and S. sichuanica (Fig. S4B). In contrast, D_XY patterns differed, with the highest mean absolute differentiation between S. septentrionalis and S. obvallata, and the lowest between S. kawakarpo and S. sikkimensis (Fig. S4C and D; Table S5). LD decay patterns differed among taxa within the complex. For several small-sample species (e.g., KS, SS, SR), r² did not decay below ~0.5 at the maximum distance; consistent with low n (Table S1), these curves are shown in Fig. S2C and not interpreted further due to small-sample bias. We therefore focus on the four core species with adequate sampling (see Table S1), for which LD decays toward background levels (Fig. 1E). TreeMix inferred minimal migration, with one gene flow event between S. septentrionalis and S. glandulosissima (Fig. 2B), consistent with predominantly allopatric divergence.

3.4. Genotype-environmental association

Given pervasive structure (reflected in elevated long-range LD for some taxa; Fig. S2C), we complemented our baseline RDA with partial RDA (pRDA) that conditions on neutral population structure and space. Baseline RDA identified least correlated environmental predictors (Fig. S5) that collectively explained up to 50% of genetic variation depending on variable groups (Fig. 3A, B, and S6A–H), confirming model robustness (Tables 1 and S6). Baseline RDA was significant for all predictors sets, but as expected, included shared-environment signal (adj.R² ~ 0.27, 0.21, 0.15, 0.19 and 0.21 for climate, cloud cover, habitat heterogeneity, land cover and topographic respectively; Table 1). RDA identified 1476 candidate SNPs associated with environmental predictors, with some SNPs linked to multiple predictors (533 with climatic, 547 with cloud cover, 753 with habitat heterogeneity, 235 with topographic, and 594 with land cover; Fig. S7). The largest numbers of candidate SNPs were associated with solar radiation (srd, 458 candidate SNPs) and spatial variance (varian, 365 candidate SNPs), followed by June cloud cover (jun, 290 candidate SNPs) and potential evapotranspiration (pet, 288 candidate SNPs) (Fig. 3E). These values reflect the number of loci detected along each environmental gradient, rather than implying stronger effect sizes (Fig. 3A). Species-specific responses highlight complex selective landscape across the HHM region (Fig. S6).

Specifically, we derived neutral PCs from a structure-only SNP panel, including up to ten positive dbMEMs, and fitted Y ~ Environment + Condition (PCs + dbMEMs) with restricted permutations within populations (999 perms). After formal conditioning, global pRDA models remained significant (all p = 0.001) with smaller, structure-free effects (Fig. S8; Table 1): climatic adj. R² = 0.05 (six constrained axes), cloud = 0.04 (five axes), habitat heterogeneity = 0.03 (three axes), land cover = 0.03 (five axes), and topographic = 0.05 (six axes). The corresponding numbers of SNPs exceeding |z| ≥ 3.5 on significant axes were 411 (climatic), 504 (cloud cover), 110 (habitat heterogeneity), 335 (land cover), and 444 (topographic) (counts are per set and may overlap among sets) (Table 1). Variance partitioning showed that pure environment accounts for ~3–7% of allele-frequency variation after removing structure/space, with a comparable fraction attributed to pure structure and the remainder shared or residual (Fig. 3G). Null env-label permutations within populations produced adj.R² near zero, whereas observed values lay in the extreme right tail for every set (Fig. S9), indicating that the retained signal is unlikely to arise from demographic confounding alone. Across pRDA and Gradient Forest, temperature seasonality (bio4), metrics of cloud regime (e.g., interannual intensity and angular scatter, theta/rgb), and local relief (slope/roughness; slp/vrm) consistently emerged as the strongest environmental correlates of allele-frequency turnover (Figs. 3E, F, and S8F–I). These drivers align with the HHM's steep climatic gradients and complex topographic, and they remained influential after rigorous control for neutral structure and spatial autocorrelation.

We next asked whether single-locus tests supported these multivariate patterns. LFMM, which models latent structure directly, identified 1296 adaptive SNPs significantly associated with June cloud cover (|Z| > 2, FDR < 0.01; Figs. S10 and S11), corroborating the pRDA emphasis on cloud regime. Although PCQ was not retained as a top pRDA predictor, high LFMM z-values (7–11) for SNPs correlated with June cloud cover suggests an indirect link (Fig. 3D). Negative z-values for annual aridity and temperature seasonality indicate reduced allele frequencies under harsh conditions, consistent with stress driven selection, whereas positive associations with habitat heterogeneity and evergreen/deciduous needleleaf cover point to adaptive response in structurally complex or stable environments (Fig. S7). We interpret these as signals within the detectable portion of the genome sampled by RAD-seq and emphasize that polygenic adaptation with many small-effect loci may be under-detected.

Mantel and partial Mantel tests revealed the role of both IBD and IBE in genetic divergence (Fig. S12). Genetic distance was significantly correlated with both geographic (r = 0.594, p = 0.001; Fig. S12A) and environmental distances (r = 0.476, p = 0.001; Fig. S12D). These correlations remained significant after controlling each other (r = 0.5 and r = 0.322, respectively; Fig. S12B and E). Residual regression plots further visualized the partial Mantel effects (Fig. S12C and F).

3.5. Spatial prediction of adaptive genetic variation

Gradient Forest modeling based on adaptive SNPs evaluated the relative strength of environmental variables highlighting slope, PCNM2, elevation, temperature seasonality, and PCQ as the most influential predictors of allele frequency turnover (Figs. 3E and S6I). Populations in the western and eastern Himalayas shared similar adaptive genomic compositions, distinct from those in the Hengduan Mountains (Figs. 4A, C and S13B, C). In environmentally heterogeneous regions, genomic patterns remained relatively stable and were shaped primarily by habitat uniformity, wet-season temperature, and water vapor availability (Figs. S13A and S14). These spatial trends mirrored IBD and IBE signals (Fig. S12), suggesting that both geographic isolation and environmental gradients jointly shape adaptive variation across the landscape.

Procrustes analysis revealed a spatial mismatch between neutral and adaptive SNP turnover (Fig. 4B, E), highlighting regions where local adaptation divergence departs from neutral population structure. GDM further supported this pattern, identifying temperature seasonality, PCQ, and slope as dominant predictors of genomic turnover, with sharp I-spline curves indicating strong environmental selection (Fig. S15). Turnover rates peaked along with the temperature seasonality gradient in the western Himalayas (Fig. 4C and D), and adaptive SNPs showed stronger environmental responses than neutral loci. The strong model performance (R² > 0.7) confirmed that key climatic variables effectively captured genomic differentiation. Despite distinct climatic conditions, populations in the northwestern and southeastern Himalayas exhibited similar patterns of adaptive variation, suggesting parallel adaptation across these regions. In contrast, the Hengduan Mountains harbored more distinct genetic compositions, indicating unique divergence shaped by localized environmental factors (Fig. 4C–E). Both GF and GDM analyses consistently identified temperature seasonality as the primary driver of genomic turnover in the Western Himalayas, while environmental gradients appeared weaker across the eastern Himalayas and Hengduan ranges (Figs. 4D and S13B, C).

Projected genomic offset under the 2070 (RCP4.5) highlights vulnerability hotspots in the southwestern and southern Eastern Himalayas, and parts of the Hengduan Mountains (Fig. 4F), indicating an increased risk of climate-driven maladaptation. Because four climatic variables lacked future projections and were held static, these projections are likely conservative (i.e., genomic offset may be underestimated). We interpret genomic offset as a vulnerability proxy at RAD-seq sampled loci. Higher values reflect greater predicted allele-environment turnover at sampled loci. Nevertheless, unmeasured genomic regions may harbor adaptive variants, so high offset does not preclude latent adaptive potential.

3.6. Uncovering environmental heterogeneity using species distribution modeling

Species distribution models were built with both least correlated climatic and environmental predictors (Table S7), showing that top predictors varied across species (Fig. S5). The constructed model for the complex using ensemble modeling (eight algorithms) achieved high accuracy (AUC: 0.94–1.0; TSS > 0.85; Fig. S16 and Table S8). While niche models don't reflect species-level distributions, they capture the broader habitat envelope in which adaptive divergence occurs.

Current projections matched known distributions (Fig. 5A), with highest suitability in the Hengduan Mountains and Eastern Himalayas (Fig. S17). Future projections (2070, RCP4.5) suggested northward expansion and southern contraction, particularly in the western Himalayas (Fig. 5B and C). Habitat suitability was strongly influenced by temperature, precipitation, and land cover variables (Figs. 5A and S16), underscoring the role of heterogeneous high-altitude environments in shaping adaptive potential.

4. Discussion 4.1. Unraveling the genetic consequences: evolutionary insights into the Saussurea obvallata complex

The HHM region has long been recognized as a cradle for alpine biodiversity and speciation shaped by orogeny and late Miocene–Pleistocene (Wen et al., 2014; Sun et al., 2017; Ding et al., 2020). The emergence of isolated “sky island” habitats likely promoted habitat fragmentation and allopatric divergence (Quade et al., 1989; An et al., 2001; Guo et al., 2008; Miao et al., 2012; Li et al., 2023). Reflecting this biogeographic complexity, the Saussurea obvallata complex occupies a broad range of environmental conditions across its distribution, suggesting that the formation and persistence of sky island habitats may have provided the ecological foundation for local adaptation within each lineage. Within this context, the S. obvallata complex diversified from ~8.5 to 1.6 Mya, consistent with a scenario in which tectonic uplift and glacial–interglacial cycling jointly drove isolation and ecological opportunity (Hewitt, 2004).

Nuclear SNPs reveal seven cohesive but not fully reproductively isolated genetic clusters, representing a speciation continuum rather than discrete, fully diverged species. The moderate to strong genetic differentiation among clusters, particularly between Saussurea sikkimensis and S. sichuanica, suggests the presence of historical barriers to gene flow, potentially reinforced by steep ecological gradients. In contrast, admixture signatures in S. septentrionalis and S. glandulosissima point to ongoing or recent gene flow. These results align with the hypothesis that S. obvallata and close allies represent lineages in various stages of divergence within a reticulate speciation framework.

Despite this structure, the overall heterozygosity is low (F_IS ≈ −0.014 to −0.002), consistent with self-compatibility, small effective population size, and restricted pollen and seed dispersal patterns reported in Saussurea obvallata and congeners (Semwal et al., 2019; Zhang et al., 2023) and common in alpine taxa occupying fragmented and elevation-constrained ranges (Keller and Waller, 2002; Wu et al., 2019; Semwal et al., 2020). LD decay patterns also require caution. Elevated long-range r² in several small-sample taxa (e.g., S. kawakarpo, S. sichuanica, S. septentrionalis) can reflect localized selection (long-range r² 0.5) (Chhatre and Rajora, 2014) but it is also expected under strong background structure and limited n (Kemppainen et al., 2015). Because r² is upward-biased at small sample sizes, we present full LD curves in the supplement and refrain from over-interpretation in the main text (see Table S1 for per-species n). In well-sampled lineages LD decays toward background. Given the deep divergence and pervasive structure within the complex, elevated genome-wide LD can covary with environment and potentially confound genotype–environment associations. Accordingly, we treat RDA candidates as hypotheses pending validation, emphasizing concordant evidence from structure-aware LFMM and Gradient Forest importance.

Previous phylogenetic studies based on plastid DNA nuclear ITS regions, typically treated the Saussurea obvallata complex as a single species and only sampled one representative (Raab-Straube, 2003; Wang et al., 2009; Xu et al., 2019). However, these markers are insufficient for capturing the full extent of genetic variation in rapidly diversifying lineages. In our study, nuclear SNP clustering revealed clear genetic structure among lineages, while plastome phylogenies exhibited weakly supported relationships, reflecting cytonuclear discordance. Although TreeMix identified a migration edge between S. septentrionalis and S. glandulosissima, this single signal of gene flow is insufficient to disentangle incomplete lineage sorting (ILS) from historical hybridization, both of which are plausible in recently radiated lineages (Degnan and Rosenberg, 2009). Given the shallow nuclear divergence and recent diversification times, ILS remains a parsimonious explanation alongside localized gene flow. More rigorous statistical approaches, such as D-statistics (Durand et al., 2011), HyDe (Blischak et al., 2018), or coalescent-based model testing, will be required to robustly resolve the drivers of plastome–nuclear discordance. We therefore interpret this discordance cautiously, acknowledging the joint contributions of both processes (Liu et al., 2024; Xu et al., 2025), and reinforce the necessity of integrating both nuclear and organellar data when resolving species boundaries in recently radiated alpine lineages.

4.2. Decoding the drivers of divergence: environmental selection and the speciation continuum

Environmental heterogeneity is a major organizer of genomic variation in the Saussurea obvallata complex, even after accounting for background structure. We used structure-aware GEA framework that combines LFMM and RDA with partial RDA conditioning on neutral genetic PCs and spatial dbMEMs, and we evaluated significance with permutations restricted within populations (Borcard et al., 1992; Dray et al., 2006; Legendre and Legendre, 2012; Capblancq et al., 2018; Forester et al., 2018). Across predictor sets, we detected associations with several climatic variables, notably precipitation of the coldest quarter (PCQ), temperature seasonality, and mean temperature of the wettest quarter; June cloud cover emerged as a particularly strong correlate in LFMM. In pRDA, the pure environmental fraction (partial adj.R²) was modest but significant under blocked permutations (P = 0.001), ranging ~0.03–0.07 with 3–6 significant constrained axes per set. Partial adj. R² was ~0.05–0.07 for climatic and topographic variables and ~0.03–0.05 for cloud cover, land cover, and habitat heterogeneity, indicating that, beyond neutral structure and spatial autocorrelation, environmental predictors explain a detectable component of allele-frequency variation. While loadings do not measure effect size, the concordance between LFMM and pRDA (despite structural conditioning) supports a genuine environment-linked component of genomic turnover (Capblancq et al., 2018; Forester et al., 2018). Moreover, null-label tests (permuting environments within populations) produced partial R² distributions well below observed values for each set, reinforcing that signals exceed a structure-only null. Collectively, these results point to growing-season cloud cover, temperature seasonality, and related climatic gradients as repeatable axes of localized adaptation across the HHM (Frichot et al., 2013; Fitzpatrick et al., 2021), especially in the Western Himalayas.

Spatial modeling further corroborates the ecological basis of adaptive divergence. Gradient Forest highlighted temperature seasonality, slope, elevation, and cloud metrics among the most important predictors of allele-frequency turnover, and GDM I-splines showed steep genomic turnover along seasonality, PCQ, and slope gradients (Ellis et al., 2012; Fitzpatrick and Keller, 2015; Capblancq et al., 2018; Fitzpatrick et al., 2021). The resulting east–west partition, with the Hengduan Mountains harboring distinct allelic compositions relative to the Western/Eastern Himalayas, mirrors the joint action of IBD and IBE, as supported by partial Mantel analyses and Procrustes comparisons. In the HHM's rugged terrain, steep climatic and topographic gradients likely constrain gene flow and reinforce local adaptation, promoting divergence even among geographically proximate lineages.

Crucially, our results provide empirical support for a speciation continuum within the Saussurea obvallata complex. The interplay between environmental selection and spatial isolation appears to drive divergence at varying stages across the range. Lineages in the Hengduan Mountains exhibit stronger genomic structure and more distinct adaptive turnover (later-stage divergence), whereas several Western/Eastern Himalayan lineages retain weaker structure and more shared adaptive variants (earlier-stage divergence with ongoing gene exchange). This gradient of divergence stages reflects how environmental heterogeneity and geographic structure collectively shape lineage differentiation in high-altitude plant complexes.

We emphasize, however, that reduced-representation RAD-seq sparsely samples the genome and is underpowered for highly polygenic architecture; it preferentially detects locus-specific signals and may underestimate diffuse, small-effect adaptation (Bourret et al., 2014; Andrews et al., 2016; Byers et al., 2016). Accordingly, we treat pRDA candidates as hypotheses unless supported by LFMM and GF/GDM, and we recommend whole-genome resequencing and functional validation (common-garden/reciprocal transplant) to resolve the polygenic basis of high-altitude adaptation in Saussurea obvallata. Notably, while high-elevation populations exhibit robust adaptive responses, low-elevation ecotypes show reduced genomic flexibility, suggesting greater vulnerability to rapid climate change.

4.3. Genomic vulnerability and adaptive capacity under climate change

Quantifying genomic vulnerability and adaptive capacity is essential for understanding how high-elevation taxa such as the Saussurea obvallata complex will respond to rapid climate change. We used GF to translate environment-allele relationships into genomic offset under future climate (RCP 4.5, 2070), interpreting offset as a risk, often regarded as a proxy for climatic mismatch at RAD-sampled loci rather than direct adaptive potential (Fitzpatrick and Keller, 2015; Rellstab et al., 2021; Chen et al., 2022). Resulting maps reveal strong spatial heterogeneity: high-offset areas occur in the southern Eastern Himalayas, south-western Western Himalayas, and peripheral Hengduan Mountains, whereas several central Himalayan populations exhibit minimal offset, consistent with putative evolutionary refugia where extant genotype-environment combinations may remain comparatively robust (Rellstab et al., 2021; Chen et al., 2022). Ensemble SDMs corroborate these projections, highlighting upslope and poleward shifts with contractions along southern margins, aligning with genomic offset hotspots (Swab et al., 2014; Platts et al., 2019).

Signals of adaptive capacity align with elevational and spatial context. High-altitude ecotypes exhibit greater allelic diversity and broader climatic tolerances, whereas lower-elevation lineages show reduced genomic flexibility, implying greater vulnerability. Some peripheral populations, despite small size or marginal distribution, harbor distinct allele and plausibly constitute evolutionary significant units (ESUs) relevant to long-term adaptation. As habitats fragments and shrinks, the long-term persistence will depend on in-situ buffering and access microrefugia (Xu et al., 2017; Rana et al., 2021). We stress two interpretive caveats. First, genomic offset reflects vulnerability only at RAD-sampled loci. Populations may harbor unmeasured adaptive variation in unsampled genomic regions (Rellstab et al., 2021; Chen et al., 2022). Second, realized outcomes depend on multiple biological dimensions – local adaptation, phenotypic plasticity, effective population size, dispersal ability, and demographic regime (Bush et al., 2016; Mokany et al., 2019; Thomas et al., 2023).

A central determinant of genomic vulnerability is generation time. In long-lived woody plants, slow turnover constrains allele-frequency change, so high genomic offset is more likely to manifest as persistent maladaptation (Fitzpatrick and Keller, 2015; Bay et al., 2018). In contrast, herbaceous alpine taxa such as Saussurea typically have shorter generations, exposing more cohorts to selection and enabling faster tracking of shifting climates. When effective sizes and connectivity are adequate, shorter-generation species can reduce genomic mismatch more rapidly, weakening the link between static offset and realized fitness declines under sustained selection (Exposito-Alonso et al., 2019; Feng et al., 2025). Therefore, genomic offset maps should be interpreted as probabilistic, with demographic consequences mediated by life history, not as deterministic forecasts. Generation time interacts with other factors affecting vulnerability. Low diversity and small N_e can negate the benefits of short life cycles in “sky-island” populations (Willi et al., 2006). Degraded connectivity limits the influx of pre-adapted alleles along elevational corridors (Manel et al., 2003), while seed-bank dynamics and clonality modulate genotypic turnover (Evans and Dennehy, 2005). Plasticity may buffer short-term but is often constrained in narrowly adapted ecotypes (Nicotra et al., 2010). Together, these interactions help explain why central HHM regions with low offset and robust demographics served as refugia, whereas high-offset peripheries are most concerned where populations are small, fragmented, or connectivity is degraded.

Finally, signals converge across methods on the same climatic axes (i.e., temperature seasonality, slope/elevation, and cloud metrics), both as drivers of present allele-frequency turnover (LFMM, GF, and pRDA with structure conditioning) and as gradients along which offset peaks, reinforcing a climate-selection basis for predicted genomic change. Because our pRDA conditioned on neutral PCs, spatial dbMEMs used permutations restricted within populations, and partial R² values exceeded within-population label-permutation nulls (Capblancq et al., 2018; Forester et al., 2018), the environmental component we detect is unlikely to reflect background structure alone. Even so, offset remains a RAD-limited proxy and may under-detect polygenic architectures dispersed across the genome. Consequently, the most informative management path is to pair offset mapping with demographic and genetic context (e.g., N_e, recruitment, seed-bank status, connectivity), and to validate fitness consequences with common-garden and reciprocal-transplant experiments that partition plastic and genetic responses over realistic generation times. Overall, integrating generation time with genomic offset clarifies which lineages possess the demographic capacity to track rapid climate change and which are intrinsically constrained.

4.4. Conservation insights from a genomic and spatial framework

The Saussurea obvallata complex serves as a compelling model for connecting evolutionary history with forward-looking conservation planning under climate change. Across the HHM, a legacy of orogeny and Quaternary climate dynamics has produced a mosaic of glacial refugia and dispersal barriers, with “sky island” populations retaining distinct genetic variants that likely contribute to adaptive capacity (Byers et al., 2016; Kopuchian et al., 2016). By integrating population genomic data (RAD-seq) with structure-aware GEA (LFMM; baseline partial RDA conditioned on neutral PCs and dbMEMs), GF-supported genotype-environment associations, GDM, genomic offset, and ensemble SDMs, we jointly mapped historical legacies and prospective genomic vulnerability. Under RCP 4.5 (2070) climate scenario, models indicate northward, and upslope shifts with contraction along the southern margins of the Western and Eastern Himalayas, coincident with genomic offset hotspots, emphasizing the need to maintain elevational corridors and mid–altitude transition zones as critical migration pathways (Platts et al., 2019).

From a conservation standpoint, these results support a tiered dual strategy. First, safeguard core refugia in the central HHM where high diversity co-occurs with low genomic offset, as these areas likely harbor resilient genotype–environment combinations that can function as evolutionary anchors. Second, mitigate risk in high-offset peripheries by maintaining or restoring connectivity, prioritizing microrefugia, and implementing targeted genomic monitoring to reduce climate–genome mismatch and support adaptive responses (López-Uribe et al., 2014; Fitzpatrick and Keller, 2015; Chen et al., 2022). Practically, integrating GF, GDM, and ensemble SDMs pinpoints elevational and cross-valley gene-flow corridors and fine-scale microrefugia where management can have outsized benefits. Because partial RDA detected a modest but significant pure environmental signal after conditioning on structure (blocked permutations P = 0.001; within-population label permutations below observed partial R²), and because GF importance converged on similar predictors, we view these spatial priorities as robust to major confounders.

Where projected change is likely to outpace adaptive responses, we also outline proactive interventions. For moderate-high offset populations with limited connectivity or small N_e, assisted gene flow (AGF) can introduce climate-matched alleles while maintaining ecological compatibility (climate-adjusted provenance; Prober et al., 2015) and minimizing outbreeding risk (Aitken and Whitlock, 2013; Frankham, 2015; Aitken and Bemmels, 2016). For extreme offset locales where suitable climate is projected to shift beyond dispersal distances or local habitat is lost, assisted migration (AM) may be evaluated as a last-resort option following pilot trials, pathogen screening, and adaptive monitoring (McLachlan et al., 2007; Aitken and Bemmels, 2016). Across risk tiers, “no-regrets” actions include: (a) protecting refugia and corridors, (b) reducing non-climatic stressors in high-offset zones, (c) genomic monitoring (diversity, N_e, allele-environment turnover) to track vulnerability, (d) seed banking across elevational and climatic gradients, and (e) safeguarding genetically unique populations, even when vulnerable, given their irreplaceable alleles and potential for long-term adaptive value (Frichot et al., 2013; Fitzpatrick and Keller, 2015; Chen et al., 2022).

5. Conclusion

This study shows that evolutionary dynamics in the high-altitude Saussurea obvallata complex reflect the combined effects of historical tectonics, climatic oscillations, and spatially heterogeneous environmental selection across the HHM region. By integrating plastome and RAD-seq data with a structure-aware GEA framework (baseline RDA and pRDA), spatial genomic models (Gradient Forest and GDM), and ensemble SDMs, we resolve a speciation continuum characterized by swallow nuclear divergence and geographically structured allele-environment associations. After conditioning on neutral structure and space, the pure environmental signal remains modest but significant (partial adj.R² ~0.03–0.07; 3–6 constrained axes per predictor set), with null within-population label permutations fall well below observed values. Concordance among pRDA, LFMM (notably June cloud metrics), and Gradient Forest/GDM (temperature seasonality, slope/elevation, cloud regime) supports a genuine environment-linked component of genomic turnover, while acknowledging limits of reduced-representation data and possible under-detection of polygenic effects. Future projections identify genomic offset hotspots along the southern Western/Eastern Himalayas and parts of the Hengduan Mountains, whereas the central HHM populations exhibited lower offset, suggesting putative refugial potential. Offset should be interpreted as a proxy of climatic mismatch at RAD-sampled loci, not a direct measure of total adaptive capacity. We recommend a dual conservation strategy: safeguard diversity-rich, low-offset refugia and mitigate risk in high-offset peripheries via connectivity, microrefugia, and genomic monitoring. Because generation time and demography mediate the translation of genomic offset to realized risk, integrating connectivity and population management with genome-wide resequencing, reciprocal-transplant or common-garden validation, and long-term genomic–demographic surveillance will best strengthen climate-aware conservation for alpine biodiversity.

Acknowledgements

The authors would like to acknowledge Dr. Bo Xu, Dr. Yang Yang, Dr. Xu Zhang and Dr. Yazhou Zhang for material's collections, and Miss Minshu Song for her help with the molecular lab work. This study was supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) Program (2024QZKK0200), the Key Projects of the Joint Fund of the National Natural Science Foundation of China (U23A20149), the National Natural Science Foundation of China (32070232).

CRediT authorship contribution statement

Hum Kala Rana: Writing – original draft, Writing – Review & Editing, Software, Validation, Formal analysis, Investigation, Resources, Data Curation, Visualization. Santosh Kumar Rana: Writing – Review & Editing, Methodology, Software, Formal analysis, Visualization. Jacob B. Landis: Writing – Review & Editing, Methodology, Validation. Hang Sun: Writing – Review & Editing, Conceptualization, Resources, Project administration, Funding acquisition. Dong Luo: Writing – Review & Editing, Conceptualization, Investigation, Resources, Data Curation, Supervision, Project administration, Funding acquisition.

Data availability statement

The raw reads generated by shallow genome sequencing for the annotated complete chloroplast sequences, as well as raw RAD-seq reads, have been deposited in the National Genomics Data Center (NGDC; https://ngdc.cncb.ac.cn/) under the GSA accession number CRA035159.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could 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.12.004.