1. Introduction

Relict species—lineages that have persisted as range-restricted remnants of formerly widespread clades—offer rare, valuable windows into macroevolution, historical biogeography, and extinction risk under contemporary environmental change (Grandcolas et al., 2014). They are striking because many retain conservative morphologies reminiscent of their fossil relatives, yet continue to evolve genetically and ecologically (Grandcolas et al., 2014; Shear and Werth, 2014). In contrast to their broad Tertiary distributions across the Northern Hemisphere, most extant Tertiary relicts are now confined to refugial regions in East Asia, North America, and southwestern Eurasia, typically occurring as small, isolated populations shaped by complex geology and Quaternary climate oscillations (Milne and Abbott, 2002; Milne, 2006; Qiu et al., 2011; Tang et al., 2018). Among these regions, East Asia harbors an exceptional concentration of relict trees and “living fossil” genera—such as Metasequoia glyptostroboides, Ginkgo biloba, Davidia involucrata, and Tetracentron sinense—that occur in forests notable for high biodiversity and endemism (Manchester et al., 2009; Lu et al., 2018; Zhu et al., 2020). Understanding how such relicts diversified, persisted, and currently respond to rapid climate change remains a central challenge in evolutionary biology and conservation genomics.

Over the past two decades, phylogeographic studies in China and adjacent regions have collectively shown that repeated glacial–interglacial cycles reshaped ranges and connectivity, yielding genetically structured lineages associated with long-term refugia and dispersal corridors shaped by mountainous topography and monsoon dynamics (Wang et al., 2009; Qiu et al., 2011; Sakaguchi et al., 2012; Sun et al., 2014; Zhang et al., 2016). Complementary ecological niche modeling has further highlighted “long-term stable refugia” that likely buffered relict plants during Pleistocene climate variability, leaving a legacy of divergent intraspecific lineages that persist in parapatry or secondary contact (Tang et al., 2018). However, much of this work has relied on limited molecular markers, which often lack power to resolve genome-wide patterns of diversity, selection and adaptation, and provide limited insight into genetic basis of environmental adaptation and genetic load (Ellegren, 2014; Song et al., 2023; Wang et al., 2023).

The genomics era has transformed what can be inferred about relict trees. Chromosome-level reference genomes and dense resequencing now enable precise estimation of genome-wide diversity, population structure, runs of homozygosity (ROH) and inbreeding, mutational load, and demographic history (Ellegren, 2014; Song et al., 2023; Lin et al., 2025). Whole-genome resequencing of multiple East Asian relict trees has revealed cryptic intraspecific divergences, recurrent hybridization at lineage contact zones, and heterogeneous lineage-specific demography—features with direct implications for delimiting conservation units and assessing vulnerability. In the “living fossil” dove tree Davidia involucrata, population-genomic analysis has uncovered three deeply divergent local lineages (emerging ca. 3.1–0.3 million years ago; Ma) and frequent hybridization among them, highlighting the interplay of geography, climate oscillations, and secondary contact in structuring diversity (Ren et al., 2024). Failure to recognize such cryptic lineages biases predictions of environmental responses and misallocates conservation resources (Ren et al., 2024). Comparable conclusions have been reported for Tetracentron sinense, which shows multiple divergent lineages and hybridization and lineage-specific demography. Genotype–environment association and genomic offset analyses have further delineated the populations most vulnerable under future climates (Jing et al., 2025). Broadly, Pleistocene climate fluctuations have strongly shaped divergence and demography across East Asian relicts, with repeated expansions and contractions aligned with glacial–interglacial cycles and refugial dynamics (Jing et al., 2025). Nevertheless, the evolutionary consequences of Quaternary change are far from uniform across species, and for many relicts, genome-scale assessments of demographic history, genetic diversity, inbreeding, and genomic vulnerability are still lacking—information essential for both robust evolutionary inference and conservation planning.

While most plant population-genomic inferences have historically relied on single-nucleotide polymorphisms (SNPs), mounting evidence indicates that structural variants (SVs)—including insertions, deletions, copy-number variants (CNVs), inversions, and translocations—constitute a pervasive, functionally consequential layer of standing variation that can strongly influence adaptation, demography, and genetic load (Feulner et al., 2013; Collins et al., 2020; Sang et al., 2022; Stuart et al., 2024; Zhang et al., 2024; Lin et al., 2025; Long et al., 2025). SVs can rewire gene dosage and regulation, alter recombination landscapes, and generate large-effect alleles underlying ecological differentiation and reproductive barriers (Wellenreuther et al., 2019). Empirical evidence increasingly links SVs to local adaptation and multiple components of vulnerability. For example, CNVs can surpass SNPs in detecting environment–genotype associations in marine systems and outperform them in capturing temperature-associated genomic variation (Dorant et al., 2020). In plants, SVs have been shown to both constrain and facilitate adaptation in natural populations, consistent with their dual capacity to generate beneficial novelty while also perturbing genome integrity (Hämälä et al., 2021). Pan-genome studies now underscore the role of SVs as major drivers of gene-content diversity and environmental adaptation across plant lineages, strengthening the case for routine SV interrogation in conservation genomics (Liang et al., 2025). Critically, SVs contribute to realized genetic load and fitness variation: in the inbred Scandinavian wolf population, structural genomic variation measurably contributed to realized genetic load, underscoring the need to quantify SV burden alongside deleterious SNPs when assessing population viability (Smeds et al., 2024). Because SNPs and SVs capture complementary facets of genomic variation, integrating both data types is critical for a comprehensive view of diversity, demography, adaptation, and vulnerability in relict trees.

