1. Introduction

Understanding how organisms adapt to environments and undergo speciation has been a central focus of evolutionary biology for nearly two centuries (Tobler et al., 2018). Recent studies have examined how species adapt to a variety of habitats, such as high altitudes (An et al., 2024; Long et al., 2025), low temperatures environments (Zhang et al., 2020), deep-sea high pressures ecosystems (Wang et al., 2019), and nutrient-deficient habitats (He et al., 2020; Shi et al., 2025). Adaptation to such habitats often drives divergence among closely related lineages, with accumulated genetic variants ultimately triggering speciation (Todesco et al., 2020; Nosil et al., 2021; Geng et al., 2025). Contemporary biologists widely concur that most speciation events require an allopatric phase (Bolnick and Fitzpatrick, 2007; Anderson and Weir, 2022; Feng et al., 2025), where geographic isolation or ecological divergence reduces gene flow between populations (Qian et al., 2021). Natural selection accelerates this process, as heterogeneous environments exert strong selection pressures that prompt adaptive divergence and speciation (Yona et al., 2015). Allopatric speciation is thus recognized as the dominant mode for many organisms (Bian et al., 2020; Anderson and Weir, 2022; Cao et al., 2023). Therefore, elucidating species divergence processes and their adaptations in specialized habitats is essential for understanding speciation mechanisms and the origins of Earth’s biodiversity.

Since Wallace’s seminal work in the Malay Archipelago, limestone karst habitats have played a pivotal role in understanding speciation and adaptation (Wallace, 1858). These habitats are characterized by elevated levels of calcium (Ca²⁺) and magnesium (Mg²⁺), higher pH values, reduced water retention capacity, and lower availability of essential plant nutrients (Hao et al., 2015; Geekiyanage et al., 2019). Together, these factors impose strong selective pressures on plants, fostering high species endemism, particularly in the karst regions of southern China (Chung et al., 2014). Therefore, karst environments are regarded as ‘natural laboratories’ for exploring species divergence and adaptive evolution (Clements et al., 2006; Oliver et al., 2017; Chen et al., 2024). The elevated calcium and reduced water retention of karst environments have been shown to drive adaptive evolution in several plant species (Jin et al., 2018; Zhou et al., 2023). Tao et al. (2016) demonstrated key genetic adaptations enabling plants to thrive in karst environments, including the (ROS)-responsive Ca²⁺ channel gene TPC1 contributing to Primulina adaptation. Similarly, rapid expansion of the WRKY gene family facilitated adaptation in Primulina huaijiensis (Feng et al., 2020). Adaptive genetic signals linked to karst survival have also been identified in other taxa such as Lonicera confusa (Jin et al., 2018) and Marsdenia tenacissima (Zhou et al., 2023). However, although these studies investigate speciation and adaptation through comparisons between karst and non-karst species or by examining individual species, a significant gap remains in understanding the differentiation and adaptive evolution of closely related species within karst environments. The lack of genomic investigations into such species differentiation and adaptive evolution in karst systems continues to impede our understanding of species’ evolutionary patterns and complicates conservation efforts for endemic and endangered species in karst environment.

Urophysa Ulbr. (Ranunculaceae) is an endemic genus comprising only two species, U. rockii and U. henryi, distributed exclusively in the limestone karst environment of southern China, primarily on the Yunnan-Guizhou Plateau and adjacent regions (Xie et al., 2017). The two species are distinguished by key floral morphology: Urophysa rockii develops nectariferous spurs on its petals, whereas U. henryi forms saccate structures (Fig. 1). Ecologically, U. rockii is restricted to the western Sichuan Basin, with approximately 2000 wild individuals sustained across five disjunct populations (Wang and He, 2011). Its limited population size and narrow distribution have led to its classification as Critically Endangered (CR) on both the China Biodiversity Red List and the IUCN Red List. In contrast, U. henryi inhabits a broader yet fragmented range in southern karst landscapes, with populations isolated by limestone cave systems and topographic barriers (Xie et al., 2021). Tourism development and anthropogenic disturbance threaten the species, resulting in its Vulnerable (VU) status on the IUCN Red List. Given their regional endemism, high habitat specificity, and allopatric distributions, Urophysa species serve as excellent models for studying species differentiation and evolutionary adaptation to karst limestone habitats, and for exploring the mechanisms of survival and endangerment in such unique environments.

In this study, we investigated the evolutionary divergence and speciation of Urophysa rockii and U. henryi, assessed the conservation status of U. rockii, and identified mechanisms underlying their adaptation to karst habitats. For these purposes, we generated a newly sequenced genome for U. rockii and performed whole-genome sequencing of 97 individuals from 14 populations across the genus. We also examined the U. rockii genome for signatures of positive selection and tested whether candidate gene regulation responds to calcium, a critical environmental factor in karst environments.

2. Materials and methods 2.1. Genome sample collection, sequencing and de novo assembly

Fresh leaves were collected from a wild Urophysa rockii individual in Pengzhou City, Sichuan province, China (31.194°N, 103.913°E, 1087 m). The voucher specimen was deposited in the Sichuan University Herbarium (SZ) under accession numbers DF220117-1 to DF220117-2. For genome sequencing, paired-end short reads were generated using an Illumina HiSeq X platform, and long reads were generated using the PacBio Sequel Ⅱ platform. The Hi-C data were collected for chromosome-level assemblies. Genome size was estimated from k-mer analysis, and sequencing depth was calculated based on the ratio of the total number of bases sequenced to the estimated genome size. This process utilized > 70× coverage of paired high-fidelity (HiFi) long reads and > 230× coverage of Hi-C data in sequencing depth, all with a Q30 score > 93%. Based on our preliminary experiments and previous studies, Urophysa rockii is confirmed to be a diploid species, with karyotype formula of 2n = 2x = 14 (Zhang et al., 2013, 2019b).

The genome was assembled using the following pipeline: a primary assembly generated from filtered PacBio reads with HiFiasm v.0.16.1 (Cheng et al., 2021), followed by three rounds of error correction using Illumina short reads with Pilon v.1.24 (Walker et al., 2014) to produce a draft genome. Finally, Hi-C reads were aligned to the draft genome using BWA-MEM2 v.2.0pre2 (Vasimuddin et al., 2019), and contigs were clustered into chromosomes by Lachesis based on interaction data, generating a chromosome-level reference genome of U. rockii. The completeness of the final genome assembly was assessed using Benchmarking Universal Single-copy Orthologs (BUSCO) v.3.1.0 (Simão et al., 2015). Genome gene prediction, including repeats and gene function, was assisted by RNA extracted from fresh leaves, flowers and stems from the same individual. Detailed methods can be found in Supplementary Methods Appendix S1.

