1. Introduction

Genetic diversity is critical to the long-term survival and adaptability of wild plant and animal species (Frankham, 1996; Barrett and Schluter, 2008). This importance has been formally recognized at the policy level in the Kunming-Montreal Global Biodiversity Framework (Ma, 2023). Conserving and sustainably managing genetic diversity requires a thorough understanding of the distribution and amount of genetic diversity within populations and species, as well as identifying the key factors that shape this diversity. The level and distribution of genetic diversity within a plant species can be influenced by multiple factors, such as geographic range, population isolation, effective population size, and life history traits (Hamrick et al., 1992; De Kort et al., 2021; Salgotra and Chauhan, 2023). Endangered species, often confined to small geographic ranges and small effective population sizes, typically exhibit low levels of genetic variability and increased rates of inbreeding (e.g., Arafeh et al., 2002; Cole, 2003; Vellend, 2004; Honnay and Jacquemyn, 2007; Aguilar et al., 2008; Edwards et al., 2013; Turchetto et al., 2016; Ma et al., 2021; Ren et al., 2024; Tao et al., 2024). In these cases, strong genetic drift, high levels of inbreeding, and loss of genetic variation may facilitate the expression of harmful recessive alleles, reduce the adaptive potential of the endangered species (Kohn et al., 2006), and thereby increase risk of extinction (Frankham, 2005; Kardos et al., 2018; Kleinman-Ruiz et al., 2022). However, endangered species do not always exhibit low genetic diversity; substantial genetic diversity has been documented in several rare and endangered plants (Barrett and Kohn, 1991; Ellstrand and Elam, 1993; Lutz et al., 2000). It is essential to recognize that the genetic patterns observed in a species can frequently be driven by evolutionary constraints and historical factors, rather than solely current processes (Hamrick et al., 1992; Excoffier and Ray, 2008). Consequently, comparing related species is recommended to account for evolutionary history and identify species-specific differences (Soltis and Gitzendanner, 1999; Gitzendanner and Soltis, 2000; Ouborg et al., 2006; Keller et al., 2013; Lanes et al., 2018).

Historically, genetic diversity was characterized using allozyme markers and later microsatellite loci, with metrics typically reported as observed heterozygosity, mean allelic diversity, and the proportion of polymorphic loci. These methods provided a foundational understanding of genetic variation but were limited in their resolution (Allendorf, 2017). With the advent of next-generation sequencing technologies, however, we now have the capacity to generate de novo assembled and annotated reference genomes for an ever-increasing number of non-model organisms (Primmer, 2009; Allendorf et al., 2010). This significant development has revolutionized the field by enabling accurate estimation of genome-wide diversity and inbreeding levels, as well as the assessment of genetic load, even with limited sample sizes (Ouborg et al., 2010; Dussex et al., 2021; Hohenlohe et al., 2021; Khan et al., 2021; Yi et al., 2024). Moreover, whole-genome resequencing data can facilitate the detection of population bottlenecks and past demographic events, offering a historical context that can guide future conservation efforts (Feng et al., 2019). This is of particular importance when dealing with endangered species, as they are often elusive and have dwindling populations. Specifically, genomic comparison of endangered species to their closely related widespread congeners can improve our understanding of the genetic trajectory of endangered populations and inform conservation management strategies.

