1. Introduction

Plant species with narrow habitat ranges face heightened extinction risks from climate change, environmental shifts, and human activities (Bijlsma and Loeschcke, 2012; Ralls et al., 2018). Threatened plant species with small, fragmented populations often experience genetic drift, which diminishes their genetic variability and heightens differentiation as a result of bottleneck effects (Kéry et al., 2000; Haag et al., 2010). Such genetic divergence may trigger outbreeding depression, exacerbating threats to population survival (Ellstrand and Elam, 1993; Spielman et al., 2004; Frankham et al., 2011; Brys and Jacquemyn, 2016). Collectively, these factors diminish a species' ability to adapt to habitat degradation or climatic shifts (Rodger et al., 2021). Consequently, assessing genetic diversity, population structure, demographic history, fitness, and deleterious mutations is essential for conserving rare plant species (Aguilar et al., 2008; Hill, 2010; Hoban et al., 2013; Schäfer et al., 2020; Kardos et al., 2021; Kyriazis et al., 2021; Ducrettet et al., 2023).

One group of plants at critical risk of extinction are species endemic to karst ecosystems in China (Liu et al., 2024). Research that has assembled and annotated whole genomes have provided insights into the adaptive mechanisms, population genetic structure, and potential conservation strategies for karst-adapted and threatened plant species. For instance, a genomic study of Begonia masoniana (Chen et al., 2024) revealed that isolated populations exhibit limited gene flow, high inbreeding, and an accumulation of deleterious mutations. Similarly, in Urophysa species, bottlenecks, restricted gene flow, niche specialization, and inbreeding have been shown to threaten species, while physiological and transcriptomic analyses have identified positively selected genes associated with ion regulation and cell wall maintenance, which are key adaptations for survival in limestone habitats (Xie et al., 2026). Furthermore, research has also shown that Camellia rubituberculata adapted to karst environments through gene family expansion and positive selection on genes related to photosynthesis, phytohormone signaling, and calcium/ion transport, providing genomic evidence for adaptive evolution and population divergence between wild and cultivated groups (Yan et al., 2026). Nevertheless, research on the diversity, conservation, and adaptive evolution of karst-adapted plants remains limited (Feng et al., 2020; Cao et al., 2023).

Previous studies have examined the adaptive ability of karst adapted species in the genus Oreocharis. Analysis of the genomes of two species (O. esquirolii and O. maximowiczii) distributed in Guizhou Province, southeastern China, indicated that O. esquirolii is endangered due to diminished genetic diversity, pronounced inbreeding, and reduced ability to eliminate deleterious alleles (Peng et al., 2025). However, our understanding of the adaptive ability of other karst-adapted Oreocharis species remains limited. For example, we know little about the conservation status and capacity for adaptation in species such as Oreocharis mileensis. Originally described as Paraisometrum mileense (Weitzman et al., 1997), this species was reclassified under Oreocharis s.l. (Möller et al., 2011). O. mileensis is endemic to karst forests in southwest China (1100–2150 m), restricted to Yunnan, Guizhou, and Guangxi (Xu et al., 2009; Gao and Xu, 2011; Chen et al., 2014). This resurrection plant possesses high potential ornamental value and is extremely drought resistant (Chen et al., 2014). It is listed as Endangered (EN) in China's biodiversity red lists (Qin et al., 2017) and prioritized in national PSESP rescue plans (Sun et al., 2013, 2019a). The fragmented distribution and small populations of O. mileensis make the species highly susceptible to habitat destruction, climate change, and anthropogenic pressures (Chen et al., 2014; Liu et al., 2024). For instance, surveys indicate that there remain only 10 remnant populations/locations, one of which has been lost in Shui Jing Wan due to the construction of the highway. Low genetic diversity and environmental sensitivity further threaten its survival (Chen et al., 2014; Li et al., 2014).

While traditional conservation genetics relies on neutral markers (Allendorf, 2017), advancements like whole-genome resequencing (WGR) now enable genome-wide conservation analyses, improving the accuracy of genetic assessments (Abzhanov et al., 2008; Weber et al., 2013; Poelstra et al., 2014; Jing et al., 2023). These high-throughput methods facilitate the reconstruction of species' evolutionary histories (Leitwein et al., 2020), identification of putatively adaptive loci, and estimation of critical parameters such as effective population size (Ne), inbreeding, and mutation load (Hohenlohe et al., 2010; Nielsen et al., 2011; Lamichhaney et al., 2017). Such tools are invaluable for conserving plant species with extremely small populations (PSESP), a category of highly threatened flora demanding targeted interventions (Ma et al., 2013; Sun et al., 2019b, 2024; Yang et al., 2020). The PSESP approach combines ecological, propagation, and genetic data to design conservation strategies, as demonstrated in Acer, Dipteronia, Magnolia, and Rhododendron (Yang et al., 2019, 2024; Ma et al., 2021, 2022; Cai et al., 2024; Feng et al., 2024; Li et al., 2024).

This study employs conservation genomics, based on a de novo genome assembly and annotation, to evaluate the genetic diversity of Oreocharis mileensis, identify drivers of population decline, and assess the adaptive capacity of the species under future climates.

2. Materials and methods 2.1. Collection of plant material

Young leaves for whole-genome sequencing were obtained from a cultivated individual at Kunming Botanical Garden (origin: Shilin County, Yunnan). For transcriptome analysis, we collected roots, stems, leaves, and inflorescences from wild plants of the same Shilin population. All samples were flash-frozen in liquid nitrogen and maintained at −78.5 ℃. For population resequencing, fresh leaves were collected from 10 natural populations (5–20 m sampling intervals, adjusted for population size) and preserved in silica gel. Finally, 107 individuals of Oreocharis mileensis and two samples of O. hekouensis (as an outgroup) were selected for genetic diversity and phylogenetic analyses through whole genome resequencing (Table S1).

2.2. Genome sequencing

The ploidy level of Oreocharis mileensis is 2n = 34 (Tan et al., 2011). To estimate nuclear DNA content, we followed the flow cytometry protocol of Qiu et al. (2020). Leaf tissue was chopped into small fragments (1–2 cm) in 500 μl of CyStain UV Precise P extraction buffer (Sysmex) and filtered through a 50-μl mesh. A staining solution containing propidium iodide (4.29 μl), CyStain staining buffer (0.71 ml), and RNAse A (2.14 μl) was added to stain the nuclei. Samples were kept on ice before analysis on a BD LSRII H4760 flow cytometer (BD Biosciences) using the PI laser detector at 480 V, with a minimum of 2000 nuclei counted. Data was analyzed using FACSDiva v.8.0.1. Finally, genome size was estimated as 1.74 Gb.

