1. Introduction

The limestone karst formations in southwestern China, characterized by unprecedented edaphic and topographic heterogeneity, host a global biodiversity hotspot with a diverse array of unique and endemic species (Wang et al., 2014; Hao et al., 2015; Pu et al., 2017). The soils of the area are characterized by shallow profiles, limited moisture and nitrogen, while the sedimentary rock outcrops are mainly composed of calcium carbonate, which exerts significant evolutionary pressure on karst plant populations (Zhang et al., 2011; Kang et al., 2015; Wu et al., 2022). The complex karst topography functions as a terrestrial archipelago with small and isolated populations in cave or cave-like microhabitats (Sarthou et al., 2001; Clements et al., 2006; Tseng et al., 2019). Small and isolated populations often experience reduced genetic diversity, increased inbreeding, and increased genetic load, which can lead to reduced fitness and an increased risk of extinction (Spielman et al., 2004; Jangjoo et al., 2016; Ma et al., 2021a; Rodger et al., 2021).

Recent advances in sequencing technologies and bioinformatic tools have revolutionized the field of conservation genetics, providing unprecedented opportunities for comprehensive genome-wide analysis of nucleotide variation (Davey et al., 2011; Funk et al., 2012; Bertorelle et al., 2022; Paez et al., 2022; Robinson et al., 2023). We can now investigate genetic structure, divergence trajectories, historical population fluctuations, and the impact of inbreeding at the genomic scale. Additionally, these advances allow us to assess the risk of extinction due to deleterious mutations in small populations (Liu et al., 2020; Hohenlohe et al., 2021). These breakthroughs in genomic research have already shed light on the genetics of threatened species, and provided invaluable insights into the conservation needs of several endangered animals and plants, such as Rhododendron (Ma et al., 2021a), Acer yangbiense (Ma et al., 2021b), pangolins (Hu et al., 2020; Wang et al., 2022), killer whales (Kardos et al., 2023), and scimitar-horned oryx (Humble et al., 2023). Despite significant advances in genomic research on endangered animal species, only a few plants have been subjects of conservation genomic studies (Yang et al., 2018; Li et al., 2020; Feng et al., 2024; Tao et al., 2024), and there is currently no knowledge of genome-wide inbreeding and mutation load in karst plants. Understanding genome-wide diversity, inbreeding, and the burden of accumulated deleterious mutations (i.e., mutation load) in small populations is essential for predicting and increasing population persistence and resilience.

Begonia L. (Begoniaceae, Cucurbitales), one of the largest angiosperm genera, thrives in tropical forests and has significant economic value as an ornamental plant (Goodall-Copestake et al., 2010). The limestone karst landscapes in southwestern China and northern Vietnam support a remarkable diversity of Begonia species (Peng et al., 2013), which thrive in moist, shaded areas such as karst walls, cave entrances, and watersheds. Previous studies on Begonia have primarily focused on resource surveys, morphological anatomy, cytology, and phylogenetic relationships (Brennan et al., 2012; Moonlight et al., 2017; Emelianova and Kidner, 2022; Oginuma and Peng, 2022; Tsai et al., 2022). Although a few studies have used microsatellite markers or chloroplast DNA to investigate population genetic structure in several Begonia species (Goodall-Copestake et al., 2010; Twyford et al., 2013; Chung et al., 2014; Xiao et al., 2023), genome-wide population genetic studies have yet to be reported in this mega-diverse genus. The recent release of a high-quality nuclear genome of Begonia masoniana facilitates population and conservation genetic studies in this group (Li et al., 2022).

In this study, we focus on the Begonia masoniana complex, which consists of four closely related species (B. liuyanii, B. longgangensis, B. masoniana and B. variegata). B. liuyanii, B. longgangensis and B. masoniana are endemic to southwestern Guangxi, near the border with northern Vietnam (Fig. 1a), whereas B. variegata is endemic to northern Vietnam. The four species of the B. masoniana complex are morphologically distinct (Peng et al., 2005, 2013). Geographically, these species are distributed in similar habitat with small and scattered ranges in limestone karst landscapes. Due to habitat destruction caused by human activities, B. liuyanii, B. longgangensis and B. masoniana are currently classified as vulnerable (VU) species in China (Qin et al., 2017). Understanding patterns of genetic diversity, levels of inbreeding, and genetic load will help guide future conservation efforts.

