During our field survey of X. vietnamensis in Guangxi, China from 2013 to 2020, there were eight karst mountains in Guangxi where this species was found. Five and three mountains were located in the northern and southern border of Guangxi respectively (Figure 1; Table 1). Each mountain had a small number (1–20) of wild individuals growing on or near steep slopes, and the total number of individuals was approximately 45. Fresh leaves of 33 individuals were collected from six sites (mountains), immediately dried in silica gel and stored at room temperature prior to DNA extraction. Additionally, fresh and young leaves off our mature individuals were collected from each of the four sites with more than four individuals (Table 2) and kept frozen at −80°C until analyzed for ploidy level determination.

Sampling was limited because X. vietnamensis individuals were rare, and some individuals or sites were inaccessible (Table 1). Genomic DNA was extracted from the leaf tissue of all sampled individuals following the cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1990).

Total RNA was extracted from fresh leaves of one X. vietnamensis individual for Illumina sequencing (Mao et al., 2019). The cDNA library was sequenced using the Illumina HiSeq 2000 system at Novogene (Beijing, China). Sequences were filtered, and clean reads were assembled de novo using Trinity v.2.8.5 (Grabherr et al., 2011). Short sequence repeat (SSR) of motifs were identified using the Perl script MISA (Thiel et al., 2003). Primer3 v.2.3.6 was used to design primer pairs for the detected markers (Rozen and Skaletsky, 2000; Untergasser et al., 2012). A total of 300 SSR loci were applied for pilot screening of polymorphic makers based on 16 individuals, and these polymorphic loci were then employed to detect genetic variation for all sampled individuals.

PCR amplification was performed in 25 μl reactions containing 50 ng of genomic DNA, 12.50°μl of 2× PCR buffer, 300 μM of each dNTP, 0.3 μM of primer pairs, and 1.25 U of Taq DNA polymerase (all from Vazyme Biotech, Nanjing, China). The amplification program was performed as follows: initial denaturation at 98°C for 3°min; followed by 40 cycles at 98°C for 10 s, annealing temperatures (see Supplementary Table S1) of the primer pair for 30°s, 72°C for 45°s; and a final extension step at 72°C for 10 min. Each forward primer was labeled with either FAM, HEX, and TAMRA at the 5′ extremity to allow fragment detection (Supplementary Table S1). Amplified products were inspected in a 1% agarose gel and electrophoresed on the ABI PRISM 3100 genetic analyzer (Applied Biosystems, Foster City, CA, United States). The microsatellite genotype at each locus, for each individual was determined using GeneMapper v.4.1 (Soft Genetics, State College, PA, United States). The codominant SSR data were analyzed using Cervus3.0 software (Kalinowski et al., 2007) to detect genotyping errors at each locus, such as stuttering, large allele dropouts and null alleles. When the proportion of null alleles at a locus is greater than 40% (F _(Null) > 0.4), it was deemed to contain excessive null alleles. Loci were tested for natural selection using BayeScan (Foll and Gaggiotti, 2008), and the loci under selection were excluded for further population genetic analyses.

The following measures of genetic diversity were calculated for each locus, population, and geographic group (south and north) using GenAlEx v.6.5 (Peakall and Smouse, 2012): number of different alleles (A), number of effective alleles (A _E), Shannon’s diversity index (I), observed heterozygosity (H _O), expected heterozygosity (H _E), inbreeding coefficient (F _IS), and genetic differentiation coefficients between groups (F _ST) per locus. The statistical significance of the deviation of loci from Hardy-Weinberg equilibrium (HWE) and the linkage disequilibrium test between two loci were assessed by the Markov chain method in Genepop v.4.7 (Rousset, 2008). Sequential Bonferroni correction was used to determine significance level at p < 0.05 (Rice, 1989). We also implemented hierarchical analysis of molecular variance (AMOVA) (Excoffier et al., 1992) in GenAlEx v.6.5 to calculate the level of genetic variation and differentiation among geographic groups, among populations within groups, and among individuals within populations based on F-statistics (Wright, 1978), using data of the four populations with more than one sample. The significance of the differences was tested using permutation procedures (999 replicates).