Genomic DNA and transcriptome sequencing were performed on various platforms to ensure accuracy and efficiency. For DNA extraction, library construction, and PacBio sequencing, we first extracted genomic DNA from leaf samples using the CTAB method (Doyle and Doyle, 1987) and DNeasy Plant Mini Kit. DNA purity, quantity, and integrity were assessed using a NanoDrop One spectrophotometer, Qubit 3.0 Fluorometer, and Femto Pulse. DNA (8 μg) was fragmented using Megaruptor 3, purified with AMPure PB beads, and used to construct SMRT bell libraries with the Pacific Biosciences SMRT bell Express Template Prep Kit 2.0. Libraries were size-selected on the BluePippin^TM system (15 Kb insert size), primer-annealed, and bound to polymerase using the DNA/Polymerase Binding Kit. Sequencing was performed on the PacBio Sequel Ⅱ platform, generating HiFi reads.

For BGI-Next Generation DNA Library construction and sequencing, we used genomic DNA (1 μg) fragmented with a Covaris ultrasonicator, filtered (200–400 bp), PCR-amplified, purified, and circularized. Libraries were sequenced on the DNBSEQ platform, with raw reads filtered using fastp (v.0.21.0) to discard low-quality reads (Chen et al., 2018). For Hi-C DNA, genomic material was fragmented, blunted, phosphorylated at the 5′ ends, and dA-tailed at the 3′ ends. Adapters were ligated, biotin-labeled regions captured, and Hi-C DNA amplified via PCR. After quality control, quantification, and circularization, libraries were sequenced on the DNBSEQ-T7 platform using paired-end sequencing. Iso-Seq transcriptome analysis was based on full-length mRNA. Total RNA (500 ng) was diluted to 9 μL, reverse transcription primer added (Table S2), and samples underwent reverse transcription, amplification, purification with AMPure beads, and adapter ligation. Libraries were sequenced on the PromethION platform (Oxford Nanopore Technologies).

2.3. Genome assembly

We generated ~91 Gb HiFi data (~46×), ~288 Gb (~144×) DNBSEQ reads, ~222 Gb Hi-C reads, and ~15 Gb Iso-Seq reads (Tables S3–S6). Contigs were generated from HiFi and Hi-C reads using Hifiasm (v.0.19.8-r602; Cheng et al., 2021), with the haplotype assembly chosen for downstream analysis. Hi-C reads were aligned with Juicer, and chromosomes were assembled using 3D-DNA (Dudchenko et al., 2017). Manual corrections in Juicebox refined chromosome boundaries (Durand et al., 2017). Final scaffolds were generated with 100 bp gaps filled by quarTeT (Lin et al., 2023). Telomeres were extended by realigning HiFi reads and assembling contigs with Hifiasm. Genome polishing used NextPolish2 with HiFi and second-generation reads (Hu et al., 2023). Redundans (Pryszcz and Toni, 2016) was employed to eliminate redundant sequences, while low-coverage regions and rDNA were detected either manually or using Barrnap.

2.4. Assessment of genome assembly

Assembly completeness was evaluated using BUSCO (Simão et al., 2015) and LAI analyses. Sequencing reads were aligned with BWA and Minimap2 (Li, 2013, 2018), followed by removal of non-primary alignments. Genome coverage and mapping efficiency were then calculated from filtered alignments. Coverage depth was checked for a Poisson distribution expected for redundancy-free single- and multi-copy core genes (Fig. S1). HiFi and second-generation sequencing depth across GC content was evaluated for bias (Fig. S2). Minimap2 confirmed chromosome order consistency in collinearity analysis (Fig. S3). Telomeres, rDNA, tandem repeats, and other sequences were mapped to chromosomes (Fig. S4).

2.5. Genome annotation

EDTA identified transposable elements de novo (parameters: –sensitive 1 –anno 1), generating a TE library (Ou et al., 2019). RepeatMasker identified repetitive regions in the genome. A non-redundant set of 384,287 protein sequences from species like Antirrhinum majus, Boea hygrometrica, Coffea canephora, and Primulina huaijiensis was reference for annotation. Iso-Seq reads were aligned to the genome using Minimap2 (Li, 2018) and assembled with StringTie (Pertea et al., 2015), providing transcript evidence for annotation and model training. PASA pipeline was used to annotate gene structure (Haas et al., 2003), and full-length genes were identified by aligning transcripts with reference proteins. The gene prediction tools AUGUSTUS and SNAP underwent five iterative rounds of training and optimization utilizing the complete set of full-length genes (Korf, 2004; Stanke et al., 2008).

The MAKER2 pipeline performed ab initio prediction, integrating transcript and homologous protein evidence (Cantarel et al., 2008). Due to MAKER's lower accuracy (Cook et al., 2018), MAKER and PASA annotations were merged using EvidenceModeler (EVM) for consistent gene annotation (Haas et al., 2008). TE protein domains were identified with TEsorter and masked in EVM (Zhang et al., 2022). PASA updated EVM annotations by adding UTRs and alternative splicing information (Haas et al., 2003). Annotations with abnormal coding frames or short lengths (< 50 amino acids) were removed. tRNA was annotated with tRNAscan-SE (Lowe and Eddy, 1997), rRNA with Barrnap, and ncRNAs with RfamScan. Redundancy was removed, and BUSCO assessed the integrated protein annotations.

Protein-coding gene functions were annotated using (1) eggNOG-mapper for GO and KEGG annotations (Huerta-Cepas et al., 2017); (2) DIAMOND alignment with Swiss-Prot, TrEMBL, NR, and Arabidopsis thaliana databases (Buchfink and Huson, 2015), requiring > 30% sequence identity and E-value < 1e-5; and (3) InterProScan for domain analysis against PRINTS, Pfam, SMART, and CDD databases (Jones et al., 2014).

2.6. Gene family identification and phylogenetic analysis

To determine orthology and orthogroups, we retrieved genome data from 18 plant species from public databases. These species were selected to represent a broad phylogenetic range across Lamiales and related angiosperm lineages, allowing robust inference of orthogroups, gene family evolution, and paleopolyploidy events. Close relatives within Lamiales (e.g., Streptocarpus, Primulina, Boea and Oreocharis) provided resolution within the clade, while more distant eudicot species (e.g., Vitis, Coffea, Jasminum, Tectona) served as outgroups for tree rooting. OrthoFinder2 analyzed orthology and orthogroups (OGs) with "-M msa" (Emms and Kelly, 2019). Protein alignments from 1273 OGs (≥ 84.2% species with single-copy genes) were used to build an ML tree in IQ-TREE (Nguyen et al., 2015) under the JTT + F + R4 model (1000 bootstraps). A species tree was inferred with ASTRAL (Yin et al., 2019) using 4416 single-copy gene trees (≥ 70% species), differing slightly from the ML tree. ASTRAL's result was used for further analysis.

Divergence times were estimated with MCMCTree in PAML (Yang, 2007) using 199 non-missing, single-copy orthologous genes, calibrated with TimeTree. CAFÉ was used to identify gene family expansion, contraction, and rapid evolution using the time tree and 154,621 orthologous families (Hahn et al., 2005).

