Genetic structure of an endangered species Ormosia henryi in southern China, and implications for conservation
BMC Plant Biology volume 23, Article number: 220 (2023)
The evergreen broadleaved forest (EBLF) is an iconic vegetation type of East Asia, and it contributes fundamentally to biodiversity-based ecosystem functioning and services. However, the native habitat of EBLFs keeps on decreasing due to anthropogenic activities. Ormosia henryi is a valuable rare woody species in EBLFs that is particularly sensitive to habitat loss. In this study, ten natural populations of O. henryi in southern China were sampled, and then genotyping by sequencing (GBS) was applied to elucidate the standing genetic variation and population structure of this endangered species.
In ten O. henryi populations, 64,158 high-quality SNPs were generated by GBS. Based on these markers, a relatively low level of genetic diversity was found with the expected heterozygosity (He) ranging from 0.2371 to 0.2901. Pairwise FST between populations varied from 0.0213 to 0.1652, indicating a moderate level of genetic differentiation. However, contemporary gene flow between populations were rare. Assignment test and principal component analysis (PCA) both supported that O. henryi populations in southern China could be divided into four genetic groups, and prominent genetic admixture was found in those populations located in southern Jiangxi Province. Mantel tests and multiple matrix regression with randomization (MMRR) analyses suggested that isolation by distance (IBD) could be the possible reason for describing the current population genetic structure. In addition, the effective population size (Ne) of O. henryi was extremely small, and showed a continuous declining trend since the Last Glacial Period.
Our results indicate that the endangered status of O. henryi is seriously underestimated. Artificial conservation measures should be applied as soon as possible to prevent O. henryi from the fate of extinction. Further studies are needed to elucidate the mechanism that leading to the continuous loss of genetic diversity in O. henryi and help to develop a better conservation strategy.
The evergreen broadleaved forest (EBLF) is the typical vegetation biome in the subtropical regions of China with a range from 22 to 34°N latitude and 99 to 123°E longitude. Considering that the other areas in the same latitude were covered by desert or semi-desert, this biome type is unique owing to its dependence on East Asian monsoon . As this region was not covered by large ice sheets during the Last Glacial Maximum (LGM), many rare or relict plant species have been preserved, making it a hotspot for biodiversity conservation and research [2, 3]. However, natural EBLFs have greatly diminished due to anthropogenic deforestation and the unpredictable extreme climates in recent years seem to increase the potential risk of biodiversity loss. The degradation of habitats caused by these factors has led to the endangerment of many rare species in EBLFs . Although the evergreen broadleaved forests are predicted to expand in the next 100 years according to the representative concentration pathways (RCPs) scenarios under the framework of the Coupled Model Intercomparison Project Phase 5 (CMIP5) , the fate of these rare and endangered plant species in EBLFs is still under debate. As an important component of EBLFs, rare and endangered plants are more vulnerable to population loss due to their small population size and low genetic diversity . Therefore, there is an urgent need for a robust assessment of the standing genetic variation and evolutionary potential of the endangered EBLF components to help develop better conservation strategies, especially for those species with high utility value.
Single nucleotide polymorphisms (SNPs) are the most abundant and universal variations across the genome, thus making them very useful markers for population genetic analysis [7, 8]. With the rapid development in next generation sequencing technologies, a variety of reduced-representation genome sequencing (RRGS) methods have been developed, making it possible for non-model species to generate thousands of genome-wide SNP markers without a reference genome . Among these RRGS methods, genotyping by sequencing (GBS) is a flexible, rapid and cost-effective genomic approach combining both marker discovery and genotyping and has been widely used to investigate genetic diversity and population structure of non-model species [10,11,12]. GBS better detect subtle changes in population structure than traditional markers do, which allows a more comprehensive description of the patterns of genetic variation in non-model species [10,11,12].
The genus Ormosia Jacks. (Fabaceae, Papilionoideae) comprises about 150 species of trees and exhibits a striking disjunct geographical distribution between the New World- and Asian and Australasian wet tropics and subtropics . It is a moderately large genus of particular importance in tropical or subtropical broadleaved forests, with most species generally growing as evergreen (or rarely deciduous) trees [14, 15]. China is the northernmost distribution of the genus Ormosia, and there are approximately 35 Ormosia species distributed in southern China. Among these species, O. henryi is a rare and endangered evergreen tree species with multiple usages. It is highly prized not only for the durable and attractive wood, but also as an excellent ornamental plant [16, 17]. O. henryi was recently found to be a new bioactive resource with potential antidepressant activity . However, during the past several decades, the natural populations of O. henryi in China have dramatically declined due to over-logging, habitat fragmentation and difficulties in regeneration. O. henryi has been classified as a national second-class protected wild plant in China and was evaluated as Vulnerable (VU) by IUCN. Previous studies of O. henryi focused on its bioactive ingredients, medical usages and tissue culture [18,19,20]. By far, we still know little about its genetic variation, population structure or demographical history, which would lead to poor efficiency in conservation, restoration and utilization. Given the lack of a reference genome and molecular markers for O. henryi, genotyping-by-sequencing (GBS) method was adopted in this study. We identified genome-wide SNPs from ten natural populations of O. henryi in southern China to address the following questions: (1) to assess the genetic diversity and population structure of O. henryi, (2) to evaluate the impact of geography and environment on genetic differentiation of O. henryi and (3) to infer the demographic history of O. henryi. All these results will provide a molecular basis for the conservation of O. henryi natural populations and provide a technical reference for their rational utilization.
Reference survey genome and SNP discovery
After removing reads with adapters and reads with low quality, a total of 190,771,231 pairs of reads were obtained, reaching 57.23 Gb of clean data. A total of 4,066,523 contigs were assembled by SOAPdenovo software, and the maximum and N50 contig lengths were 54,197 bp and 656 bp, respectively (Additional file 1: Table S2). The GC content of contigs was 35.27%. The Further assembly generated 3,866,097 scaffolds, and the maximum and N50 scaffold lengths were 59,305 bp and 762 bp, respectively. The assembled survey genome of O. henryi was 1,364,689,951 bp in length (including N) and served as the reference genome for GBS data analysis.
A total of 1764 million reads, 264.67 Gb clean data were generated by GBS. The mapping rates of all samples ranged from 34.20% to 41.55%, and the effective sequencing depth ranged from 4.19 to 7.97 × (Additional file 1: Table S3). With our parameter settings, 331,206 SNPs were identified from 129 samples by bcftools. After filtering for minor allele frequency (MAF ≥ 0.05), missing rate (≤ 0.2), and read depth (DP ≥ 4), 64,158 high-quality SNPs were obtained for subsequent data analysis. In the outlier detection analysis, only six SNPs were identified as putative selected sites by BayeScan with a q-value threshold of 0.05, while no outlier SNPs were detected with a q-value threshold of 0.01 (Additional file 2: Fig. S1). Therefore, most SNPs were putatively neutral and used for population genetic analyses after filtering.
Genetic diversity, population differentiation and contemporary gene flow
At the species level, the observed heterozygosity (Ho), expected heterozygosity (He), and nucleotide diversity (π) were 0.1460, 0.1594 and 0.1132, respectively, suggesting a low level of genetic diversity in O. henryi populations. At the population level, genetic diversity based on the three statistics showed different trends of ranking. The Ho values ranged from 0.2282 (XF) to 0.2874 (SX) with an average of 0.2675, while the He values ranged from 0.2371 (XF) to 0.2901 (AY) with an average of 0.2669. The nucleotide diversity (π) varied between 0.1216 (YK) and 0.1429 (JH) with an average of 0.1312. The results showed that PIC values varied from 0.2023 to 0.2335 with an average of 0.2213. The inbreeding coefficient (FIS) values among populations were close to zero, ranging from -0.0232 to 0.0224 (Table 1).
Pairwise FST values between O. henryi populations varied from 0.0213 to 0.1652, with an average value of 0.0918, suggesting a moderate level of genetic differentiation among these O. henryi populations (Table 2). The geographically most distant JH and YK populations had the highest FST value, while the geographically closest AY and XF populations had the lowest FST value.
Contemporary gene flow analysis detected significant migration rates from XF to AY, from MTS to NH, and from SC to SX (Table 3). All other migration rates were very low (< 0.05) and not significant from zero, suggesting that overall contemporary gene flow between O. henryi populations was limited.
Population genetic structure
Based on the putative neutral SNP dataset, ADMIXTURE identified K = 4 as the most appropriate number according to the CV error (Additional file 2: Fig. S2). Under K = 2, the 129 O. henryi individuals from 10 populations were clustered into three groups. The first group was comprised of all the individuals from JH, the second group was comprised of all individuals from XF and AY, and the third group consisted of the remaining individuals. Under K = 3 and K = 4, all 129 individuals were divided into four groups, the first and second group were maintained unchanged, but the third group was split into one large group with all the individuals from NH, SC, MTS, SX and JO, and one small group with all the individuals from LQ and YK. Regardless of the K values, clear genetic admixtures were observed in AY and XF populations from southern Jiangxi province (Fig. 1).
To better understand the population structure of O. henryi, PCA and phylogenetic tree construction were performed. PCA results supported the division of four genetic groups by ADMIXTURE analysis under K = 3 or 4. Two principal components PC1 and PC2 accounted for 13.61% and 9.40% of total variability, respectively. At PC1, individuals from JH were clearly separated from the other individuals, and individuals from AY and XF were slightly separated from the remaining individuals. At PC2, individuals from NH, SC, MTS, SX and JO were mixed but distinct from individuals from LQ and YK (Fig. 2). The unrooted NJ tree of 129 O. henryi individuals also supported the division of four genetic groups. There were clear boundaries among the four genetic groups, as individuals from the same genetic group clustered together on the NJ tree (Fig. 3). AMOVA analysis revealed that 89.92% of the total variation was found within populations, and only 4.22% and 5.86% were found among groups and among populations within groups, respectively (Table 4).
Effects of IBD and IBE on genetic differentiation
Mantel tests showed a significant correlation between pairwise FST/(1-FST) and geographic distance for the 10 O. henryi populations (r = 0.7985, p = 0.001; Fig. 4), but no significant correlation was detected either between pairwise FST/(1-FST) and environmental distance (r = 0.2742, p = 0.157) or between geographic distance and environmental distance (r = 0.1511, p = 0.271). The results of partial Mantel tests were similar. Pairwise FST/(1-FST) displayed a significant correlation with geographic distance after controlling for environmental distance (r = 0.7095, p = 0.001), but no significant correlation was observed between pairwise FST/(1-FST) and environmental distance after controlling for geographical distance (r = 0.2177, p = 0.177). The multiple matrix regression with randomization (MMRR) analysis revealed the relative contributions of IBD and IBE, and further confirmed strong and significant effects of IBD (βD = 0.5698, p = 0.0001) but weak and non-significant effects of IBE (βE = 0.0587, p = 0.3120) on population genetic differentiation.
Demographic history of O. henryi
We estimated a contemporary effective population size (Ne) of 3.4 (95% CI: 2.6—13.6) for all the O. henryi individuals. The Ne estimates of JH, XF and AY populations were larger than 40, while the Ne estimates of other populations were below 20 (Table 1). Demographic analysis by Stairway Plot 2 showed that O. henryi experienced a population expansion of about 0.2–0.6 Mya and began to decline during the Last Glacial Period. About 50, 000 years ago, the effective population size of O. henryi declined sharply, and no recovery trend of Ne was observed (Fig. 5, Additional file 2: Fig. S3).
Genetic diversity and differentiation
Generally, species or populations with higher genetic diversity are more capable of adapting to environmental changes [21, 22]. In this study, a relatively low level of genetic diversity was found in O. henryi at both species and population levels. The observed heterozygosity (Ho) and expected heterozygosity (He) of O. henryi were 0.1460 and 0.1594 at the species level respectively, which were lower than those of some other tree species revealed by SNP markers, for example, Cinnamomum camphora (Ho = 0.319, He = 0.322, ), Quercus acutissima (Ho = 0.2319, He = 0.2272, ) and Eucalyptus regnans (Ho = 0.22, He = 0.232, ). Notably, the genetic diversity of O. henryi (nSSR, Ho = 0.340, He = 0.492, unpublished data) is much lower than that of O. hosiei (nSSR, Ho: 0.560–0.712, He: 0.723–0.835, ), which is also classified as a national second-class protected wild plant in China. As one of the most northerly distributed and cold-tolerant species in Ormosia, O. hosiei occurs in our sampling area as well. O. hosiei was evaluated as Endangered (EN) by IUCN, which is supposed to be rarer than O. henryi (Vulnerable, VU). The endangered status of O. hosiei with high genetic diversity is likely due to the recent population declines caused by human activities, and the high level of genetic diversity observed in O. hosiei could be a transitory state [17, 26, 27]. Our demographic analysis indicated that the effective population size (Ne) of O. henry began to decline during the Last Glacial Period (LGP) and the contemporary Ne estimate was very low at the species level. Thus, the relatively low genetic diversity of O. henryi appears to be a result of both evolutionary history and human-induced population declines.
Our results showed that the rankings of genetic diversity of different populations based on Ho and He were similar but quite different from that based on nucleotide diversity (π). Considering the three statistics (Ho, He and π) together, the AY population had relatively high genetic diversity despite only eight individuals being found in our field investigation. Contemporary Ne is usually smaller than the real population size . A relatively large contemporary Ne was estimated in AY population, suggesting that there may be more undiscovered individuals of O. henryi in AY population.
Many forest tree species show a high level of genetic variation within populations, with relatively low genetic differentiation among populations . As expected, AMOVA results showed that the genetic variation of O. henryi was mainly maintained within populations, with a moderate and significant genetic differentiation detected among populations (FST = 0.1008, p < 0.001). It is widely accepted that IBD and IBE are two main models of population genetic differentiation, but IBD and IBE are not mutually exclusive and can jointly contribute to population divergence [30, 31]. In this study, IBD was strongly supported by the Mantel tests and MMRR analysis, suggesting that gene flow is geographically limited in O. henryi. Seeds of O. henryi are dispersed over short distances mainly by gravity . However, due to the low germination rate under natural conditions, O. henryi does not regenerate well in the understory and seedlings are rarely seen . Seeds of O. henryi may also attract arboreal frugivorous birds and rodents for long-distance seed dispersal like the other tree species of Ormosia [33, 34], but effective gene flow between O. henryi populations has not developed, resulting in isolation. Populations with low genetic diversity are very sensitive to genetic drift. The factors above may enhance the effect of IBD. In addition, physical barriers such as mountains, rivers or man-made landscapes can facilitate genetic divergence by obstructing dispersal among populations [35, 36]. In southern China, there are several major mountains (e.g., Wuyi Mountains, Nanling Mountains) that might serve as physical barriers to prevent the gene exchange between populations. Strong IBD patterns have been found in many other tree species in this region, such as Castanopsis fargesii , Schima superba  and Sassafras tzumu . On the other hand, no significant effects of IBE were detected by either Mantel tests or MMRR analysis. A strong pattern of IBE is typically interpreted as evidence that divergent selection associated with environmental heterogeneity is maintaining population differentiation . In consist with the insignificant effects of IBE, we found only a few SNPs possibly under selection in the outlier detection analysis, suggesting that the genetic differentiation of O. henryi populations was mainly driven by neutral processes like genetic drift, but not local adaptation to different climate conditions.
Population structure and gene flow
Understanding the population structure of a threatened species can help identify conservation units and develop effective conservation strategies . Significant population structure of O. henryi was identified by the ADMIXTURE analysis, PCA and phylogenetic analysis, supporting the division of four groups. Echoing the IBD pattern, the population structure of O. henryi also exhibited variation with geographic distance. The JH population located in the range of Nanling Mountains is geographically isolated from other populations. The genetic composition of JH population was distinct from other populations regardless of the K value. The NH, SC, MTS, SX and JO populations from the third group located in the adjacent areas of Wuyi Mountains were closely related, and some individuals even clustered with individuals from other populations on the phylogenetic tree (e.g., SC11, SC13, JO12 were clustered with some individuals from SX population). According to the results of PCA and ADMIXTURE analysis, the LQ population from the fourth group was genetically more similar to the populations in the third group than the YK population. The AY and XF populations in southern Jiangxi province were observed with clear genetic admixtures, suggesting that southern Jiangxi province may be a transition zone between eastern and western populations.
Gene flow plays an important role in shaping the spatial genetic structure . Small, isolated populations often suffer from genetic drift and reduced gene flow, which is more likely to generate a prominent population structure . Contemporary gene flow analysis showed that significant recent migrate rates (m) were detected between only a few nearby populations (from XF to AY, from MTS to NH, and from SC to SX). Therefore, the genetic admixtures observed in the AY and XF populations were more likely due to historical gene flow rather than contemporary gene flow. Given the poor dispersal ability of O. henryi, it probably relies primarily on short-distance, gradual dispersal. Without habitat destruction, there can be high levels of gene flow among O. henryi populations. Once the population connectivity is broken by the processes such as habitat fragmentation and human activities, gene flow will be restricted and a significant population structure will generate.
Implications for conservation of O. henryi
The exploration of genetic diversity and population structure can provide vital information for the conservation and management of threatened species . Our population genetic analysis suggests the possible worsening of the O. henryi endangerment situation. The contemporary Ne estimates of O. henryi were much lower than the numbers generally suggested in the 50/500 rule. A Ne of 50 is recommended to avoid short-term inbreeding depression, and a Ne of 500 to maintain the long-term evolutionary potential of a species [45, 46]. The demographic analysis also revealed the sustained declining trend of Ne since the LGP. Accordingly, only a small number of individuals were found in most natural populations during our field investigation. Many isolated trees were found distributed in fragmented habitats such as roadsides, riversides and beside villages. The low Ne estimates, along with the low level of genetic diversity and the observed small population size are warning signs that O. henryi is at increased risk of genetic erosion and should therefore be subject to conservation intervention urgently. Moreover, a moderate level of genetic differentiation and strong IBD pattern were detected in the ten O. henryi populations from southern China. If the conservation and management of O. henryi is not strengthened in time, genetic differentiation among O. henryi populations will further increase, eventually leading to a more endangered status.
O. henryi has attracted much attention for its valuable timber, high ornamental value, and potential medicinal value. In recent years, O. henryi has been gradually used in the transformation of low-quality and inefficient forests such as pure fir forests and pure pine forests in southern China. However, there are still some problems in the conservation and utilization of O. henryi such as inadequate ex-situ conservation, arbitrary selection of provenances in the reforestation process and lack of good varieties. Considering the underestimated endangered status of O. henryi, urgent actions should be taken for both in-situ and ex-situ conservation. In addition, genetic assessment of natural populations outside the highly suitable areas should be conducted to determine a more comprehensive conservation strategy. Based on the findings of this study, several conservation measures are suggested here: (i) In situ conservation should be further strengthened. Priority should be given to the protection of large populations with high genetic diversity in natural forest stands. Maintaining their existing natural habitats is one of the most important ways to protect them. Our results have divided the ten populations in highly suitable areas into four genetic groups, which can be protected as separate evolutionary units. (ii) Ex-situ conservation should be carried out for the long-term survival of O. henryi. The establishment of germplasm resource nurseries and low-temperature seed banks is recommended. Individuals from different genetic groups especially superior trees should be collected to ensure genetic representativeness. (iii) Reintroduction and restoration programs are necessary for increasing population size and preventing further isolation. Seed germination of O. henryi is limited by its hard, dense seed coat and inhibitory substances [32, 47], while seedlings are slow-growing and poorly resistant to drought and frost . Therefore, artificial propagation is recommended. The seedlings can be transplanted to the land near the original populations as they grow up, to help the population recover. Furthermore, the selection of provenances should follow the principle of 'right tree for right place', taking into account factors such as genotypes, forest stands, pests, soil texture and plant growth-promoting rhizobacteria in order to improve the survival rate of seedlings.
In this study, we performed a genome-wide population genetic investigation on natural populations of O. henryi in southern China. Based on the GBS method, we found a low level of genetic diversity, a moderate level of genetic differentiation and significant genetic structures in natural O. henryi populations. Isolation by distance (IBD) is considered as the main reason for the current population genetic structure. Continuous declining trend of the effective population size was suggested, which occurred during the Last Glacial Period. The results from the current genomic data combined with our field investigation suggest that the endangered status of O. henryi is seriously underestimated. Both ex situ and in situ conservation measures should be applied as soon as possible to prevent O. henryi from the fate of extinction. To gain more complete scenario, we call for more studies to decipher the mechanism that leads to the continuous loss of genetic diversity in O. henryi. Our population genomic study provides valuable information for the conservation and management of O. henryi and contributes to the further understanding of preserving the biodiversity in EBLFs. As our sampling area did not cover the entire range of O. henryi, the genetic diversity and differentiation of O. henryi may have been underestimated in this study. In addition, we lacked an in-depth understanding of the important biological characteristics of O. henryi, such as mating system, seed dispersal and natural population regeneration. Therefore, further studies are needed to clarify the endangerment mechanism of O. henryi.
Sample collection and DNA extraction
Results from our preliminary ecological niche modeling of O. henryi suggest that, under current climatic conditions, the core distribution of O. henryi is mainly concentrated in Jiangxi, Hunan, Fujian, Zhejiang provinces in southern China . To understand the genetic structure of O. henryi in its core distribution range, 129 O. henryi individuals were collected from 10 populations in the four provinces of southern China (Fig. 6, Additional file 1: Table S1). In each population, only mature plants were sampled with a minimum distance of 50 m. The samples were carefully identified by Dr. Faxin Yu from Institute of Biological Resources, Jiangxi Academy of Sciences based on the descriptions in Flora of China, a voucher specimen was deposited in the Herbarium of Jiangxi Academy of Forestry with an accession number ZGSZOHE0006. Fresh young leaves from each healthy O. henryi individual were collected and quickly frozen with liquid nitrogen. Samples were taken back to the lab and frozen at -80℃ until DNA extraction.
Total genomic DNA was extracted by the modified CTAB method . DNA purity and concentration were checked by the NanoDrop 1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). DNA interity was then evaluated on 1% (w/v) agarose gel by electrophoresis.
Reference genome construction and genotyping by sequencing (GBS)
One O. henryi sample from XF population was selected for genome survey sequencing. Qualified genomic DNA was randomly sheared into fragments using an ultrasonicator (Covaris, Woburn, MA, USA). A paired-end DNA library was prepared and sequenced with Illumina HiSeq X-ten platform. After removing reads with adapters and reads with low quality, the survey genome assembly was performed using SOAPdenovo  software with the K-mer set at 41.
The GBS libraries were prepared following the method developed by Elshire et al. . Two restriction enzymes (MseI and NlaIII) were used to digest previously extracted DNA. Barcoded adapters were ligated to the digested fragments, and multiplex amplification by PCR was performed. All the samples were mixed, and the desired fragments were selected for GBS library construction in Novogene Co., Ltd (Beijing, China). Paired-end sequencing was performed on an Illumina HiSeq sequencing platform.
Processing of sequencing data and SNP calling
The barcoded adapter sequences of raw reads were removed by the process_radtags program in the STACKS v2.1 software . Raw reads were further filtered and trimmed with the FASTX toolkit v0.0.01 (http://hannonlab.cshl.edu/fastx_toolkit/). The obtained clean reads were aligned to the previously assembled survey genome of O. henryi using bwa v0.7.8  with the following settings: mem -t 4 -k 32 -M. The resulting sam files were converted to bam files by samtools v 1.7 . Picard v1.119 (https://broadinstitute.github.io/picard/) was used to sort bam files by coordinates and mark PCR duplicates. SNP calling was performed using bcftools v1.7  with default parameters. The raw SNP dataset was further filtered with a minor allele frequency of at least 0.05 (–maf 0.05), a genotyping rate of at least 80% (–max-missing 0.8) and a minimum depth of 4 (–minDP 4) by vcftools v0.1.13 .
Population genetic analysis
To fulfil neutrality assumptions in population genetic analyses, we detected potential SNPs under selection using the Bayesian method as implemented in BayeScan v2.1 . To assesss genetic diversity of O. henryi, observed heterozygosity (Ho) and expected heterozygosity (He) for each population was estimated by Arlequin v3.5 , and nucleotide diversity (π) for each population was calculated by vcftools v0.1.13. The inbreeding coefficient (FIS) and polymorphism information content (PIC) was calculated according to the formula [60, 61]. Pairwise genetic differentiation (FST) between populations were calculated by Arlequin v3.5. The contemporary gene flow among populations was estimated by BA3-SNPs program  with 10 million iterations, a burn-in of 1 × 106, and a sampling frequency of 1000. Mixing parameters for allele frequencies, inbreeding coefficients and migration rates were adjusted to achieve acceptance rates between 0.2 and 0.6 as recommended by the software's manual. Five independent runs with different random seeds were performed. Migration rates were considered significant if the 95% credible sets did not include zero.
The population genetic structure was investigated using phylogenetic tree construction, principal component analysis (PCA) and the model-based evolutionary clustering approach as implemented in ADMIXTURE v1.2 . An unrooted neighbor-joining (NJ) tree of all the 129 samples based on all the neutral SNPs was inferred using megacc software with 1000 bootstraps . Principal component analysis (PCA) was performed by the –pca function in PLINK v1.9 software . In the ADMIXTURE analysis, the pre-defined genetic clusters (K) were set from 2 to 10, and different random seeds were set for 10 repeated runs. The most appropriate K value was selected according to the lowest value of cross-validation (CV) error rate, and the results were visualized by R package Pophelper v2.3.1 .
AMOVA analysis was carried out by Arlequin v3.5 to quantify the genetic variation among and within the genetic groups inferred by the population structure analyses. The significance levels for the variance components were tested using 1,000 permutation steps.
Isolation by distance (IBD) and isolation by environment (IBE)
To investigate the roles of geographic and environmental factors in shaping the spatial genetic differentiation, we calculated: (a) the correlation between environmental and geographic distance, (b) isolation by distance (IBD) and (c) isolation by environment (IBE). The geographic distance matrix was calculated based on the GPS coordinates of sampling sites using GenAlex v6.51 . A total of 19 bioclimate variables for each sampling site were obtained from the WorldClim database averaged for the years 1970–2000 at 30 s spatial resolution (~ 1 km) (https://worldclim.org/, ). To avoid multicollinearity, Pearson's correlation coefficients between bioclimate variables were calculated by "cor" function in R. Highly correlated variables (|r|> 0.8) were removed and 10 bioclimate variables were kept for IBE analysis (Additional file 1: Table S4). The environmental distance matrix was calculated by first scaling and centering the 10 selected bioclimate variables to account for differences in magnitude, then calculating pairwise Euclidean distances between sampling sites using R package vegan . Pairwise genetic distance FST /(1- FST) between populations was calculated by Arlequin v3.5. Mantel tests were performed to assess the associations between genetic distance, geographic distance and environmental distance with 9999 permutations in R package vegan. The relative role of geography (IBD) and environment (IBE) on population divergence was also evaluated using multiple matrix regression with randomization (MMRR) implemented in the "MMRR" function of R . MMRR allows quantification of how a dependent variable (genetic distance) responds to changes in explanatory variables (geographic and environmental distance). Regression coefficients of IBD (βD) and IBE (βE) were determined after 9999 permutations.
Demographic history inference
Contemporary effective population size (Ne) for each population was estimated using the linkage disequilibrium method in NeEstimator v.2.1  with the random mating model and the minor allele frequency cutoff of 1/(2n) (n = sample size). The demographic history of O. henryi was inferred using Stairway Plot 2  based on site frequency spectrum (SFS). The folded SFS was constructed by a custom python script easySFS.py with a down projection method (https://github.com/isaacovercast/easySFS). The mutation rate for O. henryi was taken from literature-based estimates for genome-wide mutation rates and was eventually set at 4e-9 per site per year [73,74,75]. The generation time was set as 7 years, which was estimated from our long-term field observation of the time from seeds to flowering in Jiangxi Academy of Forestry. Other parameters were set following the manual's instructions, and bootstrap iterations of 200 were implemented.
Availability of data and materials
The plant materials were collected from natural populations in the geographic distribution of Ormosia henryi. The raw reads files can be accessed on NCBI Sequence Read Archive (SRA) under the BioProject number of PRJNA909598.
Genotyping by sequencing
Principal component analysis
Multiple matrix regression with randomization
Reduced representation genome sequencing
- Ho :
- He :
Polymorphism information content
- F IS :
Single nucleotide polymorphism
Analysis of molecular variance
Isolation by distance
Isolation by environment
Song YC, Da LJ. Evergreen broad-leaved forest of East Asia. In: Box EO, editor. Vegetation structure and function at multiple spatial, temporal and conceptual scales. Switzerland: Springer International; 2016. p. 101–28.
Tang CQ. Evergreen broad-leaved forests. In: Tang CQ, editor. The subtropical vegetation of southwestern China: plant distribution, diversity and ecology. Dordrecht: Springer; 2015. p. 49–112.
Mi XC, Feng G, Hu YB, Zhang J, Chen L, Corlett RT, et al. The global significance of biodiversity science in China: an overview. Natl Sci Rev. 2021;8(7):nwab032.
Dong M, Zang RG. Research in population and conservation biology of rare and endangered plant species. Beijing: Science Press; 2022.
Dai EF, Wu Z, Ge Q, Xi W, Wang X. Predicting the responses of forest distribution and aboveground biomass to climate change under RCP scenarios in southern China. Glob Chang Biol. 2016;22(11):3642–61.
Matocq MD, Villablanca FX. Low genetic diversity in an endangered species: recent or historic pattern? Biol Conserv. 2001;98(1):61–8.
Feng L, Du FK. Landscape genomics in tree conservation under a changing environment. Front Plant Sci. 2022;13: 822217.
Li LF, Cushman SA, He YX, Ma XF, Ge XJ, Li JX, et al. Landscape genomics reveals genetic evidence of local adaptation in a widespread tree, the Chinese wingnut (Pterocarya stenoptera). J Syst Evol. 2022;60(2):386–97.
Christiansen H, Heindler FM, Hellemans B, Jossart Q, Pasotti F, Robert H, et al. Facilitating population genomics of non-model organisms through optimized experimental design for reduced representation sequencing. BMC Genomics. 2021;22(1):625.
Jia KH, Zhao W, Maier PA, Hu XG, Jin YQ, Zhou SS, et al. Landscape genomics predicts climate change-related genetic offset for the widespread Platycladus orientalis (Cupressaceae). Evol Appl. 2020;13(4):665–76.
Liu WZ, Xie JG, Zhou H, Kong HH, Hao G, Fritsch PW, et al. Population dynamics linked to glacial cycles in Cercis chuniana F. P. Metcalf (Fabaceae) endemic to the montane regions of subtropical China. Evol Appl. 2021;14(11):2647–63.
Ji R, Yu X, Ren T, Chang Y, Li Z, Xia X, et al. Genetic diversity and population structure of Caryopteris mongholica revealed by reduced representation sequencing. BMC Plant Biol. 2022;22(1):297.
Torke BM, Cardoso D, Chang H, Li SJ, Niu M, Pennington RT, et al. A dated molecular phylogeny and biogeographical analysis reveals the evolutionary history of the trans-pacifically disjunct tropical tree genus Ormosia (Fabaceae). Mol Phylogenet Evol. 2022;166: 107329.
Li XC, Manchester SR, Xiao L, Wang Q, Hu Y, Sun BN. Ormosia (Fabaceae: Faboideae) from the Miocene of southeastern China support historical expansion of the tropical genus in East Asia. Hist Biol. 2021;33(12):3561–78.
Wang ZX, Shi GL, Sun BN, Yin SX. A new species of Ormosia (Leguminosae) from the middle Miocene of Fujian, Southeast China and its biogeography. Rev Palaeobot Palyno. 2019;270:40–7.
Ma F, Zhang JL. Investigation on key protected wild plant resources in China. Beijing: China Forestry Publishing House; 2009.
Zhang R, Zhou Z, Du K. Genetic diversity of natural populations of endangered Ormosia hosiei, endemic to China. Biochem Syst Ecol. 2012;40:13–8.
Lu Y, Zhu S, He Y, Peng C, Wang Z, Tang Q. Phytochemical profile and antidepressant effect of Ormosia henryi Prain leaf ethanol extract. Int J Mol Sci. 2019;20(14):3396.
Feng S, Hao J, Xu Z, Chen T, Qiu SX. Polyprenylated isoflavanone and isoflavonoids from Ormosia henryi and their cytotoxicity and anti-oxidation activity. Fitoterapia. 2012;83(1):161–5.
Wu GY, Wei XL, Wang X, Wei Y. Induction of somatic embryogenesis in different explants from Ormosia henryi Prain. Plant Cell Tiss Organ Cult. 2020;142:229–40.
Reed DH, Frankham R. Correlation between fitness and genetic diversity. Conserv Biol. 2003;17(1):230–7.
Hughes AR, Inouye BD, Johnson MTJ, Underwood N, Vellend M. Ecological consequences of genetic diversity. Ecol Lett. 2008;11(6):609–23.
Gong X, Yang A, Wu Z, Chen C, Li H, Liu Q, et al. Employing genome-wide SNP discovery to characterize the genetic diversity in Cinnamomum camphora using genotyping by sequencing. Forests. 2021;12(11):1511.
Gao J, Liu ZL, Zhao W, Tomlinson KW, Xia SW, Zeng QY, et al. Combined genotype and phenotype analyses reveal patterns of genomic adaptation to local environments in the subtropical oak Quercus acutissima. J Syst Evol. 2021;59(3):541–56.
von Takach B, Ahrens CW, Lindenmayer DB, Banks SC. Scale-dependent signatures of local adaptation in a foundation tree species. Mol Ecol. 2021;30(10):2248–61.
Li F, Zhou Z, Xie Y. Genetic diversity and genetic differentiation of different populations of Ormosia hosiei in three watersheds. Mol Plant Breeding. 2017;15(10):4263–74.
Li F, Chen H, Liu S, Zhang H, Zhou Z. Mating systems of single families and population genetic diversity of endangered Ormosia hosiei in South China. Genes. 2022;13(11):2117.
Zhang X, Sun Y, Landis JB, Zhang J, Yang L, Lin N, et al. Genomic insights into adaptation to heterogeneous environments for the ancient relictual Circaeaster agrestis (Circaeasteraceae, Ranunculales). New Phytol. 2020;228(1):285–301.
Porth I, El-Kassaby YA. Assessment of the genetic diversity in forest tree populations using molecular markers. Diversity. 2014;6(2):283–95.
Wright S. Isolation by distance under diverse systems of mating. Genetics. 1946;31(1):39.
Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23(23):5649–62.
Wei X, Meng X, Deng Z. Relation between being endangered and seed reproductive ecology of a rare species Ormosia henryi. Seed. 2014;33(01):82–6.
Foster MS, Delay LS. Dispersal of mimetic seeds of three species of Ormosia (Leguminosae). J Trop Ecol. 1998;14(4):389–411.
Liu X, Wang Z, Xiao Z. Patterns of seed predation and dispersal of an endangered rare plant Ormosia hosiei by Edward’s long-tailed rats and Chinese white-bellied rats. Biodivers Sci. 2011;19(01):93–6.
Su H, Qu LJ, He K, Zhang Z, Wang J, Chen Z, et al. The Great Wall of China: a physical barrier to gene flow? Heredity. 2003;90(3):212–9.
Bai WN, Liao WJ, Zhang DY. Nuclear and chloroplast DNA phylogeography reveal two refuge areas with asymmetrical gene flow in a temperate walnut tree from East Asia. New Phytol. 2010;188(3):892–901.
Sun Y, Hu H, Huang H, Vargas-Mendoza CF. Chloroplast diversity and population differentiation of Castanopsis fargesii (Fagaceae): a dominant tree species in evergreen broad-leaved forest of subtropical China. Tree Genet Genomes. 2014;10(6):1531–9.
Hahn CZ, Michalski SG, Fischer M, Durka W. Genetic diversity and differentiation follow secondary succession in a multi-species study on woody plants from subtropical China. J Plant Ecol. 2017;10(1):213–21.
Wang S, Wang Y, Zhou J, Li P, Lin H, Peng Y, et al. Genetic diversity and population structure of an arctic tertiary relict tree endemic to China (Sassafras tzumu) revealed by novel nuclear microsatellite (nSSR) markers. Plants (Basel). 2022;11(20):2706.
Nosil P, Funk DJ, Ortiz-Barrientos D. Divergent selection and heterogeneous genomic divergence. Mol Ecol. 2009;18(3):375–402.
Funk WC, McKay JK, Hohenlohe PA, Allendorf FW. Harnessing genomics for delineating conservation units. Trends Ecol Evol. 2012;27(9):489–96.
Slatkin M. Gene flow and the geographic structure of natural populations. Science. 1987;236(4803):787–92.
Ellstrand NC, Elam DR. Population genetic consequences of small population size: implications for plant conservation. Annu Rev Ecol Syst. 1993;24:217–42.
Laikre L, Allendorf FW, Aroner LC, Baker CS, Gregovich DP, Hansen MM, et al. Neglect of genetic diversity in implementation of the convention on biological diversity. Conserv Biol. 2010;24(1):86–8.
Jamieson IG, Allendorf FW. How does the 50/500 rule apply to MVPs? Trends Ecol Evol. 2012;27(10):578–84.
Frankham R, Bradshaw CJA, Brook BW. Genetics in conservation management: revised recommendations for the 50/500 rules, Red List criteria and population viability analyses. Biol Conserv. 2014;170:56–63.
Deng Z, Wei X. Study on method of breaking seed dormancy for rare tree species Ormosia henryi. Seed. 2016;35(11):1–4.
Zhou J, Zhu T, Wen H, Wu H, Duan X, Feng J, et al. Research progress of artificial reproduction techniques of Ormosia henryi. J Sichuan For Sci Tech. 2019;40(5):104–7.
Song Y, Zhang GL, Jia QQ, Wen Q, Mo XY, Liu LT. Prediction of potential distribution of Ormosia henryi in China under climate change. J Northwest For Univ. 2021;36(6):108–15+273.
Allen GC, Flores-Vergara MA, Krasynanski S, Kumar S, Thompson WF. A modified protocol for rapid DNA isolation from plant tissues using cetyltrimethylammonium bromide. Nat Protoc. 2006;1(5):2320–5.
Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience. 2012;1(1):18.
Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE. 2011;6(5): e19379.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22(11):3124–40.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011;27(21):2987–93.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.
Foll M, Gaggiotti O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics. 2008;180(2):977–93.
Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.
Wright S. Coefficients of inbreeding and relationship. Am Nat. 1922;56(645):330–8.
Botstein D, White RL, Skolnick M, Davis RW. Construction of a genetic linkage map in man using restriction fragment length polymorphisms. Am J Hum Genet. 1980;32(3):314–31.
Mussmann SM, Douglas MR, Chafin TK, Douglas ME. BA3-SNPs: Contemporary migration reconfigured in BayesAss for next-generation sequence data. Methods Ecol Evol. 2019;10(10):1808–13.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19(9):1655–64.
Kumar S, Stecher G, Peterson D, Tamura K. MEGA-CC: computing core of molecular evolutionary genetics analysis program for automated and iterative data analysis. Bioinformatics. 2012;28(20):2685–6.
Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4:7.
Francis RM. Pophelper: an R package and web app to analyse and visualize population structure. Mol Ecol Resour. 2017;17(1):27–32.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics. 2012;28(19):2537–9.
Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37(12):4302–15.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Vegan: Community Ecology Package. R Package version 2.5–2. Available at: https://CRAN.R-project.org/package=vegan (Published: November 28, 2020).
Wang IJ. Examining the full effects of landscape heterogeneity on spatial genetic variation: a multiple matrix regression approach for quantifying geographic and ecological isolation. Evolution. 2013;67(12):3403–11.
Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol Ecol Resour. 2014;14(1):209–14.
Liu X, Fu YX. Stairway Plot 2: demographic history inference with folded SNP frequency spectra. Genome Biol. 2020;21(1):280.
Shalev TJ, Gamal El-Dien O, Yuen MMS, Shengqiang S, Jackman SD, Warren RL, et al. The western redcedar genome reveals low genetic diversity in a self-compatible conifer. Genome Res. 2022;32(10):1952–64.
Yang F, Cai L, Dao Z, Sun W. Genomic Data Reveals Population Genetic and Demographic History of Magnolia fistulosa (Magnoliaceae), a Plant Species With Extremely Small Populations in Yunnan Province. China Front Plant Sci. 2022;13: 811312.
He Z, Feng X, Chen Q, Li L, Li S, Han K, et al. Evolution of coastal forests based on a full set of mangrove genomes. Nat Ecol Evol. 2022;6(6):738–49.
We would like to thank Dr. Yao Zhao at Center for Watershed Ecology, Institute of Life Science and School of Life Sciences, Nanchang University, for his help in using the software and valuable discussions.
This work was supported by Major Programs of Scientific Research and Achievement Popularization of Jiangxi Academy of Forestry (2017511201), Forestry Science and Technology Innovation Project of Jiangxi Forestry Bureau (20236) and Special Fund of Fundamental Research and Talents of Jiangxi Academy of Forestry (2023521202).
Ethics approval and consent to participate
The authors have complied with all relevant institutional, national and international guidelines and the collection permissions for plant materials were granted by the local Provincial Forestry Bureau.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sampling information of Ormosia henryi populations. Table S2. Statistics of the O. henryi genome assembly. Table S3. Summary of GBS data. Table S4. Ten bioclimate variables used in this study.
Potential SNPs under selection identified by BayeScan. The vertical solid line represents the threshold for being under selection after correction with false discovery rate (0.05). Six SNPs colored in red were probably under selection. Fig S2. Cross-validation (CV) errors of 10 repeat runs under different K values. Fig. S3. Demographic history of three genetic groups of Ormosia henryi inferred by Stairway Plot 2. The x-axis indicates the time before the present, and the y-axis represents the historical effective population size. The first genetic group were not used in the demographic history analysis due to the small number of individuals.
About this article
Cite this article
Zhou, C., Xia, S., Wen, Q. et al. Genetic structure of an endangered species Ormosia henryi in southern China, and implications for conservation. BMC Plant Biol 23, 220 (2023). https://doi.org/10.1186/s12870-023-04231-w