Miocene diversification of a golden‐thread nanmu tree species (Phoebe zhennan, Lauraceae) around the Sichuan Basin shaped by the East Asian monsoon

Abstract Understanding the role of climate changes and geography as drivers of population divergence and speciation is a long‐standing goal of evolutionary biology and can inform conservation. In this study, we used restriction site‐associated DNA sequencing (RAD‐seq) to evaluate genetic diversity, population structure, and infer demographic history of the endangered tree, Phoebe zhennan which is distributed around the Sichuan Basin. Genomic patterns revealed two distinct clusters, each largely confined to the West and East. Despite sympatry of the two genomic clusters at some sites, individuals show little or no evidence of genomic introgression. Demographic modeling supported an initial divergence time between the West and East lineages at ~15.08 Ma with further diversification within the West lineage at ~7.12 Ma. These times largely coincide with the two independent intensifications of the East Asian monsoon that were initiated during the middle (Langhian) and late Miocene (Messinian), respectively. These results suggest that the Miocene intensification phases of the East Asian monsoon played a pivotal role in shaping the current landscape‐level patterns of genetic diversity within P. zhennan, as has been found for the interspecific divergence of other subtropical Chinese plants. Based on isolation‐by‐distance and species distribution modeling, we hypothesize that P. zhennan followed a ring diversification which was facilitated by the Sichuan Basin acting as barrier to gene flow. In situ and ex situ conservation management plans should consider the results obtained in this study to help secure the future of this beautiful and culturally significant endangered tree.


| INTRODUC TI ON
It is well known that global climate change and geography have profound impacts on population divergence and speciation (Schluter & Pennell, 2017;Zou, Yang, Doyle, & Ge, 2013) and the climate of Asia has experienced a series of drastic changes since the Miocene.
During the Pleistocene, the monsoon system interacted with glacial-interglacial cycles, producing a more variable monsoon climate Sun & Wang, 2005). Previous studies have shown that the dynamic history of the East Asian monsoon played a major role in organism evolution in the QTP and adjacent regions (see review in Favre et al., 2015), largely by controlling the hydrological cycles (Tada et al., 2016). Most of these earlier studies focused on interspecific evolutionary process, and few have emphasized the role of past historical events and evolutionary factors on population divergence and speciation since the Miocene; a period that witnessed some of the most significant climate changes and orogenesis in Asia Guo et al., 2008;Zachos, Dickens, & Zeebe, 2008).
Geographic isolation can facilitate divergence and speciation events, because barriers such as distance, water bodies, and mountains can impede gene flow and drive genetic differentiation through allopatric or parapatric speciation (Pritchard, Stephens, & Donnelly, 2000;Winger & Bates, 2015). Southwest China, particularly the mountains surrounding the Sichuan Basin (e.g., the Hengduan Mountains to the west, the Daba Mountains to the north and northeast, and the Wuling and Dalou Mountains to the south and southeast), has been identified as a global diversity hotspot, due in part to its mild monsoon climate and complex topography (Myers, Mittermeier, Mittermeier, Da Fonseca, & Kent, 2000). This region has been recognized as a "plant museum" for relictual angiosperms and other seed plants by providing long-term ecological and environmental stability (López-Pujol, Zhang, Sun, Ying, & Ge, 2011). Several phylogeographic studies suggested that geographic and climatic events during the late Miocene (11.2-5.3 Ma) played important roles in the development of intraspecific lineages in southwest China (Qi et al., 2012;Sun et al., 2014;Zhang et al., 2013). However, only a few of these studies have focused specifically on the role of the Sichuan Basin as driver of speciation and for most taxa and the interactions between geography and climate as far as back as the late Miocene remain understudied.
Phoebe zhennan (Lauraceae) is a diploid (2n = 24), dioecious tree up to 30 m tall which produces purple-brown drupes and grows in subtropical evergreen broad-leaved forests (EBLFs) of China (Li, Li, & van der Werff, 2008). Chinese Phoebe L. species are insect-pollinated, and the seeds are dispersed primarily by frugivorous birds (Li et al., 2008;Li, Liu, Ma, Zhang, & Xu, 2018). Phoebe zhennan is the major source of the well-known wood "Golden-thread nanmu," which is extremely valuable due to its durability, unique special fragrance, and attractive golden color. This timber has been used for the manufacture of coffins and the construction of palaces and furniture for over a thousand years (Lan, 1994;Zhang, 2010). As a result of its desirability, the current sporadic distribution of the around the Sichuan Basin is largely due to habitat loss and logging and most extant populations consist of fewer than 70 individuals (Ding, Xiao, Huang, & Li, 2015;Fu, 1992) and the species is listed as endangered in the Chinese Plant Red Book (Fu, 1992).
Although previous studies of P. zhennan mainly focused on its habitat characteristics , seedling biology (Xie et al., 2017), and species delimitation from other "Nanmu" tree species (Ding, Xiao, Li, Conran, & Li, 2019), little is known about its evolutionary history. An earlier investigation of genetic diversity using amplified fragment length polymorphisms (AFLPs) in P. zhennan showed low genetic diversity and significant divergence, apparently linked with topography in southwestern China (Gao et al., 2016).
However, that study was inconclusive about the demographic history of P. zhennan, or whether past climate changes may have had a strong impact on its evolutionary history, especially in relation to the observed "ring-species" pattern of variation.
Restriction site-associated DNA sequencing (RAD-seq) is a popular approach for ecology, evolutionary population genomic and conservation genomics of nonmodel species (with or without a prior genomic resource), which can generate genome-wide single nucleotide polymorphisms (SNPs) cost-effectively and quickly (Baird et al., 2008;Andrews, Good, Miller, Luikart, & Hohenlohe, 2016;Parchman, Jahner, Uckele, Galland, & Eckert, 2018). Additionally, RAD-seq leverages the site frequency spectrum (SFS), providing a powerful tool for demographic history inference (Excoffier, Dupanloup, Huerta-Sánchez, Sousa, & Foll, 2013;Gutenkunst, Hernandez, Williamson, & Bustamante, 2009;Hotaling et al., 2018;Sousa & Hey, 2013). Therefore, we applied SNPs markers generated through RAD-seq to investigate genetic differentiation patterns of P. zhennan around the Sichuan Basin and infer its demographic history to understand processes that may have led to extant patterns. We addressed the following two questions in regard to the population genomics and the demographic history of P. zhennan: (a) How is genetic diversity distributed and what is the population structure? (b) Can we identify the factors that impacted on the current distribution of genetic diversity within the species?