Proteins were aligned using DIAMOND (Buchfink et al., 2015), and then WGDI was used to analyze collinearity and calculate Ks values (Sun et al., 2022). Orthologous collinear blocks were identified using the Orthology Index (Zhang et al., 2024, 2025). Paleopolyploidy events and ploidy levels were determined via collinearity dot plots, with Ks values estimating event timing and frequency.

2.7. Genome mapping and SNP calling

Whole-genome resequencing generated approximately 10.95 billion reads, totaling 2.24 terabases of raw sequence data. The raw reads were trimmed for adapters and ensured for high quality (Chen et al., 2018). Further, they were mapped to the reference genome with BWA-MEM (Li, 2013), followed by the removal of duplicate reads via SAMtools (Li et al., 2009). To ensure data quality, only bases with Phred quality scores above 20 and mapping quality scores exceeding 30 were retained, yielding a final set of 206,951,631 high-confidence genomic sites (dataset 1). Initial variant calling with GATK identified 133,071,832 raw SNPs in Oreocharis mileensis. These variants were further filtered using VCFtools with stringent criteria (Danecek et al., 2011): sites with genotype quality < 20 or read depth < 5 were marked as missing; non-biallelic sites and non-SNP variants were excluded; SNPs with missing data rates exceeding 20% were removed; and variants with minor allele frequencies below 0.05 were discarded. This filtering process resulted in a refined dataset of 7,415,406 high-quality SNPs (dataset 2) suitable for downstream analyses.

2.8. Inference of population structure based on all loci, putatively adaptive loci, and neutral loci

We used PLINK 1.90b4.1 (Chang et al., 2015) to remove loci in linkage disequilibrium (–indep-pairwise 50 5 0.35), from dataset 2 and retained 752,569 SNPs (dataset 3, all loci). This dataset was further processed to determine environment associated SNPs. The environmental data, which included 19 bioclimatic variables at 30-s resolution, was downloaded from WorldClim v.2.1 database (Fick and Hijmans, 2017). Each environmental factor was extracted through the coordinates of sampling points. We selected the top climatic factors based on weighted R² importance (Fig. S5) using a machine learning gradient forest (GF) model in the R package gradientforest v.0.1–37 (Ellis et al., 2012) by modeling the relationship of climatic variables and 50,000 randomly selected SNPs with 500 regression trees. To reduce multicollinearity among environmental variables, we retained only those with an absolute Pearson correlation coefficient below 0.7, calculated using the R package corrplot v.0.92 (Friendly, 2002; Wei et al., 2017). From this analysis, we selected the four most informative and uncorrelated variables: Mean Temperature of Driest Quarter (BIO9), Mean Temperature of Warmest Quarter (BIO10), Precipitation of Wettest Month (BIO13), and Precipitation of the Driest Month (BIO14). We then applied BayeScEnv (Villemereuil and Gaggiotti, 2015) to detect SNPs associated with environmental factors. BayeScEnv is a univariate method for genotype–environment association. For this analysis, environmental variables were standardized by mean variance, and genotype data in codominant format (derived from dataset 3) were prepared using PGDSpider v.2.1.1.5 (Lischer and Excoffier, 2012). RDA, conducted with vegan v.2.6.2 (Dixon, 2003), was used with VCF files converted via PLINK v.1.9 (Purcell et al., 2007), with missing values replaced by the most frequent non-missing value (Forester et al., 2018). A 999th ANOVA assessed significance, with candidate SNPs identified at ±3 SD (p = 0.0027). SNPs determined by both methods were classified as putatively adaptive SNPs (dataset 4; 1737 SNPs), other SNPs were classified as neutral SNPs (dataset 5; 750,832 SNPs). Gene Ontology enrichment via eggNOG-mapper (Cantalapiedra et al., 2021) was conducted to determine those adaptive SNPs functions.

Population structure was then inferred using ADMIXTURE 1.3.0 with the optimal number of populations (K) determining the minimum cross-validation error. We performed PCA in R for each dataset and constructed a Maximum Likelihood (ML) phylogenetic tree using RAxML (Stamatakis, 2014) specifically for dataset 3. This dataset was chosen for tree construction as previous studies have shown that different tree-building methods produce similar results for comparable genomic datasets (Shen et al., 2024).

Because ADMIXTURE does not capture recent gene flow and generates highly structured populations (sometimes even between close populations), we utilized fineSTRUCTURE (Lawson et al., 2012) on dataset 2. It is a haplotype-based clustering method that identifies groups of individuals with shared ancestry through recombination-aware modeling. Phase files and recombination maps were provided as input, and coancestry matrices were generated using Chromopainter. The fineSTRUCTURE algorithm was then run using a minimum SNP threshold of 5000 per chromosome (-s1minsnps 5000), 10,000 MCMC iterations for burn-in (-s3iters 10,000), and 10,000 iterations for tree building (-s4iters 10,000). To support fineSTRUCTURE results, we used PopVAE, a variational autoencoder framework and unsupervised deep learning method tailored to population genomics. It projects individuals into a low-dimensional latent space, capturing both linear and nonlinear population structure. We used popvae.py with the following parameters: –seed 42 to ensure reproducibility and –patience 300 to prevent overfitting by using early stopping. The input file was dataset 2, containing representative individuals and high-quality biallelic SNPs. Output latent space coordinates were used to infer population clusters and visualize patterns of divergence and admixture.

2.9. Population genetic parameters and runs of homozygosity

We analyzed population genetic parameters using dataset 2 in VCFtools v.0.1.17 (Danecek et al., 2011), computing nucleotide diversity (π), pairwise genetic differentiation (F_ST), and Tajima's D statistics across 10 kb non-overlapping windows. For π calculations, we included both polymorphic and monomorphic sites. Additional metrics including observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficients (F_IS) were generated using VCFtools' –het function (Danecek et al., 2011). Linkage disequilibrium (LD) patterns were evaluated across populations using PopLDdecay v.3.4.1 (Zhang et al., 2018), which quantified the rate of LD decay. For inbreeding analysis, we detected runs of homozygosity (ROH) in dataset 2 using VCFtools v.0.1.15 (Danecek et al., 2011) with the –LROH parameter, considering only ROH segments exceeding 100 kb. We calculated the genomic inbreeding coefficient (FROH) as the fraction of the genome covered by these long ROH segments (≥ 100 kb) relative to the effective genome size (Yang et al., 2018).

2.10. Identification of deleterious mutations

We identified deleterious variants in Oreocharis mileensis using SIFT4G (Vaser et al., 2016) with plant orthologs from TrEMBL (Boeckmann et al., 2003). Coding SNPs in dataset 1 (low-frequency variants) were categorized as Deleterious (SIFT score < 0.05), Tolerated (SIFT score ≥ 0.05), or Synonymous. Low-confidence and "NA" sites were filtered prior to analysis (Sim et al., 2012).

2.11. Demographic history inference