Here we investigate genomic divergence and mutation load in the Begonia masoniana complex using whole-genome re-sequencing data. We aimed to (1) clarify population structure and phylogenetic relationships, with a particular focus on understanding the contribution of gene flow, especially in species pairs with sympatric distributions; (2) infer demographic history, including population divergence, fluctuations in effective population size (N_e), and their potential correlations with climatic fluctuations; and (3) assess levels of inbreeding and mutation load accumulation in these species. Our results have important implications for the development of more targeted and effective conservation strategies for small and isolated populations in karst regions.

2. Materials and methods 2.1. Population sampling, genome sequencing and variant calling

We collected 71 samples from different locations in the karst landscapes of southern China (Fig. 1a). These samples included 62 individuals of Begonia liuyanii (21), B. longgangensis (16), B. masoniana (22) and B. variegata (3) (Table S1). Nine individuals of B. ferox were used as the outgroup for phylogenetic analyses. The distribution ranges of these species were highly restricted, with B. liuyanii and B. longgangensis displaying geographic sympatry, while individuals of B. masoniana are allopatric to the other species (Fig. 1a). The three individuals of B. variegata were collected from northern Vietnam and cultivated at Kunming Institute of Botany, Chinese Academy of Sciences, however, the specific location of these samples is unknown.

Genomic DNA was extracted from dried leaves using a customized cetyltrimethylammonium bromide (CTAB) method (Doyle, 1987). Each sample was re-sequenced on Illumina Novaseq 6000 PE150 platform using paired-end sequencing libraries. Quality control was performed using FastQC v.0.11.8 followed by trimming of raw reads using TRIMMOMATIC v.0.39 (Bolger et al., 2014; Chen et al., 2018). Clean reads were aligned to the Begonia masoniana reference genome using the BWA-MEM algorithm within BWA v.0.7.15 (Li, 2013). Then, SAMTOOLS v.0.1.19 was used to convert sequence alignment files in sequence alignment map (SAM) format into binary alignment map (BAM) format and sort them by chromosome number, and calculated depth and mapping rate of each sequence (Li et al., 2009). Picard v.2.20.2 (DePristo et al., 2011) was then used to mark and remove PCR duplicates. The resulted BAM files constitute Dataset 1 for this research. Single nucleotide polymorphisms (SNPs) and genotypes were called using the Genome Analysis Toolkit v.4.1.2.0 (GATK) (McKenna et al., 2010). Only bases with a quality score of ≥ 30 were included and all indels and variants containing more than two alleles were excluded. To minimize the influence of mapping bias, we further discarded the following sites: (i) sites with anomalous depth (mean depth < 5 or > 100), and (ii) sites with more than 30% missing data (genotyped in fewer than 70% of all individuals). A total of 102,645,147 sites were finally obtained (Dataset 2). SNPs with a minor allele frequency < 0.05 were then also removed using VCFTOOLS v.0.1.13 (Danecek et al., 2011), leaving 2,711,146 SNPs (Dataset 3) for downstream analysis. Next, a moderate linkage disequilibrium (LD) pruning was performed using PLINK v.1.90, removing SNPs within a window of 50 SNPs with r² > 0.5, resulting in 574,987 SNPs (Dataset 4).

2.2. Population genetic structure and genetic parameters