2.2. Gene family cluster, phylogenetic reconstruction, and divergence time estimation

To perform gene family cluster analysis, the protein-coding sequences of 11 other species were downloaded from Phytozome v.12, NCBI, and the Beijing Institute of Genomics (BIG) Data Center. We employed two strategies for taxon selection (see Supplementary Methods Appendix S2). OrthoFinder v.2.4.1 (Emms and Kelly, 2019) was employed to identify orthologues, and 1009 single-copy orthologous genes (SCGs) were detected and used in both concatenation and coalescent-based phylogenetic reconstruction. Divergence times were estimated with three calibration constraints based on the single-copy orthologous genes identified. WGDI (Sun et al., 2022) was applied to detect whole-genome duplication (WGD) events. CAFE v.4.2.1 (Han et al., 2013) was employed to deduce gene family expansion and contraction along all phylogenetic branches. Genes that had undergone significant expansion were further annotated. Detailed information on orthologues identification, phylogenetic reconstruction, divergence time estimation, WGD events, gene family expansion and contraction, and candidate selection regions detection is provided in Supplementary Methods Appendix S2.

2.3. Population genome resequencing and SNP calling in Urophysa species

To better understand how Urophysa rockii populations diverged from U. henryi, we first re-sequenced and characterized the genomes of 97 individuals from 14 natural populations across the entire distribution range of the two Urophysa species, including 45 samples from 5 populations of U. rockii, and 52 samples from 9 populations of U. henryi. To minimize genetic overlap given restricted population size, individuals were sampled at least 10 m apart. For outgroups, we used Semiaquilegia adoxoides (genomic data from one collected individual) and four Aquilegia specimens downloaded from the National Genomics Data Center (NGDC) (Table S1).