We reconstructed the demographic history of Oreocharis mileensis through two complementary approaches. First, we applied the Pairwise Sequentially Markovian Coalescent (PSMC) model to analyze a consensus genome. The input data comprised BAM files (dataset 1) converted to 'psmcfa' format with depth thresholds of 10 × (minimum) and 100 × (maximum). PSMC analysis employed standard parameters: time intervals set to '4 + 25*2 + 4 + 6' and 30 optimization cycles (-N). For each sample, we analyzed 500 kb genomic segments with 10 bootstrap replicates. Second, we performed independent demographic inference using Stairway Plot 2 (Liu and Fu, 2020). The folded site frequency spectrum (SFS) was generated with ANGSD v.0.921 (Korneliussen et al., 2014) from dataset 1 using quality filters "-minMapQ 30 -minQ 20". Both analyses incorporated identical evolutionary parameters: a mutation rate (μ) of 3.7 × 10⁻⁹ substitutions/site/year (calculated as Ks/2T, where Ks = 0.12 represents the synonymous substitution rate and T = 16.2 million years reflects the divergence time between O. mileensis and Primulina huaijiensis; Fig. S6). We assumed a generation time of 3 years based on field observations of O. mileensis life history.

2.12. Isolation by distance (IBD), by environment (IBE) and drivers of genetic variation

To evaluate spatial genetic patterns, we examined relationships between genetic differentiation (F_ST/(1-F_ST)) and both geographic and environmental (BIO9, BIO10, BIO13, BIO14) distances. Geographic and environmental distance matrices were computed as Euclidean distances from coordinates. We performed IBD (environmental distance was controlled) and IBE (geographic distance was controlled) analyses using partial Mantel test with statistical significance assessed through 999 permutations. These analyses were conducted separately for putatively adaptive (dataset 4) and neutral (dataset 5) SNP sets. To evaluate the relative contributions of climate, geography, and population structure to genomic variation, we performed partial redundancy analyses (pRDA) on putatively adaptive loci and neutral loci. The genotype matrix (coded 0, 1, 2; missing data as 9) was imputed using most frequent genotypes and standardized prior to analysis. Those four selected climate variables, geographic coordinates (latitude, longitude), and ADMIXTURE-derived population structure were used for this analysis. We constructed a full model incorporating all three predictor sets, along with partial models assessing each component's unique contribution while conditioning on the others. Model significance was evaluated through 999 permutation tests, with variance partitioning used to separate shared and independent effects. All analyses were conducted in R using the vegan package (Dixon, 2003).

2.13. Gene flow, introgression and genetic offset

We utilized TreeMix (Pickrell and Pritchard, 2012) to infer population segregation and migration events, testing models with 1–15 migration edges (m) over 20 iterations. SNP dataset 3 (Rellstab et al., 2015) was used with "-global", "-se", and "-noss" parameters for this analysis. The optimum of migration was determined using OptM v.0.1.5, based on 99.8% variance explained and maximized Δm statistic.

To investigate historical introgression among populations, we applied Patterson's D-statistic and the f4 admixture ratio using Oreocharis hekouensis as the outgroup, following the approach of Patterson et al. (2012). These analyses were conducted using the Dtrios module from the Dsuite package (Malinsky et al., 2021), which allows for inclusion of sites with non-missing genotypes in the outgroup taxa. Statistical significance of the D values was assessed through block jackknife resampling, providing a robust measure of introgression signals.

The R package gradientforest v.0.1–37 was used to predict genomic vulnerability by GF model (Ellis et al., 2012). Dataset 4 of 1737 SNPs from Bayescenv and RDA was compared to 1737 randomly selected SNPs (dataset 3). SNPs with MAF > 10% were converted to population-level MAFs. GF, with 500 regression trees per SNP, modeled BIO9, BIO10, BIO13, and BIO14. GO was calculated as Euclidean distance between current (1970–2000) and future (2081–2100) climates under BCC-CSM2-MR and SSP126/SSP585 scenarios. Local, forward, and reverse offsets (Gougherty et al., 2021) were calculated using bioclimatic variables to assess vulnerability, gene flow potential, and preadaptation.

3. Results 3.1. Genome sequencing and assembly

The libraries sequenced on the PacBio Sequel Ⅱ platform obtained a total of 6.24 million reads with ~91.15 Gb HiFi data (~46×), and an average read length of 15 kb (the longest read was 47 kb, and N50 was 14 kb) (Table S1). A total of 1480 million reads (~221.97 Gb) of Hi-C sequencing data and 1918 million reads (~287.67 Gb, ~144×) of whole-genome short-read sequencing data were produced (Tables S4 and S5). A total of ~11 million reads (~14.65 Gb) of Iso-Seq sequencing data were generated using the ONT platform for third-generation transcriptome (Table S6). Based on PacBio HiFi sequencing data and Hi-C-assisted assembly, we obtained a haplotype-resolved genome.

Oreocharis mileensis genome assembly spans ~3.99 Gb, with exceptional continuity demonstrated by contig and scaffold N50 values of ~124 Mb each (Table S7). The phased assembly comprises two haplotypes: haplotype A (2.01 Gb; contig N50 123 Mb, scaffold N50 124 Mb) and haplotype B (1.98 Gb; contig/scaffold N50 124 Mb) (Tables S8 and S9). This assembly size shows strong agreement with independent k-mer-based estimates from FindGSE and Genomescope (Fig. S7). A total of 99.97% of the sequences from the two haplotypes were anchored onto 34 chromosomes (Fig. 1A and Table S10). Most of the chromosomes were assembled with complete telomeric sequences, and the number of gaps was only 1, indicating extremely high continuity (Table 1). The mitochondrial and chloroplast genomes were assembled into circular DNA molecules of 515,902 bp and 154,542 bp, respectively (Table S10). Evaluation of sequencing depth demonstrated comprehensive genome coverage, with strong mapping efficiency metrics (Table S11). The assembly demonstrated strong completeness with 1588 (98.3%) complete BUSCO groups out of 1614 groups, including 12 (0.7%) single-copy and 1576 (97.6%) duplicated genes. Few genes were fragmented (0.7%) or missing (1.0%) (Table S5; Figs. 1B, S2 and S3). A total of 56,834 genes and 130,476 transcripts were also assembled from the Iso-seq reads, with average lengths of 2606 and 2554, respectively (Table S12).

3.2. Genome annotation

The Oreocharis mileensis genome contained 3.46 million repetitive sequences spanning ~3.16 Gb (79.05% of the genome). Long terminal repeats (LTRs) dominated these repeats, comprising ~2.90 Gb (72.69%), with Gypsy (41.84%) and Copia (18.98%) elements being the most abundant LTR families (Table S13). After the integration of annotations, 89,009 protein-coding genes, 114,466 transcripts and 114,466 coding sequences were identified, with mean lengths of 3719, 1727 and 1135, respectively (Table S14). In addition, 3111 rRNA, 2211 tRNA, and 7324 ncRNA sequences were identified (Table S15). Protein annotation completeness was evaluated using BUSCO, showing 98.5% core gene coverage (97.1% duplicated, 1.4% single-copy). Only 0.4% were fragmented and 1.1% missing (Table S7).

