
==== Front
Gigascience
Gigascience
gigascience
GigaScience
2047-217X
Oxford University Press

39311762
10.1093/gigascience/giae070
giae070
Research
AcademicSubjects/SCI00960
AcademicSubjects/SCI02254
Genomic insights into endangerment and conservation of the garlic-fruit tree (Malania oleifera), a plant species with extremely small populations
Shen Yuanting Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
University of Chinese Academy of Sciences, Beijing 100049, China
State Key Laboratory of Plant Diversity and Specialty Crops, Institute of Botany, Chinese Academy of Sciences, Beijing 100093, China

https://orcid.org/0000-0002-1396-0524
Tao Lidan Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
University of Chinese Academy of Sciences, Beijing 100049, China

https://orcid.org/0000-0002-8028-9208
Zhang Rengang Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
University of Chinese Academy of Sciences, Beijing 100049, China

https://orcid.org/0000-0002-3628-7088
Yao Gang Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China

Zhou Minjie Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
University of Chinese Academy of Sciences, Beijing 100049, China

https://orcid.org/0009-0009-5246-9226
Sun Weibang Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China

https://orcid.org/0000-0002-7725-3677
Ma Yongpeng Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China
Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China

Correspondence address: Yongpeng Ma, NO. 132 Lanhei Road, Panlong Street, Kunming City, Yunnan Province, China. E-mail: mayongpeng@mail.kib.ac.cn
Correspondence address: Weibang Sun, NO. 132 Lanhei Road, Panlong Street, Kunming City, Yunnan Province, China. E-mail: wbsun@mail.kib.ac.cn
These authors contributed equally to this work.

23 9 2024
2024
23 9 2024
13 giae07010 5 2024
17 7 2024
22 8 2024
© The Author(s) 2024. Published by Oxford University Press GigaScience.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Advanced whole-genome sequencing techniques enable covering nearly all genome nucleotide variations and thus can provide deep insights into protecting endangered species. However, the use of genomic data to make conservation strategies is still rare, particularly for endangered plants. Here we performed comprehensive conservation genomic analysis for Malania oleifera, an endangered tree species with a high amount of nervonic acid. We used whole-genome resequencing data of 165 samples, covering 16 populations across the entire distribution range, to investigate the formation reasons of its extremely small population sizes and to evaluate the possible genomic offsets and changes of ecology niche suitability under future climate change.

Results

Although M. oleifera maintains relatively high genetic diversity among endangered woody plants (θπ = 3.87 × 10−3), high levels of inbreeding have been observed, which have reduced genetic diversity in 3 populations (JM, NP, and BM2) and caused the accumulation of deleterious mutations. Repeated bottleneck events, recent inbreeding (∼490 years ago), and anthropogenic disturbance to wild habitats have aggravated the fragmentation of M. oleifera and made it endangered. Due to the significant effect of higher average annual temperature, populations distributed in low altitude exhibit a greater genomic offset. Furthermore, ecological niche modeling shows the suitable habitats for M. oleifera will decrease by 71.15% and 98.79% in 2100 under scenarios SSP126 and SSP585, respectively.

Conclusions

The basic realizations concerning the threats to M. oleifera provide scientific foundation for defining management and adaptive units, as well as prioritizing populations for genetic rescue. Meanwhile, we highlight the importance of integrating genomic offset and ecological niche modeling to make targeted conservation actions under future climate change. Overall, our study provides a paradigm for genomics-directed conservation.

recent inbreeding
deleterious mutation
demographic history
genomic offset
ecological niche modeling
conservation genomics
Natural Science Foundation of Yunnan Province 10.13039/501100005273 202001AS070019
==== Body
pmcIntroduction

Historical climate disturbances and frequent human activity have caused many species that were once widespread with continuous distributions to become small, fragmented populations [1]. High levels of inbreeding continually occur in these populations, leading to the accumulation of deleterious mutations and low species adaptability, ultimately increasing the risk of extinction [2, 3]. The genome contains evolutionary footprints that can be used to estimate inbreeding levels of species even without detailed pedigrees [4, 5]. For example, runs of homozygosity (ROH), genome regions with a certain length that is identical by descent, have been widely used as an indicator of inbreeding [6, 7]. The long ROH indicate a closer relationship to the most recent common ancestor, implying a higher level of inbreeding. For small and isolated populations with high inbreeding levels, genetic rescue is necessary to introduce beneficial mutations by establishing gene flow between populations [8]. However, caution should be taken when making decisions regarding genetic rescue, and comprehensive exploration of the genetic background of these small populations must be done in advance [9, 10].

Rapid climate change in the future is a widely recognized threat to global biodiversity [11–13]. The threatened degree of species depends on how they respond to climate change. There are 2 main mechanisms: migrating to new habitats or adapting to the changing environments through phenotypic plasticity or de novo mutations [14–16]. However, in the case of long-lived forest trees, individual organisms are nearly incapable of migrating to keep pace with changing climate and may experience maladaptation throughout their lifetimes [17–19]. Thus, if species’ responses to future climate change can be predicted, their extinction risks can be estimated, and targeted conservation guidelines and management strategies can be developed in advance.

