A total of 89 individuals from 16 populations of Malania oleifera were collected, including five populations from Yunnan and eleven from Guangxi (Figure 1). Genomic DNA was extracted from silica-dried leaf tissue using a modified CTAB protocol (Doyle and Doyle, 1987). DNA concentration and purity were initially measured using a Nanodrop 1000 spectrophotometer (Nanodrop, MA, USA), followed by accurate quantification using a Qubit fluorometer (Thermo Fisher Scientific, CA, USA). DNA integrity and potential RNA contamination were assessed via 0.8% agarose gel electrophoresis.

High-quality DNA samples were digested with two restriction enzymes, EcoRI and MspI, and adapters were ligated to both ends of DNA fragments using T4 DNA ligase. The ligated products were amplified by PCR, and 400–600 bp fragments were recovered using 1% agarose gel electrophoresis and purified with AMPure XP magnetic beads. Purified products were sequenced on the Illumina HiSeq 4000 platform (Illumina, CA, USA) with paired-end reads of 150 bp.

Raw reads were filtered using Fastp (Chen et al., 2018) to remove adapter sequences, abnormal bases, and low-quality reads (reads with ≥10% unidentified bases or ≥50% bases with Phred quality score ≤10). Clean reads were aligned to the M. oleifera reference genome (GCA_029873635.1) using BWA (Li and Durbin, 2009) (parameters: mem -t 4 -k 32 -M). The alignment results were sorted and indexed using SAMtools (Weeks and Luecke, 2017), and the average sequencing depth was calculated with the samtools depth command. SNP calling was performed using bcftools (Weeks and Luecke, 2017), and the following filtering criteria were applied: minimum sequencing depth ≥3×, missing genotype rate ≤30%, and minor allele frequency (MAF) ≥0.05. The remaining high-quality SNPs were used for downstream population genetic analyses.

Genetic diversity indices, including nucleotide diversity (π), observed heterozygosity (Ho), expected heterozygosity (He), Tajima’s D, pairwise genetic differentiation (F_ST), and gene flow (Nm), were calculated using VCFtools (Danecek et al., 2011).