A total of 100,655 genes in both haplotypes were annotated using different software tools. The annotation results based on eggNOG-mapper, including GO (33,188, 37.29%), KEGG_KO (31,159, 35.01%), EC (12,998, 14.60%), KEGG_Pathway (19,457, 21.86%), eggNOG (62,202, 69.88%), and COG (66,381, 74.58%). The annotation results based on DIAMONG included Swiss-Prot (48,892, 54.93%), TrEMBL (76,802, 86.29%), NR (67,152, 75.44%) and Arabidopsis thaliana (58,924, 66.20%). The annotation results based on InterProScan included Pfam (55,121, 61.93%), SMART (21,007, 23.60%), Gene3D (45,672, 51.31%), and CDD (21,622, 24.29%) and others (Table S16).

3.3. Analysis of phylogeny and WGD

We used haplotype A of Oreocharis mileensis, consisting of 44,569 genes, for phylogeny and WGD events. OrthoFinder2 analysis of 19 genomes identified 65,638 genes grouped into 22,933 orthologous families, with O. mileensis showing significant gene family evolution (2519 expanded and 1539 contracted; Fig. S8). Phylogenetic reconstruction using 1273 conserved orthogroups revealed divergence times within Gesneriaceae: O. mileensis split from O. esquirolii and O. maximowiczii at ~8 (6.3–9.9) Mya, Oreocharis split from Primulina at ~14.7 (12.1–17.8) Mya and from Streptocarpus at ~19.4 (16.2–22.9) Mya, with all three genera diverging from Boea earlier at ~22.9 (19.3–26.8) Mya (Fig. 1C).

A total of 23,126 post-filter collinear gene pairs and 1224 post-filter collinear blocks were inferred within the Oreocharis mileensis genome (Table S17). Unlike most core eudicots, Vitis vinifera and Coffea canephora appear to have retained their diploid state following the ancestral γ hexaploidization, showing no evidence of lineage-specific whole-genome duplications (WGD). Therefore, these two species are used as references to investigate the relative ploidy levels and inferring the WGD events in other species. O. mileensis and the closely related species Primulina huaijiensis, and Streptocarpus rexii show an obvious 1:8 relationship in collinearity compared to the reference species (Coffea canephora), which may be the result of three lineage-specific tetraploidization events (Fig. S9). The synteny ratio between O. mileensis and O. esquirolii, O. mileensis and O. maximowiczii, O. mileensis and P. huaijiensis retained its diploid state, indicating that these species share all WGD events (only O. mileensis and P. huaijiensis synteny plot is provided in Fig. S6). The synteny ratio of O. mileensis/P. huaijiensis to Streptocarpus rexii is 2:2, indicating that each species may have experienced a lineage-specific tetraploidization event (a deduction that also matches the peaks of Ks values), meaning that they share the first two tetraploidization events (Fig. S10). When comparing O. mileensis with itself or with P. huaijiensis, there is first a recent (Ks ≈ 0.2) duplicated collinear block, which can be confirmed as a tetraploidization event (Fig. S11). Each chromosome roughly corresponds to two older collinear blocks (with Ks ≈ 0.8) and approximately four even older collinear blocks (with Ks ≈ 1.2), which correspond to two more ancient tetraploidization events.

3.4. Population structure and phylogeny

Whole-genome resequencing generated ~11.12 Gb (15.68× coverage) of high-quality data per sample, with 97.25% of filtered reads (Q20: 96.68%; Q30: 94.82%) mapping to the reference (Table S18). Population analyses using two SNP datasets (all and neutral loci) revealed optimal genetic structure at K = 8 through ADMIXTURE (Fig. S12A and B). The results from these two datasets consistently indicate a highly structured population, with similar patterns of genetic differentiation across the groups (Figs. 2A, S13A and S14A). DHB1 and DHB2 exhibited identical genetic compositions, suggesting they constitute a single population. Furthermore, two populations (XY and DHB) shared genetic profiles, with only a slight genetic admixture observed from an additional population (QT). Minor genetic mixing was also detected between two other populations (MXJ and YL), whereas the remaining populations had a significant genetic structure. For putatively adaptive loci, the best K was 7, merging DHB1, DHB2 and XY as the same component of ADMIXTURE and MXJ and YL as well (Fig. S15A). Due to this pronounced genetic differentiation, PCA was conducted and visualized for components 1 to 6 across all three datasets (Figs. 2B, S13B, S14B and S15B). The PCA results revealed a clear separation into four distinct clusters in PC1 and PC2 axes with 45.66%, 45.61%, and 44.57% of genome covariance for all loci, neutral loci, and putatively adaptive loci, respectively. These clusters were identified as follows: Cluster 1 – DHB1, DHB2, QT, and XY; Cluster 2 – GB, GS, MXJ, and YL; Cluster 3 – EJIA; and Cluster 4 – WS. Further analysis demonstrated that GB and GS were separated from each other in PC1 and PC3 in all datasets. Similarly, MXJ and YL were distinguished in PCA 1 and PCA 5 in the all and neutral loci datasets (Fig. S14B), while EJIA and WS were differentiated in PC 1 vs. PC 4, 5, and 6. The results of fineSTRUCTURE were consistent with these findings (Fig. S15C and D). The GS/GB/YL/MXJ, DHB/QT/XY, WS and EJIA populations formed separate clusters, however, co-ancestry within these clusters was found as weak, highlighting ancient not recent geneflow. popVAE result was similar to previous analyses, however, QT was separated from DHB and XY while GS was separated from GB/MXJ/YL populations (Fig. S15E).

Maximum-likelihood phylogeny showed strong concordance with ADMIXTURE results, where populations sharing genetic ancestry clustered together phylogenetically (Fig. 2C). In brief, DHB, XY and QT, MXJ and YL, GB and GS were clustered together, while EJIA and WS were both distinct from all other populations.

3.5. Parameters of population genetics

Nucleotide diversity (π), Tajima's D, inbreeding coefficient (F_IS), observed heterozygosity (Ho), and expected heterozygosity (He) were calculated using VCFtools v.0.1.17 (Table S19 and Fig. S16). π values ranged from 0.000285 (GS) to 0.000422 (WS), indicating low diversity across populations. Tajima's D ranged from 0.0829 (MXJ) to 0.53 (WS), while F_IS values were high (0.59 in WS to 0.86 in GS). He values were almost consistent, ranging from 0.2538 (GS, XY) to 0.2560 (WS). Genetic differentiation (F_ST) ranged from 0.20 (MXJ-YL) to 0.77 (EJIA-GS), showing high differentiation (Fig. 3A and Table S20).

We next assessed genome-wide runs of homozygosity (ROH), which reflects recent inbreeding through extended homozygous segments. Following the method outlined by Robinson et al. (2022), we estimated the number of generations since the most recent common ancestor (g) using the formula g = 100/(2*L), where L represents the mean ROH length in megabases (Mb). Inbreeding events were estimated to occur between 250 (GS) and 1322 (WS) generations ago (Fig. S17). FROH was highest in GS (median = 0.1038) and lowest in WS (median = 0.0008) (Fig. 3B). Populations with severe inbreeding had longer ROH (> 1 Mb), with GS showing the highest proportion (3.86% of the genome), while WS and XY had none (Fig. 3C and Table S21).