Perkinsiodendron macgregorii (Styracaceae) is a monotypic East Asian relict that exemplifies these needs. The species was originally described within Halesia, a classic East Asia–North America disjunct genus comprising two species in the eastern United States (H. carolina, H. diptera) and a single species in China (H. macgregorii) (Fritsch et al., 2016). However, combined morphological evidence and multilocus phylogenies show that the East Asian taxon formerly known as H. macgregorii is sister to Rehderodendron, whereas the two North American Halesia species cluster with Pterostyrax hispidus. Accordingly, Fritsch et al. (2016) erected the monotypic genus Perkinsiodendron to reflect these evolutionary relationships, and subsequent family-level treatments and plastid phylogenies in Styracaceae clarified generic boundaries and corroborated this revision (Fritsch, 2001; Yan et al., 2018). Perkinsiodendron macgregorii occurs in montane forests and forest margins at elevations of ~700–1200 m in southeastern China (Chen and Chen, 1996; Fritsch et al., 2016) and is valued horticulturally for its straight trunk, abundant and fragrant white blossoms, uniquely shaped fruits, and striking red autumn foliage (Fig. 1a). Despite its apparently broad regional distribution, populations are typically small and narrowly localized, with limited regeneration linked to defective embryo development, low germination rates, and recalcitrant seeds (Liao et al., 2024). Although prior work has addressed seed biology and germination (Liao et al., 2024), phylogenetic placement within Styracaceae (Fritsch, 2001; Fritsch et al., 2016; Yan et al., 2018), and ecological niche modeling (Yan et al., 2022), genomic resources for P. macgregorii remain limited. No genome-wide studies have examined the genetic diversity, demographic history, or genomic vulnerability of P. macgregorii—information essential for evaluating genetic health and predicting adaptive potential under future climates.

Here, we combine a chromosome-level reference genome with population resequencing to address these questions across the species' range. We generate a high-quality genome assembly for Perkinsiodendron macgregorii and analyze whole-genome resequencing data from 167 individuals representing 30 populations spanning its distribution. Specifically, we aim to: (ⅰ) produce and annotate a reference genome suitable for joint SNP and SV discovery; (ⅱ) resolve lineage divergence and demographic history to contextualize Perkinsiodendron within the context of East Asian relict diversification and potential secondary contact among lineages; (ⅲ) quantify genome-wide diversity, inbreeding, and mutational load at both SNPs and SVs to assess the current genetic health at lineage and population scales; and (ⅳ) evaluate SNP- and SV-based signals of local adaptation and project genomic vulnerability under future climates. Our results provide an integrative population-genomic baseline for P. macgregorii and demonstrate that integrating SNPs and SVs sharpens evolutionary inference and improves conservation assessments in East Asia's relict flora.

2. Materials and methods 2.1. Genome assembly and quality assessment

Fresh young leaves were collected from a wild Perkinsiodendron macgregorii individual in Shaoguan City, Guangdong Province, China (113°2′3″E, 24°52′54″N). For PacBio long-read sequencing, HiFi libraries were sequenced on the PacBio Sequel Ⅱ platform (Pacific Biosciences, Menlo Park, CA, USA). Short-read libraries (350-bp insert) were constructed and sequenced on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA), generating 150 bp paired-end reads. To achieve chromosome-scale assembly, a Hi-C library was prepared using the DpnII restriction enzyme and sequenced on the Illumina NovaSeq 6000 platform. To facilitate gene annotation, total RNA was extracted from four tissues (flower, branch, leaf, and bud), and subsequent transcriptome libraries were sequenced on the Illumina NovaSeq 6000 platform.

Genome characteristics (size, heterozygosity, and repeat content) were surveyed using GCE v.1.0.2 (Liu et al., 2013). HiFi reads were assembled using hifiasm v.0.20.0 (Cheng et al., 2021), and redundant sequences were removed with purge_dups v.1.2.5 (Guan et al., 2020). Chromosome-scale scaffolding used the HapHiC v.1.0.6 (Zeng et al., 2024) pipeline, and scaffolds were then adjusted for scaffolding errors in Juicebox v.2.15 (Durand et al., 2016). To evaluate the quality of the genome assembly, we estimated the Benchmarking Universal Single-Copy Orthologs (BUSCO) score using BUSCO v.5.8.0 (Manni et al., 2021) with the embryophyta_odb10 database, and calculated the Illumina short-read mapping rate and coverage using BWA v.0.7.18 (Li and Durbin, 2009) and SAMtools v.1.20 (Li et al., 2009). Additionally, to assess the quality of repetitive genome regions, LTR_FINDER v.1.0.7 (Xu and Wang, 2007) and GenomeTools v.1.6.5 (Gremme et al., 2013) were used to predict LTR-RTs, and LTR_retriever v.3.0.1 (Ou and Jiang, 2018) was then used to calculate LTR Assembly Index (LAI).

2.2. Genome annotation

Prior to gene prediction, repetitive elements were identified using EDTA v.2.2.2 (Ou et al., 2019). Protein-coding genes were predicted utilizing BRAKER3 (Gabriel et al., 2024) by integrating homolog-based, ab initio, and RNA-seq-assisted evidence. Input data included protein sequences from five plant species (Actinidia chinensis, Arabidopsis thaliana, Nyssa sinensis, Rhododendron simsii, Sinojackia xylocarpa) and RNA-seq reads generated herein. Predicted gene models encoding proteins <100 amino acids were removed. The longest transcript per gene was retained, and annotation completeness was assessed using BUSCO with the embryophyta_odb10 database. Gene functions were inferred using eggNOG-mapper, Pannzer2, and SwissProt/TrEMBL.

2.3. Comparative genomic analyses

