Species endangerment represents one of the most pressing challenges to global ecological stability. The latest assessments reveal that many plant species face imminent extinction threats due to factors such as climate change, habitat loss, invasive species, and overexploitation (Thuiller et al., 2005; Humphreys et al., 2019). Consequently, understanding the mechanisms driving plant endangerment is urgently required to inform effective conservation strategies and mitigate further losses.

Myristica yunnanensis is a member of the ancient and ecologically significant Myristicaceae family within the Magnoliales order (Gentry, 1982; Pitman et al., 2001; Valencia et al., 2004). The Myristicaceae is under substantial conservation pressure, with 162 species assessed as endangered—well above the average per family (IUCN, 2024). As the sole representative of the Myristicaceae and the Myristica genus in mainland China (Li, 1976; Mao et al., 2019), M. yunnanensis possesses significantly economical and medicinal potential. Its seeds are rich in oils with diverse compositions (Zhang et al., 2019), and its tissues contain bioactive compounds such as β-sitosterol, ergosterol, afzelin, and quercitrin, known for anti-inflammatory and anti-tumor properties (Li and Ding, 2001). Despite these attributes, the species was listed as Critically Endangered (CR) within two decades of its discovery in 1976 (Li, 1976; Fu, 1992; World Conservation Monitoring Centre, 1998). Although found in parts of Southeast Asia, such as Vietnam and Thailand (Hoàng Văn Sâm, 2008), its wild distribution in mainland China has dramatically declined, now limited to a few isolated populations in southern Yunnan Province, which show clear signs of further decline in wild (Fig. S1). Historical field records from Xishuangbanna in the late 1990s documented ~20 wild trees; subsequent surveys continue to report only scattered, highly fragmented individuals with only a few localities in mainland China (Ma et al., 2017; Xu et al., 2017a, 2017b). In recognition of its precarious status, the species has been designated in Yunnan Province as a Plant Species with Extremely Small Populations (PSESP; Sun, 2021). Recent genomic advances have enabled investigations into genetic mechanisms behind population decline and offer valuable tools for conservation (Yang et al., 2018, 2024; Ma et al., 2022; Theissinger et al., 2023; Liu et al., 2024; Tao et al., 2024).

To date, no genome has been reported from the Myristicaceae, limiting our understanding of its phylogenetic position and evolutionary patterns. Here, we present the first chromosome-level genome assembly of Myristica yunnanensis, providing a foundational resource for investigating the genomic evolution of Myristicaceae. We further analyzed its genomic features and conducted population-level analyses using whole-genome resequencing of 20 wild individuals from mainland China, offering insights into the species’ lineage-specific dynamics and potential causes of its rarity of Chinese population.

The genome size of Myristica yunnanensis was estimated to be 628.82 Mb based on k-mer frequency analysis of 345.03 Gb BGI clean data (Fig. S2; Tables S1 and S2), and 638.06 Mb as determined by flow cytometry (Fig. S3 and Table S3). The de novo genome assembly of M. yunnanensis was performed with the NextDenovo program using 74.84 Gb (~122 ×) of long reads generated by Oxford Nanopore Technologies (ONT), followed by polishing with both ONT and 345.03 Gb (~566 ×) of BGI Next-generation sequencing (NGS) data. Haplotigs and redundant overlaps were removed using the purge_dups package, and chromosome-level scaffolds were generated using Hi-C contact maps derived from 227.84 Gb (~374 ×) of Hi-C data (Table S1). The initial assembled M. yunnanensis genome was 604.82 Mb (accounting for 96.18% of the estimated genome size), containing 104 contigs after long-read assembly and short-read correction (contig N50: 23.36 Mb) (Fig. 1a and b; Table S4). Around 97.73% of the sequences (591.11 Mb) were anchored to 19 pseudochromosomes, a number consistent with the chromosome numbers of several published species of Magnoliaceae (Fig. S4; Tables S5 and S6). The Bench-marking Universal Single-Copy Orthologs (BUSCO, v.5.2.2, embryophyta_odb10) analysis of the final assembly indicated an assembly completeness of 91.82% (Table S7), and the long terminal repeat (LTR) assembly index (LAI) score of M. yunnanensis was 20.40 (Fig. S5), indicating a high quality of genome assembly.