Genomic DNA was extracted from leaf samples using the Qiagen DNeasy Plant Kit. Whole-genome paired-end sequencing was then performed on the Illumina NovaSeq 6000 platform, targeting an average coverage of 25× for each individual. We used Trimmomatic v.0.367 (Bolger et al., 2014) to remove adapter sequences and trim bases from the start or end of reads where the quality score was below 20. Reads shorter than 36 bases in length were discarded. The resulting high-quality reads were aligned to our de novo assembled U. rockii genome using the BWA-MEM algorithm from BWA v.0.7.17 with default parameters (Li and Durbin, 2009). The alignment results were processed with SAMtools v.1.9 and Picard v.2.18.11 (http://broadinstitute.github.io/picard/) for sorting and marking PCR duplicates.

The Genome Analysis Toolkit (GATK v.3.7) (DePristo et al., 2011) was used for single nucleotide polymorphism (SNP) calling, and ANNOVAR (Wang et al., 2010) was used for variant annotation and the prediction of the effect of variants on gene function. We applied a strict filtering process to reduce low-quality SNPs, as follows: (ⅰ) Only those genotypes with a Phred-scaled likelihood of 20 or above were retained, ensuring a genotyping accuracy rate of at least 99%; genotypes with lower quality scores were treated as missing. (ⅱ) To maintain uniformity, only biallelic SNPs were retained for population-level analysis. (ⅲ) To reduce the incidence of false positives, SNPs found within INDEL regions, including a 5-base pair flanking area, were excluded. (ⅳ) SNPs with a non-reference allele frequency below 1% were removed to concentrate on more common variants. (ⅴ) Additionally, sites with a depth of coverage (DP) below 50% were filtered out to guarantee adequate sequencing depth for accurate genotyping. After filtering, a high-quality SNP dataset was preserved for further analysis.

2.4. Population structure and phylogenetic analyses of Urophysa rockii and U. henryi

To investigate population structure of Urophysa rockii and U. henryi populations, we first mitigated the effects of linkage disequilibrium using PLINK v.1.90 (Purcell et al., 2007) with the parameters –indep-pairwise 50 5 0.2. Population structure was analyzed with ADMIXTURE v.1.23 (Alexander et al., 2009). The ADMIXTURE was run with predefined numbers of clusters (K) ranging from 1 to 10, each run repeated 20 times. The K value with the lowest cross-validation error was chosen as the most likely number of putative genetic groups. Principal component analysis (PCA) was conducted with PLINK v.1.90.

Phylogenetic reconstruction of Urophysa rockii and U. henryi populations were based on our SNP database. Specifically, maximum-likelihood (ML) phylogenetic trees were generated with IQ-TREE v.1.6.8 (Nguyen et al., 2015) using 1000 bootstrap replicates and the GTR + G model, with Semiaquilegia adoxoides and four additional Aquilegia species serving as outgroups. In addition, coalescent-based species trees were estimated using multispecies coalescent methods implemented in MP-EST v.3.0 (Liu et al., 2010) and ASTRAL-Ⅲ v.5.6.3 (Zhang et al., 2018). Detailed analysis procedures and parameters setting are documented in Supplementary Methods Appendix S3.

2.5. Genetic diversity, linkage disequilibrium and genomic inbreeding

We employed the software VCFTOOLS v.0.1.14 (Danecek et al., 2011) and the genomics_general tool (https://github.com/simonhmartin/genomics_general) to compute the nucleotide diversity (π), relative differentiation (F_ST), and absolute differentiation (d_XY) between populations and clades of Urophysa rockii and U. henryi. To detect potential deviations from neutrality in each population of the two Urophysa species, we used Tajima’s D statistics with VCFTOOLS v.0.1.14 (Danecek et al., 2011). All statistics were estimated across 10-kb non-overlapping sliding windows. To discern and compare the linkage disequilibrium (LD) patterns across different Urophysa clades, we calculated the squared correlation coefficient (r²) between pairwise SNPs with a 300-kb window using PopLDdecay v.3.40 (Zhang et al., 2019a). Linkage disequilibrium decay measured the distance at which the Pearson’s correlation efficient (r²) dropped to half of the maximum.

We also determined the extent of inbreeding in Urophysa rockii and U. henryi populations. Runs of homozygosity (ROHs) and genomic regions devoid of heterozygous loci, were identified in all individuals using PLINK’s homozygosity function. To specifically detect ROHs that exceed 10 Kb in length, we meticulously applied a 50-SNP sliding window analysis, rigorously excluding any heterozygous loci. ROHs were then categorized by length: those ranging from 10 to 100 Kb were labeled as short, those from 100 to 200 Kb as medium, and those exceeding 200 Kb as long. Long ROHs are typically indicative of recent inbreeding events, whereas short ROHs are linked to ancient inbreeding (Meyermans et al., 2020). We also quantified the extent of genome-wide inbreeding by calculating the F_ROH, which is the ratio of the cumulative length of all ROHs (greater than 10 Kb) within an individual to the total genome length. Concurrently, we determined the classical individual-based inbreeding coefficient, F_IS, using Plink’s het parameter. A high F_IS value is indicative of reduced heterozygosity among individuals.

2.6. Demographic history analysis

Historical changes in the effective population size (Ne) of Urophysa rockii and U. henryi were inferred using the PSMC program (Li and Durbin, 2011) with default parameters. All individuals of each population were selected to run the PSMC analysis, respectively, with 100 bootstrap replicates per individual. Demographic modeling outcomes were calibrated based on an estimated rate of synonymous substitutions per site per year of approximately 0.5 × 10⁻⁸ for Ranunculales (Guo et al., 2018), with a generation time of three years assumed for Urophysa species.

The demographic history of Urophysa rockii and U. henryi were further examined in fastsimcoal2 v.2.6 (Excoffier et al., 2013). Briefly, various hypotheses were tested regarding species divergence order and time, the presence or absence of asymmetric gene flow between each lineage, and recent fluctuations in effective population sizes. Our models of demographic history were based on our finding that U. rockii and U. henryi could be divided into four phylogenetic clades (two clades per species; see Results, section 3.2). We assigned all populations to four groups based on these clades, facilitating demographic history analysis for each geographical unit. Initially, a multidimensional joint site frequency spectrum (SFS) was constructed from the neutral SNPs using a tailored Python script easySFS.py (https://github.com/isaacovercast/easySFS). The two species diverged at different times without gene flow (models 1 to 3) (see Results part, section 3.3), confirming the species divergence order. Models 4–9 were designed on the basis of models 1–3, incorporating asymmetric gene flow between lineages during different periods. Then, based on the best model detected from models 4–9, historical and recent (75 thousand years ago; kya) gene flows were evaluated (models 10–12). Considering that Urophysa species are unable to flower and bear seeds in the first two years, the generation time was assumed to be three years. The mutation rate was 0.5 × 10⁻⁸ for Ranunculales (Guo et al., 2018). Each model underwent 100 independent runs, with 100,000 coalescent simulations and 50 conditional maximization algorithm cycles to ensure robust global maximum likelihood (ML) estimation. The optimal demographic model was identified using Akaike Information Criterion (AIC). The model with the lowest Akaike value was determined to be the best-fitting model. Finally, parametric bootstrapping (based on the parameter estimates from the best-fitting model) was employed to establish 95% confidence intervals, with 100 independent replications to reinforce the statistical reliability of the results.

2.7. Ecological niche modelling and mantel tests

Current potential distributions of Urophysa rockii and U. henryi were predicted by ecological niche modelling in MAXENT v.3.3.4 with 19 bioclimatic variables at a 2.5 arc-minute resolution (Table S2) from the WorldClim database (v.1.4, accessible at http://www.worldclim.org). We also performed niche overlap and identity tests using ENMTOOLS v.1.4.3 (Warren et al., 2008). The Schoener’s D and the standardized Hellinger distance (I) were used to quantify niche differences between species and clades. Species distribution dynamics were projected for historical periods (LIG, LGM) and future scenarios using ecological niche modelling. Mantel tests were performed to explore the influence of geography and environment on species and population divergence in Urophysa. Regression analyses of F_ST/(1-F_ST) against pairwise geographical and environmental distance matrices were conducted to test for “isolation-by-distance” (IBD) and “isolation-by-environment” (IBE) using the ‘vegan’ R package. For detailed methodology on distribution modelling and Mantel tests, see Supplementary Methods Appendix S3.

2.8. Genomic signatures

We used F_ST analysis to identify genomic regions of karst-adapted Urophysa species that have undergone differential selection. Briefly, we compared F_ST and π ratio of karst-adapted Urophysa species (97 individuals) to that of non-karst group (all outgroups). Genomic regions under selection were indicated by significantly elevated genetic differentiation (top 5% F_ST) as well as notable reductions in genetic diversity within each group (top and tail 5% π ratio). Comparative analysis was first used to assess whether such categorization strategy would influence the detection of candidate genes. A 50-kb window and 10-kb step size were over a total of 306 windows to identify 257 genes potentially under selection. To validate the accuracy of these analyses, we performed positive selection analysis using the branch-site model (Yang and Dos, 2011) and Bayesian empirical Bayes (BEB) methods (Yang et al., 2005) with the PAML software. A detailed description of these analyses is included in Supplementary Methods Appendix S3. After validation, we considered the 257 identified genes as positively selected genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on these candidate genes was performed in R package topGO (Alexa and Rahnenfuhrer, 2022), with p-values adjusted for false discovery rate (FDR).

2.9. Transcriptome analysis of candidate genes positively selected for karst environments

Several genes identified as under positive selection for karst environments are related to calcium transport and signaling. Thus, we determined whether candidate genes positively selected for karst environments were differentially expressed in Urophysa rockii grown in distinct calcium concentrations. U. rockii seedlings were cultivated from seeds as suggested by Hu et al. (2017), and this experimental procedure lasted for 92 days (see Fig. S1 for the detailed processes). Uniform seedlings were cultured in Hoagland hydroponics solution. Plants were then treated daily with Ca²⁺-free (0 mmol/L), medium Ca²⁺ (15 mmol/L), or high Ca²⁺ (30 mmol/L) nutrient solution. Each treatment consisted of three biological replicates, with one seedling per replicate. After 5 days, the leaves and roots of each seedling were collected, flash-frozen in liquid nitrogen, and stored at −80 ℃.

Total RNA was extracted using the Plant RNA Kit (Tiangen, Beijing, China). RNA-seq libraries were constructed and sequenced on the Illumina NovaSeq 6000 platform. Raw reads containing adapters, poly-N sequences, or low quality were filtered out using Trimmomatic v.0.367 (Bolger et al., 2014). The remaining high-quality reads were then aligned to the genome of Urophysa rockii using HISAT2 (Kim, 2019) with default parameters. Gene expression levels were quantified and normalized using the TMM (Trimmed Mean of M-values) method (Robinson and Oshlack, 2010). DESeq2 R package v.1.10.1 (Love et al., 2014) was used to perform differential expression analysis, and genes with threshold of |log₂FC| > 1 and P < 0.05 were regarded as differentially expressed genes (DEGs). DEGs were subjected to GO and KEGG enrichment analysis to determine their functions. Finally, we integrated analyses of DEGs, PSGs, unique genes, and expanded genes. We considered genes present in each of these categories as strong candidates for genes adapted to karst environments in Urophysa.

To determine whether these strong candidate genes had undergone selection across the genome, we used a 10-kb window (as described above) to scan for various molecular signatures indicative of selection. Specifically, we calculated ⅰ) genetic diversity (π) to assess the level of genetic polymorphism within a region; ⅱ) Tajima’s D statistics to detect potential deviations from neutrality (Tajima, 1989); ⅲ) the composite likelihood ratio (CLR) to identify signatures of natural selection in a region (Nielsen et al., 2005); and ⅳ) the integrated haplotype score (iHS) to identify regions under positive selection (Voight et al., 2006). Finally, we used the DCMS method to integrate these four statistics into one DCMS statistic (Ma et al., 2015).