To infer whole-genome duplication (WGD) events, we employed WGDI (Sun et al., 2022) to calculate the distribution of synonymous substitutions per synonymous site (Ks). Genome synteny between Perkinsiodendron macgregorii and Sinojackia xylocarpa was identified using MCScanX (Wang et al., 2012) and visualized using JCVI (Tang et al., 2024). OrthoFinder v.3.1.0 (Emms and Kelly, 2019) was applied to identify single-copy and low-copy orthologous (≤ 2 copies) orthogroups across P. macgregorii and 14 other angiosperms. Multiple sequence alignments of low-copy orthologous proteins were generated with MAFFT v.7.525 (Katoh and Standley, 2013), and poorly aligned regions were filtered out using trimAL v.1.5.0 (Capella-Gutiérrez et al., 2009). IQ-Tree v.3.0.1 (Wong et al., 2025) was used to construct a maximum-likelihood (ML) phylogenetic tree. Divergence times were estimated using MCMCTree in PAML v.4.10.9 (Yang, 2007). Three calibration points were chosen from the TimeTree website (http://www.timetree.org/), including Oryza sativa–Tetracentron sinense (142.1–163.5 million years ago; Ma), V. vinifera–Dipteronia sinensis (109.8–124.4 Ma), and S. xylocarpa–Rhododendron simsii (82.8–106.0 Ma). Gene family expansion and contraction were analyzed using CAFE v.5.1.0 (Mendes et al., 2021).

2.4. Genome resequencing and variant calling

Leaf samples from 167 individuals were collected across 30 populations of Perkinsiodendron macgregorii. Genomic DNA was extracted from silica-dried leaves using the CTAB method (Doyle and Doyle, 1987) and sequenced on the DNBSEQ-T7 platform to generate 150-bp paired-end reads. Raw reads were filtered using fastp v.0.23.4 (Chen et al., 2018) to remove adaptors and low-quality bases, and then the clean reads were mapped to the reference genome with BWA v.0.7.18 and sorted using SAMtools v.1.20. The MarkDuplicate and addOrReplaceReadGroup tools from Picard v.2.18.29 (https://github.com/broadinstitute/picard) were used to mark PCR duplicates and assign read group information, respectively. The resulting BAM files were used for downstream variant calling. SNP calling and joint genotype calling were performed using HaplotypeCaller and CombineGVCFs in the Genome Analysis Toolkit (GATK) v.4.5 (McKenna et al., 2010), respectively. Hard filters were applied, and variants were further refined with VCFtools v.0.1.16 (Danecek et al., 2011) by removing sites with > 20% missing data and mean depth < 5 or > 100, yielding 767,554,148 variant sites (Dataset 2). Subsequently, SNPs with a minor allele frequency (MAF) > 0.05 were retained using VCFtools, resulting in 22,282,414 sites (Dataset 3). Finally, linkage disequilibrium (LD) pruning was performed using PLINK (Chang et al., 2015) (r² > 0.2), generating 2,563,601 SNPs (Dataset 4).

Structural variants (SVs) were identified using SMOOVE v.0.2.8 (https://github.com/brentp/smoove), DELLY v.1.3.2 (Rausch et al., 2012), and MANTA v.1.6.0 (Chen et al., 2016). For DELLY and SMOOVE, SVs were first called independently per individual and then merged into one a cross-sample, nonredundant SV set, which was subsequently used for joint genotyping across all samples. BCF files generated by DELLY were converted to VCF format using bcftools. Due to computational constraints, MANTA was run per individual using the default workflow (configManta.py; runWorkflow.py), and inversions were additionally extracted from breakend calls using convertInversion.py. Following Jeffares et al. (2017), calls from the three methods were then merged with SURVIVOR v.1.0.7 (parameters: 1000 2 1 1 0 50), harmonizing 0/0, 0/1, and 1/1 classes; retained variants were then merged into a joint sample file using SURVIVOR with the same parameters. To ensure high-quality SVs, we filtered sites with quality scores ≤ 30 and a missing rate >20% using VCFtools. After filtering, a total of 10,357 SVs were retained. Then, VCFtools was used to filter SVs with MAF < 0.02, generating 8549 SVs. LD pruning was performed using PLINK v.1.9 (r² > 0.2), resulting in 5575 SVs for downstream analyses.

2.5. Population structure and genetic diversity

After LD pruning, we excluded the core adaptive variants (i.e., SNPs and SVs identified by genotype–environment association analyses) from the pruned SNP and SV VCFs using VCFtools. The resulting “neutral” SNP and SV datasets were then analyzed separately to infer population structure within the species. We investigated population structure using ADMIXTURE v.1.3.0 (Alexander and Lange, 2011) with the number of genetic clusters (K) ranging from 1 to 10. Principal component analysis (PCA) was performed on all individuals using GCTA v.1.26.0 (Yang et al., 2011). Subsequently, VCF2Dis (Xu et al., 2025) was used to calculate pairwise genetic distances based on Dataset 3, and a neighbor-joining (NJ) tree was then constructed using FastME (Lefort et al., 2015). Next, the identity-by-state (IBS) kinship matrix between individuals was analyzed using VCF2PCACluster v.1.41 (He et al., 2024) based on Dataset 2. Additionally, the D-statistic was estimated using Dsuite v.0.5 r58 (Malinsky et al., 2021) to perform the ABBA–BABA test among populations.

To evaluate genetic diversity, we calculated unbiased nucleotide diversity (π), absolute diversity (D_XY), and fixation statistics (F_ST) with PIXY v.1.2.11 (Korunes and Samuk, 2021) using Dataset 2 in a sliding window of 200 kb. Whole-genome heterozygosity (H_e) was estimated using realSFS implemented in ANGSD v.1.9 (Korneliussen et al., 2014) for each individual. Linkage disequilibrium (LD) decay was analyzed using PopLDdecay v.3.43 (Zhang et al., 2019).

2.6. Demographic history

We used three complementary approaches to reconstruct the demographic history of Perkinsiodendron macgregorii. First, the pairwise sequentially Markovian coalescent (PSMC) model (Li and Durbin, 2011) was applied to infer effective population size (N_e) trajectories from genome-wide heterozygosity, using the highest-coverage individuals from representative populations of each lineage. Due to no direct data are available on age at first reproduction or generation time for P. macgregorii, generation time was set to 10 years and the mutation rate to 2.5 × 10⁻⁸ per site per generation, following estimates of Sinojackia xylocarpa (Zhu et al., 2024). To capture more recent dynamics, we additionally inferred N_e fluctuations using Stairway Plot v.2.1.2 (Liu and Fu, 2020) and SMC++ (Terhorst et al., 2017), applying the same generation time and mutation rate.

2.7. Inbreeding and mutation load

Runs of homozygosity (ROHs) were identified using PLINK. ROHs were defined as segments > 10 kb and classified as short (10 < ROH ≤ 100 kb), medium (100 < ROH ≤ 200 kb), and long (> 200 kb). For each individual, the inbreeding coefficient based on ROHs (F_ROH; Kardos et al., 2015) was calculated as the proportion of the genome covered by ROHs > 10 kb.

To estimate genetic load, SnpEff v.5.2e (Cingolani et al., 2012) was used to annotate SNPs and SVs for each individual. An individual of Sinojackia xylocarpa was used as the outgroup to polarize alleles (ancestral vs. derived). For SNPs, variants annotated as synonymous were defined as synonymous (SYN) mutations. Nonsynonymous mutations were categorized by Grantham score (GS; Grantham, 1974) as tolerated (TOL; 5 < GS < 150) or deleterious (DEL; GS ≥ 150). Putative loss-of-function (LOF) mutations were identified with SnpSift using the criterion “LOF[*].PERC = 1.00”. For SVs, SnpEff impact categories were used to classify variants by predicted effect, and we summarized per-individual proportions across HIGH, MODERATE, and LOW impact classes. For all SNP classes (TOL, DEL, LOF), we tabulated per-individual proportions of derived homozygotes and heterozygotes. We refer to the proportion of derived homozygotes as “realized load” and to the proportion of derived heterozygotes as “masked load”.

2.8. Genomic vulnerability analysis

We identified climate-associated variants using two complementary methods applied to SV and LD-filtered SNP datasets, with SVs phased and imputed using Beagle v.5.4 (Browning and Browning, 2009). Nineteen bioclimatic variables were downloaded from WorldClim 2.1 at 30-arc-second (~1 km) resolution (Fick and Hijmans, 2017). First, latent factor mixed models (LFMMs) were run in the R package LEA (Frichot and François, 2015), with the number of latent factors set to the K inferred with ADMIXTURE; associations were retained at a Benjamini–Hochberg FDR ≤ 0.05. Second, redundancy analysis (RDA) was used to test multilocus genotype–environment associations with the R package vegan (Dixon, 2003). To reduce collinearity, predictors were ranked with gradientForest (Ellis et al., 2012), and five variables with pairwise correlation |r| < 0.6 were retained (Table S1). RDA outliers were defined as variants with loading beyond the mean ± 3 SD. Variants detected by both LFMM and RDA were retained as “core adaptive variants” for genomic vulnerability analysis. Core adaptive variants were annotated using SnpEff, and gene ontology (GO) enrichment was performed with the R package clusterProfile (Wu et al., 2021).

To partition the roles of geography and environment in shaping genetic variation of core adaptive and neutral variants, we performed Mantel and partial Mantel tests to assess isolation by distance (IBD) and isolation by environment (IBE), with significance determined from 999 permutations in the R package vegan. For partial Mantel tests, we tested IBE by correlating genetic distance with environmental distance while controlling for geographic distance, and tested IBD by correlating genetic distance with geographic distance while controlling for environmental distance.

The potential distribution range of Perkinsiodendron macgregorii was inferred using ecological niche modeling (ENM) with MaxEnt v.3.4.4 (Phillips and Dudík, 2008). For genomic offset under future climate change, we considered two Shared Socioeconomic Pathway (SSP) scenarios (SSP126, SSP585) for two periods (2021–2040, 2081–2100) under the BCC-CSM2_MR model. We used the R package gradientForest to estimate genetic offset under the different future scenarios with five environmental variables, and then plotted them in ArcGIS v.10.8 (ESRI, Redlands). Genetic offset was calculated as the Euclidean distance in gradientForest-transformed genomic composition between present and future climates. Risk of non-adaptation (RONA) was computed as a weighted average across core adaptive variants and the five climate variables, following the methods outlined by Sang et al. (2022). Local, forward, and reverse offsets were quantified with gradientForest following Gougherty et al. (2021): local offset is the minimal displacement cost within the current range; forward offset is the minimal future location matching each population's current genomic–environmental composition; reverse offset is the set of present-day locations matching projected future composition.

3. Results 3.1. Genome assembly and comparative genomics

We generated a chromosome-level assembly of the Perkinsiodendron macgregorii genome using 73.54 Gb of Illumina short reads (65 ×), 48.5 Gb of PacBio HiFi reads (43 ×), and 207.14 Gb of Hi-C data (182 ×) (Table S2). A 17-mer analysis of the Illumina reads estimated a genome size of ~1.20 Gb (Table S3 and Fig. S1). The final assembly spans ~1.15 Gb (contig N50 = 7 Mb; scaffold N50 = 86 Mb), with 92.33% of the assembly anchored to 12 pseudochromosomes (Table S4; Figs. 1b and S1). Quality metrics indicate a reference-grade genome assembly (Illumina read-mapping rate = 91.5%; LAI = 17.9; BUSCO completeness = 98.1%; Table S5). Repeats account for 60.65% of the genome, including Gypsy (18.37%) and Copia (3.69%) elements (Table S6). Gene prediction identified 31,851 protein-coding genes, of which 30,767 (96.60%) were functionally annotated; BUSCO recovered 1572 complete genes (97.4%) in the annotation set (Tables S7 and S8), indicating a high level of annotation completeness.

Comparative analyses support conserved collinearity within Styracaceae. Synteny with Sinojackia xylocarpa revealed extensive chromosome-scale collinearity (Figs. 1d and S2). Ks distributions showed a peak at ~1.35 that is shared with Vitis vinifera, consistent with the core eudicot whole-genome triplication (WGT), and a second peak at ~0.47 shared by Perkinsiodendron macgregorii and S. xylocarpa, indicating that there is no lineage-specific WGD in P. macgregorii (Fig. 1e). Orthologous syntenic blocks between P. macgregorii and S. xylocarpa showed an orthologous Ks peak at ~0.04, suggesting relatively recent divergence between the two species (Fig. 1e). Phylogenetic analyses based on 655 low-copy orthologs from 15 angiosperms recovered P. macgregorii and S. xylocarpa as sister lineages that diverged ~10.0 Ma (95% HPD: 17.0–4.1 Ma) (Figs. 1c and S3).

Gene family analysis identified 1305 expansions and 526 contractions in Perkinsiodendron macgregorii (Tables S9 and S10). Expanded families were enriched for epigenetic regulation of gene expression (GO:0040029), post-embryonic plant morphogenesis (GO:0090698), nuclear speck (GO:0016607), cellular hormone metabolic process (GO:0034754). These functional categories plausibly contribute to survival in heterogeneous and stressful environments by promoting regulatory plasticity and developmental flexibility under fluctuating conditions (Liu and Zhong, 2024; Srikant and Drost, 2026). Consistent with this interpretation, several members within expanded families have documented roles in flowering and seasonal regulation (VIL1; Sung et al., 2006), RNA-mediated regulation and immunity (PRP16; Qiao et al., 2015), DNA demethylation and genotoxic-stress responses (ROS1; Qüesta et al., 2013), nutrient recycling via autophagy under limitation (ATG11; Li et al., 2014), and specialized metabolism linked to defense (C71Z1; Krieger et al., 2018). In contrast, contracted families were enriched for anthocyanin and flavonoid biosynthesis (GO:0009718; GO:0009813), nucleotidyltransferase activity (GO:0016779), and auxin polar transport (GO:0009926), which may reflect lineage-specific shifts in secondary metabolism and growth regulation, although their ecological consequences require functional validation.

3.2. Population structure, gene flow and genomic diversity

Whole-genome resequencing of 167 individuals yielded an average mapping rate of 92.6% and a mean coverage of 31 × (Table S11). ADMIXTURE (K = 2) resolved two primary lineages that were consistent with their sampling geography—East and West—with finer substructure within each lineage in both SNP and SV datasets (Figs. 2a, b and S4). This clustering pattern was also corroborated by PCA, a maximum-likelihood phylogeny, and identity-by-state (IBS) analysis (Figs. 3a and S4). Notably, the JXGS and JXYJF populations were strongly differentiated across methods (Figs. 3a and S4). Mantel (r = 0.388, P = 0.001) and partial Mantel (r = 0.359, P = 0.001) tests supported isolation by distance (IBD) across the species’ range (Fig. 2c and d), with significant IBD within the West lineage (r = 0.490, P = 0.012) but not within the East lineage (r = 0.033, P = 0.407) (Fig. S5). D-statistics (Dsuite) indicated inter-lineage introgression between East and West, as well as more frequent intra-lineage gene flow among West populations than among East populations (Fig. S6).

Genome-wide nucleotide diversity (π) for Perkinsiodendron macgregorii was 0.0102, with higher diversity in the East lineage (0.0103) than the West lineage (0.0086). Population-level π ranged from 0.0048 (JXJGS) to 0.0098 (FJYX), and per-individual heterozygosity ranged from 0.0073 (JXJGS) to 0.0132 (FJXY) (Table S12). Pairwise differentiation among populations ranged from F_ST = 0.0284 to 0.5222 and from D_XY = 0.0076 to 0.0113 (Fig. S7), reflecting pronounced East–West divergence (mean F_ST = 0.1171; mean D_XY = 0.0107). Linkage disequilibrium (LD) decay was faster in the East than in the West, suggesting a larger effective population size (N_e) and/or higher recombination in the East lineage (Fig. 3b).

3.3. Demographic history

PSMC analyses for both lineages indicated modest growth during the late Miocene followed by a long-term decline (Fig. 3c). The East lineage exhibited a slight expansion from ~0.1 to 0.07 Ma, followed by a bottleneck at ~0.07 Ma and a subsequent expansion centered at ~0.03 Ma before stabilization (Fig. 3c). By contrast, the West lineage continued declining until ~20 thousand years ago (Ka), followed by an expansion around ~10 Ka, and then stabilized (Fig. 3c). At the onset of the Naynayxungla Glaciation (NG), the East lineage had a larger N_e than the West. SMC++ dated East–West divergence to ~0.16 Ma (Fig. S8). Stairway Plot 2 inferred post-LGM (26.5–19 Ka) stability followed by a continuous decline beginning ~3 Ka, with a more pronounced decline in the West (Fig. 3d).

3.4. Inbreeding and genetic load

Short ROHs were more frequent than medium and long ROHs across the species (Fig. S9), indicating predominantly ancient inbreeding with additional evidence for recent inbreeding. Across lineages, the West lineage harbored significantly more short, medium, and long ROHs than the East lineage (Fig. S9). The genome-wide inbreeding coefficient based on ROHs (F_ROH) ranged from 0.0116 to 0.1937; lineage means were higher in the West (0.054) than in the East (0.030), with JXJGS showing the highest values, followed by HNSHS (Fig. S9). Consistent with the expectation, F_ROH and per-individual heterozygosity were negatively correlated (Fig. S10).

To estimate genetic load, SNPs were partitioned into synonymous (SYN), tolerated (TOL; 5 < GS < 150), deleterious (DEL; GS ≥ 150), and loss-of-function (LOF) classes, and SVs were summarized by predicted impact (LOW/MODERATE/HIGH). Relative to the East, the West lineage showed higher proportions of homozygous TOL and DEL sites and higher realized (homozygous) load overall, despite smaller heterozygous burdens (Fig. 4a). JXJGS (followed by HNSHS) had the highest proportions across mutation categories and exhibited the highest realized load with the lowest masked load, consistent with conversion of masked to realized load under strong inbreeding (Fig. 4b). Both JXJGS and HNSHS had elevated proportions of homozygous LOF mutations within ROHs, implicating recent inbreeding in elevating LOF homozygosity (Fig. 4c). SVs showed a higher proportion of HIGH-impact sites than SNPs (Fig. S11), highlighting their potential to disproportionately affect gene function and fitness.

3.5. Genotype–environment associations and genomic vulnerability

LFMM analysis identified 1884 SNPs and 317 SVs associated with one or more of 19 bioclimatic variables (FDR ≤ 0.05) (Fig. S12 and S13). RDA using five uncorrelated variables (Bio3, Bio8, Bio9, Bio18, Bio19) recovered 9801 SNPs and 170 SVs (Fig. S14). The intersection of the two approaches yielded 755 SNPs and 63 SVs as “core” adaptive variants (Fig. S12 and S13). Most adaptive SNPs and SVs were located in non-coding regions (96.29% and 96.83%, respectively). Among the coding SNPs (3.71%), non-synonymous changes accounted for 35.71% (1.32% of all adaptive variants). Genes associated with core adaptive variants were enriched, for SNP-linked genes (within gene bodies or ±2 kb flanks region), in post-transcriptional regulation of gene expression, response to nematodes, and changes in DNA conformation (Fig. S12), and for SV-linked genes (overlapping gene bodies), in immune response, innate immune response, thylakoid components, and magnesium ion binding (Fig. S13).

Patterns of spatial genetic structure differed between neutral and adaptive markers. Neutral variants (r = 0.388, P = 0.001) exhibited slightly stronger isolation-by-distance (IBD) than adaptive variants (r = 0.191, P = 0.049). For adaptive variants, IBD disappeared after controlling for environmental distance (Fig. 2c and d). Within the West lineage, adaptive variants (r = 0.618, P = 0.001) showed stronger IBD than neutral variants (r = 0.49, P = 0.012), while neither IBD nor isolation-by-environment (IBE) was significant in the East lineage, likely reflecting the broader geographic range and more complex topography of the West lineage (Fig. S5).

Gradient Forest projected increasing genetic offset under higher-emission scenarios and later time windows, with high-offset areas expanding markedly (Fig. 5a–d). High-risk regions were concentrated in western Jiangxi, western Fujian, and southwestern Hunan, with broadly concordant patterns between SNP and SV datasets (Figs. 5 and S15). Risk of non-adaptation (RONA) likewise increased for most environmental variables under more severe climate change scenarios (Fig. S16 and S17). Using two representative variables (Bio9 and Bio19), we found that populations in regions with greater environmental change were expected to show higher RONA values (Fig. S18). A decrease in precipitation during the coldest quarter increased RONA in the GDDDS population. Forward-offset analyses highlighted western Jiangxi, the Hunan-Guangxi border region, and the Fujian-Zhejiang border region, mirroring patterns of local genetic offset (Fig. 5e). Reverse-offset analyses indicated that southwestern portions of the current suitable area would likely no longer support population survival under the SSP585 scenario in 2081–2100 (Fig. 5f). Consistent results were obtained from both SNP and SV datasets (Figs. 5 and S15).

4. Discussion 4.1. Lineage divergence and demographic history

Mountain regions are widely recognized as crucial refugia and centers of endemism across central and southern China, where long-term environmental stability can promote persistence, divergence, and speciation (López-Pujol et al., 2011; Ren et al., 2025). In Perkinsiodendron macgregorii, we detected a strong pattern of isolation by distance (IBD), indicating that geographic continuity contributes to broad-scale population structure. Superimposed on this clinal signal is a discrete East–West split that aligns with the Wuyi Mountains. Similar East–West phylogeographic breaks have been observed in the Opsariichthys acutipinnis-evolans complex (Gao et al., 2023) and in Deinagkistrodon acutus (Huang et al., 2007), supporting a general role of the Wuyi Mountains as a barrier to dispersal. The Wuyi Mountains (~27°–28°N and 117°–118°E) extend northeast–southwest across southeastern China and have likely impeded gene flow historically, coincident with the East–West lineage split observed in P. macgregorii. Previous studies similarly identified the Wuyi Mountains as both a glacial refugium and a biogeographic barrier for many organisms (Huang et al., 2007; Guo et al., 2025; Ren et al., 2025). At the same time, alternative processes could reinforce this divergence. First, the Wuyi-associated boundary may also track contemporary environmental contrasts (e.g., monsoon-related precipitation and seasonality gradients), raising the possibility of isolation by environment (IBE), whereby environmentally mediated selection reduces effective gene flow even across moderate geographic distances. Second, ongoing habitat loss and fragmentation can reduce landscape connectivity, erode genetic diversity, and amplify population differentiation (Johnson et al., 2017; Exposito-Alonso et al., 2022); thus, recent anthropogenic disturbance—especially in low-elevation corridors between mountain blocks—may have further weakened contemporary dispersal and strengthened the apparent East–West split. Together, these patterns suggest that the divergence across the Wuyi Mountains likely reflects the combined effects of historical topographic barriers, environmentally structured selection, and human-driven fragmentation, rather than a single driver acting in isolation.

By contrast, the Nanling Mountains connect eastward to the Wuyi Mountains and northward to the Luoxiao Mountains, and together they function as major east–west (Nanling) and north–south (Luoxiao) dispersal corridors (Qiu et al., 2011; Zhang et al., 2015; Tian et al., 2018; Yang et al., 2019; Li et al., 2023; Lin et al., 2024). These corridors are consistent with the admixture we observed between western Jiangxi and Guangdong populations of the West lineage and populations in the East lineage. Notably, the predominantly N–S orientation of the Luoxiao Mountains may simultaneously act as a partial barrier between the East and West lineages—an effect also inferred in Oreocharis auricula (Sun et al., 2024). Chloroplast data further suggest that both the Wuyi and Nanling Mountains harbored glacial refugia for P. macgregorii (Yan et al., 2025), a view supported by elevated genetic diversity within Wuyi populations. Together, our results suggest that the three mountain ranges (Wuyi, Nanling, and Luoxiao) have jointly shaped the phylogeographic pattern of P. macgregorii, highlighting the profound influence of mountainous terrain on lineage divergence.

Demographic history analyses revealed a long-term decline in effective population size (N_e) since the Late Pliocene, with East–West divergence dated to ~0.16 Ma. This split likely reflects intensified climatic oscillations during the Middle Pleistocene (Pisias and Moore, 1981; Hofreiter and Stewart, 2009). Notably, this divergence coincides with pronounced climatic instability and monsoon variability in East Asia, which could have increased habitat fragmentation and promoted lineage isolation in this relict tree. Since the onset of the Naynayxungla Glaciation (NG), the West lineage has experienced more severe population contraction and weaker postglacial expansion than the East lineage. Consistent with this, the results of Stairway Plot 2 showed a more pronounced post-LGM (26.5–19 Ka) decline in the West lineage. Continued post-LGM contraction may further reflect limited natural regeneration of P. macgregorii and ongoing anthropogenic disturbance (Cao et al., 2015). Although PSMC suggests the East lineage maintained larger N_e around the onset of the NG, SMC++ infers a more recent divergence (~0.16 Ma). We interpret this difference as primarily methodological rather than truly contradictory: PSMC is a single-genome, within-lineage coalescent-based method with limited resolution for recent events and can inflate ancient N_e in recently diverged lineages by treating between-lineage coalescence as within-population events (Mazet et al., 2016). By contrast, SMC++ leverages multiple genomes and site-frequency information to better resolve more recent demographic changes and split times (Terhorst et al., 2017). In addition, absolute timing in both frameworks is sensitive to assumed mutation rates and generation times, so divergence estimates should be viewed as approximate; however, the qualitative inference of long-term decline and stronger recent contraction in the West lineage is consistent across methods.

4.2. Complementary insights from sequence and structural variation

Structural variants (SVs) represent a major—and often underappreciated—fraction of genomic variation and can exert disproportionate functional effects relative to SNPs (Wellenreuther et al., 2019; Stuart et al., 2024). Analyses of SNP and SV datasets yielded congruent population structure, jointly supporting a robust East–West split. In the SV dataset, one individual from HNCZ displayed a unique genetic signature, potentially reflecting localized rearrangements or calling artifacts and thereby warranting targeted validation. Compared with SNPs, SVs showed a significantly higher proportion of high-impact sites (HIGH) and fewer moderate- and low-impact sites (MODERATE and LOW), consistent with their larger mutational targets (deletions, duplications, inversions) and higher potential deleteriousness (Catanach et al., 2019; Collins et al., 2020; Hämälä et al., 2021).

Genotype–environment association results were complementary across marker classes: 755 candidate SNPs and 63 candidate SVs associated with bioclimatic variation mapped to largely distinct functional processes, implying partially distinct molecular routes to environmental response. SNP-linked genes were enriched for post-transcriptional regulation, DNA conformation, and defense processes, whereas SV-linked genes were enriched for immune functions, thylakoid components, and magnesium ion binding. Similar SNP–SV complementarity has been reported in Atlantic salmon (Lecomte et al., 2024). One plausible implication is that SNPs may often contribute to fine-scale tuning of gene expression and RNA metabolism—consistent with enrichment for post-transcriptional regulation and chromatin/DNA-conformation terms—thereby modulating developmental and stress-response programs in a more incremental manner. In contrast, because SVs can encompass larger genomic segments and alter gene dosage, exon structure, or local regulatory architecture, they may disproportionately affect stress-response and organellar/metabolic functions that are central to abiotic tolerance. In Perkinsiodendron macgregorii, the enrichment of SV-linked signals in immune pathways and chloroplast/thylakoid-related functions (and magnesium ion binding) is consistent with the idea that structural changes could influence photosynthetic performance, redox homeostasis, and stress-associated signaling, traits that are often sensitive to drought and temperature extremes and therefore relevant to the persistence and projected vulnerability of relict lineages. Representative SV-linked candidates include stress- and defense-related regulators such as ERF4 and ABA-associated signaling components (e.g., MKK3), along with immune-associated genes (e.g., NDR1, PAT1, CDC5) and jasmonate-related pathways (e.g., LPA1), supporting a potential role for SVs in coordinating broad stress-defense networks. Importantly, the widespread occurrence of non-coding adaptive variants highlights the importance of regulatory elements in shaping gene expression responses to environmental stress. By influencing key regulatory networks, they may enable phenotypic plasticity and rapid adjustment to heterogeneous or changing habitats (Romero et al., 2012; Marand et al., 2023). Overall, SNPs and SVs may contribute to adaptation via distinct molecular mechanisms, with SNPs potentially fine-tuning gene expression and regulatory processes, while SVs may exert larger-scale effects on metabolic and immune-related pathways, collectively enhancing the species’ capacity to cope with heterogeneous environments.

Despite these mechanistic differences, projections of genomic vulnerability (genetic offset and RONA) exhibited the same qualitative pattern for both datasets: maladaptation risk increases under more severe emission scenarios and over longer time horizons. This trend has also been observed in many other species (Sang et al., 2022; Jiang et al., 2025; Lin et al., 2025). However, high-offset regions were more extensive for SNPs than for SVs in Perkinsiodendron macgregorii, suggesting that SNP-based signals may be more spatially diffuse. Notably, some taxa show tight SNP–SV concordance (e.g., Oncorhynchus nerka; Tigano et al., 2024), suggesting that marker behavior can be system-specific.

4.3. Genetic load, diversity and vulnerability

As a relict tree species, Perkinsiodendron macgregorii has comparatively high genome-wide diversity (mean π = 0.0102; Fig. S19), exceeding values reported for several East Asian relicts—including Sinojackia xylocarpa (π = 0.0050; S. Zhu et al., 2024), the Chinese lineage of Cercidiphyllum japonicum (π = 0.0011; Zhu et al., 2020), Liquidambar formosana (π = 0.0044; Xu et al., 2024), and Davidia involucrata (π = 0.0047; Ren et al., 2024)—and comparable to Tetracentron sinense (π = 0.0115; Jing et al., 2025), though lower than Bretschneidera sinensis (π = 0.0636; X.-L. Zhu et al., 2024) and Euptelea pleiosperma (π = 0.0890; Cao et al., 2020). While cross-study comparisons should be interpreted cautiously (e.g., due to sequencing depth, filtering, reference genomes), these benchmarks suggest that P. macgregorii is not characterized by extremely low diversity, which is consistent with its current Vulnerable (VU) status rather than imminent extinction risk. Gene flow among populations likely contributes to maintaining genetic diversity and adaptive potential (Weeks et al., 2016, 2017; Clarke et al., 2024).

At the same time, small, isolated populations are prone to converting masked load into realized load via drift and inbreeding, while purifying selection can purge strongly deleterious alleles (Yang et al., 2018; van der Valk et al., 2019; Dussex et al., 2021; Khan et al., 2021; Kardos et al., 2023; Mooney et al., 2023; Robinson et al., 2023; Feng et al., 2024; Yi et al., 2024; Bemmels et al., 2025). Compared with the East lineage, the West lineage exhibited higher inbreeding, fewer heterozygous deleterious sites (DEL and TOL), and more homozygous deleterious sites. LOF burdens do not differ markedly between lineages in homozygotes, but the West retained fewer heterozygous LOF variants than the East. Together with severe population contraction since the onset of the Naynayxungla Glaciation (NG), these patterns suggest that bottleneck-associated inbreeding may have exposed (and partially purged) strongly deleterious LOF alleles while allowing weakly deleterious mutations to accumulate. Consistently, realized load, F_ROH, and the proportion of homozygous LOF sites within ROH were higher in the JXJGS and HNSHS populations, indicating active conversion of masked load to realized load via inbreeding and heightened vulnerability in these populations.

Genetic offset and RONA analyses have been used to predict genomic vulnerability in many species (Sang et al., 2022; Jiang et al., 2025; Long et al., 2025). Multiple offset metrics consistently pinpoint the same regions as genomic vulnerability hotspots. When considered together with elevated realized load and extensive runs of homozygosity (ROHs), these signals indicate that certain western populations face the greatest near-term risk, whereas several eastern populations exhibit distinctive variation consistent with incipient divergence. Overall, the concordance between climate-associated variation and patterns of genetic load highlights spatial heterogeneity in both risk and evolutionary potential across the species’ range, providing a results-based foundation for the management recommendations.

4.4. Conservation implications

Genomic evidence from SNP and SV datasets, demographic inference, and IBD jointly support the East and West lineages. These two evolutionarily cohesive lineages should be designated as the primary conservation units (CUs). Within each lineage, we further delineate secondary CUs to capture contrasting levels of risk and evolutionary distinctiveness. For example, JXJGS and HNSHS show elevated inbreeding, evidence of conversion of masked load to realized load, and occur within future genomic-offset hotspots in western Jiangxi and southwestern Hunan. These populations should be prioritized for habitat protection, demographic reinforcement, and cautious within-lineage assisted gene flow to reduce inbreeding while avoiding outbreeding depression. In contrast, JXGS and JXYJF exhibit consistent genetic differentiation and harbor unique variants indicative of incipient divergence; they should be managed as distinct secondary CUs to preserve evolutionary potential. Reverse-offset analyses indicate that parts of the southwestern range may become unsuitable under late-century high-emission scenarios, arguing for ex situ safeguards alongside in situ actions.

Acknowledgements

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

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 Genomics Data Center (NGDC; https://bigd.big.ac.cn/bioproject) under BioProject accession PRJCA050730 (genome assembly) and PRJCA051433 (whole-genome resequencing).

CRediT authorship contribution statement

Jiaxin Li: Writing–original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. Lihua Yang: Investigation, Writing–review & editing. Danqi Li: Investigation. Chen Feng: Supervision, Resources, Investigation, Writing–review & editing. Ming Kang: Supervision, Resources, Project administration, Funding acquisition, Conceptualization, Writing – review & editing.

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