Linkage disequilibrium (LD) decay varied across populations, with MXJ showing the slowest decay and EJIA, QT, and WS the fastest (Fig. 3D). Based on SIFT prediction approach, we detected a total of 327,684 deleterious mutations, 624,179 tolerated mutations and 859,658 synonymous mutations (Table S22). Our results showed that the number of homozygous deleterious sites to total deleterious mutations in the MXJ and GS populations were significantly higher than in other populations, suggesting that these populations had a higher homozygous genetic load (Fig. 3E). YL and GB populations also had more significant genetic load compared to other populations.

3.6. IBD, IBE, drivers of genetic variation, and gene flow

We identified candidate genomic variants associated with climate adaptation. BayeScEnv analysis detected 6348 loci associated with four selected environmental variables, whereas RDA identified 5368 loci across five axes. To minimize false positives, 1737 loci overlapping in both analyses were classified as environment-associated (Table S23). Of those 1737 loci, 201 correlated with Mean Temperature of Driest Quarter, 1332 correlated with Mean Temperature of Warmest Quarter, 184 correlated with Precipitation of Wettest Month, and 189 correlated with Precipitation of the Driest Month. The divergence pattern was tested by IBD and IBE analyses using partial Mantel test. According to IBD and IBE analyses, geographic distance did not explain genetic variation (r = 0.173, p = 0.26 for neutral loci and r = 0.21, p = 0.13 for putatively adaptive loci), but environmental distance was significant (r = 0.65, p = 0.002 for neutral loci and r = 0.21, p = 0.01 for putatively adaptive loci; Fig. 4A and B).

To assess the relative contributions of climate, geography, and population structure to genetic variation, we performed pRDA analyses on both putatively adaptive and neutral loci (Table S24). For putatively adaptive loci, the full model (climate + geography + structure) explained 87.87% of genomic variation, with population structure alone accounting for 87.18%, climate 72.02%, and geography 31.13% when controlling for other factors. The first two RDA axes captured 79.01% of constrained variation, and permutation tests confirmed all predictors were significant (p = 0.001). Neutral loci showed weaker overall associations, with the full model explaining 63.85% of variation. Population structure explained most variation (61.02%), while climate (40.48%) and geography (18.94%) had reduced but still significant effect (p = 0.01). TreeMix analysis identified three migration events (99.986% variance explained), with the optimal model (iteration 11) showing minimal residual error (Tables S25 and S26; Fig. S18). TreeMix indicated limited and weak gene flow or admixture among O. mileensis populations (Figs. 4C and S19). Three weak gene flow events were predicted: between YL and GB, XY and the ancestors of YL, MXJ, GB, and GS, and between GS and the ancestors of WS and EJIA. WS and EJIA were genetically distinct from other populations.

The D-statistics analysis revealed distinct patterns of gene flow for EJIA and WS, indicating that each has undergone introgression with different genetic clusters (Fig. 4D). EJIA showed significant but minor introgression with GB (p < 0.001, Z > 3, Table 2). Also, gene flows from YL to GB (p < 0.001, Z > 3), and from GB to MXJ (p < 0.001, Z > 3), were determined in f branch. In addition, WS exhibited gene flow signals with the DHB/QT/XY cluster. We found significant but modest introgression between WS and QT, with WS likely receiving gene flow from QT or a QT-related lineage. In addition, comparisons such as XY-DHB-QT showed significant but minor introgression.

3.7. Genomic vulnerability to climate change

GF modeling showed higher genomic vulnerability for all populations (Fig. 5A–D). Putatively adaptive SNPs exhibited higher GO than reference SNPs under both scenarios (ssp585 and ssp126), indicating greater sensitivity to climate change. Analysis of local, forward, and reverse offsets using an RGB model indicated that northeastern regions showed were more vulnerable despite potential migration (Fig. 5E and F). However, for other populations, local offset was higher than forward and reverse offsets suggesting urgent need for assisted gene flow or habitat restoration. Functional annotation of outlier SNPs predicted functions for 276 genes (p < 0.05; Table S27). These genes were assigned to categories, including 25 biological processes, such as ion balance maintenance, oxidative damage reduction, signal transduction, and metabolism. Further, Cell Protection & Repair were related to stress tolerance and Cell Expansion & Hydration Recovery, Metabolic Reactivation and Hormonal & Transcriptional Reactivation were related to rapid growth after hydration.

3.8. Demographic history of Oreocharis mileensis

Because PSMC is primarily effective at reconstructing ancient demographic histories and Stairway Plot is particularly effective at resolving recent demographic changes, we decided to use both approaches. PSMC analysis identified two severe population declines, with continuous reductions in effective population size (Ne) and no recovery (Fig. 6A). The first decline occurred ~0.6 Mya during the Mid–Pleistocene Transition, reducing Ne to ~15.4 × 10⁴. After a small stable period, Ne declined (only one population had a recovery followed by a sharp decline) again between 0.67 and 0.2 Mya. Stairway Plot analysis revealed a continuous decline in Ne across all populations, punctuated by some sharper reductions (Figs. 6B and S20). Notably, the XY population exhibited a stable Ne trend ~100,000 years ago (Fig. S20), though this anomaly may reflect bias rather than true demographic stability. The first major decline (~0.86 Mya), reducing Ne to 1.08 × 10⁶, aligned with the Mid–Pleistocene Transition, a period of intensified glacial cycles and climatic instability (Fig. 6B). The second decline (~50,000 years ago) reduced Ne to ~1.42 × 10⁵, aligning with the Late Pleistocene and Last Glacial Maximum. Further declines occurred ~6000 years ago (Ne ~ 29,600) at the end of the Holocene Climate Optimum and ~3000 years ago (Ne ~9400) during the Late Holocene Neoglacial period. Finally, Ne dropped to ~950 by 114 years ago, coinciding with the early Industrial Era and human-driven environmental changes.

4. Discussion 4.1. High-quality genome assembly of Oreocharis mileensis advances conservation genomics and evolution

Conservation genomics depends on high-quality genetic data and whole-genome sequencing provides unparalleled resolution for studying evolutionary processes and population dynamics (Barrick and Lenski, 2013; Jing et al., 2023; Song et al., 2023). Despite its potential, application remains limited for many plant groups, particularly karst-adapted species (Feng et al., 2020; Cao et al., 2023; Peng et al., 2025). In this study, we assembled the genome of Oreocharis mileensis, a resurrection plant endemic to karst forests in southern China. The high-contiguity assembly exhibits excellent gene completeness and sequencing coverage, establishing a foundation for conservation prioritization. The Oreocharis assembly also provides a reliable high-quality reference genome for the study of the drought resistance and resurrection processes of Oreocharis species in karst regions in the near future.