3. Results 3.1. Genome evolution of Urophysa rockii

Urophysa rockii genome size was estimated to be approximately 289.03 Mb and the heterozygosity level 0.26% (Fig. S2 and Table S3). We generated 21.47 Gb (74×) paired high-fidelity (HiFi) long reads and 66.53 Gb (230×) Hi-C data for genome assembly (Fig. S3 and Table S4). The assembled genome was 303.21 Mb, with a contig N50 length of 20.75 Mb (Table 1). The assembled contigs were assembled into 7 pseudo-chromosomes, which is consistent with the previous cytological investigation (Zhang et al., 2019b), ranging in size from 28.23 to 43.56 Mb (Fig. 1e and Table S5). The genome’s completeness was determined to be 97.60% (1575) based on Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis (Table S6), indicating a complete and high-quality genome.

The Urophysa rockii genome contains 32,080 protein-coding genes, with an average coding sequence length of 3690 bp (Table 1). Repeat sequences comprised 48.80% of the genome assembly, totaling 147.98 Mb (Table S7). The genome component and annotation information, see Tables S8 and S9 for more details. The U. rockii genome contains 824 gene families unique within Urophysa (Fig. S4). GO enrichment annotation revealed that the unique gene families of U. rockii are mainly enriched in DNA integration, transposition, transmembrane transport, mitochondrial functions, and ribonuclease activity (Table S10).

A total of 1009 single-copy orthologous genes were extracted from 12 angiosperm genomes, and concatenation and coalescent methods yielded congruent trees with full support values (Figs. S5a and S6). The lineages of monocots, Ceratophyllales and eudicots were all recovered as monophyletic groups, and Ranunculales formed a sister clade to the core eudicots with 100% support, consistent with previous studies. Papaver somniferum was the earliest-diverging species in Ranunculales, and U. rockii was sister to Aquilegia species. Divergence time estimates indicate that the crown-group eudicots arose approximately 141 million years ago (131.3–151.4 Ma, 95% highest posterior density (HPD)). Papaver somniferum diverged from other Ranunculales approximately 124.4 Ma (116.9–132.2 Ma, 95% HPD) and U. rockii diverged from Aquilegia species at 6.0 Ma (3.67–8.7 Ma, 95% HPD). The slightly wide 95% HPD intervals for divergence times may result from limited fossil calibrations and relatively restricted taxon sampling (only 12 species representing angiosperms).

The synonymous substitution rate (Ks) indicates that the Urophysa rockii genome shares the eudicot whole-genome duplication (WGD) event, slightly prior to its divergence from Aquilegia coerulea (Ks peak: ~0.1.0767) but postdating the divergence from V. vinifera (Ks peak: ~1.3391) (Fig. S7a). No additional duplication was detected. Intergenomic analyses revealed a clear 1: 2 syntenic depth ratio in large collinear blocks between U. rockii and A. trichopoda, whereas a 3:2 synteny ratio was observed between U. rockii and V. vinifera (Fig. S7b). CAFE analysis of the U. rockii genome revealed that 754 gene families expanded and 896 contracted (Fig. S5b and Table S11). Gene Ontology (GO) enrichment analysis indicated that the expanded gene families in U. rockii are primarily enriched for ion channel functions and binding (e.g., heme binding; magnesium ion binding; iron ion binding; and calcium ion binding), metabolism-related enzymes activities (e.g., monooxygenase activity; catalytic activity; RNA-directed DNA polymerase activity), membrane-bounded organelle organization; cellular response to stimulus; nutrient reservoir activity; and cell surface receptor signaling pathways (Fig. S5b and Table S11). Functional analysis also revealed that the contracted gene families were mainly enriched for O-methyltransferase activity, ATP metabolic process, methylation, and proton-associated pathways (Table S12).

3.2. Population structure and genetic diversity

Whole-genome resequencing of 98 individuals from two Urophysa species and the outgroup Semiaquilegia adoxoides generated 1134.269 Gb of clean data (Table S1). On average, approximately 95% of the clean reads were successfully aligned to the reference genome, achieving an average sequencing depth of 40.1× and a coverage rate of 94.73%. Using this dataset, we detected 12,925,384 high-quality SNPs for subsequent analyses. To investigate the population structure, we first performed ADMIXTURE analysis, K = 2 was the optimal K-value with the lowest cross-validation error (0.365) (Figs. S8 and S9), which separated all individuals of U. rockii and U. henryi into two distinct genetic groups (Fig. 2). When K = 3, U. henryi populations subdivided into two geographically distinct clades, designated as the West-clade (HW) and East-clade (HE). At K = 4, U. rockii populations clustered into two additional geographically aligned clades, the West-clade (RW) and East-clade (RE), with no detectable admixture. PCA results largely corroborated the ADMIXTURE findings (Fig. 2): the first principal component distinguished U. rockii from U. henryi, the second principal component separating the two clades of U. henryi (HW and HE). These patterns were further confirmed by phylogenetic results inferred through concatenation-based (PP/IQ-TREE) (Fig. 2d) and coalescent-based methods (BS/MP-EST, LPP/ASTRAL-Ⅲ) (Fig. S10). Both sets of analyses collectively resulted in the division of all individuals into two major lineages, which were subdivided into four clades. This division is consistent with the inter- and intra-specific relationships.