Transposable elements constituted 56.68% of the Myristica yunnanensis genome (Tables S8 and S9). Among the repetitive elements, LTR retrotransposons were the most abundant repetitive elements in the M. yunnanensis genome, constituting 36.10% of the total sequence, in which the Copia subfamily was the most dominant (23.73%), followed by the Gypsy subfamily (3.79%). Notably, DNA transposons of the Helitron subfamily represented a significant portion, making up 17.80% of the genome (Tables S8–S10). The majority of these repetitive sequences were located in intergenic regions, with a smaller fraction present within intronic regions (Fig. S6). Analysis of insertion timing suggests that DNA transposons are relatively ancient in the M. yunnanensis genome, whereas LTR transposons were inserted more recently (Fig. S7). For protein-coding genes annotation, a total of 22,002 protein-coding genes were predicted, averaging 9336.57 bp in gene length of and 1173.06 bp in CDS length (Table S11). The BUSCO completeness assessment resulted in a score of 91.40% (Table S7). Additionally, 93.42% of these genes were assigned putative functional annotations (Table S12). In summary, these data suggest high quality genomic data of M. yunnanensis, adequate for downstream analysis.

To explore the underlying causes of population decline, we used genomic data to infer the origin of Myristica yunnanensis. Possibly as an effect of the WGD, we were only able to identify 173 single-copy orthologous families from 12 species. Therefore, we opted to employ BUSCO genes for the extraction of the gene set used for phylogenetic analysis. A phylogenetic tree was constructed using 4DTv sites of 1152 single-copy BUSCO genes from 12 species, with Arabidopsis thaliana as the outgroup. Divergence times were estimated using MCMCtree, with four fossil calibration points sourced from Timetree (Figs. 1e and S8). Consistent with previous studies (Hu et al., 2019; Qin et al., 2021; Zheng et al., 2024), the phylogenetic tree revealed that Piperales, an early diverging group within the Magnoliids, diverged from Magnoliales and Laurales around 128.4 million years ago (Mya). Our data suggest the Myristicaceae is the basal clade in Magnoliales, diverging from the common ancestor of Magnoliaceae and Annonaceae approximately 101 Mya.

To assess the occurrence of WGD events in Myristica yunnanensis, we conducted a comparative genomic analysis using the following species: Aristolochia fimbriata (no WGD) (Qin et al., 2021), Liriodendron chinense (WGD event ~116 Mya) (Chen et al., 2019), and Nymphaea colorata (WGD event ~133 Mya) (Zhang et al., 2020). Synteny dot plot analysis revealed fragmented self-syntenic blocks with a corresponding syntenic depth ratio of 2:2 (Fig. S9). Genomic synteny comparisons of M. yunnanensis genome with the genomes of Aristolochia fimbriata, Liriodendron chinense, and Nymphaea colorata yielded corresponding collinearity depth ratios of 2:1, 2:2, and 2:2, respectively (Figs. 1c, S9 and S10). Furthermore, synonymous substitution rate (Ks) distribution of gene pairs in syntenic blocks peaked at 1.01, indicating that the WGD event occurred after the divergence from Piper nigrum and Arabidopsis thaliana, and before the divergence from Liriodendron chinense and Cinnamomum kanehirae (Fig. 1d). Integrating WGD estimations from other studies (Chaw et al., 2019; Chen et al., 2019; Dong et al., 2021; Yin et al., 2021), we infer that the Myristica WGD event was shared by Magnoliales and Laurales, but not by Piperales.