Input files for ADMIXTURE were generated using PLINK (Purcell et al., 2007). Population structure and ancestral components were inferred using ADMIXTURE (Harris and DeGiorgio, 2017) with predefined clusters (K = 1–10). The optimal K value was determined based on the minimum cross-validation (CV) error. The results were visualized using the R package ggplot2 (https://github.com/tidyverse/ggplot2). Principal component analysis (PCA) was also conducted using PLINK (Purcell et al., 2007), and the results were visualized in three dimensions using the R package plotly (https://github.com/plotly/plotly.R). Pairwise genetic distances among individuals were calculated in PLINK (Purcell et al., 2007), and a neighbor-joining (NJ) tree was constructed using the R package ape (https://github.com/emmanuelparadis/ape). The NJ tree was visualized and annotated using iTOL (Letunic and Bork, 2024).

Gene flow among populations was inferred using Treemix (Pickrell and Pritchard, 2012). Migration parameters (m) were set from 1 to 10, and each m value was repeated ten times. The optimal number of migration edges was determined using the R package optM (Fitak, 2021). Visualization of migration events was performed with the plotting_funcs.R script provided in Treemix (Pickrell and Pritchard, 2012).

Pairwise genetic distances were estimated based on F_ST/(1−F_ST) using VCFtools (Danecek et al., 2011). To quantify isolation by environment (IBE), environmental distance was calculated using a subset of key variables. Based on the major environmental drivers identified by Zhang et al. (2025), we performed a Variance Inflation Factor (VIF) analysis to minimize multicollinearity. The retained variables (Bio07, Bio15, Bio16, srad, and t_usada_tex) were then used for subsequent analysis. Geographic distance, representing isolation by distance (IBD), was computed from population coordinates using the R package geosphere. Partial Mantel tests were performed using the R package vegan (https://github.com/vegandevs/vegan) to disentangle the effects of geography and environment. Specifically, IBD was tested while controlling for environmental distance, and IBE was tested while controlling for geographic distance. The significance of correlations was assessed using 999 permutations.

Redundancy analysis (RDA) was conducted using the R package vegan to evaluate the influence of environmental variables on genetic variation among SNPs. SNPs significantly associated with environmental factors (p < 0.05) were mapped to the M. oleifera reference genome to identify candidate genes.

Historical changes in effective population size (Ne) were inferred using Stairway Plot 2 (Liu and Fu, 2020). The site frequency spectrum (SFS) was generated from SNP data using the Python script easySFS (https://github.com/isaacovercast/easySFS), with the folded SFS configuration to account for minor allele frequencies. The generation time was set to 10 years, and the neutral mutation rate to 2.5 × 10^-8 per site per generation (Shen et al., 2024). Projection values maximizing the number of segregating sites were used to output SFS data, which were then provided as input for Stairway Plot analyses. Sixty-seven percent of the sites were randomly selected, and 200 bootstrap replicates were performed to estimate the median Ne and its 95% pseudo-confidence intervals (CIs).

A total of 16 populations of Malania oleifera were genotyped using the GBS approach. In total, 3,278,540 raw SNPs were identified, and after stringent quality filtering and validation, 332,551 high-quality single nucleotide polymorphisms (SNPs) were retained for downstream analyses (Figure 2). The raw reads generated in this study have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under BioProject accession PRJNA1418834 with SRA run accessions SRR37114155-SRR37114243.

Genetic diversity indices, including nucleotide diversity (π), observed heterozygosity (Ho), expected heterozygosity (He), and Tajima’s D, were calculated for each of the 16 populations (Table 1). At the provincial level, populations from Guangxi exhibited significantly higher nucleotide diversity compared to those from Yunnan. Conversely, no significant differences were observed between the two regions for mean Ho, He, and Tajima’s D.

Nucleotide diversity (π) at the population level ranged from 0.109 in MLS to 0.264 in BM, indicating that the BM population possessed the highest genetic variability, while MLS exhibited the lowest. Observed heterozygosity (Ho) ranged from 0.333 (SG) to 0.458 (LG), whereas expected heterozygosity (He) ranged from 0.304 (SG) to 0.506 (MLS). Most populations showed similar Ho and He values; however, in XL, BB, MLS, and LG, Ho was markedly lower than He, suggesting a deficiency of heterozygotes. This pattern was most pronounced in MLS (Ho = 0.434, He = 0.506) and XL (Ho = 0.431, He = 0.484), implying a potential inbreeding effect within these populations.

Tajima’s D values were positive in all 15 populations with sufficient data (ranging from 0.317 in BB to 0.670 in NL), suggesting a deficit of rare alleles and a predominance of intermediate-frequency alleles. Because the JX population included only a single sample (N = 1), Ho, He, and Tajima’s D could not be estimated.

Population structure was assessed using ADMIXTURE, principal component analysis (PCA), and neighbor-joining (NJ) phylogeny based on SNP data. ADMIXTURE cross-validation analysis indicated that K = 6 was the optimal number of clusters with the lowest CV error (Supplementary Figure 1). Based on this model, the 16 populations were divided into six major genetic clusters: Clade 1: JX, BM, and LA; Clade 2: LY and FS; Clade 3: BS and TL; Clade 4: NL; Clade 5: GL, AD, and XL; Clade 6: MLS, LG, SG, BB, and FN (Figure 3A). PCA results further supported this structure (Figure 3B). The first three principal components clearly separated individuals from different geographical regions, consistent with the ADMIXTURE clustering. Similarly, the NJ tree (Figure 3C) grouped individuals into six branches, largely corresponding to their geographic origins. Collectively, these analyses indicate that M. oleifera exhibits strong geographic structuring and significant population differentiation.

Pairwise genetic differentiation (F_ST) and gene flow (Nm) analyses demonstrated widespread yet variable levels of differentiation among populations (Figure 4). F_ST values ranged from 0.070 (JX–BM) to 0.495 (JX–MLS). According to Wright’s guidelines (Wright, 1978), values above 0.25 indicate high differentiation; most population pairs exceeded this threshold, reflecting pronounced population structure.

Nm values ranged from 0.255 (JX–MLS) to 3.426 (FS–LY). The majority of population pairs showed Nm < 1.0, suggesting restricted gene flow and a strong effect of genetic drift. The JX–MLS pair exhibited the highest F_ST (0.495) and lowest Nm (0.255), indicating strong geographic or reproductive isolation. In contrast, the JX–BM pair displayed the lowest F_ST (0.070) and relatively high Nm (3.321), reflecting frequent genetic exchange. Several other pairs also exhibited elevated gene flow (Nm > 2.0), including FS–LY (3.426), JX–GL (3.083), LG–SG (2.762), JX–LA (2.177), and GL–AD (2.131).

The optM analysis indicates that m = 4 is the optimal number of migration edges (Supplementary Figure 2). Treemix analysis with two migration events (m = 4) revealed a stepwise divergence pattern along the drift axis (Figure 5). Four directional migration events were detected: from FN to JX, from FN to Clade 1/Clade 2/Clade 3, from LA to BM, and from LA to LY/FS/BM, mainly among geographically adjacent populations. These findings suggest that localized gene flow still occurs along geographic contact zones despite overall strong differentiation.

To identify the main drivers of genetic differentiation, Mantel tests were conducted to evaluate the relationships between genetic differentiation (F_ST), geographic distance (IBD), and environmental distance (IBE). The results revealed a highly significant positive correlation between genetic and geographic distances (p < 0.001; Figure 6A), supporting a clear isolation-by-distance (IBD) pattern. This indicates that gene flow is strongly constrained by geographic barriers, leading to greater differentiation among distant populations. In contrast, the correlation between genetic and environmental distances was not significant (p = 0.19; Figure 6B), suggesting that environmental heterogeneity had a limited role in shaping population structure at this spatial scale. Overall, geographic isolation rather than environmental adaptation appears to be the primary driver of population differentiation in M. oleifera.

RDA identified 4,361 SNPs significantly associated with environmental variables, of which 738 (16.9%) were successfully annotated to genes or nearby regulatory regions, while 3,623 (83.1%) were located in noncoding or unannotated regions (Figure 7). The main environmental variables included precipitation and soil characteristics. The first two RDA axes cumulatively explained 36.24% of the total genetic variation. These results indicate that precipitation and soil factors may contribute to local adaptation in M. oleifera, providing important candidate loci for further Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses and functional validation.

Historical demographic reconstruction using Stairway Plot 2 revealed distinct fluctuations in the effective population size (Ne) of M. oleifera. In the early period, Ne remained large and relatively stable. Approximately 30 kilo years ago (kya), population size sharply declined, which was likely associated with the onset of and climatic fluctuations during the Last Glacial Maximum (LGM, 21 kya). Thereafter, even with the warming climate during the early Holocene (10–5 kya), the downward trend of the population did not reverse, showing only a slight deceleration in the rate of decline. Subsequently, the population continued to undergo a progressive decline (Figure 8).

Based on GBS data, M. oleifera exhibited a moderate level of nucleotide diversity, whereas genome resequencing revealed a relatively higher level of nucleotide polymorphism (Shen et al., 2024). This discrepancy primarily arises from the methodological differences between the two approaches: GBS captures SNPs only near restriction enzyme cut sites and does not provide uniform coverage across the entire genome, which may lead to an underestimation of genome-wide genetic diversity (Davey et al., 2011; Elshire et al., 2011). Nevertheless, both methods yielded consistent assessments of relative diversity patterns among populations. Regional comparisons indicated that the Yunnan populations exhibited lower genetic diversity than those from Guangxi, while the BM population showed the highest level of genetic diversity, suggesting that it may have retained more ancestral variation or maintained a larger effective population size. In MLS and XL populations, the Ho was lower than He, suggesting heterozygote deficiency and potential inbreeding (Keller, 2002). Such patterns are common in small, geographically isolated plant populations and are usually driven by genetic drift and limited gene flow (Ellstrand and Elam, 1993). As a root hemiparasitic tree species dependent on specific hosts and microhabitats, M. oleifera may have intensified inbreeding due to habitat fragmentation and ecological specialization (Li et al., 2019; Wang et al., 2025).

Population structure analyses consistently revealed M. oleifera into six major genetic clusters, which correspond closely with geographic distribution, indicating that geographic isolation plays a dominant role in shaping the genetic structure. The unique geomorphological features of karst regions, including deep valleys, limestone peaks, and narrow gorges, may act as natural barriers that restrict long-distance pollen and seed dispersal, thereby promoting inter-population genetic divergence (Shen et al., 2019). The Yujiang River and its tributary network form a strong dispersal barrier within the karst context, limiting across-river gene flow and accelerating genetic differentiation among river-separated populations (Wu et al., 2024). The river barrier effect has been documented in other mountain and karst species, such as Parrotia subaequalis (Hung T. Chang) R.M. Hao & H.T. Wei (Geng et al., 2015), Primulina tabacum Hance (Ni et al., 2006) and Quercus fabri Hance (Xiong et al., 2020), where inter-river populations exhibit significantly restricted gene flow. Therefore, the genetic structure of M. oleifera is shaped not only by mountain isolation but also by the configuration of regional river systems.

High pairwise F_ST values indicate strong genetic differentiation, suggesting that gene flow among populations has been largely absent over time. The JX–BM and FS–LY population pairs showed relatively low F_ST and higher Nm values, implying localized connectivity. Migration events inferred by Treemix also support the occurrence of historical gene flow among geographically adjacent populations, suggesting that local connectivity may have partially mitigated complete genetic isolation (Perez-Garcia et al., 2025).

A significant positive correlation between genetic and geographic distances, but no correlation with environmental distance, suggests that the genetic structure of M. oleifera is primarily driven by IBD rather than IBE. This pattern aligns with findings in Fagus crenata Blume (Hiraoka and Tomaru, 2009) and Euryodendron excelsum Hung T. Chang (Shen et al., 2009; Su et al., 2009). M. oleifera relies partly on animal-mediated seed dispersal, but due to gravity-dependent seed fall and limited pollination range, gene flow mainly occurs over short distances, resulting in a pronounced spatial genetic gradient (Wu et al., 2004; Huang et al., 2008). Furthermore, geographic barriers such as valleys and rivers intensify isolation effect at the regional scale. Although environmental factors may influence genetic variation locally, they do not appear to dominate the pattern of population differentiation at the spatial scale examined in this study.

RDA analysis identified 4,361 SNPs significantly associated with environmental variables, among which 738 were mapped to known genes or regulatory regions. These candidate loci were primarily related to precipitation patterns and soil properties, suggesting potential signals of local adaptation in M. oleifera. Similar environment-driven genetic variation has been reported in other plant species, such as moisture-related adaptive genes in Cunninghamia lanceolata (Lamb.) Hook. (Gao et al., 2021) and temperature and precipitation associated loci in Ulmus elongata L.K. Fu & C.S. Ding (Lin et al., 2025). Although the overall effect of environmental isolation was not significant, microhabitat variation may still impose selection pressures that shape adaptive differentiation in specific genomic regions, warranting further functional validation.

Demographic inference revealed that M. oleifera once maintained a large and stable Ne but experienced a sharp and persistent contraction beginning around 30 kya. This downward trajectory aligns with the climatic fluctuations leading into the LGM (Bai et al., 2018), Unlike many species that exhibited a robust rebound during the postglacial warming of the early Holocene, the population of M. oleifera only showed a slowed rate of decline, without a clear recovery. The subsequent and more recent decline in Ne may be further exacerbated by anthropogenic disturbances and increasing habitat fragmentation. This pattern of sustained decline despite climatic amelioration has also been documented in other endangered species, such as Begonia masoniana Irmsch. ex Ziesenh. complex (Chen et al., 2024) and Ostrya rehderiana Chun (Yang et al., 2018), highlighting the species’ vulnerability to long-term environmental and human-induced pressures.

From a conservation perspective, M. oleifera should be regarded as a geographically isolated species with small and fragmented populations. Conservation strategies, both in situ and ex situ, should be implemented based on genetic differentiation units, with emphasis on maintaining potential gene flow corridors across rivers and mountain ranges. In addition, conservation and restoration programs should consider local adaptation differences to preserve both overall genetic diversity and adaptive evolutionary potential.