We used a Python script (https://github.com/edgardomortiz/vcf2phylip) to convert the SNPs from Dataset 3 into concatenated sequences for all individuals. Then, we constructed a maximum likelihood (ML) phylogenetic tree using IQ-TREE v.1.6.12 based on these concatenated sequences, using Begonia ferox as an outgroup (Nguyen et al., 2015). Bootstrap values were calculated from 1000 repetitions and subsequently visualized with FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). Next, the identity-by-state (IBS) kinship matrix between individuals was analyzed based on Dataset 2 using VCF2PCACluster v.1.1.40 from github (https://github.com/hewm2008/VCF2PCACluster).

Subsequent analyses were based on Dataset 4. For population structure analysis, we used ADMIXTURE v.1.3.0 (Alexander and Lange, 2011). For each K value ranging from 1 to 10, we conducted 10 replicates and determined the optimal K based on the cross-validation error. To confirm the result, we performed principal component analysis (PCA) on all individuals using PLINK v.1.90 to examine genetic structure and relatedness. In addition, we used the R package ggplot2 to plot the results (Villanueva and Chen, 2019).

We also calculated the unbiased nucleotide diversity (π), absolute divergence (d_xy) and fixation statistics (F_ST) by PIXY v.1.2.6 (Korunes and Samuk, 2021) using Dataset 2 through 200-kb windows. The values of Tajima's D were then calculated using VCFTOOLS v.0.1.13. Whole genome heterozygosity (H_e) represents the proportion of heterozygous sites in the whole genome, excluding any gaps, and we used ANGSD v.1.9 software to calculate the H_e of each individual (Korneliussen et al., 2014).

2.3. Inferring demographic history

To reconstruct the demographic history using the pairwise sequential Markovian coalescent (PSMC) model (Li and Durbin, 2011), we selected the individual with the highest sequencing depth and three additional individuals from different populations each of the species (Figs. 2a and S3). Parameter -p (which specifies the atomic time intervals) was set as “4 + 25*2 + 4 + 6” and other parameters were set as default. The generation time was set at 2 years, and the mutation rate (μ) was set to 6.5 × 10^-9 per site per year, following the estimates in rice (Molina et al., 2011). A bootstrapping approach with 100 replicates was performed to evaluate the variation in the inferred N_e trajectories (Fig. 2a).

In addition, we aimed to confirm the specific demographic history using dadi based on site frequency spectrum (SFS) calculations. However, for an accurate estimation of population history based on SFS calculation, it is advisable to have a sample size more than 10 individuals (Gutenkunst et al., 2009; Beichman et al., 2017). Therefore, B. variegata was excluded from this analysis due to small sample size (N = 3). We used the Python script easySFS.py (https://github.com/isaacovercast/easySFS) to estimate the folded SFS based on LD-pruned SNPs data without missing sites (Gutenkunst et al., 2009). We then analyzed the SFSs of the three species using the diffusion approximation method of dadi and fitted 20 demographic models (Fig. S4), following the demographic modeling workflow (dadipipeline) described by Portik et al. (2017). First, we performed a comprehensive optimization process (dadi_Run_Optimizations.py), which included model fitting with specific configurations for a given number of replicates. Next, we used the parameters derived from the best-performing repetition to initiate a subsequent round of model fitting with updated configurations. This optimization process consisted of four rounds, with each round consisting of 10, 20, 30, and 40 repetitions, respectively. The four-round optimization process contributed significantly to reducing discrepancies in the likelihood values obtained at the end of the runs. Finally, to establish confidence intervals (CIs) for the parameters derived from the demographic models with the highest likelihood, we used a procedure that involved the generation of 100 bootstrap replicates of the spectra.

2.4. Estimating gene flow and inbreeding

To identify historical introgression events, Patterson's D-statistic (Durand et al., 2011) was performed on the whole genome of the four species using Dtrios implemented in Dsuite v.0.5 r48 (Malinsky et al., 2021). The D-statistic is derived from the assessment of the proportional occurrence of two different topological configurations in variable genetic sites, referred to as “ABBAs” and “BABAs”. In summary, in the context of a given relationship, represented as (((P1, P2) P3) P4), the genetic alleles hosted by P4 (the outgroup) are designated as the ancestral state A, while the altered state is denoted as B (Durand et al., 2011). To estimate contemporary migration rates, we used the Bayesian MCMC method of BayesAss3-SNPs (BA3-SNPs) (Mussmann et al., 2019; Wilson and Rannala, 2003). Delta values (-m, -a, and -f) were each set to 0.5, and we ran the analysis for 10,000,000 iterations with a burn-in of 1,000,000 steps, and 1000 sampling intervals. These analyses were performed on Dataset 3.

Individual-based F_IS values were calculated with Dataset 2 by using VCFTOOLS v.0.1.13. Population scale inbreeding coefficients were calculated using PLINK v.1.90 for SNPs across different lineages. To further investigate inbreeding, we used PLINK v.1.90 to identify runs of homozygosity (ROH), retaining those longer than 100 kb (Yang et al., 2018; Ma et al., 2021a). The inbreeding coefficient based on ROH (F_ROH) was calculated for each individual as the total length of all its ROH divided by the total length of the autosomal genome (herein 799,826,415 bp).

2.5. Characterization of deleterious mutations and gene annotations

To estimate the genetic load in four species of the Begonia masoniana complex, we annotated SNPs across the entire genome using SNPEFF v.4.3t (Cingolani et al., 2012). Individuals of B. ferox were used as an outgroup to determine the ancestral state for each SNP. The ancestral state was defined as the nucleotide that was homozygous in the outgroup individuals and matched at least one of the alleles within the ingroup (Chen et al., 2020b). We then classified each species within the ingroup into heterozygous mutation sites and homozygous mutation sites based on the ancestral and derived states from the outgroup. We followed the Grantham score (GS) rating criteria (Grantham, 1974), which quantifies the occurrence of three distinct mutation types: synonymous mutations, tolerated mutations, and deleterious mutations. Sites with scores ranging from 5 to 150 were designated as “tolerated” mutation sites, while SNPs with a GS ≥ 150 were categorized as “deleterious mutation” sites (Li et al., 1984). However, mutations that significantly affect the transcription of the entire gene are considered highly deleterious, and from this point on, we referred to them as loss-of-function (LOF) mutations. We then calculated the ratios of three types of harmful mutations to synonymous mutations and all derived sites, respectively.

Finally, we characterized the detected LOF mutations by gene ontology (GO) analysis and Kyoto encyclopedia of genes and genomes (KEGG) according to the gene annotation results of the previous study (Li et al., 2022) using TBtools (Chen et al., 2020a). Analyses were performed for our selected candidate genes, focusing on those genes associated with statistically significant terms (P < 0.05).

3. Results 3.1. Sampling and sequencing

Raw reads were trimmed and mapped onto the Begonia masoniana reference genome, resulting in a mean depth of 15.90× (Table S1) and generating 71 BAM files (e.g., Dataset 1). Then SNPs were called using HaplotypeCaller in GATK and filtered using stringent criteria (hard filter in GATK and VCFTOOLS), which finally generated Dataset 2 (102,645,147 sites). After using maf 0.05 and mac 1, we retained 2,711,146 high quality nuclear SNPs as Dataset 3, which were further reduced to 574,987 SNPs through LD pruning (Dataset 4).

3.2. Population genetic structure and differentiation

Admixture identified K = 3 as the optimal number of clusters (Fig. S1). At K = 3, the populations of Begonia longgangensis and B. masoniana formed into two distinct clusters, while the third cluster includes both B. liuyanii and B. variegata, with individuals of B. variegata showing substantial admixture (Fig. 1b). However, when K = 4, four distinct clusters were observed for the corresponding four species (Fig. S1). Nevertheless, we observed genetic admixture in several individuals, which may be indicative of potential hybridization between species (Fig. 1b). The inferred phylogenetic relationships based on nuclear SNPs revealed four genetic clusters with strong support (Fig. 1b), which is consistent with the pattern identified by structure analysis. Such a genetic pattern was also supported by PCA analysis (Figs. 1c and S2).

Pairwise estimates of genetic differentiation between species ranged from 0.1347 (Begonia masoniana vs B. liuyanii) to 0.3026 (B. variegata vs B. liuyanii) for F_ST, and from 0.0113 (B. longgangensis vs B. liuyanii/B. masoniana) to 0.0124 (B. longgangensis vs B. variegata) for d_xy (Fig. 1d). These results suggest pronounced genomic divergence among these species. In addition, intraspecific differentiation (F_ST) is higher in B. masoniana (0.1437) compared to B. liuyanii (0.0946) and B. longgangensis (0.0827) (Table S3).

3.3. Demographic history and gene flow

The PSMC analysis suggested that all four species had experienced a severe bottleneck in response to global temperature fluctuations throughout a series of glacial and interglacial periods (Figs. 2a and S3). Specifically, Begonia liuyanii and B. masoniana exhibited a significant expansion prior to Naynayxungla glaciation (NG, 0.5–0.72 million years ago [Mya]) (Zheng et al., 2001), followed by a sharp decline in N_e until the end of the last glacial maximum (LGM, 10–20 thousand years ago [kya]) (Bai et al., 2018), with the most substantial decline occurring in B. masoniana (Figs. 2a and S3). Additionally, the positive values of Tajima's D in each group may also indicate a population bottleneck in their evolutionary history (Fig. 1d and Table S2). This long-term bottleneck corresponded with global temperature fluctuations in the karst region (Zachos et al., 2008). Moreover, around ten thousand years ago, the effective population size (N_e) of B. masoniana reached its minimum among the four species, with B. masoniana having a N_e less than twenty thousand (Figs. 2a and S3).

We used coalescent-based simulations in dadi to estimate gene flow and divergence events amongst the three species (Begonia liuyanii, B. longgangensis, B. masoniana) (Tables S4 and S5). The best fitted model fit the actual data relatively well (Figs. 2b and S5). These results suggest that B. liuyanii diverged from ancestor A approximately 88.7 kya (95% CIs: 66.1–11.1 kya) with high gene flow with ancestor B. About 72.4 kya (95% CIs: 57.2–87.6 kya), B. longgangensis and B. masoniana also diverged, with only B. longgangensis maintaining a more pronounced historical gene flow with B. liuyanii. Then about 11,000 years ago, three species became isolated from each other (Fig. 2b). We also found that the N_e of B. masoniana was the largest during the early stages of divergence, around 70,000 years ago (Fig. 2b; Table S5), a pattern also observed in the PSMC analysis (Fig. 2a).

Gene flow is present among four closely distributed lineages, albeit with low contemporary migration rates (Fig. 2c and d). Dsuite detected three instances of gene flow (Fig. 2c and Table S6). Notably, gene flow between Begonia liuyanii and B. longgangensis reinforces the evidence for genetic admixture in the LLG2 population observed in the admixture analysis, and the IBS analysis also confirms the close genetic relationship between LLG2 of B. liuyanii and B. longgangensis (Figs. 1b and S6; Table S7). Moreover, the BA3-SNP analysis reveals diverse and asymmetric contemporary migration rates among the four species, although they all remain relatively low (Fig. 2d).

3.4. Inbreeding and genetic diversity

Amongst the four species, we observed the highest levels of inbreeding in Begonia masoniana (Fig. 3 and Table S8). The F_IS for B. masoniana is 0.212, while for the other species it is less than 0.1 (Table S8). Similar trends were also observed in F_ROH values, which were found to be 0.116 for B. liuyanii, 0.106 for B. longgangensis, 0.263 for B. masoniana, and 0.030 for B. variegata, respectively (Fig. 3a; Table S8). Additionally, it is evident that individuals of B. masoniana have accumulated a greater number of long ROH segments. As a result, B. masoniana exhibits the lowest heterozygosity among the four species (mean heterozygosity rate = 0.01). However, all of the four species show a similar level of nucleotide diversity (π ranging from 0.005 to 0.007) (Fig. 3b and c; Tables S2 and S8).

3.5. Mutation load

The accumulation of genetic load varies among species, with Begonia masoniana notably carrying a relatively higher homozygous genetic load (Fig. 4 and Table S9). When compared with synonymous mutations, tolerated mutations predominate in both homozygous and heterozygous states across all species, with a lower occurrence of strongly deleterious mutations and a relatively small proportion of LOF mutations (Fig. 4a and b). Moreover, in B. masoniana, the average proportion of observed missense homozygosity is higher than in the other three species for three types of missense mutations, while its missense heterozygous sites were notably the least frequent when compared with the other three species (Fig. 4a and b; Table S9).

We further assessed mutation load by comparing the proportions of homozygous sites to all derived sites for three mutation types. We found that Begonia masoniana showed average proportions exceeding 0.6 for tolerated and deleterious types, with LOF mutations being more than twice as prevalent as those in other species (Fig. 4c).

Furthermore, we found a positive correlation between F_ROH levels and the total count of genome-wide heterozygous deleterious mutations in all species except Begonia masoniana (Fig. S7). In contrast, when assessing the association between F_ROH and the total number of homozygous deleterious mutations, only B. masoniana had a slight positive correlation without statistical significance (P = 0.59) (Fig. 4d). In addition, B. masoniana showed the highest proportion of homozygous deleterious mutations within the ROH segments compared to the other three species (Fig. 4e).

3.6. Functional enrichment of genes with loss-of-function mutations

A total of 435 genes in Begonia liuyanii, 467 genes in B. longgangensis, 668 genes in B. masoniana, and 96 genes in B. variegata were identified as carrying homozygous LOF mutations, in some cases the same genes with LOF mutations were found in multiple species. B. masoniana not only had the largest number of homozygous LOF mutant genes, but also had 533 unique homozygous LOF mutant genes (Fig. 4f and Table S10). In B. masoniana, homozygous LOF mutant genes were enriched in “oxidoreductase activity” (GO:0016627, P < 0.05) and pathways related to leaf, phyllome, shoot system, tissue and plant organ development (Fig. S8 and Table S11). GO pathway analysis identified “ATP-dependent activity (GO:0140657, P < 0.05)” as significantly enriched in homozygous LOF mutant genes in B. liuyanii, B. masoniana, and B. longgangensis. Moreover, the results indicated that these LOF mutation genes were significantly associated with phosphate metabolism and DNA repair in the four species (Figs. S8 and S9; Tables S11 and S12).

4. Discussion 4.1. Population structure of the Begonia masoniana complex

Increased geographic distance leads to decreased gene flow and consequently increased genetic differentiation between populations due to increased genetic drift (Bohonak, 2002). Karst landscapes in southern China, characterized by ubiquitous local-scale soil type mosaics, harbor a majority of species that are narrow endemics. Previous studies suggested that species divergence of karst plants has been influenced by uniquely heterogeneous habitats (Wang et al., 2017a; Zhu et al., 2022). Recent population genetic studies in Begonia species have revealed significant population structure due to geographic isolation (Matolweni et al., 2000; Forrest et al., 2005; Hughes and Hollingsworth, 2008; Thomas et al., 2012; Chung et al., 2014; Twyford et al., 2014). However, no study to date has attempted to assess whole-genome differentiation in this mega-diverse genus, nor has it been linked to edaphic or landscape factors.

Our admixture analysis, together with PCA and phylogenetic tree construction, revealed that the Begonia masoniana complex exhibits four distinct genetic clusters (Fig. 1). Furthermore, despite the close phylogenetic relationship between B. masoniana and B. longgangensis (Fig. 1b), they exhibited minimal gene flow (Fig. 2b and c) and no recent geographic overlap (Fig. 1a). It appears that the significant genetic structure among species can be attributed to the distinctive karst topography, leading to species isolation and subsequent speciation. These findings are consistent with previous observations documenting strong reproductive barriers among certain Begonia species in nature (Twyford et al., 2015).

We found genomic admixture between the close neighbors Begonia liuyanii and B. longgangensis (Fig. 1b), likely due to their close geographic proximity (Fig. 1a). Additional tests confirmed the occurrence of gene flow between these two species (Fig. 2). IBS analysis suggests a relatively close genetic relationship between the LLG2 population of B. liuyanii and B. longgangensis, further indicating hybridization between these species (Fig. S6 and Table S7). Although hybridization has been widely observed in Begonia species, it has been suggested that the majority of hybrid individuals belong to the F1 generation (Li et al., 2015; Tian et al., 2017, 2018). Therefore, despite the high genetic differentiation caused by karst habitat heterogeneity, natural hybridization can still occur between species in sympatry and parapatry.

4.2. Demographic history of the Begonia masoniana complex following climate-induced divergence and isolation

Our reconstructed demographic history suggests a consistent decline in population size after NG for the four species. Remarkably, the population size of Begonia masoniana appeared not to stabilize during interglacial periods, but to continue to decline throughout the LGM and Holocene (Figs. 2a and S3). This pattern likely reflects the influence of Quaternary glaciations on the B. masoniana complex. Previous studies suggested that historical climate temperature fluctuations may have influenced the speciation of karst plants and imposed constraints on their distribution (Kong et al., 2017; Ke et al., 2022). We used dadi to infer an optimal model of speciation events. However, previous studies have shown that a substantial sample size (N > 10) is essential for a reliable analysis of the SFS-based modeling (Gutenkunst et al., 2009; Beichman et al., 2017). Due to the unknown sampling locations and a small sample size (N = 3) for B. variegata, our analysis focused solely on the other three species. Notably, these divergence events are consistent with the timing of the last glacial period (LGP, 11.7–115 kya) (Chen et al., 2011). Furthermore, B. longgangensis and B. liuyanii are characterized by higher gene exchange compared to other species pairs. In addition, dadi analysis revealed that these three species continued to be isolated until approximately 11,000 years ago (Fig. 2b), which ultimately strengthened the formation of distinct new species.

Another notable finding is that after differentiation, Begonia masoniana virtually ceased gene flow with other species (Fig. 2b). The species was also observed to have a higher level of inbreeding than the other species, which may underlie the lower heterozygosity and more extended ROH in B. masoniana (Fig. 3). These results indicate that B. masoniana, which experienced the highest level of isolation and had the smallest recent N_e (Figs. 2a and S3), also suffered the most severe inbreeding. Small and isolated populations are at higher risk of experiencing inbreeding depression due to the increased likelihood of mating between closely related individuals (Bortoluzzi et al., 2020). It is reasonable to speculate that these patterns could be due to the prolonged isolation experienced by B. masoniana, potentially leading to increased inbreeding and the accumulation of deleterious mutations. These findings highlight the impact of Quaternary glacial climate fluctuations on the genetic differentiation and population structure of karst endemics.

Previous studies have found that endangered species with small populations have high levels of inbreeding, resulting in reduced N_e and low genetic diversity (Willi et al., 2006; Provan and Bennett, 2008; Jangjoo et al., 2016). However, contrary to our expectations, despite the more pronounced isolation and lower N_e in Begonia masoniana (Fig. 1, Fig. 2a), the genetic diversity of this species is similar to that of the other species (Fig. 3c). We speculate that the high genetic diversity may have resulted from large ancestral populations of B. masoniana accumulating a high level of genetic polymorphisms prior to Naynayxungla glaciation (Figs. 2a and S3). Relatively high genetic diversity has been reported for many endangered plants with small populations, such as Rhododendron meddianum (Zhang et al., 2021), Cupressus chengiana (Li et al., 2020), Tricyrtis ishiiana (Setoguchi et al., 2010), and Rhododendron protistum var. gigantum (Wu et al., 2014). Therefore, assessing the conservation status of species requires multifaceted approaches that consider not only genetic diversity within populations, but also demographic history.

4.3. The accumulation of deleterious mutations due to geographic isolation

Small and isolated populations that have experienced long periods of decline are particularly vulnerable to the influence of genetic drift (LaBar and Adami, 2017; Rodger et al., 2021; Humble et al., 2023). This, in turn, compromises the efficiency of purifying selection to purge deleterious variants (Whitlock, 2000; Wang et al., 2021). Consequently, these populations are susceptible to the accumulation of deleterious alleles, i.e., mutation load, which can lead to species extinction, also known as “mutation meltdown” (Poon and Otto, 2000; Charlesworth, 2009; Agrawal and Whitlock, 2012; Robinson et al., 2019, 2023).

In our study, we observed the accumulation of deleterious mutations in the Begonia masoniana complex, particularly in the case of B. masoniana. We found that compared to all synonymous mutations or to all derived alleles, B. masoniana had the highest accumulated proportion of three types of homozygous mutations (Fig. 4b and c). Under evolutionary constraints, homozygous derived allele mutations represent realized genetic load and are predicted to have fitness consequences (Marsden et al., 2016; Feng et al., 2019; van der Valk et al., 2019). This suggests that B. masoniana may face a higher risk of extinction compared to other species. Notably, our analysis shows a positive correlation between the accumulation of deleterious type homozygous mutations and F_ROH in B. masoniana (Fig. 4d), although the correlation is not statistically significant. In addition, we found that the proportion of deleterious homozygous mutations within ROH segments in B. masoniana is the highest among all four species (Fig. 4e). This suggests that an increased inbreeding within this species may contribute to the accumulation of genetic load, potentially affecting the fitness and long-term survival of the population. Empirical and simulation studies have shown that long-term small effective population sizes (N_e) could be conducive to the reduction of deleterious variation through the purging of strongly deleterious mutations, which could help the long-term survival of endangered species (Yang et al., 2018; Feng et al., 2019; van der Valk et al., 2019; Khan et al., 2021), but this may not hold true for B. masoniana in karst landscapes. However, the genetic differentiation (F_ST) within B. masoniana is higher than that of B. liuyanii and B. longgangensis (Table S3). It seems that B. masoniana populations may have experienced recent population size decline, which led to an increase of genetic drift and decrease of efficacy for purging deleterious mutations.

4.4. The functional enrichment of loss-of-function genes and conservation implications

We analyzed the genes affected by the homozygous putative deleterious LOF mutations in the four species and found that a number of GO and KEGG categories were involved in processes related to DNA repair and phosphate metabolism (Figs. S8 and S9; Tables S11 and S12). Notably, there is the marked enrichment of nitrogen metabolism pathways (KEGG, 00910) in Begonia longgangensis (Fig. S9 and Table S12), which may be intricately linked to the distinctive N, P-deficient karst soil environment (Wang et al., 2017b). These findings suggest a potential adaptive evolutionary strategy adopted by the B. masoniana complex in response to the unique karst soil conditions. In B. masoniana, we also observed a significant enrichment of LOF genes associated with oxidative reduction processes (GO, GO0016627). Moreover, we found that genes associated with GO pathways related to leaf, phyllome, shoot system, tissue and plant organ development were highly enriched in the LOF mutant genes (Fig. S8 and Table S11). This observation suggests a higher burden of deleterious mutations compared to other species, posing a significant threat to its long-term survival. Taken together with previous analyses (Grossen et al., 2020; Ma et al., 2021b; Wang et al., 2021), we propose that B. masoniana may be at higher risk of extinction due to low N_e, a higher genetic load, and prolonged isolation. Therefore, conservation efforts should prioritize B. masoniana. Introducing hybridization via exogenous breeding programs could mitigate the genetic erosion resulting from recent demographic events.

However, although hybridization between Begonia liuyanii and B. longgangensis can mitigate the accumulation of homozygous deleterious mutations, both species still carry a substantial burden of heterozygous deleterious mutations (Fig. 4b). In the event of further isolation, their genetic integrity could be rapidly compromised. The fragile karst ecosystem in Guangxi has witnessed an increasing influx of human activities, leading to a degree of anthropogenic interference with karst flora in recent years (Chen and Xu, 2023). Therefore, it is imperative to protect the ecological habitats of these species and mitigate the risk of anthropogenic habitat degradation. The results of this study highlight the importance of not relying solely on genetic diversity when assessing the conservation status of these plants, but also considering factors such as population size and accumulation of deleterious mutations for a comprehensive evaluation of conservation strategies.

Acknowledgment

We are grateful to Ms. Jifeng Long from Nonggang National Nature Reserve, Guangxi Zhuang Autonomous Region and Ms. Jingxiu Li from Kunming Institute of Botany, the Chinese Academy of Sciences for their kind help on sample collecting. This work was supported by Key-Area Research and Development Program of Guangdong Province (Grant No. 2022B1111230001) and National Natural Science Foundation of China (31860048).

Data availability statement

The genome re-sequencing data has been deposited to SRA database at NCBI under BioProject PRJNA1043235.

CRediT authorship contribution statement

Yiqing Chen: Writing – original draft, Formal analysis, Data curation. Lina Dong: Resources, Investigation, Funding acquisition. Huiqin Yi: Methodology. Catherine Kidner: Writing – review & editing. Ming Kang: Writing – review & editing, Project administration, Methodology, Funding acquisition, Conceptualization.

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