| Sample collection and RAD-seq sequencing
A total of 72 P. zhennan trees were sampled from across 12 sites around the Sichuan Basin (Table 1), with voucher specimens processed and deposited in the Herbarium of Xishuangbanna Tropical Botanical Garden (HITBC). Fresh leaves were dried in silica gel and stored at −20°C before processing for DNA extraction. Genomic DNA was extracted using a modified CTAB protocol (Rogers & Bendich, 1989), and quality checked by electrophoresis on 1% agarose gels and normalized to a DNA concentration of 30 ng −mL . RADseq libraries were prepared following the protocol described in Baird et al. (2008). Briefly, genomic DNA was digested using the restriction enzyme EcoRI (5′-G*AATTC-3′). Fragmented DNA was then ligated to Illumina sequencing adaptors (Illumina Biotechnology Company) containing sample-specific barcode sequences, followed by PCR.
Library preparation and pair-end sequencing (85 bp) were carried out by the Beijing Genomics Institute using an Illumina Hiseq 2500.

| De novo assembly and SNPs calling
We performed de novo assembly and single nucleotide polymorphisms (SNPs) using the STACKS 1.35 pipeline (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011;Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). Prior to analysis, poor quality reads (phred scores less than 30), reads with possible adapter contamination, and those lacking restriction sites (using process_ RADtags) were removed. Remaining reads were then assembled into loci in "ustacks" with a maximum distance between stacks of M and minimum read depth of m. The loci were clustered further into "cstacks," with mismatches allowed between samples. To increase coverage depth and maximize loci, a minimum read depth of 6 (m = 6) and a maximum distance between stacks of 2 in "ustacks" (M = 2), and 6 mismatches in "cstacks" (n = 3) in STACKS was set. After finishing de novo assembly, SNPs calling and genotyping were performed using "populations" model in STACKS. Loci with 100% presence among populations (p = 12) were retained, whereas loci with >25% missing data within populations were removed (r = 75%) to maximize shared SNPs across individuals. The first SNP per locus and SNPs with minor allele frequencies (MAF) ≥ 0.05 were retained. After completing SNP calling, those remaining with a Hardy-Weinberg equilibrium (HWE) p-value less than .05 were also removed using VCFtools v4.0 (Danecek et al., 2011). PGDspider v2.02 (Lischer & Excoffier, 2012) was used subsequently for file conversion to program-specific formats.