Genetic diversity is higher in U. henryi than in U. rockii, at both the species level (π = 0.0055 ± 0.00157 for U. henryi, and π = 0.0027 ± 0.0008 for U. rockii) and the clade and population levels. Genetic diversity was lowest in two western populations of U. rockii (CXS: π = 0.0000876l; PZ: π = 0.0000963) (Fig. 2e and Fig. S11; Table S13). Tajima’s D distributions were significantly skewed toward middle-frequency variants in both species (D_{U. rockii} = 0.9809 ± 0.5135, D_{U. henryi} = 0.3151 ± 0.2691) and within clades (Fig. 2f and Fig. S12; Table S13). Tajima’s D values were generally higher in U. henryi populations than in U. rockii populations, while the two populations that formed the western U. rockii clade (i.e., CXS and PZ) displayed a pronounced shift toward low-frequency variants. Notably, the number of medium and long ROHs were higher in U. rockii individuals, especially in the two western populations (Fig. 2g and h; for more details, see Table S14). Furthermore, both F_ROH and F_IS coefficients were significantly higher in U. rockii than in U. henryi, a pattern evident at the population level, particularly in population PZ (F_ROH = 0.4878 ± 0.0006, F_IS = 0.4885 ± 0.1601) and CXS (F_ROH = 0.4748 ± 0.0005, F_IS = 0.4870 ± 0.1592) (Table S13). These results indicate significantly higher inbreeding in U. rockii, with populations PZ and CXS (RW clade) showing the highest inbreeding coefficients.

3.3. Demographic history of Urophysa species

Significant genetic differentiation was observed between Urophysa rockii and U henryi (F_ST = 0.247, d_xy = 0.256). Similarly, significant genetic differentiation was detected among different clades and among populations of these species (Figs. 3a, S13 and S14). LD decay distance was shorter in U. henryi (600 bp at r² = 0.0501) than in U. rockii (2200 bp at r² = 0.1605). In addition, LD decay distances were shorter in the two clades of U. henryi (900 bp at r² = 0.0742 in HW, 900 bp at r² = 0.1041 in HE) than in the two clades of U. rockii (1200 bp at r² = 0.3959 in RW, 1600 bp at r² = 0.1835 in RE) (Fig. 3b).

Using the PSMC approach, we reconstructed historical effective population size (N_e) trajectories for four clades of two Urophysa species. Each showed different demography (Fig. 3c). Both the western (RW) and eastern clades (RE) of U. rockii had a stable N_e from 5.0 to 8.0 million years ago (Ma), followed by a rapid expansion peaking at 4.0–5.0 Ma. The effective populations of these clades then diverged. The western clade of U. rockii (RW) decreased continually until reaching its present level. In contrast, the eastern clade of U. rockii (RE) experienced a significant bottleneck between 0.4 and 0.7 Ma, followed by a dramatic recovery at ~0.25 Ma and subsequent stabilization. The effective populations of both the eastern and western clades of U. henryi (HE and HW) expanded between 4.0 and 7.0 Ma, reaching their maximum at 4.0 Ma. The effective populations of the two U. henryi clades diverged markedly after 0.6 Ma. The western clade of U. henryi (HW) expanded (0.3–0.6 Ma), then decreased, with a brief bottleneck (~0.6 Ma); in contrast, the eastern clade (HE) declined continuously until being interrupted by a minor expansion around 0.1 Ma, and a bottleneck (0.3–0.4 Ma) before stabilizing.

We determine the best-fit model (model 10) based on Akaike’s weight values (Tables S15 and S16), which revealed no gene flow since approximately 75 thousand years ago (kya) (Fig. 3d), and the two Urophysa species diverged about 5.95 Ma (95% HPD: 5.18–6.29 Ma). The eastern and western clades of U. henryi diverged approximately 5.35 Ma (95% HPD: 4.75–5.91 Ma). The U. rockii clades diverged about 5.31 Ma (95% HPD: 4.34–6.14 Ma), with both experiencing significant contraction, particularly for the western clade of U. rockii (changed from 45,813 to 5189), although a slight population expansion occurred in the eastern clade (from 8613 to 11,533).