The bottleneck test for population was estimated using the M-ratio method (Garza and Williamsion, 2001), which is the ratio of the total number of alleles to the overall range of allele size. We assumed a microsatellite mutation rate per locus per generation of 10⁻⁴ and pre-bottleneck effective population size of 400 (θ = 4N _e μ = 0.16) and 4,000 (θ = 1.6) to estimate the M-ratio. The mean size of non-single-step mutations and the percentage of mutations larger than a single step were set to 3.5 and 0.1, respectively (Garza and Williamsion, 2001). The M-ratio (M ₀) estimated using the program M_P_Val was compared to the 95% critical M-ratio (M _c) estimated using the Critical_M program. We assumed that the population experienced a significant bottleneck if M ₀ < M _c when N _e = 4,000, and a moderate bottleneck if M ₀ < M _c when N _e = 400.

The Bayesian clustering method in STRUCTURE v.2.3.4 (Pritchard et al., 2000) was used to assess the genetic clustering of sampled X. vietnamensis populations. This analysis was run for 10 independent runs per K value (1–10) with a burn-in period of 50,000 iterations and 100,000 Markov chain Monte Carlo (MCMC) iterations, using the admixture (allele frequencies correlated) model. Structure Harvester (Earl, 2012) was used to visualize the best K value based on delta K and maximum log likelihood L(K) (Evanno et al., 2005). Discriminant analysis of principal components (DAPC) was performed to identify the differentiated subpopulations using the R package Adegenet (Jombart, 2008; Jombart et al., 2010). The number of clusters was set at 10. The microsatellite data was transform into principal components (PCs) using the function find.clusters, retaining all the PCs in the analysis. The Bayesian information criterion (BIC) was used to infer the optimal number of clusters (Jombart et al., 2010). The functions a.score and optim.a.score were applied to identify the optimal number of PCs to be retained. The function dapc performed discriminant analysis using the optimal number of PCs and four eigenvalues were retained.

MIGRATE-n (version 3.3.1) was used to estimate the mutation-scaled effective population size Θ (4N _e μ; N _e, effective population size; μ, mutation rate per locus per generation), mutation-scaled migration rate M (m/μ; m, historical migration rate per generation), and the number of migrants N _m (N _e m) between the two geographic groups (Beerli, 2006). We used a uniform prior distribution to estimate Θ (rang 0–80), M (range 0–300), and reported the mean and 97.5% confidence intervals for Θ and M. The starting values for Θ and M were estimated using the F _ST. Monmonier’s maximum difference algorithm implemented in BARRIER v.2.2 (Manni et al., 2004) was used to identify potential barriers to gene flow among populations. Nei’s distance matrices, used in determining the barriers, were produced by bootstrapping over loci with Microsatellite Analyser v.4.05 (Dieringer and Schlötterer, 2003).

Genome size was assessed for 16 individuals of X.vietnamensis (Table 2). Relative nuclear DNA content of the individuals was estimated by propidium iodide (PI) flow cytometry using fresh leaves (Bourge et al., 2018; Farhat et al., 2021). We adopted the nuclei isolation buffer “Galbraith’s” (Galbraith et al., 1983), and samples of Cupressus funebris (2n = 22 and 2C = 21.52 pg DNA) (Hizume et al., 2001) and Pseudolarix amabilis (2n = 44 and 2C = 52.2 pg DNA) (Zonneveld, 2012) were selected as the internal calibration standards. For each sample, about 1 g of fresh leaves were placed into a glass Petri dish. Then leaves were chopped using a razor blade in 4 ml of 4°C cold buffer. The nuclear suspension was filtered through 37 or 48 μm nylon mesh. RNase A was added to the suspension to prevent staining of double-stranded RNA. Cell nuclei were stained with 100 μg/ml PI, and incubated in the dark for 30 min at 4°C before analysis by flow cytometry.