| Population genomic analyses
Three datasets were generated for the population genetic analyses: (a) the full matrix with all polymorphic loci for genetic sum-  (Pritchard et al., 2000). Because there are often difficulties in resolving genomic datasets using STRUCTURE, Rodriguez-Ezpeleta et al. (2016) suggested that a dataset with up to 5,000 SNPs is an optimal analysis size to capture the ancestral group. Note: For RAD-seq sequencing, we collected 12 populations. However, AX, PW, CK, and WX sites were only used for species distribution modeling (SDM) because we cannot collect samples in these collection areas and adjacent regions. The R1 point is used to estimate "ring distances" among populations for isolation-by-distance analysis (IBD).

TA B L E 1
Site characteristics, sample sizes (n), and habitat information for 12 populations of P. zhennan across its distribution in Sichuan Basin and adjacent areas Genetic summary statistics for RAD-seq genomic data, including percentage of polymorphic sites (%P), observed and expected heterozygosity (H obs and H exp ), nucleotide diversity (π), and inbreeding coefficients (F IS ), were estimated using populations in STACKS for all populations (10) with more than six samples. Pairwise F ST values were calculated using Arlequin v3.5.2 (Excoffier & Lischer, 2010), with 10,000 permutations. Hierarchical AMOVA analysis of molecular variance (AMOVA) was implemented based on our assessment the hierarchical population structure (K = 2 and K = 3, see results and discussion) in order to quantify genetic variation partitioning across the different sampling levels.
Genetic structure was inferred using Bayesian clustering in STRUCTURE based on the minimum dataset (e.g., the randomly selected 5,000 SNPs dataset). An admixture model with no priori information for sampling location was utilized to determine whether the number of clusters (K) ranged from 1 to 13. For each value of K, 10 independent simulations were conducted with a burn-in length 100,000, followed by 300,000 iterations. The best K value was identified based on the successive change of LnP(D), Evanno's delta K, and individual assignment probabilities (Evanno, Regnaut, & Goudet, 2005), as implemented in STRUCTURE HARVESTER v0.6.93 (Earl & Vonholdt, 2012). Results were plotted using DISTRUCT (Rosenberg, 2004), and population structure investigated further with PCA using the R package adegenet (Jombart, 2008).
We used TreeMix v1.13 (Pickrell & Pritchard, 2012) to address historical relationships between populations. The VCF file was converted to a frequency file that could be transformed to a TreeMix file using the plink2treemix script. The number of migration events was tested by starting at zero and adding one by one until the residual plot stopped improving. The resulting maximum-likelihood (ML) tree and residual plot were visualized in R and illustrated in combination with the STRUCTURE results.