Eleven sequenced Magnoliids species were included for the gene family evolution analysis. A total of 27,773 ortholog groups were identified, of which 6370 orthogroups were present across all species. In Myristica yunnanensis, 19,795 genes (89.9% of the total annotated genes) were clustered into 13,021 gene families, among which 797 genes were clustered into 258 M. yunnanensis-specific gene families (Fig. 1f and Table S13). The evolution analysis of M. yunnanensis identified 126 expanded and 695 contracted gene families (p < 0.05) (Fig. 1e). Compared to species from the Annonaceae and Magnoliaceae families, M. yunnanensis possessed fewer expanded gene families, which might be attributed to its relatively reduced gene content (Fig. 1e). Functional enrichment analysis of KEGG and GO revealed that expanded gene families were significantly enriched in pathways related to “membrane” and “microtubule” (Figs. 1g and S11; Tables S14 and S15). The contracted gene families were significantly enriched in pathways related to modification processes such as “protein phosphorylation”, “protein modification process” and “macromolecule modification” (adjusted p < 0.01) (Fig. S12; Tables S16 and S17). Given the key role of gene duplication in plant evolution, we identified duplicate genes in M. yunnanensis. They were divided into four types: 3682 WGD, 2531 tandem (TD), 1534 proximal (PD), and 9569 dispersed duplicates (DD), with DD being the most abundant at 44% (Fig. S13 and Table S18). Since tandem duplicates often help plants adapt to changing environments (Hanada et al., 2008; Lex and Jonathan, 2009; Schaper and Anisimova, 2015), KEGG and GO enrichment analyses were performed on the 2531 TD genes. The results revealed that these genes were significantly enriched in certain pathways such as “cellular response to alcohol”, “abscisic acid-activated signaling pathway”, and “metabolism of terpenoids and polyketides” (Tables S19 and S20).

To investigate the genetic status and evolutionary history of wild Myristica yunnanensis in mainland China, we resequenced 20 individuals from a highly restricted region in Xishuangbanna at an average depth of ~95 × (Fig. S14 and Table S21). After variant calling, over two million high-quality single nucleotide polymorphisms (SNPs) were obtained (Fig. S15). Population structure analyses (phylogeny, principal component analysis (PCA), and admixture analysis) revealed two closely related genetic clusters, within which individuals exhibited close kinship (Figs. 1h, 1i and S16). Genome-wide estimates of genetic differentiation and nucleotide diversity support this pattern: the weighted average F_ST between the two clusters was 0.044, and the mean D_XY was 0.002, calculated based on all-site data and SNP-only data (Fig. S17). Within-population nucleotide diversity (π, per site) was also similar between the clusters (0.00168 vs. 0.00174), and close to the overall population-level diversity (0.00174), which is lower than most other endangered trees (Table S22). Inbreeding levels were low (average FROH = 0.028; Table S23), and no significant difference in inbreeding was observed between clusters (0.025 vs. 0.032, p = 0.08186, Welch two-sample t-test). Genome-wide diversity was extremely low among wild M. yunnanensis in China, with heterozygosity of 1.487 per kb (Table S23). A pronounced Tajima's D peak near 1.5 suggests potential signatures of balancing selection, although it may also result from population structure or recent demographic contraction. The secondary negative peak could be related to historical bottlenecks or purifying selection (Fig. S18) (Tajima, 1989; Nielsen, 2005). Thus, we further inferred the demographic history using PSMC and SMC++ software and found two pronounced declines in effective population size: one ~8–12 million years ago and another ~80,000 years ago. These declines broadly coincide with global cooling events in the late Miocene and the Last Glacial period (Fig. 1j) (Clark et al., 2009; Burke et al., 2018; Holbourn et al., 2018). Although previous analyses indicated that the two groups were highly similar, we nonetheless treated them as distinct conservation units and inferred their demographic histories separately. Substructure analysis showed that only the northern subpopulation may have experienced a recent recovery (Fig. S19). These demographic patterns may have been influenced, at least in part, by historical climate fluctuations.