Malania oleifera Chun & S. K. Lee (NCBI:txid397392), the single species in the genus Malania (Olacaceae), is an endemic, semi-parasitic evergreen tree (diploid) that naturally scattered in the west Guangxi (a.s.l. 300–1,000 m) and southeast Yunnan province, China (a.s.l. 300–1,640 m) [20]. It adapts well to rocky desert habitats and can be used as an afforestation tree in karst landscapes [21]. Moreover, it has extremely high economic and medicinal value due to the large amount of lipids in its seed. The main lipid component is nervonic acid (cis-tetracos-15-enoic acid, >60%), which is essential for human nervous health [22]. However, mainly due to overexploitation, wild resources of M. oleifera have decreased by approximately 25,000 individuals between 2000 and 2017 solely in Guangnan County (Yunnan, China) [23, 24]. Additionally, physiological factors of M. oleifera, including large seed size, short seed life span, difficulty in natural seed germination, low rate of pollen germination, and susceptibility to root rot, have made natural propagation and regeneration difficult [25–28]. Therefore, M. oleifera has been categorized as Vulnerable (VU) on the IUCN Red List (https://www.iucnredlist.org/species/32361/9701100) and recorded in the Class II Key Protected Wild Plant List in China [29], and it has also been listed as a plant species with an extremely small population size in China [30], highlighting the urgent need for its conservation.

Previous studies related to M. oleifera mainly focused on the biosynthetic pathway of nervonic acid [22, 31] and exploring the optimal condition of growing artificial seedlings for better utilization [32, 33]. However, the conservation process of M. oleifera is limited to in situ protection of existing wild resources [34]. Recent advancements in whole-genome sequencing techniques enable covering nearly all the nucleotide variations of a genome and can provide deep insights into protecting endangered species [35]. However, there is a gap in using genomic data to guide conservation strategies, particularly for plants. Therefore, we utilize M. oleifera as a study case to investigate the above issues. We aimed to provide a comprehensive framework for the M. oleifera conservation through the perspective of conservation genomics.

Materials and Methods

Sample collection and whole-genome resequencing

A total of 165 leaf samples were collected from 16 wild populations across the entire distribution of M. oleifera from Yunnan and Guangxi provinces, China (Fig. 1A; Supplementary Table S1). Among them, 76 samples were collected based on our field investigation, and the remaining 89 samples were obtained from the Germplasm Bank of Wild Species in Southwest China. The number of individuals collected per population (including the samples from germplasm bank) varied from 5 to 17, and if populations contained fewer than 10 individuals, all individuals were sampled.

Figure 1: Population genomics of M. oleifera. (A) Geographic distribution and sampled populations of M. oleifera. Different colors in the pie chart represent the genetic groups identified by ADMIXTURE based on adaptive loci, and the size of the pie corresponds to the level of heterozygosity. The optimal population genetic structure of M. oleifera with K = 10 (B) and a neighbor-joining (NJ) phylogenetic tree (C) based on adaptive loci. Samples in the STRUCTURE and phylogenetic tree results correspond. Node bootstrap values below 0.8 are not shown. (D) Results of PCA based on adaptive loci, with the first 2 PCs explaining 25.2% of the genome covariance. Populations are defined as BB1 = Banbeng, BB2 = Babao, BL = Banlun, BM1 = Bama, BM2 = Bamei, DX = Daxin, FS = Fengshan, GL = Gaolong, JM = Jiumo, LY1 = Leye, LY2 = Linyun, ML = Mulun, NP = Nanping, SG = Shuguang, ZL = Zhemiao, and ZS = Zhesang.

Genomic DNA was extracted from silica-dried leaf tissues using a modified CTAB method [36], and the concentration and quality of the DNA was determined using a NanoDrop2000 Spectrophotometer (Thermo Fisher Scientific). Samples were sent to Beijing Ori-Gene Science and Technology Co., Ltd for Illumina sequencing library preparation according to the manufacturer’s specifications. Paired-end raw reads (150 bp) were generated on the Illumina HiSeq platform.

Read mapping and single nucleotide polymorphism calling

The raw data were filtered using Fastp (RRID:SCR_016962) v. 0.19.3 [37]. Paired-end clean reads were mapped to the chromosome-level genome of M. oleifera (∼1.5 Gb) [31] using BWA-MEM (RRID:SCR_022192) v. 2.1 [38]. SAMtools (RRID:SCR_002105) v. 1.9 [39] was used to convert sequence alignment map (SAM) format files to sorted binary alignment map (BAM) format files. Sambamba (RRID:SCR_024328) v.0.7.1 [40] was used to mark and remove duplicate reads. Freebayes (RRID:SCR_010761) v. 1.3.6 [41] was employed to call variants, and only bases with a quality score ≥20 and reads with a mapping quality score ≥30 were included. This produced a total of 43,413,408 initial variant sites. We then employed VCFtools (RRID:SCR_001235) v. 0.1.15 [42] to filter sites with the following criteria: (i) sites with coverage depth below 1/2 * average site coverage and above 2 * average site coverage were discarded after investigating the coverage distribution; (ii) sites located on the organelle genomes or contigs that were not anchored on the chromosomes were excluded; (iii) single nucleotide polymorphisms (SNPs) with a depth below 3× or a genotype quality score <20 were redefined as missing; (iv) only biallelic SNPs were reserved; (v) SNPs with a missing rate >20% were removed, leaving 2,144,506 SNPs (dataset 1); and (vi) SNPs with a minor allele frequency <0.05 were all excluded. Finally, 250,362 SNPs remained (dataset 2), which distributed on 13 pseudochromosomes (Supplementary Fig. S2).

Population genetic diversity and runs of homozygosity

Based on dataset 2, we detected a genome-wide linkage disequilibrium (LD) decay among 16 populations using PopLDdecay (RRID:SCR_022509) v. 3.4.0 [43]. Nucleotide diversity (θπ), Watterson’s θ (θw), and heterozygosity rate were calculated using ANGSD (RRID:SCR_021865) v. 0.921 [44] based on bam files, which removed duplicates (see above). In addition, we calculated the values of the 3 parameters in more specific genomic regions (intergenic, CDS, intron, fold-0, and fold-4). To examine inbreeding depression, we detected ROH using vcftools v. 0.1.15 [42] based on dataset 2 with the key parameters “–LROH”, and only ROH longer than 100 kb were kept. Moreover, we calculated the frequency of runs of homozygosity (FROH), which is equal to the sum of all ROH lengths longer than 100 kb divided by genome effective length [45].

Inference of population structure based on all loci, adaptive loci, and neutral loci

We employed 2 software to detect outlier SNPs potentially related to adaptive evolution. First, we used the sparse nonnegative matrix factorization (snmf) function applied in the R package LEA (RRID:SCR_009090) v. 3.1.4 [46] to estimate the most likely number of ancestral populations based on dataset 2. To reduce the number of false positives, we reserved SNPs with the false discovery rate (FDR) less than 0.01. Second, we applied a principal component analysis (PCA) method using R package Pcadapt (RRID:SCR_022019) v. 4.3.3 [47] with a 0.01 cutoff of FDR to identify SNPs that highly influenced the formation of observed differentiation. SNPs detected by both methods were identified as potential adaptive SNPs; otherwise, they were considered as neutral SNPs. Finally, we employed PLINK (RRID:SCR_001757) v. 1.90b4.1 [48] to filter out linkage disequilibrium sites and ultimately obtained 33,971 (dataset 3, all loci), 1,515 (dataset 4, adaptive loci), and 32,930 (dataset 5, neutral loci) SNPs for downstream analysis (Supplementary Fig. S1).

Based on the 3 datasets described above, we employed ADMIXTURE (RRID:SCR_001263) v. 1.3.0 [49] to infer the population structure. The most likely population number of K was determined by the minimizing cross-validation error. PCA was conducted in GCTA v1.94.1 [50] and MEGA (RRID:SCR_000667) v. 7.0 [51] was used to construct NJ trees. Pairwise fixation statistics (Fst) among the 16 populations were calculated using vcftools v. 0.1.15 [42].

Estimation of demographic history

We employed Stairway Plot v.2 [52] and MSMC (RRID:SCR_023677) v.2 [53] to infer the population demographic history of M. oleifera. For the analysis of the stairway plot, we first performed ancestral sequence reconstruction to infer ancestral status (see details in Supplementary Note S1 and Supplementary Table S16). To mitigate the effects of selection, we excluded the upstream and downstream 5-kb regions of genes to infer folded site frequency spectrum (SFS) and unfolded SFS using ANGSD v. 0.921 [44]. We set the average generation time of M. oleifera as 10 years because it takes about 10 years for a seed to grow into a seed-producing plant according to our field observations. The mutation rate was set to 2.5 × 10−8 per site per generation (see details in Supplementary Note S2 and Supplementary Table S17). MSMC is a multiple sequentially Markovian coalescent approach that uses the density of heterozygous sites to estimate the effective population size (Ne) through time. Different individual numbers or haplotypes provide distinct resolutions for the analysis of demographic histories. Therefore, we employed MSMC to separately estimate the coalescence rate within 2, 4, and 8 haplotypes, referring to the simulation results of Schiffels and Durbin [53]. A total of 100 random combinations of individuals were used for the 3 haplotype analyses to estimate the medians and 95% confidence interval values. The average generation time and mutation rate values were set to the same as stairway plot analysis.

Detection of deleterious mutations

Deleterious mutations in M. oleifera were predicted using the Sorting Intolerant From Tolerant (SIFT) algorithm [54]. We used a modified approach to perform the SIFT prediction (see Supplementary Note S3). The TrEMBL plant database [55] was used to search for orthologous genes, and the SIFT scores were calculated based on the degree of conservation among loci. Based on dataset 1 (included low-frequency variants), SNPs in the coding regions were categorized as deleterious (SIFT score <0.05), tolerated (SIFT score ≥0.05), or synonymous using SIFT4G (RRID:SCR_021850) [56]. The low confidence sites and “NA” sites were not considered.

To provide accurate and direct genetic rescue guidance for M. oleifera populations with a high genetic load, we selected 5 populations, including 1 as the potentially threatened population and 4 as the candidate pollen donors. We drew a Venn diagram to explore the distribution of shared or unique homozygous deleterious mutations among the 5 populations. The 4-candidate pollen donors were characterized by (i) having low genetic load, low levels of inbreeding, high genetic diversity, and high heterozygosity; (ii) low genetic differentiation with the rescued population; or (iii) sharing the same genetic lineage as the rescued population based on adaptive loci.

Identification of environment-associated adaptive variants

The environmental data, which included 19 bioclimatic variables at 2.5-minute resolution (5 km), were downloaded from WorldClim v.2.1 database (RRID:SCR_010244) (Supplementary Table S9). Each environmental factor was extracted through the coordinates of sampling points. We selected the top climatic factors based on weighted R2 importance (Supplementary Fig. S14) using a machine learning gradient forest (GF) model in the R package gradientforest v.0.1–37 [57] by modeling the relationship of climatic variables and SNPs with 500 regression trees. To avoid multicollinearity, we kept variables with |correlation coefficient| <0.7 by calculating Pearson’s correlation coefficient in the R package corrplot v.0.92 (RRID:SCR_024683) [58] and ultimately retained the 4 most important and uncorrelated factors, including BIO3 (Isothermality), BIO7 (Temperature Annual Range), BIO14 (Precipitation of the Driest Month), and BIO15 (Precipitation Seasonality). Next, we used BayeScEnv [59] and redundancy analysis (RDA) [60] to identify environment-associated SNPs. BayeScEnv represents a univariate genotype–environment association approach. For BayeScEnv method, the input files included an environmental factor standardized by the mean variance and contained codominant data, which converted by PGDSpider v. 2.1.1.5 [61] based on dataset 2. RDA is a multivariate linear regression-based method [62]. We ran RDA analysis in the R package vegan v.2.5–7 [63], using function “anova.cca” to check the significance of the RDA model and function “outliers” to identify local adaptation–associated SNPs that loaded in the tails of the ±3 standard deviation cutoff (2-tailed P = 0.0027). To recognize the gene functions of the candidate SNPs obtained from BayeScEnv and RDA, we employed a Gene Ontology enrichment analysis using the eggNOG-mapper (RRID:SCR_021165) v.2 [64].

Genomic offset modeling with gradientforest

We used the GF model in the R package gradientforest v.0.1–37 [57] to predict genomic vulnerability to future climate change. The environmental-associated SNPs detected by both BayeScEnv and RDA were denoted as the candidate SNP dataset. In addition, 500 randomly selected SNPs based on dataset 2 were recorded as a reference SNP dataset to match the magnitude of candidate SNP dataset. The SNPs data with minor allele frequencies (MAFs) >10% were converted into MAFs per population. To ameliorate the linkage effect, we only kept 1 SNP per 100,000 bp range and finally obtained 326 reference SNPs and 213 candidate SNPs. We employed the GF model with 500 regression trees per SNP to build a function for the 4 most important environmental factors (BIO3, BIO7, BIO14, BIO15). Genomic offset (GO) was defined by Euclidean distance between current (1970–2000) and future (2081–2100) climate, which used the current condition as the baseline [16]. To explore possible future climate conditions and predict GO for M. oleifera, we employed 3 widely used global climate models (BCC-CSM2-MR, CNRM-CM6-1, and CNRM-ESM2-1) and 2 emission scenarios (SSP126 and SSP585) that represent the mild and extreme future carbon emissions. The predicted GO results of 3 global climate models for each grid were averaged with assigned weights.

Ecological niche modeling

The distribution records of M. oleifera were collected from online databases, published academic articles [65, 66], and field investigation, and all records were manually verified using an online map. To reduce sampling bias, we kept only 1 record within 5 km using the rarefy function of the R package Humboldt [67], with 87 records remaining (Supplementary Table S10). The 19 bioclimatic variables (Supplementary Table S9) were also used in ecological niche modeling. Since background or pseudo-absence data of ecological niche models were sampled from the entire modeling map, we recalculated the correlation coefficients of 19 bioclimatic variables based on the whole distribution area. We removed autocorrelated variables (|Pearson’s r| >0.7 and variance inflation factor >10) using the R package usdm v.2.1 [68] and kept 6 uncorrelated variables (BIO1, BIO2, BIO7, BIO12, BIO14, and BIO18) for ecological niche modeling.

The ecological niche model (ENM) was built using an ensemble modeling method that combined outputs of 5 single models with high performance: GAM (generalized additive model by the R package mgcv v.1.9) [69], MaxEnt (tuned MaxEnt model by the R package dismo v.1.3) [70], RF (random forest with downsampling by the R package randomForest v.4.7) [71], Lasso (by the R package glmnet v.4.1) [72], and BRT (boosted regression trees by the R package dismo) [73]. Each model was run for 10 replicates, with pseudo-absence data of 10,000 points randomly generated using the R package Biomod2 [74] for 3 replicates, resulting in a total of 5 * 10 * 3 = 150 single models. However, only models with positive Somer’s D values were employed to create the final ensemble prediction, which was weighted by the true skill statistic (TSS) value of each model. The evaluation of the single model and ensemble model was performed by the R packages Ecospat v.4.0.0 [75] and prg v.0.5.1 [76].

The niche suitability of M. oleifera under future (2081–2100) carbon emission scenarios (SSP126 and SSP585) was predicted using the same climate models (BCC-CSM2-MR, CNRM-CM6-1, and CNRM-ESM2-1) as GO analysis. We used R package PresenceAbsence v. 1.1.11 [77] to calculate the threshold of ecological niche suitability, and grids with values higher than the threshold were defined as suitable habitats. Furthermore, referring to the method of Chen et al. [16], we defined niche suitability change (NSC) as the niche suitability index in the current climate minus the niche suitability index in the future climate. A positive value implies that niche suitability will decrease in the future compared to the present condition, while a negative value means increasing niche suitability. The NSC results of a single model for each emission scenario were averaged.

Results

Population structure and phylogeny of M. oleifera

Whole-genome resequencing generated an average of ∼4.71 Gb raw data and 65,007,259 paired-end reads for each sample, and the average sequencing depth was 6.5-fold. After filtering, the average Q20 and Q30 rates of paired-end reads were 97.51% and 92.68%, respectively, with an average mapping rate of 99.39% (Supplementary Tables S2 and S3). Neutral and adaptive genomic variations have inconsistent evolutionary patterns [78] and provide different types of information when determining optimal conservation measures [79]. To disentangle these discrepancies, we used 3 SNP datasets, including all loci (dataset 3), adaptive loci (dataset 4), and neutral loci (dataset 5), to decipher the genetic relationships within M. oleifera by constructing population structure, PCA, and phylogenetic trees.

The ADMIXTURE analysis results from all loci and neutral loci both revealed the optimal number of cluster (K) was 14 (Supplementary Fig. S7). Samples in most populations were relatively pure with no or only a mild genetic mixture with other populations, except for SG, ML, FS, and LY2 populations (SupplementaryFigs. S8a and S9a). However, ADMIXTURE analysis based on adaptive loci indicated that K = 10 was optimal (Supplementary Fig. S7), with ML-BB2, BB1-BM2, and ZL-LY1 paired populations having the same genetic composition, implying the paired population has similar adaptability, respectively (Fig. 1B). Notably, DX and GL populations were found to be 100% pure based on ADMIXTURE analysis of all datasets. Measures of PCA based on all datasets revealed clear separation of the DX population from other populations by PC1 and PC2, which explained 29.6%, 68.4%, and 25.2% of the genome covariance based on the results of all loci, neutral loci, and adaptive loci, respectively (Fig. 1D, Supplementary Figs. S8b and S9b). The NJ trees based on all datasets were consistent with the corresponding ADMIXTURE analysis, showing populations with similar genetic components had closer phylogenetic relationships (Fig. 1C, Supplementary Figs. S8c and S9c).

Genetic diversity, heterozygosity, and genetic differentiation

The average whole genomic genetic diversity of M. oleifera was 3.87 × 10−3 ± 1.34 × 10−3 for pairwise nucleotide differences (θπ) and 3.46 × 10−3 ± 1.23 × 10−3 for Watterson"s θ (θw) (Table 1, Supplementary Table S4 ). The BM1 population showed the highest θπ and θw compared with other populations, while NP, JM, and BM2 populations had lower θπ and θw (Supplementary Fig. S4) and showed a more obvious sawtooth-like distribution pattern of genetic diversity across the genome (Supplementary Fig. S5). When we divided the genome into 5 specific genomic regions, the values of θπ and θw showed the trend of intergenic > fold-4 > intron > CDS > fold-0, which was highly consistent among all 16 populations (Table 1, Supplementary Table S4). The mean heterozygosity rate in M. oleifera was 0.50% ± 0.14%, and it varied among populations, with the BM1 (0.56% ± 0.16%) population showing the highest values and with the lowest heterozygosity rate seen in JM (0.30% ± 0.04%), NP (0.31% ± 0.04%), and BM2 (0.35% ± 0.02%) populations (Supplementary Table S5). As expected, the results of the heterozygosity rate across more specific genomic regions showed the same trend as the genetic diversity (Supplementary Fig. S6). Moreover, the genome-wide LD decay analysis revealed that the level of LD varied greatly between populations, with the BM2 population showing the slowest decay of LD, whereas the BM1 population had the fastest LD decay (Supplementary Fig. S3).

Table 1: Sample sizes and nucleotide diversity (θπ) in M. oleifera populations within the whole genome, CDS, fold-0, fold-4, intergenic, and intron regions.

Population	Sample size	Number of SNPs	θπ_whole (× 10−3)	θπ_CDS (× 10−3)	θπ_fold-0 (× 10−3)	θπ_fold-4 (× 10−3)	θπ_intergenic (× 10−3)	θπ_intron (× 10−3)	
BB1	10	141,725	3.10 ± 0.42	1.19 ± 0.15	0.99 ± 0.13	1.93 ± 0.23	3.52 ± 0.49	1.83 ± 0.18	
BB2	10	132,480	2.79 ± 0.41	1.08 ± 0.16	0.90 ± 0.13	1.75 ± 0.25	3.17 ± 0.47	1.68 ± 0.24	
BL	10	146,763	3.56 ± 0.47	1.38 ± 0.17	1.15 ± 0.14	2.19 ± 0.27	4.07 ± 0.57	2.12 ± 0.25	
BM1	17	218,460	6.13 ± 0.58	2.41 ± 0.18	2.00 ± 0.15	3.92 ± 0.26	7.56 ± 0.44	3.87 ± 0.20	
BM2	5	101,436	2.15 ± 0.25	0.86 ± 0.13	0.73 ± 0.11	1.34 ± 0.17	2.45 ± 0.29	1.30 ± 0.13	
DX	10	155,019	4.57 ± 0.57	1.73 ± 0.27	1.44 ± 0.22	2.77 ± 0.43	5.20 ± 0.62	2.71 ± 0.39	
FS	9	189,667	5.13 ± 0.41	1.94 ± 0.22	1.61 ± 0.18	3.15 ± 0.31	5.86 ± 0.45	3.04 ± 0.26	
GL	10	137,835	3.56 ± 0.45	1.35 ± 0.18	1.13 ± 0.15	2.13 ± 0.26	4.11 ± 0.52	2.05 ± 0.26	
JM	10	98,031	2.07 ± 0.40	0.83 ± 0.20	0.69 ± 0.17	1.34 ± 0.30	2.36 ± 0.46	1.25 ± 0.22	
LY1	10	157,624	3.89 ± 0.47	1.47 ± 0.22	1.22 ± 0.19	2.39 ± 0.33	4.47 ± 0.52	2.25 ± 0.31	
LY2	14	200,098	5.50 ± 0.27	2.03 ± 0.19	1.68 ± 0.16	3.28 ± 0.26	6.28 ± 0.27	3.20 ± 0.21	
ML	10	151,563	3.53 ± 0.40	1.34 ± 0.18	1.12 ± 0.15	2.13 ± 0.26	4.05 ± 0.47	2.08 ± 0.21	
NP	10	106,945	2.03 ± 0.48	0.82 ± 0.25	0.69 ± 0.20	1.30± 0.38	2.31 ± 0.54	1.25 ± 0.30	
SG	10	133,980	2.75 ± 0.34	1.10 ± 0.20	0.91 ± 0.17	1.80 ± 0.30	3.12 ± 0.40	1.65 ± 0.18	
ZL	10	169,320	4.48 ± 0.41	1.67 ± 0.21	1.41 ± 0.17	2.61 ± 0.32	5.15 ± 0.43	2.59 ± 0.29	
ZS	10	163,672	4.35 ± 0.32	1.63 ± 0.14	1.36 ± 0.12	2.60 ± 0.21	4.99 ± 0.36	2.55 ± 0.20	

The values of pairwise Fst based on adaptive loci (dataset 4) were significantly higher than those of Fst based on all loci (dataset 3) and neutral loci (dataset 5), showing that adaptive loci > all loci > neutral loci in all paired populations (Supplementary Table S6). This was particularly prominent between the DX population and other populations, with a significant high pairwise Fst ranging from 0.87 to 0.91 based on adaptive loci, compared to 0.26–0.46 and 0.20–0.41 based on all loci and neutral loci, respectively (Supplementary Table S6). Moreover, BM2, JM, and NP populations, which have the lowest genetic diversity, showed high genetic differentiation from other populations (average pairwise Fst = 0.54 based on adaptive loci). In contrast, the BM1 population, with the highest genetic diversity, showed relatively low genetic differentiation from other populations (average pairwise Fst = 0.45 based on adaptive loci) (Supplementary Table S6).

Demographic history of M. oleifera

The stairway plot detected 2 severe population declines of M. oleifera based on unfolded SFS. The first occurred around 0.5–0.22 million years ago, corresponding to the Middle Pleistocene with climate upheaval, and the Ne was reduced to ∼8,230 (Fig. 2A). Subsequently, all the populations quickly recovered to ∼2.4 × 105 and remained stable until a recent bottleneck at around 10 Ka during the last glacial maximum (LGM), where there was a sharp population contraction to its lowest level (∼1,676). The result based on folded SFS also showed 2 bottleneck events at the corresponding time (Supplementary Fig. S10). MSMC tracked the more recent demographic trajectory of M. oleifera, especially within the past 10,000 years (Fig. 2B). Based on the results of 2, 4, and 8 haplotypes, the Ne of M. oleifera experienced a significant decline over time, reaching a nadir (below 75) around 400 to 500 years ago, followed by a slight population expansion. It is worth mentioning that both programs detected a population decline in M. oleifera during the LGM, which strengthened the reliability of the results.

Figure 2: Demographic history of M. oleifera inferred by Stairway Plot v.2 based on unfolded SFS (A) and MSMC v.2 within 2, 4, and 8 haplotypes (B). The light blue lines correspond to the upper and lower bounds of the 95% confidence intervals. The severe effective population size (Ne) declines observed during the last glacial maximum (LGM) and the Middle Pleistocene are highlighted with gray vertical bars.

Characterization of runs of homozygosity and deleterious mutations

We investigated whether M. oleifera showed signs of recent inbreeding by calculating the ROH. Referring to the method of Robinson et al. [10], we used the physical length of ROH to estimate the number of generations to the common ancestor (g) as g = 100/(2 * L), where L is the mean length of ROH in megabases (Mb). Here, the L of all 16 populations of M. oleifera ranged from 0.46 Mb (BM1) to 1.03 Mb (JM) (Supplementary Fig. S11). Our results indicated that inbreeding occurred about 49 to 112 generations ago. Specially, the effects of inbreeding varied greatly among populations (Fig. 3A). We found the FROH was significantly higher in the JM (45.37%–70.95%) population than in other populations, whereas it was lower in LY2 (3.52%–9.82%), BM1 (4.31%–13.92%), FS (7.11%–13.39%), and DX (8.51%–15.60%) populations. Moreover, populations with severe inbreeding would be predicted to have larger numbers of long ROH (>1 Mb) than short ROH (100 Kb–1 Mb) (Fig. 3B). Specifically, the JM population harbored maximum number of ROH >1 Mb, which represented 37.07% of the total genome. In contrast, the BM1 population had a minimum number of long ROH, with only 1.20% of ROH being longer than 1 Mb (Supplementary Table S7).

Figure 3: Levels of inbreeding and genetic load in different M. oleifera populations. (A) Fractions of the runs of homozygosity (FROH) show discrepancies in inbreeding levels in M. oleifera populations. (B) Distributions of long (>1 Mb) and medium (100 kb–1 Mb) runs of homozygosity (ROH) among 16 populations of M. oleifera. Solid dots represent ROH >1 Mb and hollow dots represent 1 Mb > ROH > 100 Kb. (C) Ratios of homozygous-derived deleterious mutations show the discrepancies in genetic load in different M. oleifera populations. Populations marked with the same letters in (A) and (C) are not significantly different. (D) Distributions of the ratio of 0- to 4-fold heterozygosity versus the intergenic heterozygosity across the 16 populations of M. oleifera. The dark line represents the significant negative correlation between these populations, with R = −0.21 and P = 0.006. Each dot represents an individual, which is colored by population.

Based on the modified SIFT prediction approach, we detected a total of 2,404 deleterious mutations, 5,040 tolerated mutations, and 6,172 synonymous mutations (Supplementary Table S8 and Supplementary Fig. S13a). In particular, the frequency of deleterious mutations of homozygous-derived alleles reflects genetic load and adaptability of species, and it varies greatly among populations even within the same species [80]. Our results showed that the number of homozygous deleterious sites to the total deleterious mutations in the JM population was significantly higher than in most other populations except NP and SG, suggesting that the JM population had a higher genetic load (Fig. 3C). Interestingly, the ratio of 0-fold to 4-fold degenerate site heterozygosity showed a significant negative correlation with intergenic (neutral) heterozygosity (Fig. 3D). It implied that more severely deleterious mutations can be effectively purged by purifying selection, which may be the maintaining mechanism of M. oleifera with a small population size [45]. To provide accurate and direct genetic rescue guidance for the JM population, we selected 4 populations as potential pollen donors to construct a Venn diagram of homozygous deleterious variants (Supplementary Fig. S12a). Our results showed that the most homozygous deleterious variants were shared among all the 5 populations, and the JM population had the least shared homozygous deleterious variants with the BM1 population (172 variants).

Signals of genomic offset to future climate change

Potential genomic variants related to climate adaptation were detected using BayeScEnv and RDA. For BayeScEnv analysis, with a q-value cutoff of 0.05, we identified 589 SNPs related to climate adaptation. Of these, 491 SNPs were associated with BIO7 and BIO14, respectively, followed by BIO3 (471 SNPs) and BIO15 (168 SNPs) (Supplementary Table S11). For RDA, 694 SNPs were detected along 5 significant RDA axes, of which 459 SNPs were correlated most to BIO14, 156 SNPs to BIO3, 40 SNPs to BIO7, and 39 SNPs to BIO15. To avoid false positives, we assigned the 380 SNPs detected by both BayeScEnv and RDA as environment-associated SNPs (Supplementary Table S11). To figure out the potential function of genomic variants associated with climate adaptation, we conducted a functional annotation of outlier SNPs. Gene Ontology enrichment analysis assigned a total of 258 Gene Ontology categories (P < 0.05), of which 158 categories belonged to biological processes and abundant genes were associated with metabolism, transmembrane transport, methylation, cell development, flowering, and telomere maintenance (Supplementary Table S12).

To assess which population of M. oleifera will be most likely disrupted in the future (2081–2100) under 2 greenhouse gas scenarios (SSP126 and SSP585), we employed the GF method to investigate the GO using integrated results of BCC-CSM2-MR, CNRM-CM6-1, and CNRM-ESM2-1 climate models. The GO is measured by the Euclidean distance of future climate condition compared to current climate status. Higher GO means greater allele frequency changes are required to adapt to the changing climate [81]. GF modeling showed that the degree of GO of all populations increased under scenario SSP585 compared to scenario SSP126, suggesting that extreme future climate change will cause severe genomic vulnerability to M. oleifera (Fig. 4). Compared with all SNPs (reference), we found that adaptive SNPs (candidate) exhibited higher GO under the same scenarios, implying that adaptive variants were more sensitive to climate change (Fig. 4). In addition, we found a strong negative correlation between GO and altitude, with populations at higher altitudes generally having lower GO values (Fig. 4).

Figure 4: Predicted genetic offset of M. oleifera in the year 2100 under the SSP126 and SSP585 scenarios based on all SNPs (A, B) and adaptive SNPs (C, D), with higher values (red) representing more severe genomic vulnerability to future climate change. The inner mini plot represents the correlation between altitude and GO value in the corresponding scenario.

Ecological niche modeling predicted niche suitability change

We integrated results from 5 models to perform ecological niche modeling for M. oleifera (Supplementary Table S13). The area under the curve (AUC) value, Somer’s D value, TSS value, Boyce value, and the area under the precision–recall gain curve (AUCprg) were about 0.99, 0.99, 0.98, 0.74, and 0.97, respectively, indicating high performance of the ecological niche models (Supplementary Table S14). Compared to the current state, the potential suitable region in 2100 will reduce by 71.15% and 98.79% under scenarios SSP126 and SSP585, respectively, with the threshold of ecological niche suitability equal to 0.56 (Supplementary Fig. S15). Further, we calculated NSC between current and future climate for each grid using the following equation: NSC = niche suitability index in the current climate − niche suitability index in the future climate. A positive value indicates that niche suitability will be decreased under future climate change. Higher positive value means more severe degree of unsuitability. Our results showed slightly higher NSC in the southern part of the distribution range under scenario SSP126 (Fig. 5A). However, under scenario SSP585, the NSC increased to a much higher extent in the northern part of the distribution range (Fig. 5B), which is the Karst basin harboring the Nanpan, Beipan, and Tuoniang rivers. This suggest that drastic climate change will exacerbate ecological vulnerability in karst landforms by affecting hydrological processes [82, 83].

Figure 5: Predicted niche suitability change (NSC) of M. oleifera in the year 2100 under the SSP126 (A) and SSP585 (B) scenarios. Higher positive value suggests more severe degree of unsuitability.

Discussion

The genome harbors valuable evolutionary information of a species and provides deep insights into genetic diversity and evolutionary dynamics, but the full implementation of conservation genomics in practice is still limited [84]. In this study, we conducted a conservation genomics study on M. oleifera based on population-wide genome resequencing data, including 165 individuals. We aim to distinguish the potential factors that affect genetic diversity of M. oleifera, reveal the causes for the formation of its extremely small population patterns, and assess its adaptability under future climate change. It is our hope that based on the comprehensive results of conservation genomics, it will be possible to provide practical and meaningful suggestions for the conservation actions of this ecologically and economically important species.

Recent inbreeding affected genetic diversity

M. oleifera has relatively high genetic diversity among endangered woody plants (Supplementary Table S15). This is confirmed by a population structure result, which shows K = 14 is optimal based on all loci (Supplementary Fig. S8a), suggesting M. oleifera has complex ancestral components despite occupying a narrow distribution. It seems optimistic, but genetic diversity cannot determine the endangered status of a species even though it is an important criterion for species conservation [85]. Genetic diversity is affected by many factors such as inbreeding, gene flow, life form, distribution, and rarity [86]. The key to perform conservation actions for endangered species is to realize the pivotal factor affecting genetic diversity. Our results showed that populations with low genetic diversity (JM, NP, and BM2) have a severe degree of recent inbreeding and displayed more pronounced sawtooth-like distribution patterns of nucleotide diversity across the genome due to long ROH (Fig.   3A, Supplementary Figs. S4 and S5). Furthermore, high levels of inbreeding in JM, NP, and BM2 populations have led to the accumulation of deleterious mutations (Fig. 3C), and they showed greater differentiation from other populations (Supplementary Table S6), which may result in a vicious circle of inbreeding depression without intervention [87].

Causes for the formation of small and isolated population

Historical climate disturbances were one of the reasons for the formation of currently observed small and isolated M. oleifera populations. Based on the estimations of Stairway Plot v.2 and MSMC v.2, we observed a bottleneck event of M. oleifera during the LGM, resulting in a swift decline in Ne (Fig. 2). Although the stairway plot suggested the Ne recovered to its historical peak at the end of the LGM (Fig. 2A), this inference is deemed unreliable of very recent demographic events based on SFS [88]. In contrast, the MSMC results suggested that the Ne of M. oleifera underwent a protracted decline after the LGM, with a slight recovery occurring approximately 500 years ago (Fig. 2B). Moreover, we utilized the mean length of ROH as a metric to estimate the generations of inbreeding, referring to the method of Robinson et al. [10]. Our findings revealed that the JM population had experienced the most recent inbreeding, approximately 490 years ago (Supplementary Fig. S11). Intriguingly, the demographic history inferred by MSMC showed that the Ne of M. oleifera reached its nadir (below 75) about 400–500 years ago (Fig. 2B), which may have contributed to the extensive inbreeding. At the same time, human overexploitation and destruction of wild resources have exerted unbearable demographic pressures and resulted in further population fragmentation, as it is hard to find wild individuals again according to a large number of previous distribution records like Mashan, Pingguo, Tiandong, Tianyang, Youjiang, and Longzhou counties in Guangxi province [34]. Overall, the combined effect of historical bottleneck events, recent inbreeding, and excessive human disturbance may have led to the formation of small and isolated populations of M. oleifera.

Local adaptation-related alleles lead to climate change-driven genomic vulnerability

The climate is currently shifting, and many species face the challenge of keeping pace with ongoing climate changes [89, 90]. Therefore, supporting species able to adapt to the variable climate is a key but tough task to future conservation action. For M. oleifera, we detected a total of 380 SNPs that are related to climate adaptation. The associated genes are also significantly enriched in key processes of metabolism, transmembrane transport, flower development, and telomere maintenance (Supplementary Table S12). Moreover, each climate factor is associated with dozens to hundreds of SNPs accordingly (Supplementary Table S11). This is consistent with the polygenic effects underlying local adaptability, meaning that organisms can adapt to rapid climate change through small polygenic allele frequency shifts [18, 91].

Inequivalent response to climate change exists within populations of the same species due to local adaptation to heterogeneous environments [92]. Analyzing local adaptation pattern can help us better understand how species respond to future climate change. Here, we used the Euclidean distance between future and current climate environments to measure GO. By incorporating intraspecific variations into the predictive GF model, our results showed M. oleifera populations distributed in the low elevation exhibit higher GO under both future scenarios (Fig. 4). Elevation also showed a strong negative correlation with annual mean temperature (R = –0.91, P < 2.2 × 10−16) across the distribution range of M. oleifera (Supplementary Fig. S16). This suggests populations in low altitude have a more significant adaptive lag in response to rapid climate change (especially temperature), indicating a greater risk of local extinction if appropriate conservation measures are not taken.

Ecological niche modeling provides insights into ex situ conservation

Ecological niche modeling can predict potential current distribution ranges and suitable habitats under future climate change by linking observed species distribution and abundance to selected environmental variables [93]. Our results showed the suitable habitats for M. oleifera will be decreased and ecological niche suitability will be further reduced in the future (Supplementary Fig. S15). However, the degree of niche suitability change is varied under different climate scenarios. The extreme climate (SSP585) is likely to have direct impact on hydrological processes in karst landforms, making the northern part across the distribution range that harbors the Nanpan, Beipan, and Tuoniang rivers most unsuitable for living (Fig. 5B). BM2, SG, ZL, and LY1 populations located in the area deserve the highest priority for ex situ conservation when future climate becomes extremely severe. Ecological niche modeling reveals niche suitability change, while GO provides information about genomic inadaptation to future climate change [16]. The 2 methods provide disparate views to estimate climate-driven vulnerability. It is necessary to combine the methods of ecological niche modeling and genomic offset to make conservation decisions.

Implications for conservation guidelines and management strategies

Currently, a wide range of field investigations and in situ protection of existing M. oleifera resources have been implemented [94, 95]. For example, in 2017, Guangnan County labeled 7,941 wild individuals and recorded their growth state [34]. But these measures only have a limited impact on guiding future conservation actions. The conservation guidelines and management strategies should be made under demarcating reasonable management units (MUs) and adaptive units (AUs) [96]. Based on the result of the population structure and phylogenetic tree (neutral loci), we suggest delineating 14 MUs of M. oleifera, with the most single population being separate MUs (Supplementary Fig. S9). Maintaining multiple MUs ensures long-term persistence of the species. Based on adaptive loci, we identified 10 AUs, including JM, SG-NP, ML-BB2, BB1-BM2, GL, ZS, ZL-LY1, LY2-FS, BM1, and DX AUs (Fig. 1B). Different AUs represent varied evolutionary potential. Understanding the patterns of adaptive differentiation is crucial when considering conservation priorities, assisted gene flow, migration, and supplementation [79].

For populations with recent inbreeding, genetic rescue is necessary through assisted gene flow [97]. The JM population has the lowest genetic diversity and the highest inbreeding and genetic load (Fig. 3, Supplementary Fig. S4), and thus it needs urgent genetic rescue. In some cases, a high Fst combination among parents generated offspring with a high heterozygosity, which may produce a heterozygous advantage, but in most cases, outbreeding depression has occurred [98]. It may occur between interpopulation crossing even in the same species [99]. So, to minimize the risk of outbreeding depression, we do not consider the DX population as a potential pollen donor because it has the highest Fst (0.90) with the JM population based on adaptive SNPs (Supplementary Table S6), suggesting severely adaptive differentiation has occurred. Under comprehensive evaluation, we propose the BM1 population as the prior pollen donor for the JM population, because (i) it displays moderate genetic differentiation (Fst = 0.51) with the JM population based on adaptive SNPs (Supplementary Table S6), (ii) it has the highest genetic diversity and heterozygosity (Supplementary Figs. S4 and S6), and (iii) it has the lowest inbreeding level and fewest shared homozygous deleterious mutations with the JM population (Supplementary Fig. S12a), positively reducing the impact of genetic load on hybrid offspring. Moreover, previous research has shown that M. oleifera seeds have low germination rates under natural conditions [24], so it is better to keep the seeds for artificial germination after implementing assisted gene flow measures and introduce robust seedlings back to the natural population later.

GO analysis predicts that populations located at lower altitudes require a greater change in adaptive allele frequencies to adapt to extreme climates (Fig. 4). Low-altitude areas are more susceptible to extreme temperatures than higher elevations (Supplementary Fig. S16). Therefore, we suggest cultivating heat-resistant individuals and screening a preadapted genotype under controlled conditions in the laboratory, and then regression experiments can be conducted. Future work should prioritize the conservation of M. oleifera because its lasting existence is a prerequisite to exploit resources.

Supplementary Material

giae070_GIGA-D-24-00159_Original_Submission

giae070_GIGA-D-24-00159_Revision_1

giae070_Response_to_Reviewer_Comments_Original_Submission

giae070_Reviewer_1_Report_Original_Submission Yongzhi Yang, Ph.D. -- 6/24/2024 Reviewed

giae070_Reviewer_2_Report_Original_Submission Scott Ferguson -- 6/26/2024 Reviewed

giae070_Reviewer_3_Report_Original_Submission Wei Zhao -- 6/30/2024 Reviewed

giae070_Supplemental_Files

Additional Files

Supplementary Note S1. Ancestral sequence reconstruction.

Supplementary Note S2. Estimation of mutation rate.

Supplementary Note S3. Detection of deleterious mutations based on REF-ALT strategy.

Supplementary Fig. S1. Resequencing data processing workflow of Malania oleifera.

Supplementary Fig. S2. The distribution of SNPs (dataset 2) across 13 pseudochromosomes under 1-Mb windows.

Supplementary Fig. S3. Genome-wide linkage disequilibrium (LD) decay of Malania oleifera. (a) Considering 16 populations separately and (b) as a whole.

Supplementary Fig. S4. The comparison of mean θπ and θw among 16 populations of Malania oleifera in whole-genome (a), intergenic (b), CDS (c), intron (d), fold-0 (e), and fold-4 (f) regions.

Supplementary Fig. S5. Distributions of nucleotide diversity (θπ) across the genome with (a–c) representing NP, JM, and BM2 populations (lowest average θπ), respectively, and (d) representing the BM1 population (highest average θπ).

Supplementary Fig. S6. The comparison of heterozygosity rate among 16 populations of Malania oleifera in whole-genome (a), intergenic (b), CDS (c), intron (d), fold-0 (e), and fold-4 (f) regions.

Supplementary Fig. S7. Cross-validation error curve based on all loci (a), neutral loci (b), and adaptive loci (c) for the 16 populations of Malania oleifera inferred by ADMIXTURE.

Supplementary Fig. S8. The inference of population structure (a), principal component analysis (b), and NJ tree (c) of Malania oleifera based on all loci.

Supplementary Fig. S9. The inference of population structure (a), principal component analysis (b), and NJ tree (c) of Malania oleifera based on neutral loci.

Supplementary Fig. S10. Demographic history of Malania oleifera inferred by Stairway Plot v.2 based on folded SFS.

Supplementary Fig. S11. Runs of homozygosity (ROH) frequency differences among 16 populations of Malania oleifera. The dashed lines correspond to the mean ROH lengths for each population.

Supplementary Fig. S12. The Venn diagrams of share and private homozygous deleterious mutations of JM (a), SG (b), and BM2 (c) populations with candidate pollen donors.

Supplementary Fig. S13. Differences in the number of deleterious mutations detected by REF-ALT strategy and ancestral status-based strategy of Malania oleifera (a) and Acer yangbiense (b).

Supplementary Fig. S14. The importance of environmental variables inferred by gradient forest modeling. *Top uncorrelated environment variables (Pearson’s |r| < 0.7) used in BayeScEnv, RDA, and GF analysis.

Supplementary Fig. S15. Integrated results of ecological niche modeling based on 5 models in the current (a) and future SSP126 (b) and SSP585 (c) scenarios. Higher values represent higher suitability.

Supplementary Fig. S16. Correlation of altitude with BIO1, BIO2, BIO4, BIO13, and BIO14 climate variables used in GO analysis.

Supplementary Table S1. Geographical location of all sampled individuals of Malania oleifera.

Supplementary Table S2. Statistical analysis of resequencing data before and after filtering with Fastp.

Supplementary Table S3. Statistical analysis of each individual mapping to the Malania oleifera reference genome.

Supplementary Table S4. Statistics on Watterson’s θ (θw) of Malania oleifera populations within whole-genome, CDS, fold-0, fold-4, intergenic, and intron regions.

Supplementary Table S5. Statistics of heterozygosity and homozygosity rates of all individuals used in the whole-genome resequencing.

Supplementary Table S6. Weighted fixation statistics (Fst) between populations of all loci, adaptive loci, and neutral loci and geographical distances between populations of Malania oleifera.

Supplementary Table S7. Estimation of runs of homozygosity (ROH) and frequency of runs of homozygosity (FROH).

Supplementary Table S8. Summary of deleterious, tolerated, and synonymous mutations of Malania oleifere based on REF-ALT strategy.

Supplementary Table S9. A list of 19 bioclimatic variables used in this study.

Supplementary Table S10. The projection coordinates of M. oleifera distribution records for ecological niche modeling.

Supplementary Table S11. Environment-associated SNPs and corresponding genes detected by BayeScEnv and redundancy analysis (RDA).

Supplementary Table S12. Gene ontology enrichment analysis of environment-associated genetic variants.

Supplementary Table S13. Predicted GO and NSC values of populations under future SSP126 and SSP585 scenarios using BCC-CSM2-MR, CNRM-CM6-1, and CNRM-ESM2-1 climate models.

Supplementary Table S14. The evaluation of ecological niche models.

Supplementary Table S15. Comparison of nucleic acid diversity of endangered species.

Supplementary Table S16. Information of 17 individuals for ancestral sequence reconstruction.

Supplementary Table S17. The genome list of 17 published species used in our research to estimate mutation rate for Malania oleifera.

Abbreviations

a.s.l: above sea level; AUs: adaptive units; CDS: coding sequences; ENM: ecological niche model; FDR: false discovery rate; FROH: frequency of runs of homozygosity; Fst: pairwise fixation statistics; GF: gradient forest; GO: genomic offset; IUCN: International Union for Conservation of Nature; LD: linkage disequilibrium; LGM: last glacial maximum; MAF: minor allele frequency; MUs: management units; Ne: effective population size; NJ: neighbor-joining; NSC: niche suitability change; PCA: principal component analysis; RDA: redundancy analysis; ROH: run of homozygosity; SFS: site frequency spectrum; SIFT: sorting intolerant from tolerant; SNP: single nucleotide polymorphism; SSP: Shared Socioeconomic Pathways; θw: Watterson’s θ; θπ: nucleotide diversity.

Competing Interests

The authors declare that they have no competing interests.

Ethics Statement

All plant molecular materials and specimens were collected with permission.

Author Contributions

Yongpeng Ma (Conceptualized, Supervised, Reviewed and edited the manuscript). Weibang Sun (Conceptualized, Supervised, Reviewed and edited the manuscript). Yuanting Shen (Data analysis, Interpretation, Visualization, Wrote the manuscript). Lidan Tao (Data analysis, Interpretation, Visualization). Rengang Zhang (Data analysis, Interpretation, Visualization). Gang Yao (Investigated). Minjie Zhou (Visualization). All authors reviewed and approved the final manuscript.

Funding

This work was supported by the Key Project of Natural Science Foundation of Yunnan Province (Grant No. 202001AS070019), the CAS “Light of West China” Program, the Ten Thousand Talent Program of Yunnan Province(Grant No. YNWRQNBJ-2018-174) and the Yunnan Provincial Science and Technology Mission (Grant No. 202404BI090014).

Data Availability

Raw resequencing data are available at the NCBI Sequence Read Archive under BioProject PRJNA978997. All additional supporting data are available in the GigaScience repository, GigaDB [100].
==== Refs
References

1. Miraldo A , LiS, Borregaard​​​MK, et al. An Anthropocene map of genetic diversity. Science. 2016;353 (6307 ):1532–35.. 10.1126/science.aaf4381.27708102
2. Lynch M , ConeryJ, BurgerR. Mutation accumulation and the extinction of small populations. Am Soc Naturalists. 1995;146 (4 ):489–518. 10.1086/285812.
3. Charlesworth D , WillisJH. The genetics of inbreeding depression. Nat Rev Genet. 2009;10 (11 ):783–96. 10.1038/nrg2664.19834483
4. Xue Y , Prado-MartinezJ, SudmantPH, et al. Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding. Science. 2015;348 (6231 ):242–45. 10.1126/science.aaa3952.25859046
5. Feng S , FangQ, BarnettR, et al. The genomic footprints of the fall and recovery of the crested ibis. Curr Biol. 2019;29 (2 ):340–49. e7. 10.1016/j.cub.2018.12.008.30639104
6. Keller MC , VisscherPM, GoddardME. Quantification of inbreeding due to distant ancestors and its detection using dense single nucleotide polymorphism data. Genetics. 2011;189 (1 ):237–49. 10.1534/genetics.111.130922.21705750
7. Ma Y , LiuD, WarissHM, et al. Demographic history and identification of threats revealed by population genomic analysis provide insights into conservation for an endangered maple. Mol Ecol. 2022;31 (3 ):767–79. 10.1111/mec.16289.34826164
8. Hedrick PW , Garcia-DoradoA. Understanding inbreeding depression, purging, and genetic rescue. Trends Ecol Evol. 2016;31 (12 ):940–52. 10.1016/j.tree.2016.09.005.27743611
9. Caballero A , BravoI, WangJ. Inbreeding load and purging: implications for the short-term survival and the conservation management of small populations. Heredity (Edinb). 2017;118 (2 ):177–85. 10.1038/hdy.2016.80.27624114
10. Robinson JA , KyriazisCC, Nigenda-MoralesSF, et al. The critically endangered vaquita is not doomed to extinction by inbreeding depression. Science. 2022;376 (6593 ):635–39. 10.1126/science.abm1742.35511971
11. Malcolm JR , LiuC, NeilsonRP, et al. Global warming and extinctions of endemic species from biodiversity hotspots. Conserv Biol. 2006;20 (2 ):538–48. 10.1111/j.1523-1739.2006.00364.x.16903114
12. Wiens JJ . Climate-related local extinctions are already widespread among plant and animal species. PLoS Biol. 2016;14 (12 ):e2001104. 10.1371/journal.pbio.2001104.27930674
13. Bay RA , HarriganRJ, UnderwoodVL, et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science. 2018;359 (6371 ):83–86. 10.1126/science.aan4380.29302012
14. Aitken SN , YeamanS, HollidayJA, et al. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol Appl. 2008;1 (1 ):95–111. 10.1111/j.1752-4571.2007.00013.x.25567494
15. Aguirre-Liguori JA , Ramirez-BarahonaS, GautBS. The evolutionary genomics of species' responses to climate change. Nat Ecol Evol. 2021;5 (10 ):1350–60. 10.1038/s41559-021-01526-9.34373621
16. Chen Y , JiangZ, FanP, et al. The combination of genomic offset and niche modelling provides insights into climate change-driven vulnerability. Nat Commun. 2022;13 (1 ):4821. 10.1038/s41467-022-32546-z.35974023
17. Jia KH , ZhaoW, MaierPA, et al. Landscape genomics predicts climate change-related genetic offset for the widespread Platycladus orientalis (Cupressaceae). Evol Appl. 2020;13 (4 ):665–76. 10.1111/eva.12891.32211059
18. Sang Y , LongZ, DanX, et al. Genomic insights into local adaptation and future climate-induced vulnerability of a keystone forest tree in East Asia. Nat Commun. 2022;13 (1 ):6541. 10.1038/s41467-022-34206-8.36319648
19. Yang H , LiJ, MilneRI, et al. Genomic insights into the genotype-environment mismatch and conservation units of a Qinghai-Tibet Plateau endemic cypress under climate change. Evol Appl. 2022;15 (6 ):919–33. 10.1111/eva.13377.35782009
20. Li SG . Malania, a new genus of oil-yielding plant. Bull Bot Lab N E Forest Inst. 1980;1 :67–72.. https://bbr.nefu.edu.cn/CN/Y1980/V0/I1/67.
21. Lv SH , WeiCQ, HuangFZ, et al. Fruit and seed traits and adaptability to rocky desertification mountain of rare tree species Malania oleifera. Chin J Ecol. 2016;35 (1 ):57–62. 10.13292/1.1000-4890.201601.008
22. Yang T , YuQ, XuW, et al. Transcriptome analysis reveals crucial genes involved in the biosynthesis of nervonic acid in woody Malania oleifera oilseeds. BMC Plant Biol. 2018;18 (1 ):247. 10.1186/s12870-018-1463-6.30340521
23. Lu SG , LeiLB, YangQS, et al. The current status and the cause of the endangerment of Malania oleifera Chun et Lee in southeast Yunnan. In: ChenYY, ed. Biodiversity conservation and regional sustainable development: the 4th Biodiversity Conservation and Sustainable Use Conference. Beijing: China Forestry Press; 2000:169–72.. https://www.doc88.com/p-0478327154099.html.
24. Xu DB , ChenF, GuoXC, et al. Research on the bottleneck of resource protection and industrialization development of rarely endangered Malania oleifera. Issues Forest Econ. 2018;38 (3 ):13–20. 10.16832/j.cnki.1005-9709.2018.03.003
25. Wu YQ , LiXD, HuYJ. Reproductive biology of Malania oleifera. Acta Sci Nat Univ Sunyatseni. 2004;43 (2 ):81–83.
26. Lai JY , ShiHM, PanCL, et al. Pollination biology of rare and endangered species Malania oleifera Chun et Lee. J Beijing Forest Univ. 2008;30 (2 ):59–64. 10.13332/j.1000-1522.2008.02.021
27. Li XD . Life-table analysis of Malania oleifera, a rare and endangered plant. J Central South Univ Forest Technol. 2009;29 (2 ):73–76.
28. Xu SS , KanW, KongBH, et al. First report of fusarium oxysporum and fusarium solani causing root rot on Malania oleifera in China. Plant Dis. 2020;104 (2 ):584–84. 10.1094/PDIS-07-19-1426-PDN.
29. Fu LG . Red data book of Chinese plant: the rare and endangered plants. Beijing: Science Press; 1992.
30. Ma Y , ChenG, Edward GrumbineR, et al. Conserving plant species with extremely small populations (PSESP) in China. Biodivers Conserv. 2013;22 (3 ):803–9. 10.1007/s10531-013-0434-3.
31. Yang T , ZhangR, TianX, et al. The chromosome-level genome assembly and genes involved in biosynthesis of nervonic acid of Malania oleifera. Sci Data. 2023;10 (1 ):298. 10.1038/s41597-023-02218-8.37208438
32. Chen W , WangP, PuT, et al. Symbiotic effect of co-cultivated plants on Malania oleifera seedlings. Acta Agric Univ Jiangxiensis. 2022;44 (5 ):1197–206. 10.13836/j.jjau.2022119
33. Chen Q , LiY, LiY, et al. Dynamics of tissue nutrient content in relation to declining seedling growth in Malania oleifera. Guihaia. 2024;44 (1 ):137–46. 10.11931/guihaia.gxzw202303048
34. Su C , WangG, GaoY, et al. Resource protection and development counterplants of Malania oleifera. J Anhui Agric Sci. 2023;51 (12 ):104–7. 10.3969/j.issn.0517-6611.2023.12.024
35. Supple MA , ShapiroB. Conservation of biodiversity in the genomics era. Genome Biol. 2018;19 (131 ):1–12. 10.1186/s13059-018-1520-3.29301551
36. Doyle JJ , DoyleJL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987;19 (1 ):11–15. 10.1016/0031-9422(80)85004-7.
37. Chen S , ZhouY, ChenY, et al. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34 (17 ):i884–i90. 10.1093/bioinformatics/bty560.30423086
38. Li H . Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv. 10.48550/arXiv.1303.3997. Accessed May 26 2013.
39. Danecek P , BonfieldJK, LiddleJ, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10 (2 ):giab008. 10.1093/gigascience/giab008.33590861
40. Tarasov A , VilellaAJ, CuppenE, et al. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31 (12 ):2032–34. 10.1093/bioinformatics/btv098.25697820
41. Garrison E , MarthG. Haplotype-based variant detection from short-read sequencing. arXiv. 10.48550/arXiv.1207.3907 Accessed 20 July 2012.
42. Danecek P , AutonA, AbecasisG, et al. The variant call format and VCFtools. Bioinformatics. 2011;27 (15 ):2156–58. 10.1093/bioinformatics/btr330.21653522
43. Zhang C , DongSS, XuJY, et al. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics. 2019;35 (10 ):1786–88. 10.1093/bioinformatics/bty875.30321304
44. Korneliussen TS , AlbrechtsenA, NielsenR. ANGSD: analysis of next generation sequencing data. BMC Bioinf. 2014;15 (1 ):356. 10.1186/s12859-014-0356-4.
45. Yang Y , MaT, WangZ, et al. Genomic effects of population collapse in a critically endangered ironwood tree Ostrya rehderiana. Nat Commun. 2018;9 (1 ):5449. 10.1038/s41467-018-07913-4.30575743
46. Frichot E , FrançoisO, O'MearaB. LEA: an R package for landscape and ecological association studies. Methods Ecol Evol. 2015;6 (8 ):925–29. 10.1111/2041-210x.12382.
47. Luu K , BazinE, BlumMG. pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol Ecol Resour. 2017;17 (1 ):67–77. 10.1111/1755-0998.12592.27601374
48. Chang C , ChowC, TellierL, et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4 (7 ):1–12. 10.1186/s13742-015-0047-8.25838885
49. Alexander DH , NovembreJ, LangeK. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19 (9 ):1655–64. 10.1101/gr.094052.109.19648217
50. Yang J , LeeSH, GoddardME, et al. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88 (1 ):76–82. 10.1016/j.ajhg.2010.11.011.21167468
51. Kumar S , StecherG, TamuraK. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33 (7 ):1870–74. 10.1093/molbev/msw054.27004904
52. Liu X , FuYX. Stairway Plot 2: demographic history inference with folded SNP frequency spectra. Genome Biol. 2020;21 (1 ):280. 10.1186/s13059-020-02196-9.33203475
53. Schiffels S , DurbinR. Inferring human population size and separation history from multiple genome sequences. Nat Genet. 2014;46 (8 ):919–25. 10.1038/ng.3015.24952747
54. Sim NL , KumarP, HuJ, et al. SIFT web server: predicting effects of amino acid substitutions on proteins. Nucleic Acids Res. 2012;40 (Web Server issue ):W452–57. 10.1093/nar/gks539.22689647
55. Boeckmann B , BairochA, ApweilerR, et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res. 2003;31 (1 ):365–70. 10.1093/nar/gkg095.12520024
56. Vaser R , AdusumalliS, LengSN, et al. SIFT missense predictions for genomes. Nat Protoc. 2016;11 (1 ):1–9. 10.1038/nprot.2015.123.26633127
57. Ellis N , SmithSJ, PitcherCR. Gradient forests: calculating importance gradients on physical predictors. Ecology. 2012;93 (1 ):156–68. 10.1890/11-0252.1.22486096
58. Friendly M . Corrgrams: exploratory displays for correlation matrices. Am Stat. 2002;56 (4 ):316–24. 10.1198/000313002533.
59. Villemereuil P , GaggiottiOE. A new FST-based method to uncover local adaptation using environmental variables. Methods Ecol Evol. 2015;6 (11 ):1248–58. 10.1111/2041-210x.12418.
60. Legendre P , LegendreL. Numerical ecology. Amsterdam: Elsevier; 2012.;
61. Lischer HE , ExcoffierL. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012;28 (2 ):298–99. 10.1093/bioinformatics/btr642.22110245
62. Forester BR , LaskyJR, WagnerHH, et al. Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. Mol Ecol. 2018;27 (9 ):2215–33. 10.1111/mec.14584.29633402
63. Oksanen J , BlanchetFG, KindtR, et al. Package ‘vegan’: community ccology package. 2013. R package version 2.3-0. https://cran.r-project.org/web/packages/vegan/index.html. Accessed 10 Jul 2023.
64. Cantalapiedra CP , Hernandez-PlazaA, LetunicI, et al. eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale. Mol Biol Evol. 2021;38 (12 ):5825–29. 10.1093/molbev/msab293.34597405
65. Yu X , DaiM, PuT, et al. Population structure and dynamics analysis of rare and endangered plant Malania oleifera. J West China Forest Sci. 2023;52 (3 ):8–16. 10.16473/j.cnki.xblykx1972.2023.03.002
66. Gong MJ , WangJ, FuXY, et al. Suitable regions forecasting and environmental influencing factors of Malania oleifera in Yunnan and Guangxi. J Nanjing Forest Univ. 2022;46 (2 ):44–52. 10.12302/j.issn.1000-2006.202109039
67. Brown JL , CarnavalAC. A tale of two niches: methods, concepts, and evolution. Front Biogeography. 2019;11 (4 ):e44158. 10.21425/f5fbg44158.
68. Naimi B , HammNAS, GroenTA, et al. Where is positional uncertainty a problem for species distribution modelling?. Ecography. 2013;37 (2 ):191–203. 10.1111/j.1600-0587.2013.00205.x.
69. Wood SN . Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J R Statist Soc B. 2011;73 :3–36. 10.1111/j.1467-9868.2010.00749.x.
70. Hijmans RJ , PhillipsS, LeathwickJ, et al. dismo: species distribution modelling. 2023. R package version 1.3-14. https://CRAN.R-project.org/package=dismo. Accessed 1 Apr 2024.
71. Liaw A , WienerM. Classification and regression by randomForest. R news. 2002;2 :18–22.
72. Friedman J , TibshiraniR, HastieT. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33 (1 ):1–22. 10.18637/jss.v033.i01.20808728
73. Valavi R , Guillera-ArroitaG, Lahoz-MonfortJJ, et al. Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecol Monogr. 2022;92 :e01486. 10.1002/ecm.1486.
74. Wilfried T , DamienG, MayaG, et al. biomod2: ensemble platform for species distribution modeling. 2023. R package version 4.2-4. https://CRAN.R-project.org/package=biomod2. Accessed 12 Dec 2023.
75. Di Cola V , BroennimannO, PetitpierreB, et al. ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography. 2017;40 (6 ):774–87. 10.1111/ecog.02671.
76. Flach P , KullM. prg: creates the Precision-Recall-gain curve and calculates the area under the curve. 2023. R package version 0.5.1. https://github.com/meeliskull/prg. Accessed 19 Feb 2024.
77. Freeman E . PresenceAbsence: presence-absence model evaluation. 2023. R package version 1.1.11. https://CRAN.R-project.org/package=PresenceAbsence. Accessed 22 Feb 2024.
78. Guzmán S , GiudicelliGC, TurchettoC, et al. Neutral and outlier single nucleotide polymorphisms disentangle the evolutionary history of a coastal Solanaceae species. Mol Ecol. 2022;31 (10 ):2847–64. 10.1111/mec.16441.35332594
79. Funk WC , McKayJK, HohenlohePA, et al. Harnessing genomics for delineating conservation units. Trends Ecol Evol. 2012;27 (9 ):489–96. 10.1016/j.tree.2012.05.012.22727017
80. Hu JY , HaoZQ, FrantzL, et al. Genomic consequences of population decline in critically endangered pangolins and their demographic histories. Natl Sci Rev. 2020;7 (4 ):798–814. 10.1093/nsr/nwaa031.34692098
81. Fitzpatrick MC , KellerSR. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol Lett. 2015;18 (1 ):1–16. 10.1111/ele.12376.25270536
82. Rao W , ShenZ, DuanX. Spatiotemporal patterns and drivers of soil erosion in Yunnan, Southwest China: rulse assessments for recent 30 years and future predictions based on CMIP6. Catena. 2023;220 :106703. 10.1016/j.catena.2022.106703.
83. Liu H , ZhangM, LinZ, et al. Spatial heterogeneity of the relationship between vegetation dynamics and climate change and their driving forces at multiple time scales in Southwest China. Agric For Meteorol. 2018;256–257 :10–21. 10.1016/j.agrformet.2018.02.015.
84. Theissinger K , FernandesC, FormentiG, et al. How genomics can help biodiversity conservation. Trends Genet. 2023;39 (7 ):545–59. 10.1016/j.tig.2023.01.005.36801111
85. Kahilainen A , PuurtinenM, KotiahoJS. Conservation implications of species–genetic diversity correlations. Global Ecol Conserv. 2014;2 :315–23. 10.1016/j.gecco.2014.10.013.
86. He ZZ , StotzGC, LiuX, et al. A global synthesis of the patterns of genetic diversity in endangered and invasive plants. Biol Conserv. 2024;291 :110473. 10.1016/j.biocon.2024.110473.
87. Ellegren H , GaltierN. Determinants of genetic diversity. Nat Rev Genet. 2016;17 (7 ):422–33. 10.1038/nrg.2016.58.27265362
88. Chen H , HeyJ, ChenK. Inferring very recent population growth rate from population-scale sequencing data: using a large-sample coalescent estimator. Mol Biol Evol. 2015;32 (11 ):2996–3011. 10.1093/molbev/msv158.26187437
89. Dawson TP , JacksonST, HouseJI, et al. Beyond predictions: biodiversity conservation in a changing climate. Science. 2011;332 (6025 ):53–58. 10.1126/science.1200303.21454781
90. Pacifici M , FodenWB, ViscontiP, et al. Assessing species vulnerability to climate change. Nat Clim Change. 2015;5 (3 ):215–24. 10.1038/nclimate2448.
91. Fagny M , AusterlitzF. Polygenic adaptation: integrating population genetics and gene regulatory networks. Trends Genet. 2021;37 (7 ):631–38. 10.1016/j.tig.2021.03.005.33892958
92. Yuan S , ShiY, ZhouBF, et al. Genomic vulnerability to climate change in Quercus acutissima, a dominant tree species in East Asian deciduous forests. Mol Ecol. 2023;32 (7 ):1639–55. 10.1111/mec.16843.36626136
93. Thuiller W . Ecological niche modelling. Curr Biol. 2024;34 (6 ):R225–29. 10.1016/j.cub.2024.02.018.38531309
94. Jia D , MaoJ, ChenF, et al. Investigation and analysis of wild garlic fruit resources in Guangnan. Forest by-product and Speciality in China. 2017;3 :72–76. 10.13268/j.cnki.fbsic.2017.03.032.
95. Liu Y , NingS. Status and evaluation of natural resources of emphasis protective wilding plant in Guangxi. Guangxi Sci. 2002;9 (2 ):124–32. 10.13656/j.cnki,gxkx.2002.02.012.
96. Barbosa S , MestreF, WhiteTA, et al. Integrative approaches to guide conservation decisions: using genomics to define conservation units and functional corridors. Mol Ecol. 2018;27 (17 ):3452–65. 10.1111/mec.14806.30030869
97. Pavlova A , BeheregarayLB, ColemanR, et al. Severe consequences of habitat fragmentation on genetic diversity of an endangered Australian freshwater fish: a call for assisted gene flow. Evol Appl. 2017;10 (6 ):531–50. 10.1111/eva.12484.28616062
98. Frankham R , BallouJD, EldridgeMD, et al. Predicting the probability of outbreeding depression. Conserv Biol. 2011;25 (3 ):465–75. 10.1111/j.1523-1739.2011.01662.x.21486369
99. Severns PM . Precautionary hand pollination suggests outbreeding depression between potential seed donor populations for a rare wetland plant. J Torrey Bot Soc. 2013;140 (1 ):20–25. 10.3159/TORREY-D-12-00046.1.
100. Shen YT , TaoLD, ZhangRG, et al. Supporting data for “Genomic Insights into Endangerment and Conservation of the Garlic-Fruit Tree (Malania oleifera), a Plant Species with Extremely Small Populations.” GigaScience Database. 2024. 10.5524/102555.