| Isolation-by-distance (IBD) analysis and cline analysis
An isolation-by-distance pattern was examined using the R package vegan (Oksanen, Kindt, Legendre, & O'Hara, 2006). A mantel test implementing vegan with 1,000 permutations was used to detect significant correlations between genetic and geographic distances, with F ST /(1 − F ST ) values used to represent genetic distance. Geographic distances among the sites with "straight-line distances" and "ring distances" were calculated in the R package geosphere (Hijmans, Williams, & Vennes, 2011). The ring distances assume dispersal routes following the ring distribution, instead of straight lines across intervening uninhabitable areas. For example, initial analyses suggested that there is a genetic "barrier" between QCS and QL with little evidence of gene flow across the gap. Based on predicated suitable habitat and the ringshaped distribution, we therefore used one reference point (R1) to approximate the ring distance between sites. All distances between Western and eastern sites were estimated by the sum of the distances from Western sites to R1 and then R1 to eastern sites.
Geographic cline analysis was conducted to detect any sudden changes in population genetic composition or "hybrid zones" along the distribution ring using the R package hzar (Derryberry, Derryberry, Maley, & Brumfield, 2014). Clines were defined by distance from site QCS in one dimensional space, following an approximate circle around the basin and using the first principal component analysis axis (PC1) as a proxy of genetic variation. Five cline models with different fitting tails (none fitted; left only; right only, mirror tails; both tails estimated separately) and the minimum and maximum genetic variation (pMin, pMax) fixed-to-observed values were tested. To facilitate model comparison, a null cline model was also established, which assumes that genetic variation is independent of any clines. The best model was decided upon Akaike's information criterion (AIC) scores. Finally, the maximum-likelihood clines and summary statistics were extracted from the best-fit model. Twenty demographic models were tested ( Figure S1; detailed in Appendix S1), and detailed model schematics are provided in Appendix S1. Briefly, models included scenarios specifying: (a) two divergence events for all possible topologies of the three genetic clusters (models 1-9); (b) models where admixture between two existing genetic clusters created the third (models 10-18); and (c) trifurcation models where extant genetic clusters emerged simultaneously from a common ancestor (models 19-20). For all models, the potential for bidirectional gene flow was varied both historically and recently.

| Demographic history inference
The fastsimcoal2 analyses followed an initial set of model selection runs, with comparisons of maximum observed and expected likelihoods to select the best-fit model, then subsequent parameter estimation through simulation of new SFS for the best-fit model, followed by parametric bootstrapping (Method S2; detailed in Appendix S1).
Fifty independent parameter estimations were performed to achieve the maximum composite likelihood of the joint fold-SFS, in which parameterized the simulation sample from prior and the parameter estimation was optimized through 40 cycles of a conditional maximization algorithm. The point estimates were selected from the run with the highest maximum composite likelihood, with confidence intervals of parameter estimates obtained by 100 parametric bootstrapping runs from simulated SFS of the parameter estimates.

| Species distribution modeling (SDM)
SDM was carried out in Maxent v3.4.1 (Elith et al., 2011;Phillips, Anderson, & Schapire, 2006) to predict the potential and suitable distribution range of the species and also to investigate whether Sichuan basin may have served as a barrier potentially facilitating ring diversification. A total of 47 species occurrence records were obtained from herbarium collections and specimen records (Table S1; detailed in Appendix S1). Nineteen bioclimatic variables (Table S2; (Synes & Osborne, 2011). For each highly correlated variable pair (Pearson's r ≥ 0.7), the variable that gave a higher value in the regularized gain and the percent contribution to the Maxent model was retained (Method S1; detailed in Appendix S1).
Maxent was configured with 75% of species presence data for training and 25% for testing data, and sampling procedure was replicated 50 times. The area under the curve (AUC) was estimated to test the accuracy of the model prediction (Method S1; detailed in Appendix S1). The distributions of these periods were then plotted on a Chinese map using a geographic information system as implemented in the software ArcGIS 10.4 (Environmental Systems Research Institute, Inc.)

| RAD-seq dataset
After filtering reads with an average quality score below 30, the   (Table 2). In addition, the southern and southwestern populations showed the higher genomic diversity whereas the populations at the tails of the chain (QL and TL population) had lower genetic genomic diversity (%P and π, showed the same trend when using all nucleotides (Table 2).  (Table 3). Similarly, pairwise F ST values between QL and TL, the junction between the end of the Western sites (QL) and southeastern sites (TL) was particular high (F ST = 0.446) ( Table 3).

| Genomic diversity
An unexpected relatively high pairwise F ST between YC and TL sites The first and second principal components (PC1 and PC2) explained 28.93% and 14.51% of the total variation, respectively. These clusters are congruent with clusters inferred with the K = 3 results in STRUTURE (Figures 1b and 2a).
The maximum-likelihood tree for the 12 populations (Figure 2b) also largely reiterated the relationships inferred with the cluster analysis and explained 97.08% of the population variance relatedness. When migration events were added to the tree, there were three migration edges and the model explained 99.08% of the relatedness variance. These migration edges indicated three admixture events: (a) from the ancestor of the East group into EMS (w = 47%), (b) from the ancestor of the East group into XY (w = 29%), and (c) from YC into QCS (w = 12%; Figure 2b; detailed in Appendix S1). Within the broad West group, there was strong evidence for admixture, supported by visually apparent residuals ( Figure S2; detailed in Appendix S1), which is consistent with previous clustering analysis.
Despite the different sampling levels, the AMOVA revealed significant differentiation among groups identified by STRUCTURE and most the variation (K = 2:67.3% and K = 3:68.72%, respectively) occurred within populations ( Table 4). The AMOVA results showed significant genetic differentiation in the two major groups (East and West groups) and among populations (F CT = 0.174; F ST = 0.327). The AMOVA results also suggested that there is strong genetic isolation between the East, West1, and West2 groups (F CT = 0.207; Table 4).

| IBD and cline analysis
A significant IBD relationship among the populations was also de-

| Demographic history
Results of the fastsimcoal2 analysis showed that the best-supported scenario was model1, with an initial divergence between the East and West lineages and subsequent divergence between the West1 and West2 sublineages, with bidirectional gene flow ( Figure 5;

| Species distribution model (SDM)
The potential and suitable distribution ranges of P. zhennan predicted by SDM are shown in Figure 1c, with AUC values >0.9 implying that the performance of Maxent model is outstanding. The current most suitable habitat is clearly ring-shaped, encircling the Sichuan Basin,

F I G U R E 1 (a) Map of Sichuan
Basin with sampling site of P. zhennan. Sample sites are color-coded according to their genetic clusters (K = 3). The black points represent specimen records from the Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn/) and are not subjected to sequencing. The black pentagram represents reference point, which is used for estimating "ring distance" between sites for isolation-by-distance analysis (IBD). (b) Individual assignment probability barplots from the genetic clustering analysis using the STRUCTRUE. K = 2 is the best fit for our data.

| Population genomic diversity and diversification around the basin
A pattern of low level genomic diversity within populations, coupled with most variation occurring within populations and significant genetic differentiation between groups (   Note: Demographic parameters are reported for the M1 model passing both the AIC and the goodness-of-fit criteria. Confidence intervals were generated by parametric bootstrapping.

TA B L E 5 Parameter estimates from Fastsimcoal2 analyses
Second, compared to wind-pollinated species, pollen of insect-pollinated species (such as most species in Lauraceae, Li et al., 2018;Rohwer, 1993) disperses over shorter distances, resulting in limiting gene flow and increased genetic differentiation. For example, the pollinators of Neolitsea sericea (Blume) Koidz. generally only travel about 4 km (Chung, Chung, Oh, & Epperson, 2000), so short-distance pollination could similarly be a factor influencing the population structure of P. zhennan. Lastly, the patterns observed here may be a consequence of past genetic bottlenecks, or genetic drift associated with small effective population sizes and founder effects (Birky, Fuerst, & Maruyama, 1989). This bottleneck model was supported by demographic history inference, which showed that the ancestral population size for P. zhennan was larger than any of the descendant ones (Table 5;  Similar to previous AFLP results (Gao et al., 2016), the clustering analyses (PCA and STRUCTURE) were congruent and indicated that all individuals clearly group into two main groups.
These groups were defined largely by geography (West and East), but with a few exceptions, with the West and East (including south and southeast populations) groups showing limited connectivity (Figures 1b and 2). Interestingly, at two sites (YC and QCS), individuals were sampled that showed no signs of genomic introgression but which clustered with samples from the divergent geographic group (Figures 1b and 2). The PCA and ML analysis results revealed two distinct sublineages within the West region: West1 and West2 (Figures 1b and 2; also congruent with

| The dynamics of the East Asian monsoon shaped a middle Miocene diversification of Phoebe zhennan
Several processes likely contributed to the patterns of divergence inferred with genomic data. The intensification of the East Asian monsoon has previously been proposed as a process that shaped intra species diversification. Our demographic analysis suggested that the timing of divergence between the two major groups occurred during the early middle Miocene (15.08 Ma), with subsequent divergence within the West clade during the early late Miocene (7.12 Ma). of the Northeast India (e.g., P. sublanceolata; Bhattacharyya, 1983), extending there to the Pliocene (Khan & Subir, 2014). Phoebe fossils have also been reported from the early Miocene to late Pliocene of China (Huang et al., 2016;WGPC, 1978). Second, divergence timing in demographic inference is influenced strongly by the settings used for Fastsimcoal2 with SFS data. We chose the model process after 50 replications for every candidate model and also considered the influence of different mutation rates, with the results almost always identical for the origin of P. zhennan (see supplementary information).
This study also provides an opportunity to explore the impact Second, the initial intensification of the East Asian monsoon in the middle Miocene triggered the intraspecific level diversification of lineages inhabiting subtropical EBLFs in East Asia (Yu et al., 2017).
For example, divergence in Cyclocarya paliurus coincides with the first intensification of the East Asian monsoon, which provided the necessary humid climatic conditions for it to survive in southwestern of China (Kou et al., 2016). Similarly, the distribution of the fern genus Lepisorus (J. Sm.) Ching correlates with precipitation brought by the summer monsoon and its radiation matches the initial intensification of the East Asian monsoon ~15 Ma . The second divergence event seen between the West1 and West2 sublineages seems to have occurred during the late Miocene, coinciding with the second phase of the East Asian monsoon intensification at ~9-7 Ma Harrison, Copeland, Kidd, & Yin, 1992). Similar diversification patterns have also been reported for other Asian plants, with a phylogeographical study of Tetracentron sinense Oliv. (Trochodendraceae) in southwestern China revealing intraspecific species divergence ~9.6 Ma (Sun et al., 2014). Deep genetic divergences are also seen along the Eastern Margin of the Yungui Plateau, a well-known geological boundary (CASPG, 1984), and these may have been influenced by periods of rapid uplift of the QTP. Additionally, many subtropical EBLF trees (including P. zhennan) are very sensitive to winter temperatures and will die at temperatures below −10°C (Sakai, 1979;Woodward & Williams, 1987). Such a north-south distribution shift might therefore be explained in part by late Miocene global cooling (Zachos et al., 2008) (Figure 1c). Interesting, the IBD analysis also produced a very strong correlation between genetic distance and the ring distances ( Figure 3a). This correlation was only slightly Additionally, in the populations collected around the Basin, we found a decrease in genetic diversity from southern populations toward the southeast and Western populations, that is, the two ends of the ring (Figure 4).
Topographic modeling has suggested that the Sichuan Basin is an excellent candidate for the geographic facilitation of ring speciation (Monahan, Pereira, & Wake, 2012

| Conservation management
To facilitate the long-term conservation and sustainability of a species, it is widely recognized that it is crucial to protect and encourage within-species genetic diversity, as genetic diversity increases the ability of a species to evolve and adapt in the long run to a changing environment. Processes such as habitat fragmentation or loss and selective logging can have detrimental effects on both the survival of a species and its genetic diversity, and these factors are already impacting on P. zhennan. Golden-thread wood has been part of Chinese culture for many thousands of years (Li, Jin, & Xiang, 2004), and conservation is both valuable for biodiversity and cultural heritage. Here, we identified three genetically divergent population groups within the species, each with its own unique pattern of diver-  (Meeus, Honnay, & Jacquemyn, 2012). The QCS population is located in Fengshui forest that is found near an ancestral temple, and the YC population is located in a famous tree farm which not only cultivates economical trees, but also services as an ex situ conservation unit for endangered species. Human introduction is therefore a likely source of the disjunctions, but we cannot discount frugivorous bird dispersal, particularly for the individuals in QCS belonging to the Western clade. A ring diversification pattern may have also contributed toward the disjunct nature of QCS (as discussed above).
With the knowledge obtained though population genomics, the best conservation strategy for P. zhennan is in situ. Unfortunately, several populations (e.g., QCS, QL, TL, and YC) occur in regions where the human disturbance such as fengshui forestry and tree farming greatly limits available habitat in natural reserves. These reduced habitats can also severely limit gene flow and seedling establishment due to habitat fragmentation and anthropogenic disturbance, a situation that has been observed in case studies of other Chinese Lauraceae such as Phoebe bournei (Hemsl.) Yang (Ge, Liu, Shen, & Lin, 2015) and Neolitsea sericea (Blume) Koidz. .
Phoebe zhennan populations should be protected not only for their genetic diversity but also for their cultural value. Pre-and postzygotic barriers may play a vital role in selection against biparental inbreeding, and the relatively low levels of F ST within the identified genetic clusters and higher levels of observed heterozygosity suggest that it is vital to monitor factors that facilitate outcrossing in P. zhennan such as pollination, seed dispersal and predation, seedling recruitment, and herbivory (Neuschulz, Mueller, Schleuning, & Böhning-Gaese, 2016).
When ex situ conservation or augmentation is considered, it is crucial that genetic diversity is optimized within the three structure groups, with the mixing of the Western and Eastern structure groups minimized. Translocations using individuals from both Western and Eastern groups may not induce higher genetic diversity, but instead may bring together individuals from groups that cannot readily interbreed. The species will benefit from the continuation of the ban on felling, and while ancient trees are of cultural value, these trees are also likely to decrease the effects of genetic drift and help to increase within-population genetic diversity while safeguarding the valuable diversity accumulated across past climatic periods.
Finally, the markers obtained in this study could be useful for future forensic work to identify the origin of illegally felled gold-thread nanmu samples and could be developed further for the policing of illegal trade in closely related species such as P. bournei, P. hui Cheng ex Yang, and P. chekiangensis C.B. Shang.

| CON CLUS IONS
This study used population genomics to evaluate genomic diversity and genetic structure in an endangered plant species (P. zhennan) and infer its demographic history. We found that P. zhennan displays a pattern of low level of genomic diversity within populations, cou-