Our assembly shows that the Oreocharis mileensis genome (1.9 Gb) was larger than that of both O. esquirolii (1.24 Gb) and O. maximowiczii (1.78 Gb) (Peng et al., 2025). Although O. mileensis had fewer genes (44,569) than O. esquirolii (47,637), our structural analysis of genomes showed that O. mileensis has evolved more genes and lost fewer genes than relatives within the clade. Research has shown that species with more gene gains than losses frequently exhibit improved capacities in biological processes important to key adaptations (Jacquemin et al., 2014). We have yet to determine whether gains and losses in O. mileensis genes are adaptive.

The genus Oreocharis is likely to have experienced three lineage-specific tetraploidization events, with the divergence time between Oreocharis spp./Primulina huanjiangensis and Streptocarpus rexii (22.9–16.2 Mya) being consistent with the most recent genome duplication event shared by plants in the Didymocarpoid subtribe (24.2–20.6 Mya) (Feng et al., 2020).

4.2. Inbreeding and possible selective purging shape diversity in Oreocharis mileensis populations

We found that Oreocharis mileensis populations have long ROHs, suggesting that these populations are vulnerable due to recent inbreeding and substantial genetic load. ROH patterns directly reflect increased homozygosity (Kirin et al., 2010), while the resulting accumulation of homozygous deleterious mutations creates a compounding fitness burden (Bortoluzzi et al., 2020; Kyriazis et al., 2021). Thus, our findings are consistent with previous studies that indicate small populations face elevated extinction risks from inbreeding depression, where restricted mate availability increases homozygosity of deleterious alleles (Bortoluzzi et al., 2020). We also found that nucleotide diversity is lower in O. mileensis than in other threatened herbaceous plant species (including O. esquirolii and O. maximowiczii), confirming that increased homozygosity has reduced genetic diversity in these populations (Lin et al., 2024; Peng et al., 2025) (Table S28). Several factors have led to frequent inbreeding within populations of O. mileensis. Specifically, inbreeding is more common when populations are isolated, located in independent limestone hills and/or fragmented habitats. Our ROH analysis identified four core O. mileensis populations that show signs of recent inbreeding (i.e., GS, YL, GB, and MXJ). For example, the GS population exhibited high ROH, a high inbreeding coefficient (F_IS ≈ 0.86), a significant burden of homozygous deleterious mutations, and notably low genetic diversity, consistent with inbreeding events occurring approximately 250 generations ago.

Interestingly, O. mileensis populations on the periphery of karst forests (e.g., QT) displayed a unique genetic profile: a lower proportion of ROH, fewer homozygous deleterious mutations, and reduced nucleotide diversity (π = 0.000312), yet a high F_IS value (≈0.85). This pattern may result from a combination of reduced recent inbreeding, possible genetic purging of deleterious mutations, historical demographic events, and selective pressures (Cvijović et al., 2018; Grossen et al., 2020). The absence of long ROH segments in QT suggests that recombination may have fragmented these regions, rendering them undetectable by ROH analysis. The EJIA population exhibited similar genetic indicators, further supporting this interpretation.

Other peripheral populations (WS, DHB, and XY) demonstrated higher genetic diversity but similarly lower ROH, and fewer homozygous deleterious mutations. The reduced burden of deleterious mutations in these populations may reflect the possible effective purging of harmful alleles over generations, potentially driven by natural selection. These peripheral populations also occupy basal positions in the phylogeny, suggesting they diverged earlier and potentially experienced habitat fragmentation and population isolation before the core (GB/GS/MXJ/YL) populations. Such early fragmentation may have allowed more time for possible purifying selection and recombination to act, reducing the detectability of long ROH segments today.

4.3. Drivers of genetic variation and their role in formation of the extremely small and isolated populations

Partial redundancy analysis (pRDA) revealed striking differences in the drivers of genetic variation between adaptive and neutral loci. For putatively adaptive loci, the variation was overwhelmingly aligned with pre-existing population structure (87.18% of explained variation). This strongly suggests that historical demographic processes such as long-term environmental isolation, genetic drift in small populations, and barriers to gene flow have been the primary architects of the current distribution of adaptive genetic diversity, exceeding the direct influence of contemporary climate (72.02%) or geographic distance (31.13%). These findings align with studies in other plant systems (Yang et al., 2025), where neutral processes creating structure often provide the substrate upon which selection acts. For neutral loci, the overall model fit was weaker, but the pattern held the largest proportion of explained variation (61.02%) was associated with population structure. This indicates that the same historical demographic forces (isolation and drift) have also been the principal drivers of neutral genetic divergence. The significant but smaller contributions of climate and geography likely reflect isolation-by-environment and isolation-by-distance (although IBD was not significant by partial Mantel test, pRDA analysis found its effect as well) processes operating within this overarching demographic history (Mantel et al., 2003; Sexton et al., 2014). Therefore, our results demonstrate that the complex demographic history of the populations was likely driven by the complex karst landscape.

PSMC analysis identified two severe declines in Oreocharis mileensis: the first (~1–0.6 Mya) during the Mid–Pleistocene Transition reduced Ne to ~15.4 × 10⁴, and the second (0.67–0.2 Mya) further decreased Ne. Unlike PSMC, Stairway Plot revealed continuous Ne decline, with more reductions tied to major climatic and anthropogenic events: the Mid–Pleistocene Transition (~0.86 Mya), Late Pleistocene, Holocene Climate Optimum, Neoglacial period, and early the Industrial Era leading to the lowest Ne of ~950. This temporal complementarity between methods was expected, as PSMC provides robust estimates of ancient demographic changes while Stairway Plot offers finer resolution of recent population dynamics.

TreeMix results indicated that gene flow was minimal between two geographically close populations of Oreocharis mileensis (YL and GB), confirming these populations are genetically isolated regardless of distance. Dsuite analyses indicated that all O. mileensis populations experienced minor introgression. We also found that the source of introgression into the DHB/QT/XY cluster was the WS population, while the source of introgression into the GB/GS/MXJ/YL cluster was the EJIA population. These two populations (WS and EJIA) appear to have diverged first from the common ancestor, with subsequent radiation of the remaining populations occurring long ago. Additionally, the absence of recent gene flow and the effects of genetic drift have led to highly structured populations. This is supported by positive Tajima's D values across populations (Table S19), which suggest recent population contractions, reduced genetic diversity, and the formation of small, isolated populations (Tajima, 1989).

These results show that limited gene flow, historical demography, past environmental fragmentation, and recent anthropogenic pressures have together shaped the deeply structured, small, and isolated populations seen in Oreocharis mileensis today.

4.4. Local adaptation alleles as drivers of climate change-induced genomic vulnerability