Our best model indicates that ancestral migration continued between Urophysa rockii and U. henryi after their divergence (M_{4HE_N23} = 6.91e-5 and M_{N23_4HE} = 1.34e-7), although the models project that no gene flow has occurred for approximately 75 kya. Migration between the eastern and western clades of U. henryi persisted until recently, with greater migration from the western clade to the eastern (M_{N4_NHE} = 1.41e-2, M_{N1_HE} = 5.07e-5; (M_{NHE_N4} = 4.64e-4, M_{HE_N1} = 3.53e-7). During this period, the population size of the eastern clade was stable (N_HE = 56,422), whereas the population of the western clade expanded before experiencing a significant contraction (from N₄ = 86,269 to N₁ = 94,472 and current HW is 76,487). Asymmetric gene flow was also observed between the clades of U. rockii, with migration predominantly from the eastern clade to the western (M_{N2_N3} = 1.22e-4, M_{N3_N2} = 7.03e-8), although no recent gene flow is supported.

Mantel test revealed significant positive genetic–geographical correlations across Urophysa species (R = 0.022) and within U. henryi (R = 0.522) (Fig. S16a and c), alongside significant genetic–environmental correlations for the genus (R = 0.469, P < 0.05) and U. henryi (R = 0.518, P < 0.05) (Fig. S16b and d). Ecological niche modeling indicated that key environmental predictors differed between species. Bio2 (27.2%) contributed most to U. rockii, whereas Bio6 (36.5%) contributed most to U. henryi (Fig. S17). Significant niche divergence occurred between U. rockii and both U. henryi clades (HE and HW) (D/I metrics below null expectations; Fig. S18a, c, d), and between the eastern and western clade of U. henryi (Fig. S18b). Distribution projections across four climate scenarios indicated historical range maxima during the Last Interglacial (LIG), followed by progressive contractions through the Last Glacial Maximum (LGM) to the present. Future scenarios project accelerated range loss in Urophysa species (Fig. S19), particularly for U. rockii, which faces potential extirpation across its current range (Fig. S19a–d).

3.4. Candidate karst adaptation genes in Urophysa

To explore karst-adaptation candidate genes in Urophysa species, we performed a genomic scan for outlier regions with elevated genetic divergence (top 5%) and reduced genetic diversity (top and tail 5% π ratio) between karst-adapted Urophysa and non-karst outgroups. A total 306 windows with 257 candidate genes were identified (Fig. 4a). Comparative analysis confirmed minimal influence of karst-group sample size variation on candidate genes detection (Fig. S20). We further validated these 257 candidate genes via branch-site model analysis, which revealed that 252 out of the 257 genes showed significant evidence of positive selection (Table S17), confirming methodological robustness. GO enrichment showed that these 257 candidate genes are associated with membrane, cellular response to stimulus, manganese ion transmembrane transporter activity, vesicle-mediated membrane, calmodulin binding, and calcium ion binding (Fig. 4b). KEGG enrichment indicated that these genes are involved in plant–pathogen interaction, plant hormone signal transduction, and MAPK signaling pathways (Fig. 4c).

Because many of the candidate karst-adaptation genes we identified are related to calcium transport and signaling, we determined whether these genes are differentially expressed in the roots and leaves of Urophysa rockii following calcium treatment (Fig. 5a). Following calcium treatment, 205 genes were up-regulated and 208 were down-regulated in leaves (Fig. S21); in roots, 345 genes were up-regulated and 129 were down-regulated (Fig. S22). GO enrichment revealed up-regulated genes in both tissues significantly enriched in pathways related to stress or defense response (GO: 0006950, GO: 0005516, GO: 0006952), calcium ion binding (GO: 0005509), ion channel activity (GO: 0005216), and ion transmembrane transport (GO: 0006811, GO: 0015075, GO: 0098655, GO: 0098660) (Tables S18 and S19). While most down-regulated DEGs were primarily associated with metabolic processes in both tissues (Tables S20 and S21).

We integrated differentially expressed genes (DEGs) from both tissues, yielding 537 up-regulated and 328 down-regulated DEGs (Fig. S23). Among them, 15, 23 and 20 DEGs overlapped with selected, unique and expanded genes, respectively (Fig. S24 and Table S22). Functional analysis associated these DEGs with high calcium adaptation, membrane ion exchange or transport, and cell wall-associated receptors (Fig. 5b and c; Table S22). Specifically, genes ACA11, TPC1, KIC, GLR27, CAX3, SLAH1, CNGC4, CSCL4, CMTA3 and CB60C mediate Ca²⁺ sensing, transport, and metabolism (Fig. 6) (Cao et al., 2023; Zhou et al., 2023). PTR/POT, TMN11 and ABC are associated with transmembrane transport, CCX1 modulates Ca²⁺ signaling (Li et al., 2016), and MRS2F facilitates magnesium ion influx (Yu et al., 2005).

Moreover, expression analysis revealed calcium-dependent regulation of key genes: CAX3, a gene encodes a cation/proton antiporter that helps maintain cellular ion homeostasis by sequestering calcium ions from the cytoplasm into the vacuole (Jamra et al., 2024; Wang et al., 2024), was up-regulated in both tissues under medium/high Ca²⁺ versus calcium-free conditions (Fig. 5d). TPC1 (vacuolar Ca²⁺ release channel) and PTR/POT (ion/nutrient transporter) showed tissue-wide up-regulation at medium Ca²⁺ but down-regulation at high Ca²⁺ (Fig. 5d). Functional annotations align with their roles in ion homeostasis: TPC1 mediates the transport of Ca²⁺ from the vacuole to the cytoplasm, thereby contributing to the homeostatic regulation of intracellular calcium ions (Pottosin and Dobrovinskaya, 2022), while PTR/POT transports ion and nutrients (Newstead, 2015). Genomic analyses (top 1% F_ST, DCMS scores; Fig. 5e, f and Fig. S23a, b) identified strong selection signals for TPC1, PTR/POT, and CAX3, supported by elevated π, CLR, iHS, and reduced Tajima’s D. Signals were also observed for transmembrane transporters (ABC, TMN11), the cell wall receptor WAK3 (He et al., 1999; Li et al., 2023), and genes involved in calcium binding, transport, and ions homeostasis across at least one genomic parameter (Fig. 6 and Fig. S25).

4. Discussion 4.1. The great differentiation and demographic history of the two Urophysa species

Significant geographic barriers and ecological divergence often impede gene flow among populations, driving differentiation and allopatric speciation (Yona et al., 2015; Qian et al., 2021; Anderson and Weir, 2022). This speciation pattern is prevalent in biodiversity hotspots, such as the Hengduan Mountain (HDM) and Tibetan Plateau (Wu et al., 2022; Yu et al., 2023), yet remains understudied in specialized habitats such as karst ecosystems. In this study, we used genomic data to characterize the evolution and adaptation of U. rockii to karst regions in Yunnan-Guizhou Plateau. We found significant genetic divergence (high F_ST and d_XY values) between two Urophysa species and their clades. Demographic analyses indicated that the two species diverged approximately 5.95 Ma, with major lineages splits at ~5.35–5.31 Ma, aligning with the Late Miocene intense uplift of the HDM. Thus, our findings suggest that Urophysa divergence was initiated by Late Miocene orogeny. The divergent evolutionary trajectories between Urophysa species were also evidenced by disparities in effective population sizes (Ne), gene flow, genetic diversity, and niche occupancy. Both Urophysa species underwent pronounced population contractions after this time (approximately 4.0 Ma), although the western clade of U. henryi (HW) underwent a population expansion by 0.30 Ma, which slightly recovered its population size (N₁ = 94,472). These reduced population sizes indicate that both species underwent bottlenecks, potentially caused by climatic cooling during the Pleistocene glaciation. Previous studies have shown how Pleistocene glacial-interglacial oscillations amplified differentiation between other species (Hewitt, 2004; Sun et al., 2019; Jardim de Queiroz et al., 2022). In Urophysa, population divergence and allopatric speciation may thus have been accelerated by Pleistocene glacial cycling.

This scenario is further supported by demographic results, which revealed asymmetric inter-specific gene flow during initial Urophysa divergence, ceasing by the Late Miocene, coincident with intra-specific lineage differentiation (Fig. 3d). PSMC results indicate significantly reduced population sizes in both species compared to their ancestral state, suggesting a potential bottleneck possibly caused by climatic cooling during the Pleistocene glaciation. Ecological niche modeling analyses predicted highly suitable areas for Urophysa species across the Sichuan Basin and Yunnan-Guizhou Plateau during the LIG, whereas these areas fragmented during the LGM and further contracted into restricted regions, and more pronounced in U. rockii (Fig. S19). This implies a formerly widespread Urophysa ancestor across the Yunnan-Guizhou Plateau and adjacent regions, with subsequent extensive orogenic and climatic oscillations driving species divergence and shaping current distributions.

Although our analysis of isolation by distance (IBD) were non-significant, populations are known to be partitioned by major geographical barriers, such as high mountains (Daba, Xuefeng, Wushan, and Dalou Mountains) and large rivers (Yangtze, Jialing, Minjiang, and Chishui Rivers). Pronounced intra-specific isolation occurred in U. henryi (HW vs. HE, R = 0.522, P < 0.05), and geographic barriers such as Mingjiang River and Longmen Mountains significantly isolated the two lineages of U. rockii (RW and RE). Moreover, significant ecological isolation (IBE/IBS: R > 0.4, P < 0.05) and niche divergence (D/I < null expectations) further enhanced genetic differentiation between the two species and among clades within U. henryi. These patterns align with niche-driven speciation mechanisms observed in Brassica rapa, Primulina and Poaceae (Johnson et al., 2010; Ke et al., 2022; Burrill et al., 2023; Dorey and Schiestl, 2024), further supporting that genetic differentiation between U. henryi and U. rockii was likely increased by ecological isolation and niche divergence.

4.2. Current status of the endangered Urophysa rockii and implications for conservation

Biodiversity serves as the essential bedrock for human survival and development, yet escalating climate change and anthropogenic impacts have pushed numerous plant species toward extinction (Bongaarts, 2019; Ceballos et al., 2020; Tao et al., 2024). Endemic and endangered Urophysa species are distinctive winter-blooming plants with significance horticulture and conservation value. Genetic diversity was low in current populations of U. rockii (π = 0.0027), particularly in its two populations located in the Longmen Mountains region (CXS = 8.7 × 10⁻⁵ and PZ = 9.6 × 10⁻⁵). These levels of genetic diversity are lower than those of U. henryi and other endangered species (Jing et al., 2023), such as Oryza rufipogon (π = 0.003), Davidia involucrata (π = 0.00585), and Acer yangbiense (π = 0.00313) (Huang et al., 2012; Chen et al., 2020; Ma et al., 2022; Ren et al., 2024). Our analysis of ROHs, F_ROH, and F_IS values also indicate that U. rockii populations have high levels of inbreeding depression, especially in the two populations in the Longmen Mountains region (e.g., CXS and PZ).

Previous studies have shown that temporary yet severe population reductions (bottleneck) cause significant genetic diversity loss through intensified genetic drift and reduced efficacy of natural selection (Alcala and Vuilleumier, 2014; Jangjoo et al., 2016). This pattern applies to Urophysa rockii, population sizes (Ne) of the two U. rockii clades decreased continuously since 4.0 Ma, with a prolonged bottleneck driving U. rockii toward extinction since approximately 0.1 Ma (Ne of RW = 5189) (Fig. 3d). Slower LD decay rates (r²) further confirm U. rockii experienced strong selection and genetic drift, causing population contraction and drastic genetic diversity loss (Butler et al., 2022). This conclusion was further supported by the Tajima’s D values and the lower H_E (Table S13). In addition, confined to karst barriers (Xie et al., 2017), U. rockii shows high inter- and intra-specific genetic differentiation (F_ST and d_XY) (Fig. 3a). Demographic simulations confirm restricted historical gene flow between U. rockii’s primary clades (Fig. 3d). Such significant genetic differentiation and restricted gene flow among populations may result in low heterozygosity, a crucial indicator of genetic diversity (Ellegren and Galtier, 2016). This was also supported by our findings that lower H_O and H_E values were detected in U. rockii than its relative U. henryi. Furthermore, U. rockii exhibits pronounced habitat specialization, restricted to steep karstic cliffs fragmented by mountainous terrain and valleys (Xie et al., 2017). This ecological confinement severely restricts seed dispersal and population connectivity, driving elevated inbreeding rates and diminished reproductive capacity (Zhang et al., 2013), a pattern consistent with threatened species exhibiting small, isolated populations prone to inbreeding depression and extinction risks (Neaves et al., 2015). Most studies demonstrate that progressive inbreeding erodes heterozygosity and genomic recombination (Charlesworth, 2003; Willoughby et al., 2015; Hedrick and Garcia-Dorado, 2016), exacerbating genetic isolation and diversity loss. Compounding these pressures, anthropogenic disturbances, including rock climbing, wild plant harvesting, and cliff-side infrastructure development, further intensify population vulnerability. Collectively, persistent demographic decline and genetic erosion, restricted gene flow, accumulation of inbreeding, and anthropogenic impacts have depleted the genetic diversity of U. rockii and dramatically elevated its extinction risk.

The ecological environment in karst regions is fragile and exhibits a low resistance to interference (Li et al., 2017). Our simulation predicted a further contraction of Urophysa rockii’s distribution range in the future, highlighting the need for appropriate conservation and restoration measures. We recommend that conservationists collaborate with forestry authorities to establish local conservation areas across its distribution, delineating protected zones and core habitats with strict bans on excavation, rock climbing, and destructive construction. Second, priority protection should be given to the genetically depauperate Longmen Mountain populations (e.g., PZ and CXS) by implementing cross-pollination with other populations to boost their genetic diversity. Third, artificial propagation efforts should be expanded, and reintroduction experiments conducted, alongside the establishment of long-term monitoring programs to ensure population recovery and persistence.

4.3. The karst limestone habitat adaptation of the Urophysa species

Elucidating karst-specific adaptive mechanisms of Urophysa species provides critical insights for optimizing conservation strategies and artificial propagation. Given the cell wall’s role as the primary barrier against external stressors (Miedes et al., 2014; Loix et al., 2017), our differential gene expression analysis revealed significant enrichment of cell wall organization, biogenesis, and modification genes in roots and leaves. Notably, a candidate DEG encoding wall-associated receptor kinase (WAK3) exhibited strong selection signals (Figs. 5e and S25; Table S22). Previous studies demonstrate that the WAK3 gene contains an EGF-2 like domain and calcium-binding EGF, mediating plant cell wall-plasma membrane interaction, playing a crucial role in ion recognition, cell expansion control, and morphogenesis (He et al., 1999; Li et al., 2023). These findings collectively indicate that maintaining cell wall integrity constitutes a fundamental adaptation mechanism for Urophysa species inhabiting high-calcium karst limestone environments.

In addition, our genomic scan analysis (top 5% F_ST and π ratio) detected 257 genes that showed signs of selection in Urophysa versus non-karst taxa, functionally enriched in karst limestone adaptation, such as cellular response to stimulus, ion transmembrane transporter activity, calmodulin binding and calcium ion binding (Fig. 4b). Our transcriptome analysis of U. rockii seedlings identified genes under positive selection that were differentially expressed in response to calcium, including ion homeostasis genes TPC1, CAX3 and PTR/POT (Figs. 5d–f and 6 and Fig. S25; Table S22). Previous studies have indicated that TPC1, a calcium influx channel found in all land plants, maintains a conserved function across both monocots and dicots (Dadacz-Narloch et al., 2013). It contains an EF-hand domain for cytoplasmic Ca²⁺ sensing and is implicated in high-calcium stress adaptation (Cao et al., 2023). Positive selection on TPC1 in the karst-adapted herbaceous Primulina species and in the woody Platycarya longipes corroborates its adaptive significance (Tao et al., 2016; Cao et al., 2023). In Urophysa, TPC1 expression exhibited a biphasic pattern: up-regulated at moderate Ca²⁺ levels to facilitate uptake, but down-regulated during calcium overload (Fig. 5d), suggesting concentration-dependent inhibition protects against cytotoxic influx via Ca²⁺-mediated negative feedback. Similarly, the peptide transporter PTR/POT was up-regulated at intermediate calcium to enhance absorption, but suppressed under high calcium to prevent toxicity (Jørgensen et al., 2015; Parker et al., 2017). Conversely, the calcium transporter CAX3 (a Ca²⁺/H⁺ antiporter) was up-regulated across both calcium levels, sequestering excess calcium into vacuoles to maintain cytosolic homeostasis (Zhao et al., 2009; Wdowiak et al., 2024). Collectively, the coordinated regulation of TPC1 (uptake modulator), PTR/POT (absorption buffer), and CAX3 (vacuolar sequestration) constitutes a critical adaptive echanism enabling Urophysa species to thrive in calcium-saturated karst environments.

Excess intercellular ions pose significant threats to plant survival (White and Broadley, 2003; Case et al., 2007). To mitigate such stress, Urophysa species employ dual strategies. ⅰ) Ion flux regulation: they can mitigate ion accumulation by regulating absorption or excretion via positively selected or expanded genes (ACA11, GLR27, CSCL1, SLAH1, and CNGC4) involved in calcium transport (Fig. 6) (Li et al., 2007; Gilliham et al., 2011; Jin et al., 2018). ⅱ) Subcellular compartmentalization. Urophysa species can sequester ions in organelles via compartmentalization, such as in vacuoles and mitochondria (Petersen, 2002; Orrenius et al., 2003). Specifically, candidate gene ACA11 facilitates the storage of excess calcium ions in vacuoles, while gene MRS2F transfers magnesium ions to mitochondria (Yu et al., 2005). This mode of ion transport and storage likely helps Urophysa species effectively alleviate the stress of high calcium conditions. However, plant’s capacity for ion absorption, excretion, and storage is inherently limited (Jin et al., 2018). Here, we detected expression of candidate genes such as KIC, CB60C, and CMTA3 in Urophysa species (Table S22). They encode proteins that detect intracellular calcium ion concentrations, transmit calcium signals, and selectively bind to calcium ions within the cell, thereby regulating ion homeostasis and maintaining normal physiological functions. Collectively, Urophysa species have evolved an integrated defense network, combining controlled ion flux, targeted compartmentalization, and real-time sensing, equipping them to survive in specialized karst environments. Additional karst-adaptive genes (CNGC, GLR27, MRS2F) exhibiting differential expression patterns were also detected in Urophysa and other karst taxa (Jin et al., 2018; Feng et al., 2020; Cao et al., 2023; Zhou et al., 2023). Further study is needed to determine whether these genes play a role in mitigating ion-related stress in karst habitats.