A significant impact of a small effective population size is the reduced efficacy of natural selection, leading to an increase in the number of deleterious mutations (Yang et al., 2018; Dussex et al., 2023). Such accumulation imposes a genetic load, which refers to the reduction in mean fitness of a population due to the presence of deleterious alleles, relative to an idealized, mutation-free genotype (Agrawal and Whitlock, 2012; Grossen and Ramakrishnan, 2024). To evaluate the efficacy of selection and explore the patterns of genetic load, we analyzed the deleterious mutations in the wild population of Myristica yunnanensis in China. A total of 30,399 variants affecting 3953 genes were jointly identified as deleterious (DEL) or tolerated (TOL) variants (Table S24). These included 3833 DEL variants affecting 1930 genes, and 1276 loss-of-function (LOF) variants affecting 1122 genes (Fig. S20; Tables S24 and S25). Deleterious mutations were broadly distributed across the genome. A region on Chr9 exhibited a high SNP density, notably including multiple LOF variants between 8.2 and 8.6 Mb region of Chr9 (Fig. S21). GO and KEGG analyses revealed that affected genes were enriched in phosphorylation, protein modification, environmental adaptation, and plant–pathogen interaction pathways (Tables S26 and S27). In particular, six LOF genes located on Chr9 (8.2–8.6 Mb) encoded LRR receptor-like protein kinases (Fig. 1k; Tables S28 and S29), including Myu11015 and Myu11017, orthologous to AT4G08850 (encoding MIK2), AT1G35710 of Arabidopsis thaliana, and Os01t0296000, Os02t0553000 of Oryza sativa. These transmembrane receptors play roles in immunity, development, and stress signaling (Van der Does et al., 2017; Coleman et al., 2021). The presence of such deleterious mutations highlights the genetic load linked to population decline and may further threaten the survival of M. yunnanensis.

In this study, we assembled a chromosome-level genome for the critically endangered Myristica yunnanensis, providing the first genomic insight into the evolution and population dynamics of the Myristicaceae. Phylogenomic analyses showed that Myristicaceae represented the basal clade within Magnoliales, which diverged from the common ancestor of Magnoliaceae and Annonaceae approximately 101 Mya. Compared with other magnoliids, M. yunnanensis exhibited notable gene family contraction, especially in those related to post-translational modifications. Sequencing of 20 wild individuals from southwestern China revealed low genetic diversity, high genetic load, and functional loss in gene clusters linked to immunity, signaling, and water transport. Demographic inference further suggested population contraction during glacial periods. Together, our results indicate that both genomic constraints and climatic fluctuations may have contributed to the endangered status of M. yunnanensis population of China, offering important implications for tropical plant conservation.

Acknowledgments

We thank many people who provided feedback, samples and encouragement, especially Xishuangbanna Tropical Botanical Garden. This work was supported by the National Key R&D Program of China (2022YFC3400300).

Data availability statement

All raw sequencing data from Illumina, ONT, Hi-C, and RNA-seq used for genome assembly, annotation and WGS data for Myristica yunnanensis population analysis have been deposited in the National Genomics Data Center (https://ngdc.cncb.ac.cn/) database, under Bioproject PRJCA034816 and PRJCA034985 respectively. The genome assembly has also been deposited in the Genome Warehouse with accession number GWHGDHP00000000.1.

CRediT authorship contribution statement

Yongting Zhang: Writing – original draft, Visualization, Validation, Investigation, Formal analysis, Data curation, Conceptualization. Zihe Li: Writing – original draft, Methodology, Formal analysis, Data curation. Xue Liu: Validation, Methodology, Investigation. Peng Zeng: Visualization, Methodology. Chuan Peng: Resources. Botong Zhou: Writing – review & editing. Yingmei Peng: Methodology, Data curation. Wenbo Zhu: Methodology. Jian Huang: Resources. Jing Cai: Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition.

Declaration of competing interest

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

Appendix A. Supplementary data

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