Populations within the same species often exhibit differential responses to climate change due to local adaptation to heterogeneous environmental conditions (O'Neil et al., 2014; Fitzpatrick and Keller, 2015; Cortés, 2017). Understanding these patterns of local adaptation is crucial for predicting how species may respond to future climatic shifts. In this study, we quantified the GO by calculating the distance between current and future climate conditions and integrating intraspecific variation into the GF model. We found all populations O. mileensis exhibited significantly higher GO under both future climate scenarios. Our study also showed that under the ssp126 scenario, core populations (GB, GS, MXJ and YL) face significant adaptive risks, while under both climate change scenarios (ssp126 and ssp585) the local offsets were higher than forward and reverse offsets, suggesting that assisted migration or gene flow may to some extent enable adaptation to future climate change.

Shorter generation times could theoretically enhance adaptive capacity by accelerating allele frequency shifts in response to selection (Ojeda et al., 2016; Campbell-Staton et al., 2017). However, the persistence of this advantage depends on the rate and magnitude of climate change relative to the species' evolutionary potential. If climatic shifts outpace the speed of adaptive evolution, even rapidly reproducing species may face significant genomic vulnerability, as suggested by the high GO values in our study. Furthermore, while herbaceous plants may adapt more quickly than trees, their shorter lifespans and smaller effective population sizes could also increase stochasticity in adaptive trajectories (Willi et al., 2006), potentially leading to divergent outcomes from model projections. In addition, demographic costs of selection, such as elevated mortality or reproductive failure can outpace adaptive responses, causing populations to go extinction before achieving evolutionary rescue (Gomulkiewicz and Holt, 1995). Thus, while our GO models provide critical insights into near-term vulnerability, the dynamic interplay between generation time, selection, and genetic architecture warrants caution when extrapolating results over longer timescales. Future studies could integrate demographic modeling with genomic offset frameworks to better resolve how generation time mediates climate adaptation in herbaceous systems.

The physiology of resurrection plants like Oreocharis mileensis enables them to survive during severe desiccation and to regenerate when a sufficient level of water is in the environment (Li et al., 2014). Our search for adaptive loci indicated that putatively adapted SNPs are mainly located in genes responsible for ion balance maintenance, oxidative damage reduction, signal transduction, and metabolism, as well as cell protection and repair. These categories of genes correspond to complex survival strategies used by resurrection plants during dehydrated states caused by the absence of seasonal rainfall (Dinakar and Bartels, 2013). We assume that genes in these categories play key roles in the adaptation of O. mileensis to climate change. We also found that some putatively adaptive SNPs are located in genes related to cell expansion and hydration recovery, as well as metabolic reactivation and hormonal & transcriptional reactivation, which are important for rapid growth of the plant during the presence of rainfall.

4.5. Implications for conservation planning and management practices

Like many karst-dwelling species in Yunnan, Oreocharis mileensis is strongly affected by human activity (Chen et al., 2014). Our field surveys revealed that overgrazing by goats affects several populations of O. mileensis, including the DHB population, adjacent populations located within a nature reserve, and a population in the in situ conservation site established by the local government department (GS). Thus, effective conservation guidelines and management strategies must be developed based on the establishment of well-defined management units (MUs) (Barbosa et al., 2018; Yang et al., 2022; Cai et al., 2024). Our findings lead us to recommend designating eight MUs of O. mileensis. Populations that have experienced inbreeding require genetic rescue through assisted gene flow (Pavlova et al., 2017). In the case of O. mileensis, GS, GB, MXJ and YL populations have low nucleotide diversity, high inbreeding, and high genetic load. These populations appear to be the core populations over the entire geographical distribution of Oreocharis mileensis. Further, genomic offset revealed that these populations may experience higher allele frequency change due to climate change. Gene flow among populations was almost absent, even between geographically close populations. Considering the presence of homozygous deleterious mutations, we recommend facilitating assisted gene flow between MXJ (fewer shared deleterious mutations) or YL (higher genetic diversity) and GB (Fig. S21). Due to high genetic divergence (F_ST = 0.407), assisted gene flow between GB and YL should be conducted under controlled conditions. Our findings also indicate that the GB population can serve as a pollen source to rescue GS. The geographic and environmental proximity of these populations may enhance the likelihood of success (Frankham et al., 2011), although experimental trials are needed to evaluate outbreeding depression risks given their even greater divergence (F_ST = 0.485). Unfortunately, other peripheral populations (i.e., DHB, QT and XY) are isolated from the core populations (F_ST at least 0.637); thus, intercrossing these populations may be risky and lead to outbreeding depression. We therefore do not recommend intercrossing between peripheral and core populations. The QT population is located in the edge of the overall distribution of O. mileensis. The genetic make-up of XY and QT populations is similar and share fewer shared homozygous deleterious mutations. Therefore, we recommend conserving QT population and strengthening its genetic makeup using assisted gene flow from XY population to allow it to adapt to climate change. However, this also should be conducted under control as mentioned above.

The peripheral EJIA and WS populations are both isolated from all other populations and from each other as well. Due to the high risk of outbreeding depression, assisted gene flow should not be attempted here. For these populations, we first recommend protecting and restoring their natural habitats to mitigate external threats such as habitat degradation, overgrazing, or human encroachment. However, this may not be sufficient to maintain their resilience (Belovsky et al., 1994). Thus, genetic material should be preserved through collecting seeds and storing them in germplasm bank and cultivating live plants as ornamental in botanical gardens (Breman et al., 2021).

Acknowledgments

This work was equally supported by the National Key R&D Program of China (Grant No. 2024YFF1307400) and Yunnan Key R&D Program (Grant No. 202403AC100028), Yunnan Fundamental Research Projects (Grant No. 202301AT070318) and Yunnan Revitalization Talent Support Program "Young Talent" Projiect. We thank Dr. Yang Liu, Dr. Shi-Wei Guo, Dr. Ji-Dong Ya, Prof. Zhi-Kun Wu, Prof. Xin-Xiang Bai, Prof. Pan Li, Mr. Sheng-Hu Tang, Mr. Jian Xu, Mr. Bo Pan and Mr. De-Ming He, for their assistance in collecting materials and providing species distribution information, and we thank Dr. Feng-Mao Yang, Dr. Yang Liu, Dr. Bishal Gurung, Dr. Yu-Hang Chang and Ms. Na–Na Peng for their help in analyzing and providing data.

CRediT authorship contribution statement

Temur Asatulloev: Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Lei Cai: Writing – original draft, Visualization, Resources, Investigation, Funding acquisition, Formal analysis. Ziyoviddin Yusupov: Writing – review & editing. Kai-Hua Jia: Formal analysis. Ren-Gang Zhang: Formal analysis. Komiljon Sh. Tojibaev: Supervision, Writing – review & editing. Wei-Bang Sun: Supervision, Resources, Project administration, Funding acquisition, Conceptualization, Writing – review & editing.

Data availability statement

The data that supports the findings of this study are available in the supplementary material of this article and whole genome sequences and annotations are located under PRJCA037867 bioproject (https://ngdc.cncb.ac.cn/gsa/s/h9annOxq).

Declaration of competing interest

The authors Weibang Sun and Komiljon Sh. Tojibaev are editors for Plant Diversity and were not involved in the editorial review or the decision to publish this article. The other authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

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