The genus Oreocharis (Gesneriaceae) is a monophyletic lineage of about 150 species predominantly distributed in China, with several species also occurring in Japan and parts of South and Southeast Asia (Li and Wang, 2005; Möller et al., 2016, 2018). The redefined Oreocharis was recently enlarged to encompass species from 12 other monotypic or small genera of the subfamily Didymocarpoideae (Möller et al., 2011; Yang et al., 2020) and is characterized by a remarkable life history and functional trait diversity (Kong et al., 2022). Originating in the Miocene, the genus experienced an early rapid radiation and diversified during the Pliocene and Pleistocene giving rise to three current centers of diversity in Southwest and South China (Kong et al., 2022). Notably, most species currently have restricted ranges or are narrow endemic to small regions (e.g., O. esquirolii and O. pilosopetiolata), whereas few species have large geographic distribution areas (e.g., O. maximowiczii and O. auricula). In this study, we examine patterns of genetic diversity, effective population size, inbreeding, mutation load and demographic history of two Oreocharis species with contrasting distribution ranges. Oreocharis esquirolii (formerly Thamnocharis esquirolii), the primary subject of this study, is classified as a critically endangered species in the National Key Protected Wild Plants of China (2021; https://www.gov.cn/zhengce/zhengceku/2021-09/09/content_5636409.htm), with only three known populations found in Anlong and Zhenfeng Counties in Guizhou, Southwest China (Fig. 1). These populations are generally small and geographically disjunct, with two populations (ZF01 and ZF02) located in a natural reserve and one (AL) outside protected areas (Fig. 1). A recent study on photosynthetic physiological traits and anatomical structures found that O. esquirolii has the lowest overall photosynthetic performance and the weakest adaptability to environmental changes compared to its non-endangered relatives (Ou et al., 2023). O. maximowiczii is a widespread and common species distributed across southern and eastern China, ranging from Hunan, Guangdong, Jiangxi, Fujian, and Zhejiang (Wang et al., 1998). Both species typically grow on wet rocky substrates near forested slopes at elevations ranging from 200 to 800 m. Although O. esquirolii and O. maximowiczii are not sister species, O. maximowiczii represents the most closely related widespread species, allowing us to directly assess how differences in geographic distribution influence genetic diversity patterns.

Here, we undertake a comprehensive conservation genomics assessment of the endangered species Oreocharis esquirolii and draw comparisons with the more widespread O. maximowiczii. Our approach involves assembly and annotation of high-quality reference genomes for both species. We complement these genomic resources by generating whole-genome resequencing data from 28 individuals representing all known populations of O. esquirolii, alongside data from 79 individuals across 11 representative populations of O. maximowiczii. With the genomic data generated herein, we then investigate genomic features, including genetic diversity, inbreeding, and mutation load, to assess the genetic health and adaptive potential of each species. Additionally, we reconstruct the ancient and recent demographic histories of different populations, examining how these historical dynamics influenced the levels and distribution of genetic diversity and the accumulation of mutation load. Given the critically endangered status of O. esquirolii, we hypothesize that it will exhibit lower genetic diversity and higher mutation load compared to its widespread relative, O. maximowiczii. This study aims to provide critical insights into the evolutionary trajectories and conservation needs of these species, ultimately informing strategies for their preservation and management.

2. Materials and methods 2.1. Sample collection and genome sequencing

To investigate the genomic and evolutionary differences underlying the contrasting distribution patterns of Oreocharis esquirolii and O. maximowiczii, we collected fresh young leaves from one individual of each species. O. esquirolii was sampled from Anlong County (Guizhou, China), and O. maximowiczii from Nanfeng County (Jiangxi, China). Immediately after collection, leaf tissues were flash-frozen in liquid nitrogen and transported on dry ice to Novogene Corporation (Tianjin, China) for DNA extraction and sequencing. For each species, five tissues (root, stem, leaf, flower, and fruit) from the same individual were used for RNA sequencing.

We employed a combination of Pacific Biosciences (PacBio) single-molecule real-time (SMRT) long-read sequencing, Illumina short-read sequencing, and high-throughput chromosome conformation capture sequencing (Hi-C) to generate high-quality genome assemblies for both species. For PacBio sequencing, high-molecular-weight genomic DNA was used to construct 20-kb insert size libraries. Sequencing was performed on the PacBio Sequel Ⅱ platform using four SMRT cells per species to produce long-read data suitable for de novo assembly. For Illumina sequencing, short-insert libraries (~350 bp) were prepared using the NEBNext Ultra Ⅱ DNA Library Prep Kit for Illumina (New England Biolabs), following the manufacturer's protocols. Libraries were sequenced on the Illumina NovaSeq 6000 platform with paired-end reads of 150 bp to produce high-coverage short-read data for error correction and polishing of the assembly. Hi-C libraries were prepared using the Proximo Hi-C Plant Kit (Phase Genomics) to facilitate chromosome-level assembly.

2.2. Genome size estimation, assembly and annotation

To estimate genome size, heterozygosity, and repeat content, we performed k-mer frequency analysis of Illumina short reads based on the 21 k-mer depth distribution with GCE (Liu et al., 2013). The PacBio long reads and Hi-C data were integrated using the Hi-C mode of hifiasm v.0.14 (Cheng et al., 2021) to generate primary contig assemblies for each species. This approach takes advantage of the high accuracy of Hi-C data to scaffold contigs into chromosome-level assemblies. To remove redundant haplotigs and correct assembly errors due to heterozygosity, we applied Purge Haplotigs v.1.1.2 (Roach et al., 2018). The draft genomes were further refined using ALLHiC v.0.9.8 (X. Zhang et al., 2019), which utilizes Hi-C linkage information to improve assembly contiguity and accuracy. Manual curation was performed using Juicebox Assembly Tools (Robinson et al., 2018) to correct misassemblies and adjust scaffolding errors.

The quality of the genome assemblies was assessed using BUSCO v.5.6.1 (Manni et al., 2021) with the embryophyta_odb10 database to evaluate the completeness of conserved single-copy orthologous genes. To calculate coverage and average sequencing depth, Illumina short reads were mapped back to the assembled genomes using BWA-MEM v.0.7.18 (Li and Durbin, 2009). Mapping statistics were generated using SAMtools v.1.9 (Li et al., 2009). Additionally, synteny between the O. esquirolii and O. maximowiczii genomes was analyzed using the MCScanX toolkit (Wang et al., 2012) within the jcvi package (Tang et al., 2024) to evaluate the structural accuracy of the assemblies.

To annotate repetitive elements, we employed a combination of homology-based and de novo approaches. For the homology-based search, RepeatMasker v.4.1.2 (http://www.repeatmasker.org/) and RepeatProteinMask were used against the Repbase library (Bao et al., 2015) to identify known transposable elements (TEs). De novo repeat identification was performed using RepeatModeler v.2.0.1 (Flynn et al., 2020) and LTR_FINDER v.1.0.7 (Xu and Wang, 2007) to build a custom repeat library specific to each genome. Protein-coding gene prediction was conducted using a combination of ab initio, homology-based, and transcriptome-based methods. For homology-based prediction, protein sequences from eight related plant species (Arabidopsis thaliana, Henckelia pumila, Olea europaea, Primulina eburnea, P. huaijiensis, P. tabacum, Sesamum indicum, and Vitis vinifera; Table S1) were aligned to the Oreocharis genomes using BLASTp v.2.9.0+ (Camacho et al., 2009) with an E-value cutoff of 1e-5. Gene models were predicted using GeneWise v.2.4.1 (Birney et al., 2004) based on the alignments. RNA sequencing data from different tissues were mapped to the assembled genome using Cufflinks v.2.1.1 (Trapnell et al., 2012) and PASA (Haas et al., 2003).

2.3. Whole genome resequencing, mapping, and SNP calling

We collected 28 individuals from the three extant populations of Oreocharis esquirolii in Zhenfeng and Anlong Counties, covering its entire distribution range, and 79 individuals from 11 populations of O. maximowiczii from its broad geographical range (Fig. 1a). Genomic DNA was extracted from leaf tissue using the modified cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1987). Paired-end Illumina sequencing was performed using the Illumina NovaSeq 6000 platform with a read length of 150 bp. A total of 850.28 Gb and 3620.40 Gb of sequencing data were obtained for O. esquirolii and O. maximowiczii, respectively.

Raw sequencing reads were processed to remove adapter sequences and low-quality reads using Trimmomatic v.0.39 (Bolger et al., 2014). Reads with more than 10% ambiguous bases or over 50% of bases having a Phred quality score < 5 were discarded. High-quality reads were then mapped to their respective species-specific reference genomes using BWA-MEM v.0.7.18 with default parameters. SAMtools v.1.9 was used to sort and index the BAM files, and duplicate reads were marked using Picard Tools v.2.18.23 (http://broadinstitute.github.io/picard/).

Variant calling was performed using the Genome Analysis Toolkit (GATK) HaplotypeCaller v.4.1.2.0 (McKenna et al., 2010) in the GVCF mode for each individual. Joint genotyping was conducted using GATK GenotypeGVCFs to produce a raw variant set. SNPs were filtered with GATK VariantFiltration using the following criteria: Quality by Depth (QD) < 2.0, Fisher Strand Bias (FS) > 60.0, Mapping Quality (MQ) < 40.0, MQRankSum < −12.5, and ReadPosRankSum < −8.0. Additionally, indels and multi-allelic variants were excluded. SNPs with mean depth < 5 or > 100 and those with more than 50% missing data were also removed.

2.4. Population structure and gene flow analysis

We investigated the population structure of Oreocharis esquirolii and O. maximowiczii using a Bayesian clustering approach in ADMIXTURE v.1.23 (Alexander et al., 2009) and principal component analysis (PCA). For these analyses, we removed loci with minor allele frequencies (MAF) < 0.05 and used PLINK v.1.9 (Chang et al., 2015) to obtain unlinked markers with options --indep-pairwise 50 5 0.2. We ran Admixture with the number of genetic clusters (K) ranging from 1 to 10, performing 10 replicates for each K. Using the same data set, we conducted PCA with GCTA (Yang et al., 2011).

To estimate potential gene flow among populations of Oreocharis esquirolii and O. maximowiczii, we used VCFtools (Danecek et al., 2011) to generate allele frequency files, which we then converted to TreeMix format using the plink2treemix.py script (https://bitbucket.org/nygcresearch/treemix/downloads/). We performed population tree inference and modeled migration events using TreeMix (Pickrell and Pritchard, 2012), testing from 1 to 8 migration edges (m), with five independent runs for each value. The following parameters were applied: -k 500, -global, -noss, and random seeds for replicates. To determine the appropriate number of migration edges on the population tree, we inferred the optimal value of m using OptM (Fitak, 2021). To further assess potential admixture signals of O. esquirolii, we computed F₃-statistics using the threepop tool from the TreeMix package, based on the same allele frequency file. We calculated F₃ statistics in the form F₃ (A, B; C), where a significantly negative value indicates that population C is admixed between populations A and B. All possible combinations of the three groups were tested to identify signals of historical admixture events.

2.5. Linkage disequilibrium decay and genetic diversity analysis

To assess genome-wide linkage disequilibrium (LD) patterns, we calculated the squared correlation coefficient (r²) between pairs of SNPs using POPLDDECAY v.3.40 (C. Zhang et al., 2019). Genetic diversity was quantified through genome-wide heterozygosity (H), observed heterozygosity (H_O), expected heterozygosity (H_E), and nucleotide diversity (π). Individual genome-wide heterozygosity was determined by the ratio of heterozygous sites to the total number of callable sites using realSFS implemented in ANGSD (Korneliussen et al., 2014). H_O and H_E for each population and species were calculated using PLINK v.1.9. We estimated nucleotide diversity using PIXY v.1.2.5 (Korunes and Samuk, 2021), incorporating both variant and non-variant sites to account for missing data and genotype uncertainty, thereby providing unbiased estimates of π. We also calculated the fixation index (F_ST) between populations using PIXY and estimated Tajima's D within populations using ANGSD.

2.6. Inference of runs of homozygosity and inbreeding

Runs of homozygosity (ROHs) were identified to infer inbreeding events. We used PLINK v.1.9 with a sliding window approach to detect ROHs across the genome. No LD and MAF filtering was applied. The window size was set to 50 SNPs, allowing zero heterozygous calls per window. We only considered ROHs longer than 10 kb and categorized them into three classes: short ROHs (10 < ROH ≤ 100 kb), medium ROHs (100 < ROH ≤ 200 kb), and long ROHs (> 200 kb). The inbreeding coefficient (F_ROH) for each individual was calculated as the proportion of the genome contained within ROHs longer than 10 kb.

2.7. Genetic load estimation

To assess the genetic load and potential deleterious mutation burden, we categorized variants based on their predicted impact. The ancestral allele at each site was inferred as the major allele across all individuals, provided that more than 50% of individuals were homozygous for that allele. Functional annotations of SNPs were performed using SnpEff v.5.2c (Cingolani et al., 2012), classifying variants into three categories of mutations: (1) synonymous mutations, (2) missense mutations, and (3) loss-of-function (LOF) mutations.

The deleteriousness of missense mutations was evaluated using the Grantham score (Grantham, 1974), which considers the chemical differences between amino acids. Missense mutations with a Grantham score ≥ 150 were classified as deleterious, while those with scores between 5 and 150 were considered benign. LOF mutations included variants annotated as "stop_gained, " "start_lost, " and "splice_acceptor_variant." We used SnpSift (Cingolani et al., 2012) with the parameter "LOF [*].PERC = 1.00" to identify high-confidence LOF mutations.

The genetic load was quantified by calculating the ratio of deleterious mutations (missense deleterious and LOF mutations) to synonymous mutations. Functional enrichment analyses of genes affected by LOF mutations and deleterious mutations were conducted using Metascape (Zhou et al., 2019). Gene Ontology (GO) terms and KEGG pathways with enrichment factors greater than 2 and adjusted P-values less than 0.05 were considered significantly enriched.

2.8. Inference of demographic history

We used three methods to reconstruct changes in historical effective population size (Ne) for each species. We first applied pairwise sequentially Markovian coalescent (PSMC) model (Li and Durbin, 2011) to individual whole-genome sequences to infer demographic changes. To avoid bias from low coverage, we selected five individuals with the highest sequencing depth. PSMC was run with default setting (-N25 -t15 -r5 -p "4 + 25 × 2 + 4+6"), using a mutation rate of 6.5 × 10⁻⁹ per site per generation and a generation time of 2 years (Yi et al., 2022). We also used Stairway Plot v.2.0 (Liu and Fu, 2020), which infers recent population size fluctuation from site-frequency spectra (SFS). Unfolded SFS were generated from BAM files using realSFS in ANGSD v.0.933. To ensure robust demographic inference, we conducted multiple independent runs with different random seeds. Finally, we used GONE (Santiago et al., 2020) to estimate N_e over recent generations. This analysis used the observed spectrum of linkage disequilibrium (LD) of pairs of loci. We repeated the analysis 200 times for each species with default setting.

3. Results 3.1. Genome assemblies and annotations

We generated high-quality genome assemblies for the two Oreocharis species, O. esquirolii and O. maximowiczii, using a comprehensive sequencing and assembly workflow that included PacBio long reads, Illumina short reads, and Hi-C reads. For O. esquirolii, we obtained a total of 289.7 Gb of sequence data, while for O. maximowiczii, we generated 503.0 Gb of sequence data (Table S2). Based on estimated genome sizes of 1.19 Gb for O. esquirolii and 1.81 Gb for O. maximowiczii (Fig. S1 and Table S3), the sequence depth of coverage was approximately 233.63 × and 282.59 ×, respectively. The assembled genome sizes were 1.24 Gb for O. esquirolii and 1.78 Gb for O. maximowiczii, with scaffold N50 values of 69.1 Mb and 102.5 Mb, respectively (Fig. 2a and b; Table 1 and Table S4). Utilizing Hi-C data, we successfully anchored 94.81% of the O. esquirolii assembly and 98.13% of the O. maximowiczii assembly onto 17 chromosomes, respectively (Fig. S2 and Table S5), consistent with the known number of haploid chromosomes of Oreocharis (2n = 2x = 34).

BUSCO analysis identified 1580 (97.9%) conserved single-copy genes in the Oreocharis esquirolii assembly, and 1587 (98.3%) in the O. maximowiczii assembly among the total tested gene entries (Table S6), indicating high genome completeness. The reliability of the assemblies was further supported by high mapping rates of Illumina short reads to the assembled genomes, with 95.45% for O. esquirolii and 96.23% for O. maximowiczii. Repetitive elements constituted approximately 68.29% of the O. esquirolii genome and 77.72% of the O. maximowiczii genome (Table S7). For gene annotation, we integrated homologous protein mapping and RNA sequence prediction methods, resulting in the identification of 47, 637 protein-coding genes in O. esquirolii and 40, 872 in O. maximowiczii (Table S8). The large number of annotated genes and the comprehensive annotation approach contribute to the high confidence in our genome annotations. Furthermore, the two genome assemblies showed an overall high degree of collinearity, indicating that the two species exhibit relatively conserved genomic synteny (Fig. 2c).

3.2. Population structure and gene flow

We explored population structure within Oreocharis esquirolii and O. maximowiczii using Admixture based on unlinked synonymous SNPs. Admixture identified K = 3 and K = 4 as the optimal number of genetic clusters for O. esquirolii and O. maximowiczii, respectively, consistent with the PCA results (Fig. S3). O. esquirolii formed three distinct genetic clusters, each corresponding to one of its three known natural populations. The 11 populations of the widespread O. maximowiczii clustered into four groups based on geographic distribution: (1) a northern group consisting of four populations from Zhejiang, Jiangxi and Fujian; (2) a western group of four populations from Jiangxi; (3) an eastern group of two populations from Fujian; and (4) a southern group with one population in Guangdong (Fig. S3).

TreeMix analyses revealed no evidence of gene flow among Oreocharis esquirolii populations (Fig. S4). F₃-statistics for all tested population triplets yielded significantly positive F₃ values (Table S9), further confirming the absence of historical admixture among O. esquirolii populations. In contrast, TreeMix analysis indicated historical gene flow among populations of O. maximowiczii. The strongest migration signal was detected between an ancestral lineage and a population in northeastern Jiangxi Province (GX). Gene flow was also detected between two additional populations from Jiangxi (ND and NF), suggesting local admixture (Fig. S4).

3.3. Low genetic diversity in the endangered Oreocharis esquirolii

To characterize genetic diversity patterns within and between species, we performed whole-genome resequencing on 28 Oreocharis esquirolii and 79 O. maximowiczii individuals. Each sample yielded 80–173 million paired-end reads. After alignment to the reference genomes, sequencing depth ranged from 9.92 to 26.23× (Table S10). After filtering, we identified 16, 563, 591 and 27, 374, 400 high-confidence SNPs for O. esquirolii and O. maximowiczii, respectively.

As expected, the nucleotide diversity (π) of the endangered species Oreocharis esquirolii (0.0054) is about one-fifth that of the widespread O. maximowiczii (0.0259; Fig. 2d). The estimated π within populations ranged from 0.0044 to 0.0051 for O. esquirolii, which is much less than that observed in the populations of O. maximowiczii (π ranging from 0.0067 to 0.0175; Table S11). Similarly, both the observed heterozygosity (H_O = 0.184) and expected heterozygosity (H_E = 0.202) of O. esquirolii were lower than those of O. maximowiczii (H_O = 0.279, H_E = 0.217) (Table S11). Genome-wide heterozygosity (H), which is an individual-level metric rather than a population-level metric like π, showed a similar pattern. In O. esquirolii individuals, genome-wide heterozygosity ranged from 0.44% to 0.81% (average: 0.58%), markedly lower than in O. maximowiczii individuals (2.84%–3.49%; average: 3.20%; Fig. 2e). These results suggest that the endangered species O. esquirolii has substantially lower genetic diversity compared to the widespread O. maximowiczii.

The pairwise genome-wide F_ST values among the three populations of Oreocharis esquirolii ranged from 0.153 to 0.192, whereas those for O. maximowiczii ranged between 0.262 and 0.687 (Table S12), suggesting a strong level of genetic differentiation in the widespread species.

3.4. Different genomic signatures of inbreeding between Oreocharis esquirolii and O. maximowiczii

To compare inbreeding across the genomes of the two species, we quantified the extent of runs of homozygosity (ROH) grouped into three size categories using PLINK (Purcell et al., 2007). Long ROH (> 0.2 Mb) likely result from recent close inbreeding due to sharp population declines, whereas short (≤ 0.1 Mb) and medium (0.1–0.2 Mb) ROHs may reflect ancient population declines, as there would be more time for recombination to break up long segments (Ceballos et al., 2018). In both species, the majority of the ROHs were ≤ 0.1 Mb in length (Fig. 3a), while medium and long ROHs (> 0.1 Mb) accounted for only a negligible proportion of the total (Fig. 3b and c). This result suggests that most inbreeding in both species did not occur recently. Nevertheless, short and medium ROHs were in greater abundance in the endangered O. esquirolii than in the widespread O. maximowiczii (Fig. 3a and b). We evaluated the individual inbreeding coefficient (F_ROH) by assessing the proportion of ROHs longer than 10 kb. The level of inbreeding in O. esquirolii (F_ROH = 0.23; range: 0.11–0.35) was significantly higher than in O. maximowiczii (F_ROH = 0.10; range: 0.01–0.16; Fig. 3d and Table S13). Additionally, there was a significant negative correlation between the observed individual heterozygosity and F_ROH in both species (Figs. 3e, f and S5). Genome-wide linkage disequilibrium (LD) also differed significantly between O. esquirolii and O. maximowiczii. Increase in homozygosity is expected to reduce effective recombination, leading to elevated levels of LD. Indeed, LD decayed more slowly in O. esquirolii, with half the maximum r² reached at approximately 531 kb, compared to around 13 kb in O. maximowiczii (Fig. S6).

3.5. The accumulation of mutation load

To understand the effect of inbreeding and low genetic diversity on the distribution of deleterious variation, we estimated genetic load by counting the number of derived deleterious variants in coding regions and classified them into three categories: (1) loss of function (LOF) mutations, (2) deleterious mutations, and (3) benign mutations. We refer the combined total of these three types as total genetic load. The total genetic load of the endangered Oreocharis esquirolii was comparable to that of the widespread species O. maximowiczii. Consequently, there is no significant difference in total genetic load between the two species (Fig. 4a). However, the proportion of LOF and deleterious mutations was higher in the endangered O. esquirolii than in O. maximowiczii (Fig. 4b and c), indicating an accumulation of high- and moderate-impact deleterious variants. In contrast, the proportion of slightly deleterious (benign) variants were higher in O. maximowiczii than in O. esquirolii (Fig. 4d).

We further analyzed genes affected by LOF and deleterious mutations across all 28 individuals of Oreocharis esquirolii. Significantly enriched Gene Ontology (GO) pathways indicated that these genes were associated with monoatomic anion transmembrane transport (GO: 0098656; P < 0.05), chloroplast RNA processing (GO: 0031425; P < 0.05), and biological regulation (GO: 0065007; P < 0.05) (Figs. S7 and S8; Tables S14 and S15). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that "Environmental Information Processing (P < 0.05)" and "Biosynthesis of other secondary metabolites (P < 0.05)" were enriched in LOF mutations. In contrast, "Plant-pathogen interaction (P < 0.05)" and "MAPK signaling pathway - plant (P < 0.05)" were enriched in deleterious mutations (Figs. S7 and S8; Tables S16 and S17). These results suggest that genes with LOF and deleterious mutations may influence environmental adaptation and cell signaling in O. esquirolii.

3.6. Contrasting pattern of demographic history

To investigate the population dynamics of Oreocharis esquirolii and O. maximowiczii, we utilized three methods to reconstruct the demographic histories of both species. PSMC analysis of O. maximowiczii revealed a dynamic evolutionary history, with a marked decrease in effective population size (N_e) beginning at the end of the Naynayxungla Glaciation (NG, 0.7–0.5 Mya), followed by a modest growth during the Last Glacial Maximum (LGM, 26.5–19 kya), and eventually stabilization. In contrast, O. esquirolii experienced three significant declines in N_e, each corresponding to periods of temperature decreases during the three glaciation periods. After the LGM, the N_e of O. esquirolii reached critically low levels, while the population size of O. maximowiczii expanded over the last ~25, 000 years (Fig. 5). A similar pattern of population size change was detected by the Stairway Plot 2 method. Oreocharis maximowiczii underwent a population bottleneck during the Last Glacial Interglacial (LGI, 100–75 kya), after which its effective population size recovered and remained relatively stable (Fig. 5). This decline coincided with significant climatic cooling, suggesting that O. maximowiczii was affected by glacial advances but managed to maintain relatively stable populations afterward. In contrast, O. esquirolii experienced two major declines in N_e: the first during the early phase of the LGI (~150 kya), and the second persisting throughout the entire LGM (Fig. 5). These declines in O. esquirolii correspond to periods of global cooling, indicating that O. esquirolii was particularly sensitive to climatic fluctuations. The prolonged decline suggests limited capacity for population recovery between glacial periods, possibly due to its restricted distribution and specific habitat requirements. We further employed the linkage disequilibrium-based GONE method to infer recent changes in N_e. GONE analysis detected N_e fluctuations within the last 1200 years. Both species showed similar recovery during the 400–1200 years ago (200–600 generations). Over the past 400 years, the population size of O. maximowiczii remained stable and consistently higher than that of O. esquirolii, which experienced a recent population decline (Fig. S9).

4. Discussion

In the present study, we assembled and annotated chromosome-level genomes for two Oreocharis species: the endangered O. esquirolii and the widespread O. maximowiczii. The alignment between the genomes of the two species showed a high degree of synteny, suggesting that both assemblies are structurally accurate and that the two species share conserved chromosomal organization. Using genome-wide SNP datasets derived from these high-quality reference genomes, we explored the population genomics of the endangered O. esquirolii in comparison with its widespread congener, O. maximowiczii. We showed that the endangered O. esquirolii had lower genomic diversity and higher levels of inbreeding than the widespread O. maximowiczii. Despite having a lower total mutation load, O. esquirolii has accumulated a higher proportion of strongly deleterious mutations. These observed differences are likely attributable to contrasting demographic histories. Specifically, the endangered O. esquirolii experienced prolonged population decline and consistently small effective population size (N_e). This study not only contributes to the genomic resources available for Gesneriaceae but also provides critical data to guide conservation strategies for these and other related species.

4.1. Demographic histories of the two Oreocharis species

Both PSMC and Stairway Plot analyses revealed a distinct pattern of demographic history for the two Oreocharis species. In the endangered O. esquirolii, our PSMC analyses identified three bottleneck events between 0.7 and 0.4 Mya and 30–10 kya (Fig. 5). These bottlenecks coincide with major climatic changes during the Naynayxungla Glaciation (NG), the Last Glacial Interglacial (LGI), and the Last Glacial Maximum (LGM). After these events, O. esquirolii exhibited a prolonged decrease in population size, whereas O. maximowiczii experienced a slight demographic decline followed by a period of stable N_e after the LGM. Although complex topography and microclimates could theoretically offer glacial refugia (Harrison and Noss, 2017), the observed bottlenecks in O. esquirolii suggest that habitat fragmentation during extreme climatic shifts may have overridden these effects. Specifically, the abrupt cooling and aridification during the LGM likely disrupted connectivity among microhabitats, leading to population declines despite localized refugia (Chen et al., 2024; Feng et al., 2024). The allopatric distributions—with O. esquirolii in Southwest China and O. maximowiczii in East China—suggest that their contrasting demographic histories reflect region-specific responses to Quaternary climatic fluctuations. The complex topography and microclimates of Southwest China may have intensified the impact of climatic changes on O. esquirolii, leading to more severe population bottlenecks. In contrast, the more stable environments of East China may have buffered O. maximowiczii against drastic demographic shifts. These geographic contrasts in demography align with broader biogeographic patterns. In Quercus chenii, Li et al. (2019) showed that postglacial expansions into lowland regions can enhance genetic diversity via admixture, while mountainous populations retained unique lineages and exhibited less postglacial recovery. Similarly, O. maximowiczii in the East China lowlands likely benefited from interglacial range expansions that replenished its effective population size, whereas O. esquirolii in the mountainous Southwest remained constrained in fragmented refugia.

4.2. Comparison of genetic diversity

Using a comparative approach, we found that the endangered Oreocharis esquirolii had significantly lower levels of genetic diversity (in terms of π and heterozygosity) at both the species and individual population levels compared to the widespread species O. maximowiczii (Fig. 2d and e). This result is consistent with expectations when comparing rare species to their widespread congeners and has also been reported in other congeneric pairs (Gitzendanner and Soltis, 2000; Cole, 2003). Previous studies have revealed that patterns of genetic diversity vary depending on many factors related to geographic range, effective population size, life form, mating and breeding systems, and demographic processes (Hamrick and Godt, 1996; Ellegren and Galtier, 2016). However, it is often difficult to disentangle which of these factors have contributed most to current levels of genetic diversity, due to the interplay of multiple confounding variables. Comparative studies of rare and common congeneric species can help control for confounding effects of phylogeny and life history on genetic diversity (Gitzendanner and Soltis, 2000; Cole, 2003; Boyd et al., 2022). In this study, since the endangered O. esquirolii shares common life-history traits with its widespread congener O. maximowiczii, its lower genetic diversity is likely attributable to reduced population size or geographic isolation. Consistent with this, our demographic analyses revealed that O. esquirolii experienced two severe bottlenecks approximately 100 kya and 25 kya (Fig. 5), resulting in a history of persistently small effective population size.

Several studies have demonstrated that the fraction of the genome in runs of homozygosity (F_ROH) is a reliable estimator of inbreeding (Kardos et al., 2018; Alemu et al., 2021; Caballero et al., 2021; Nietlisbach et al., 2019). Genome-wide estimates indicated low inbreeding in O. maximowiczii (average F_ROH = 0.10; Fig. 3d), but substantially higher in O. esquirolii (average F_ROH = 0.23; Fig. 3d). This higher level of inbreeding observed in O. esquirolii is not surprising, given its small and isolated populations. Furthermore, a negative relationship was observed between the estimated inbreeding coefficient and individual heterozygosity in both species (Fig. 3e and f). This suggests that the observed difference in genetic diversity can be primarily driven by variation in inbreeding between the two species. In addition to quantifying inbreeding, the distribution of ROHs can shed light on past demography (Kirin et al., 2010; Pemberton et al., 2012; Ceballos et al., 2018): long ROHs are indicative of recent inbreeding, whereas short segments reflect more distant inbreeding events. The dearth of long ROHs in both O. esquirolii and O. maximowiczii suggests that they were most likely not affected by recent inbreeding (Fig. 3c). Nevertheless, the distribution of ROHs across different length categories differed significantly between the two species, likely reflecting their different inbreeding histories. In the widespread O. maximowiczii, the total length of ROHs accounts for only a negligible fraction of the genome (average F_ROH = 0.10; Fig. 3d), consistent with low levels of inbreeding. In the endangered O. esquirolii, the majority of ROH fragments were restricted to short and medium lengths (Fig. 3a and b), indicating ancient inbreeding events. This finding is consistent with the long-term small population size in O. esquirolii. Together, these findings suggest that the reduced genetic diversity in O. esquirolii is primarily the result of elevated inbreeding due to prolonged isolation and small population size.

4.3. Mutation load

Mutation load provides an important way of assessing the exposure to genetic threats in small populations (von Seth et al., 2021). Theoretical work suggests that small and isolated populations are susceptible to the accumulation of deleterious alleles, i.e., mutation load (Ohta, 1973), but this pattern does not always emerge in natural populations. For example, van der Valk et al. (2022) found that mammal populations with historically small effective population sizes carried significantly lower genetic load than larger populations, possibly due to long-term purging of deleterious variants. In agreement with this finding, we showed that the endangered Oreocharis esquirolii has a significantly lower total mutation load than the widespread O. maximowiczii (Fig. 4a–c), suggesting that purging of moderately deleterious alleles may have been an important evolutionary force. However, despite carrying a lower total mutation load, we observed that O. esquirolii has a significantly higher proportion of loss-of-function (LOF) mutations than O. maximowiczii (Fig. 4b). This suggests that purifying selection was ineffective against highly deleterious mutations in O. esquirolii. This is contrary to several recent genomic studies showing that small populations may carry relatively fewer strongly deleterious mutations than larger populations due to purging in small wild populations (Xue et al., 2015; Robinson et al., 2016, 2018, 2022; Yang et al., 2018; Grossen et al., 2020; Dussex et al., 2021; Khan et al., 2021; Ochoa and Gibbs, 2021; Kleinman-Ruiz et al., 2022; Mooney et al., 2023). The effectiveness of genetic purging in small populations of endangered species remains controversial (Hedrick and Garcia-Dorado, 2016; Xie et al., 2022). Purifying selection may become less effective at removing deleterious variants from a population during periods of small population size due to increased genetic drift (Ohta, 1973). In endangered species with extremely small population sizes like O. esquirolii, strong genetic drift may overwhelm purifying selection, resulting in the fixation of strongly deleterious mutations across the genome.

Empirical comparisons of endangered and widespread species illustrate these dynamics. For example, the critically endangered Rhododendron griersonianum has significantly more homozygous deleterious variants than its widespread congener R. delavayi, reflecting limited purging in the small population (Ma et al., 2021). Similarly, the brown eared pheasant (Crossoptilon mantchuricum) accumulated extensive runs of homozygosity and excessive loss-of-function mutations as a consequence of prolonged population decline (Wang et al., 2021). These cases illustrate that, when N_e is chronically low, genetic drift can overpower selection and permit harmful mutations to persist or even fix. In Oreocharis esquirolii, our finding of an unusually high proportion of LOF alleles thus underscores the fragility of purifying selection under persistent decline. Together, these studies emphasize the complexity of mutation accumulation in natural populations and the need for additional empirical work to understand how deleterious mutations are purged or fixed in endangered species with complex histories.

4.4. Conservation implications

Here we present a comprehensive analysis of conservation genomics of the critically endangered Oreocharis esquirolii. By assembling a chromosome-level reference genome and sequencing whole genomes from 28 individuals across all known populations, we gained valuable insights into the genetic status of the species. Our results indicate that O. esquirolii is severely threatened, exhibiting prolonged population decline, low genetic diversity, elevated inbreeding, and a diminished capacity to purge deleterious mutations. These factors substantially compromise the species' long-term evolutionary potential and its ability to respond effectively to future stochastic events and environmental changes. Given that O. esquirolii is restricted to just three populations within a limited geographic area, the preservation of its remaining habitat is critical for long-term recovery. Special attention should be given to population Anlong County, Guizhou (AL), which is located outside protected areas and is in close proximity to a roadside and a small dam, exposing it to anthropogenic disturbances. Prioritizing in situ conservation for this population is essential, and establishing mini-reserves may provide a cost-effective approach to mitigating threats. Additionally, enhancing community awareness and engagement will be vital for the long-term success of conservation efforts. Given the low genetic differentiation and geographic proximity among populations, individuals from different populations could be used interchangeably for restoration or augmentation without significantly impacting genetic integrity. In practice, seed banking and the establishment of living collections in botanical gardens can serve as genetic insurance for O. esquirolii. Overall, an integrated strategy combining strict habitat protection, active management of remaining populations, and ex situ conservation is recommended to maximize the long-term survival and evolutionary potential of this critically endangered species.

Acknowledgement

This work was supported by National Key R&D Program of China (2024YFF1307400) and Guangdong S&T Program (2022B1111230001).

CRediT authorship contribution statement

Nana Peng: Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Lihua Yang: Resources, Investigation. Xizou Shi: Investigation. Hanghui Kong: Investigation. Ming Kang: Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition, Conceptualization.

Data availability

The data supporting the findings of this work are available within the paper and Supporting Information. All of the raw sequence reads and assemblies described in this paper have been submitted to the National Center for Biotechnology Information (NCBI) under accession codes PRJNA1221333 (genome assembly of O. esquirolii), PRJNA1221335 (genome assembly of O. maximowiczii) and PRJNA1224281 (whole-genome sequencing).

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