5. Conclusion

We generated a high-quality genome assembly for the endangered cliff-dwelling plant Urophysa rockii endemic to karst limestone habitats. Comparative genomic and population re-sequencing analyses revealed that Late Miocene orogenies and climatic shifts in the HDM likely triggered the divergence of the two Urophysa species. Subsequent Pleistocene climatic fluctuations further enhanced population divergence and accelerated allopatric speciation, with niche divergence playing an important ecological role. Currently, U. rockii exhibits extremely low genetic diversity, particularly in CXS and PZ populations. Persistent demographic bottlenecks, restricted gene flow, habitat specialization, and severe inbreeding collectively contributed to its endangered status. Based on these findings, appropriate conservation and restoration measures are proposed. Moreover, by integrating seedling cultivation and physiological analyses, we identified many candidate DEGs that have experienced selection and expansion. We demonstrate that maintaining the integrity of the cell wall serves as the primary defense mechanism for Urophysa species to cope with the harsh conditions of karst limestone environments. Coordinated regulation of TPC1, PTR/POT, and CAX3 enables calcium homeostasis. Additionally, Urophysa species exhibited a significant number of selective or expanded genes associated with ion uptake/efflux, transport and cellular storage, which allow them to survive in the harsh karst environments. This study elucidates the evolutionary divergence and calciotrophic adaptations underpinning Urophysa’s persistence in karst ecosystems, offering valuable insights into endangerment mechanisms of specialized flora.