BD FACSCalibur™ Flow Cytometer (BD Biosciences, United States) equipped with a solid-state laser for PI excitation was used to detect DNA content (2C value). Three replicates were run per individual, recording at least 10,000 particles (including nuclei and fragments) per replicate. The data analysis was performed in BD CellQuest™ Pro (version 6.1). The resulting histograms were analyzed with ModFit LT 4.1.7 software. When CV values of G0/G1 peaks were below 5%, the analyses saved. The determination of ploidy level was based on the relative nuclear DNA contents and published chromosome counts (Farhat et al., 2021). The databases were accessed on August 11, 2021: Kew Plant DNA C-values database (http://data.kew.org/cvalues) and Index to Plant Chromosome Number (IPCN) (http://legacy.tropicos.org/Project/IPCN).

Twenty microsatellite markers were polymorphic in the X. vietnamensis populations analyzed in this study (Supplementary Table S1). No loci with excessive null allele and no loci under natural selection were detected. A total of 96 alleles were detected, ranging from two to ten alleles per locus, with an average of 4.8 (Table 3). The observed heterozygosity (H _O) and expected heterozygosity (H _E) per locus ranged from 0.156 to 0.940 and from 0.219 to 0.839, respectively. The inbreeding coefficient (F _IS) per locus ranged from −0.379 to 0.575. Five loci displayed significant deviations from Hardy-Weinberg equilibrium and 10% locus pairs showed evidence of linkage disequilibrium after sequential Bonferroni correction. All loci were retained for further analysis. The genetic differentiation between the geographic groups (F _ST) varied among loci, ranging from 0.008 to 0.450 (Table 3).

The genetic diversity estimates at the population and group levels are summarized in Table 4. N3 had the highest H _E (0.439), and N4 had the lowest H _E (0.373) at the population level. At the group level, the north had a slightly higher H _E (0.469) than the south (0.438). The M-ratio test for the X. vietnamensis populations detected a significant bottleneck (M ₀ = 0.7737; M _c = 0.7883 when the pre-bottleneck N _e = 4,000, p = 0.0326; M _c = 0.8693 when the pre-bottleneck N _e = 400, p = 0.0002).

The STRUCTURE results showed that the L(K) continued to increase up to K = 4, and the K = 2 model received the highest delta K value (Figure 2). Histogram showing the results of K = 2 model revealed clear differences between the northern and southern populations (Figure 3). When K = 3 and 4, the N4 population separated from other northern populations, resulting in a slight signature of separation within the northern populations (Figure 3). The BIC curve in the DAPC analysis suggested that the optimal number of clusters was five (Figure 4A). The southern populations (S1 and S2) were assigned to a single cluster while the northern populations were assigned to the other four (Figure 4C). Specifically, individuals of N3 were assigned to four different clusters, and N4 were assigned to two separate clusters (Figure 4C).

The genetic differentiation coefficient average across loci (F _ST = 0.125, SE = 0.025) suggested a moderate level of differentiation between the northern and southern population groups (Table 3). AMOVA showed more genetic variation was partitioned between the northern and southern population groups (6%) than among populations within groups (3%); most genetic variation was found within individuals (91%). Genetic differentiation between the northern and southern population groups (F _RT = 0.089) was statistically significant (p < 0.05), while the genetic differentiation among populations and among individuals were not statistically significant (Table 5).

The mean Θ values for northern and southern populations were 0.896 and 0.644, respectively. Asymmetric gene flow between north and south populations was detected, with the number of migrant individuals per generation (N _m) from south to north (0.770) being less than from north to south (2.008) (Table 6). The BARRIER analysis revealed a well-supported genetic and/or biogeographic barrier (with bootstrap support value of 75%) that separating the northern and southern populations (Figure 5).

The genome size estimates for the 16 individuals ranged 25.2–50.48 pg/2C. They could be distributed into two classes: Class 1 of lower values (25.42–26.50 pg/2C, mean = 26.08 pg/2C) corresponded to diploid with 2n = 2x = 22 for the individuals belonging to the northern populations, class 2 possessed higher values (46.59–50.48 pg/2C, mean = 48.02 pg/2C) corresponded to tetraploid with 2n = 4x = 44 for the individuals belonging to the southern populations (Table 2). Flow cytometric histograms obtained from analysis of two representative individuals, one each from the northern (N3-8) and southern (S1-2) populations, were shown in Figure 6.

Our population genetic survey of X.vietnamensis revealed that its populations from southwestern China exhibit lower diversity (mean value across all sampled populations: H _E = 0.45, H _O = 0.49) than the long-lived perennial (H _E = 0.68, H _O = 0.63), widespread (H _E = 0.62, H _O = 0.57), regional (H _E = 0.65, H _O = 0.65) and narrow range plant populations (H _E = 0.56, H _O = 0.52), but exhibit higher diversity than endemic plant populations (H _E = 0.42, H _O = 0.32) according to the mean values provided by Nybom (2004) based on microsatellite markers. Compared with other threatened species in Southern China and adjoining areas, the diversity of X. vietnamensis (H _E = 0.51, H _O = 0.45) was lower than that of C. macrolepis (H _E = 0.64, H _O = 0.71, 291 individuals) (Liao et al., 2015) and F. hodginsii (H _E = 0.64, H _O = 0.52, 427 individuals) (Yin et al., 2018), yet higher than that of Glyptostrobus pensilis (H _E = 0.27, H _O = 0.32, 343 individuals) (Wu et al., 2020), Cathaya argyrophylla (H _E = 0.41, H _O = 0.37, 49 individuals) (Wang et al., 2010) and Abies ziyuanensis (H _E = 0.44, H _O = 0.32, 139 individuals) (Tang et al., 2008). However, although the 33 individuals sampled across the six small populations in this study represent the majority of the known distribution of X. vietnamensis in southwestern China, limited sampling still may have led to an underestimation or inaccurate estimation for the genetic diversity of this species. Moreover, our loci number may be insufficient and the population indices in Table 4 should be treated with caution.

Some relict species endemic to South China were cold-adapted, for example, C. argyrophylla and A. ziyuanensis were restricted to a few isolated alpine habitats with cold microclimates. These species tend to have lower levels of genetic variation and are faced with a greater threat of losing their last refuge in South China in the face of global warming (Tang et al., 2008). In contrast, X. vietnamensis is a heliophile and subtropical or tropical plant (Averyanov et al., 2002), similar to C. macrolepis and Pinus kwangtungensis (Tian et al., 2008), there may be a wider range of potential habitats (including non-karst montane areas) in warm South China and adjacent regions in southeast Asia. But actually it is currently inhabiting isolated mountain tops, cliffs or slopes confining within karst areas and has a very small population size that might have been partly caused by historical deforestation (Zeng et al., 2007; Wang 2011; Jiang et al., 2014). At the same time, we detected a significant bottleneck event for X. vietnamensis populations in southwestern China, when assuming two different pre-bottleneck effective population size, suggesting that this species may have experienced a strong bottleneck event in the past. This provides reasonable explanation for the low-level genetic diversity of X. vietnamensis.

Our population genetic survey revealed moderate genetic differentiation (F_ST = 0.125) among X. vietnamensis populations in southwestern China. This is similar to the mean values of anemochores (0.13) (Nybom, 2004), lower than those of F.hodginsii (0.157) (Yin et al., 2018), C. macrolepis (0.163) (Liao et al., 2015), A. ziyuanensis (0.209–0.250) (Tang et al., 2008), and G. pensilis (0.452, some populations are plantations) (Wu et al., 2020) in a similar geographic region. Varied topography of mountain regions in the tropic and adjacent areas tended to harbor Tertiary relict species and subdivide some species into isolated and divergent populations, which may have evolved independently during late Tertiary climate changes and Quaternary climate fluctuations (Hewitt, 2000; Tang et al., 2018). The karst landscape in Guangxi is a residual carbonate hill protruding from a surrounding corrosion plain (Gutiérrez and Gutiérrez, 2016). The hills are generally only hundreds of meters in height, and the nearby hill tops are generally hundreds of meters apart. Even so, the karst forest ecosystem here exhibits heterogeneous geomorphology and variegated vegetation (Su et al., 2017; Guo et al., 2018). The heterogeneous geomorphology may have contributed the genetic differentiation among populations of X. vietnamensis.

At the same time, we detected significant population genetic structure in X. vietnamensis, as both STRUCTURE and DAPC analyses clearly revealed a distinct division between the northern and southern population. In addition, a relatively strong barrier to gene flow between the northern and southern populations was identified by BARRIER analysis, and the estimated N _m from south to north or vice versa (Table 6) was limited, being consistent with the moderate genetic differentiation (Slarkin, 1985). This is supported by the fact that, there are non-karst areas between or among the karst landforms, which bear different ecological environments and interrupt the continuity of karst flora (Chen et al., 2012). For example, the Youjiang Basin, which mainly consists of broad non-karst terrains (Lu et al., 2006), is located at the midpoint between the northern and southern populations of X. vietnamensis (Figure 1). The Youjiang Basin and other non-karst areas between northern and southern population of X. vietnamensis might have acted as genetic barriers and resulted in isolated and fragmented populations, and restricted gene flow and migration between the northern and southern populations of X. vietnamensis.

Generally, bottleneck effects are found in populations with reduced population size due to natural catastrophes, habitat loss, alteration and fragmentation. As such, it is likely a historical climate event in combination with the fragmented habitat in this karst landscape may have led to the significant bottleneck detected in the X. vietnamensis populations. Meanwhile, in the context of demographic bottleneck, the isolation and fragmentation could have serious negative impact on the X. vietnamensis population due to random drift and inbreeding, and may have promoted genetic differentiation among populations (Nei et al., 1975; Sękiewicz et al., 2015).

Our genome size data showed that X. vietnamensis is mixed-ploidy and populations in Guangxi possess both diploid and tetraploid cytotypes. The average tetraploid genome size of 48.02 pg/2C was similar to but slightly larger than previous study (44.60 pg/2C) using another internal standard and nuclei isolation buffers (Farhat et al., 2021); the average diploid genome size of 26.08 pg/2C is similar to that of diploid Juniperus species (Farhat et al., 2019). Although only four representative individuals per population were sampled from two populations each from northern and southern population groups, we found that all eight individuals from northern populations are diploids and these from southern populations are tetraploids.

While there was significant difference in genome size between the northern (diploid) and southern (tetraploid) populations, microsatellite diversity varied little between them. These differences may reflect the recent origin of tetraploid from the diploid (Lumaret and Barrientos, 1990; Ickert-Bond et al., 2020). However, these two cytotypes might have distinctly different genetic structure (Figures 3, 4), which was probably due to habitat differences and spatial isolation (Laport et al., 2016). Tetraploid populations appear to have greater negative inbreeding coefficients (F _IS) (Table 4). We tested the impact of small population size on F _IS (Supplementary Table S2), and found a trend of negative F _IS with greater absolute values when sample size is small, but smaller size did not cause a shift of F _IS from positive value to negative value or the other way round. These results suggest seedling production by selfing in tetraploid populations of X. vietnamensis. As reported in previous studies, chromosome doubling can lead to a breakdown of self-incompatibility (Soltis and Soltis, 2000; Sutherland et al., 2018). Thus, tetraploid populations may be more likely to establish after expansion or dispersal than diploid populations. As far as we know, many more tetraploid individuals and populations were found in the neighboring areas——Vietnam (Averyanov et al., 2013; Farhat et al., 2021). Our data revealed historical gene flow from north (diploid population) to south (tetraploid population) was higher than the reverse (Table 6). It is likely that tetraploid X.vietnamensis in southern population may be of autopolyploid origin and diploids in the northern populations may be the source populations, yet these hypotheses need to be further tested.

There is clear genetic and ploidy level distinction between the north and south populations in southwestern China according to the STRUCTURE analysis and genome size analysis. The disjunct geography and topography of the karst mountains and over 300 km of geographical distance between the northern and southern populations may have acted as barriers to gene flow. Hence the northern and southern populations should be treated as two different management units. Populations of X. vietnamensis are isolated as their habitat is confined to the mountain peaks. To reduce inbreeding and minimize bottleneck effect, artificial pollination or transplanting among fragmented populations in each unit should be implemented. Considering the rarity of wild individuals at present, we suggest in situ artificial regeneration programmes immediately for the populations, and urgent trials of ex situ plantations in Guangxi, China, through the collection, germination, cultivation of seeds and twig cuttage, similar to the strategy that has proven successful in Vietnam for this species.

In our study, sampled populations were restricted to provincial or national borders, remote from the human residential belt. The regions between the northern and southern populations are substantially impacted by humans (Zeng et al., 2007; Wang, 2011; Jiang et al., 2014). These imply that human activities led to the loss of habitats and the extremely small population size of X. vietnamensis. The northern populations are mostly located in the Guangxi Mulun National Nature Reserve, and the southern populations are located in the Guangxi Laohutiao Provincial Nature Reserve. In general, management effectiveness of nature reserves in China are relatively good, especially those with higher management level; but there are still some problems of nature reserves in China, such as the lacking information of boundaries and tenures, the poorly organized institution, the lacking of staff, fund and so on (Wang and Li, 2021). The Laohutiao Provincial Nature Reserve stretches 69 km along the border of China and Vietnam. Generally, it is not easy to effectively protect the ecology and biodiversity in the national borders (Wang and Li, 2021). Unfortunately, the distribution area of the southern populations of X. vietnamensis was accessible for agricultural exploitation, as we saw during the field survey. We recommend more attention and tighter management for the conservation of forest trees in this reserve.