Acknowledgments

We thank Dr. Minjie Li (College of Ecology, Lanzhou university) for her valuable suggestions in preparing this manuscript. This work was supported by the National Natural Science Foundation of China (Grant Nos. 32100180, 32470216, 32170209), the China Postdoctoral Science Foundation (2020M683303), the Sichuan Science and Technology Program (Grant Number 2025ZNSFSC1129), the Key project at central government level: the ability establishment of sustainable use for valuable Chinese medicine resources (2060302).

CRediT authorship contribution statement

Deng-Feng Xie: Methodology, Formal analysis, Conceptualization, Visualization, Writing review & editing, Writing original draft, Funding acquisition. Yi-Yang Zhang: Methodology, Formal analysis, Validation. Jing Cai: Methodology, Investigation, Formal analysis. Rui-Yu Cheng: Methodology, Formal analysis, Validation, Data curation. Yuan Wang: Methodology, Investigation, Formal analysis, Data curation. Jin-Bo Tan: Methodology, Investigation, Conceptualization. Yan Yu: Methodology, Formal analysis. Xing-Jin He: Conceptualization, Visualization, Writing review & editing, Writing original draft, Funding acquisition. Song-Dong Zhou: Investigation, Funding acquisition, Writing review & editing, Writing original draft.

Data availability

The assembled whole-genome has been deposited in Genome Sequence Archive (GSA) of National Genomics Data Center (NGDC) under project PRJCA035353, and all re-sequencing data of populations are deposited in Genome Sequence Archive (GSA) of National Genomics Data Center (NGDC) under project PRJCA035354; All transcriptome data have been deposited in Genome Sequence Archive (GSA) of National Genomics Data Center (NGDC) under project PRJCA035355.

Declaration of competing interest

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

Appendix A. Supplementary data

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