Zhi Yao1,2Xinyu Wang1Kailai Wang1Wenhao Yu1Purong Deng1Jinyi Dong1,2Yonghua Li1Kaifeng Cui3
Yongbo Liu1*- 1State Environmental Protection Key Laboratory of Regional Eco-Process and Function Assessment, Chinese Research Academy of Environmental Sciences, Beijing, China
- 2Hunan Provincial Collaborative Innovation Center for Field Weeds Control, Hunan University of Humanities, Science and Technology, Loudi, China
- 3Changbai Mountain Academy of Sciences, Joint Key Laboratory of Community and Biodiversity for Jilin Province and Changbai Mountain, Jilin, China
Narrow-ranged species face challenges from natural disasters and human activities, and to address why species distributes only in a limited region is of great significance. Here we investigated the genetic diversity, gene flow, and genetic differentiation in six wild and three cultivated populations of Thuja sutchuenensis, a species that survive only in the Daba mountain chain, using chloroplast simple sequence repeats (cpSSR) and nuclear restriction site-associated DNA sequencing (nRAD-seq). Wild T. sutchuenensis populations were from a common ancestral population at 203 ka, indicating they reached the Daba mountain chain before the start of population contraction at the Last Interglacial (LIG, ∼120–140 ka). T. sutchuenensis populations showed relatively high chloroplast but low nuclear genetic diversity. The genetic differentiation of nRAD-seq in any pairwise comparisons were low, while the cpSSR genetic differentiation values varied with pairwise comparisons of populations. High gene flow and low genetic differentiation resulted in a weak isolation-by-distance effect. The genetic diversity and differentiation of T. sutchuenensis explained its survival in the Daba mountain chain, while its narrow ecological niche from the relatively isolated and unique environment in the “refugia” limited its distribution.
1 Introduction
In the context of biodiversity loss worldwide, it is critical to conserve endangered species that are generally narrow ranged (Chau et al., 2013). The plant species with extremely small populations (PSESP) is proposed as a conservation concept to protect them (Yang et al., 2020). These species face extinction risk at any time because their population size is below a stable survival limit and suffers long-term interference and stress from external factors (Yang et al., 2020). The basis for species to adapt to environmental changes depends on the genetic diversity and differentiation, and even distribution area of populations (Jordan et al., 2019). Besides determining the ability to adapt to the environment, genetic diversity is the basis for maintaining long-term stability of the entire ecosystem (Hamrick et al., 1992). Species with small population size or small geographic ranges generally tend to result in low genetic diversity due to genetic drift, reduced gene flow, and increased genetic divergence (Ellstrand and Elam, 1993), which result in a low adaptability, reproductive ability, and disease resistance of species (Hamrick et al., 1992). Although there has been a debate, low genetic diversity generally increases extinction risk (Saccheri et al., 1998). Genetic variations or polymorphisms contribute to the viability and evolutionary potential of natural populations (Frankham and Briscoe, 2010); therefore, the assessment of genetic diversity is crucial in understanding the evolutionary history of endangered species and designing effective conservation and management methods (Cohen et al., 1991; Allendorf and Luikart, 2009; Ouborg, 2010).
The population demographic history is generally influenced by tectonic events and climatic change that, therefore, are key factors influencing patterns of genetic diversity (Sun et al., 2015). These tectonic and climatic events cause contraction–expansion of effective population size, thereby forming the current genetic structure (Bai et al., 2018). The extremity, variability, or stress of environmental factors can exacerbate not only population size but also latent genetic issues such as inbreeding, accumulated genetic load, and genetic variation (Fox and Reed, 2011). Thus, genetic diversity and structure of a species are indicators of life forms, reproductive systems, seed transmission mechanisms, evolutionary history, climatic factors, and human interference (Hamrick and Godt, 1996).
Here, we study the genetic diversity and structure of Thuja sutchuenensis Franch, an endangered species endemic to the Daba Mountains in northwestern Chongqing Municipality and eastern Sichuan Province, China (Xiang et al., 2002; Tang et al., 2015). Cooling and drying climate in the Quaternary forced T. sutchuenensis southward and restricted in the Daba mountains (Cui et al., 2015). When did T. sutchuenensis arrive at the Daba mountains, and can it survive only in this restricted region? Tang et al. (2015) found T. sutchuenensis populations distributed on limestone cliffs or steep slopes where there is lack of competition from other species. To address the endangered mechanisms of T. sutchuenensis and its limited distribution range in the Daba mountains as a “refugia,” wild and cultivated populations were genotyped using chloroplast simple sequence repeats (cpSSR) and nuclear restriction-site associated DNA sequence (nRAD-seq). Previous studies showed inconsistent results for T. sutchuenensis populations from different methods, high genetic diversity with ISSR (Liu et al., 2013), and low genetic diversity with six single-copy nuclear loci (Qin et al., 2020).
2 Material and Methods
2.1 Plant Materials Sampling
Thuja sutchuenensis, a monoecious and evergreen coniferous tree, distributed in the Chongqing Municipality and Sichuan Province, Southwest China (Xiang et al., 2002; Cui et al., 2015) (Figure 1). We sampled 375 plants from six wild and three cultivated T. sutchuenensis populations, and recorded longitude, latitude, and altitude of the T. sutchuenensis plants (Table 1). We sampled individuals around 5 m from each other whenever possible. A total of ∼50 g fresh leaves per plant of T. sutchuenensis was sampled and dried in allochronic silica gel.
FIGURE 1. The distribution area and sampling sites of Thuja sutchuenensis in China. (A) and (B), The locations of current T. sutchuenensis populations and T. sutchuenensis fossils in China. (C), A picture of T. sutchuenensis fossils, which were found in Zhangcun, Shanxi province (36°58′N; 112°51′E), was quoted from Cui (Cui et al., 2015). (D), Sampling populations of T. sutchuenensis in the Daba Mountains wild populations, with MZ and XY from the Chengkou County, and MY, GM, and BQ from the Kaixian County, Chongqing Municipality, as well as QS from the Xuanhan County, Sichuan Province.
TABLE 1. Wild and cultivated Thuja sutchuenensis populations in the Daba Mountains.
Thuja koraiensis (CYB) and Platycladus orientalis (ZYB), as two outgroups, were sampled from the Changbai Mountains and the Institute of Botany, CAS, respectively.
2.2 Chloroplast Simple Sequence Repeat Identification and Marker Design
Whole-chloroplast genome DNA of T. sutchuenensis was extracted from 25 mg of dry leaves using a Hi-DNA Secure Plant Kit (TIANGEN Biotech Co., Ltd.) based on a modified version of the SDS method (Yuan et al., 2012). The motifs containing the SSR loci were identified using MISA (http://pgrc.ipk-gatersleben.de/misa/). The motifs of trinucleotide, tetranucleotide, pentanucleotide, and hexanucleotide repeats included at least 5, 4, 4, and 4 repeats, respectively. SSR primers were designed using Primer 3.0 (https://sourceforge.net/projects/primer3/). The SSR primers were designed as follows: 18–24 bp in length, 40%–60% GC content, 55°C–65°C annealing temperature, and 100–300 bp PCR product. Furthermore, 400 SSR loci were selected and primers of 140 loci were screened in 20 individuals. Finally, 30 microsatellite repeats (SSR) were selected from the 140 loci to detect the polymorphic site of T. sutchuenensis populations (Supplementary Table S1), and the primers were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). Forward primers were labeled with four kinds of 5′-fluorescein bases (FAM, HEX, TRAMA, or ROX). The multiplex PCR (Supplementary Table S2) was run on a BIO-RAD T100™ Thermal Cycler. Detection of multiplex PCR products was carried out with an ABI 3730xl DNA analyzer (Sangon Biotech Co. Ltd.), setting a default range of standard length ±40 bp. The detection bands of 30 markers were scored using GeneMarker version 1.91 (size standard: GS500).
2.3 Nuclear Restriction Site-Associated DNA Sequencing
A total of 100 ng genomic DNA per sample was digested with two restriction enzymes, EcoRI, and PstI (New England Biolabs, Beverly, USA) at 37°C for 8 h. The restriction enzymes were then inactivated by heating at 65°C for 20 min. After ligation individually with barcoded EcoRI adapter and universal PstI adapter with T4 DNA ligase for each sample at 16°C for 8 h, the reaction was stopped by heating at 65°C for 20 min. The ligation products of all samples were equally pooled and size selected into 300–500 bp fragments using the agarose gel electrophoresis. After manipulating gel purification, derived fragments were used as templates (about 30 ng) for PCR amplification via 25 cycles with EcoRI and PstI adapter universal primers using PrimeStar Max DNA Polymerase (Takara, Dalian, China). Finally, the amplicons were size selected once more into 350–500 bp fragments with the method mentioned above. The resulting ddRAD library was sent to Guangzhou Jierui Biotechnology Company, LTD. (Guangzhou, China) and sequenced on the Illumina Novaseq platform using 150 nt with paired-end mode.
We used process_radtags module in the Stacks-2.4 (Catchen et al., 2011) to demultiplex the raw data, with default parameters. The reads were trimmed to 135-bp length to remove low-quality nucleotides at the 3′ end of each read. Finally, each end of the retained reads was treated as an independent locus and combined together. The ustacks module was employed to cluster the reads into exactly matching stacks, setting m = 2 as the minimum depth of coverage (m) and M = 12 as the maximum distance allowed between stacks within an individual. We then used the cstacks module to build the catalogs for all samples with n = 12 as the maximum number of mismatches allowed between individuals. The Sstacks module generated alignment results for each individual against the catalog using default parameters. In the populations module, we set p = 8 and r = 0.6 to call consensus SNP among all 137 samples (15 individuals per population).
2.4 Statistical Analysis
We employed Micro-Checker 2.2 to detect null alleles, scoring errors and allele dropout (Oosterhout et al., 2010), and no scoring errors or allele dropout were detected. Linkage disequilibrium and Hardy–Weinberg equilibrium were tested in Arlequin v. 3.5 (Excoffier and Lischer, 2010). Three of the 30 loci were null alleles, and a pattern of linkage disequilibrium was found at seven loci (p > 0.05) (Supplementary Table S3). Genetic diversity indices of cpSSR, including the number of alleles (Na), effective number of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (FIS), were calculated for all samples from nine populations by GenAlEx version 6.503 (Peakall and Smouse, 2006).
Genetic diversity indices of nRAD-seq, including the inbreeding coefficient (FIS), nucleotide diversity (Pi), observed heterozygosity (Ho), and expected heterozygosity (He), were calculated using the populations module in Stacks (Catchen et al., 2011). A variant call format (VCF) file was generated using the populations program in Stacks (Catchen et al., 2011). A maximum likelihood (ML) tree of six wild and three cultivated populations (90 individuals) was constructed with 1,000 bootstraps using the SNPhylo software (Tae-Ho et al., 2014).
We separately calculated genetic differentiation and gene flow using GenAlEx for cpSSR data (Peakall and Smouse, 2006). The genetic differentiation coefficient (FST) was calculated by analysis of molecular variance (AMOVA) (Excoffier et al., 1992). The significance of FST values was tested with 999 permutations, followed by Bonferroni correction. Gene flow (Nm) was estimated based on the formula
To test whether genetic differentiation related with geographic distance, we performed Mantel test with 999 permutations using FST/(1-FST) and logarithmic geographic distance (Smouse, 1986). The pairwise matrix of geographic distances was calculated based on geographic coordinates, which were identified by a GPS device (Garmin Oregon 450). All calculations above were performed by GenAlEx (Peakall and Smouse, 2006).
To infer population structure, Bayesian clustering analysis was conducted in Structure version 2.3.4 (Pritchard et al., 2000). We set the number of groups (K) from 1 to 9 with 10 independent runs using an admixture ancestry model and 1,000,000-step Markov chain Monte Carlo (MCMC) replicates after a 500,000-step burn-in for each run. The best K value was inferred by delta K in Structure Harvester (Earl and Vonholdt, 2012). A cluster analysis based on genetic distance was performed by the unweighted pair-group method with arithmetic mean (UPGMA) approach in MEGA 6.0 (Tamura et al., 2013). We used the TreeMix program to detect historical migration among wild populations (Pickrell and Pritchard, 2012).
BOTTLENECK was used to examine recent bottleneck events (Cornuet and Luikart, 1996), under the TPM (two phase model), which allows multiple-step mutations (Peterson et al., 1994), with a 10% infinite allele mutation, 90% stepwise mutation model, and 1,000 replicates (Wang et al., 2014). The significance of heterozygote excess was determined by Wilcoxon’s signed-rank test (p < 0.05). Departures from mutation–drift equilibrium were detected by a mode-shift test (Luikart et al., 1998).
To further investigate the demographic history of wild populations, we performed coalescent simulations of 15 evolutionary models (Supplementary Figure S1) in fastsimcoal 2 v2.5.2.2.21 (Excoffier and Foll, 2011). Based on the results of structure analysis, wild populations were clustered in three groups (best K = 3), Pop 1 (BQ and XY), Pop 2 (MZ and QS), and Pop 3 (MY and GM) (Figure 2; Supplementary Figure S2). The fastsimcoal software simulates these data following the Wright–Fisher model of evolution assuming neutrality of the genetic markers, without recombination within loci but free recombination between loci and random mating. We improved the performance of the models by reducing the number of estimated parameters (Excoffier et al., 2013), with the population parameters were calculated directly from the SNP data. Each run was performed with 200,000 simulations per likelihood estimation and 40 expectation–conditional maximization (ECM) cycles. The mutation rate was set as 1.0 × 10−8.
FIGURE 2. Genetic structure and phylogenetic tree of six wild T. sutchuenensis populations. (A) Population structure bar plots from nuclear restriction site-associated DNA sequencing (nRAD-seq) data show the clustering of samples into three clusters (best delta-K). Each vertical bar indicates an individual, and the height of each colored bar represents the proportion of assignment to that cluster. For population abbreviations, please see Figure 1. (B) Maximum-likelihood (ML) tree illustrates genetic relationships between six wild T. sutchuenensis populations based on SNPs from nRAD-seq. One color indicates one population. ZYB and CYB are outgroups.
To test the local adaptations of T. sutchuenensis, we used two genotype–environment associations analysis to detect loci related with environmental factors. One approach is OutFLANK, inferring the distribution of FST to detect loci responsible for local adaptation based on the expected Chi-square (χ2) distribution of FST in the absence of selection (Whitlock and Lotterhos, 2015). We clipped the FST distribution to 5% and 10% using the minimum expected heterozygosity of 0.10 and an error discovery rate of less than 0.05. The second approach is using latent factor mixed model (LFMM) (Frichot et al., 2013), to identify the relationship between genetic polymorphism of single nucleotide and environmental factors, using the locus GIF (genomic inflation factor) associated with the environment. The “VCF2LFMM” algorithm in the package was used to convert the nucleotide diversity files in the VCF format into the 0,1 format files of LFMM (Frichot et al., 2013). Extracted values of bioclimatic factors corresponding to species distribution points were used as environmental variables, and LFMM was used to convert environmental variables into the “env” format. Using 100,000 scans, “Burn-In” was set to 50,000 as recommended. In order to effectively control false discovery rate (FDR), we selected the K value whose estimated value of λ is close to or slightly lower than 1.0. Since the test of K values in a wide range is usually ineffective, we chose K = 1–10. Sites associated with the environment were finally picked out. The climatic data were downloaded from the WorldClim database, including 19 bioclimatic variables (bio1 to bio19) with a spatial resolution of 30 arc-seconds (http://worldclim.org) (Supplementary Table S4).
3 Results
3.1 Genetic Diversity
After all quality filters, 62,375,634 reads with no rad tag and 1,000,938 low-quality reads were discarded. We retained a total of 3,509,088,912 reads from the initial 3,609,108,036 raw reads, with an average of 2.63 Mb reads per sample. After trimming and clustering, we obtained loci number ranging at 6,933–395,804 bp and mean locus coverage across all samples was 72.64X from 32.69 to 135.05X, with an average length of 135 bp per loci and with 44.92% GC content. After the “cstacks” module processing, we obtained 4,886,198 catalogs and 6,119 SNPs.
Out of the 5,990 polymorphic SNPs loci, 2 (14063_15 and 3630303_6) and 1 (3630303_6) SNPs were considered as responsible for local adaptation by the OutFLANK and LFMM approaches, respectively. Because there is no existing reference genome for Cupressaceae, it is impossible to annotate these selected loci.
T. sutchuenensis populations showed high chloroplast genetic diversity (cpSSR) and low nuclear genetic diversity (nRAD-seq). The Ho and He of cpSSR in T. sutchuenensis populations were 0.43–0.92 and 0.23–0.53, respectively (Table 2). The Ho and He of cpSSR in BQ and QS populations were higher than in other wild populations (Table 2), while the nRAD-seq genetic diversity of T. sutchuenensis among populations were similar (Ho = 0.061–0.070, He = 0.078–0.087) (Table 3). The nucleotide diversity (pi) of T. sutchuenensis populations ranged from 0.082 to 0.092 for variant position and from 0.0053 to 0.0062 for all positions (Table 3).
TABLE 2. The characteristics of chloroplast genetic diversity (cpSSR) in Thuja sutchuenensis populations.
TABLE 3. The characteristics of nuclear genetic diversity (nRAD-seq) in T. sutchuenensis populations.
The Ho and He of cpSSR in cultivated populations BQR and GMR were comparable with wild populations BQ and QS (Table 2). The high absolute values of cpSSR FIS suggested a significant homozygous excess of T. sutchuenensis populations (Table 2), while the FIS values of nRAD-seq were low (Table 3). The tested loci were not deviated from the Hardy–Weinberg equilibrium (HWE) (Supplementary Table S5). The wild QS population and the cultivated GMR and BQR populations experienced recent bottleneck events because their loci departures from the mutation–drift equilibrium in the allele frequency distribution (Supplementary Table S6).
3.2 Genetic Differentiation and Population Structure
The cpSSR genetic differentiation varied with the pairwise comparisons of populations (FST = 0.01–0.19; Nm = 1.12–22.48) (Table 4). The genetic differentiation (FST) in pairwise comparisons of BQ-QS and MZ-MY-XY was low, with high gene flow (Nm = 7.13–10.62) (Table 4). The cultivated GMR and BQR populations and wild BQ-QS populations showed low pairwise FST values (FST = 0.01–0.03) and high gene flow (Nm = 5.56–22.48) (Table 4). The cultivated XYR population showed high cpSSR FST and low Nm from other populations (Table 4). However, the genetic differentiation of nRAD-seq in any pairwise comparisons were low (FST = 0.033–0.055) with high gene flow (Nm = 4.26–8.32) (Table 4). This indicates that T. sutchuenensis populations were genetically similar (Ewens, 1970), which was also supported by the analysis of molecular variance (AMOVA) that showed 81% of genetic variation within populations and 19% among populations (Supplementary Table S7).
TABLE 4. Matrix of pairwise FST and Nm coefficients of T. sutchuenensis populations from chloroplast simple sequence repeats (cpSSR) and nRAD-seq (values in parentheses).
The correlation between genetic and geographic distance was marginally significant (R2 = 0.029, p = 0.08) (Supplementary Figure S3), indicating a weak effect of isolation by distance.
The structure analysis of cpSSR showed two clusters (best delta K = 2) (Supplementary Figure S4), with one cluster of BQ, QS, GMR, and BQR populations and others clustered together. When K = 4 (the second top delta K) (Supplementary Figure S4), six wild populations clustered three groups and the population XYR cultivated from seeds 15 years ago differed from other populations. The nRAD-seq structure for the nine populations indicated gene flow between populations or originated from a common ancestor (Supplementary Figure S5). The structure analysis of nRAD-seq for six wild populations showed three clusters (best delta K = 3) (Figure 2A; Supplementary Figure S2). The results of UPGMA (Supplementary Figure S6), phylogenetic tree (Figure 2B), principal coordinate analysis (PCoA) (Supplementary Figure S7A), and principal component analysis (PCA) (Supplementary Figure S7B) confirmed the structure results from cpSSR and nRAD-seq, respectively.
3.3 Demographic History
Fastsimcoal results from nRAD-seq showed that the wild populations Pop 1 (BQ and XY) and Pop 3 (MY and GM) were from an ancestral population at ∼203 ka (Figure 2). An admixture event occurred at ∼12 ka between the Pop 1 and Pop 3 and formed the Pop 2 (MZ and QS) (Figure 2A, Figure 3; Supplementary Figure S2).
FIGURE 3. Demographic scenario tested in Fastsimcoal. Scenarios were demographic scenarios that were set to infer population history. Effective population size and divergence time for the three clusters (K = 3; Figure 2) for six wild populations using nRAD-seq data (see Supplementary Figure S1), i.e., BQ and XY populations (Pop 1, red line), MZ and QS populations (Pop 2, blue line), and GM and MY populations (Pop 3, green line). t2, divergence time between Pop 1 and Pop 3; t1, divergence time between Pop 1 and Pop 2.
TreeMix analysis of cpSSR revealed historic migrations between QS, BQ, and GM populations (Figure 4A; Supplementary Figures S8A, S9), while nRAD-seq data showed two migration events, from GM to MZ and from QS to XY (Figure 4B; Supplementary Figures S8B and S10).
FIGURE 4. TreeMix results showing historical mitigation between the wild T. sutchuenensis populations. (A), TreeMix results from cpSSR. (B), TreeMix results from nRAD-seq. For population abbreviations, please see Figure 1.
4 Discussion
Assessing the genetic diversity and genetic structure of rare, endemic, and endangered species, i.e., narrow-ranged species, is not only the basis for developing scientific and rational conservation strategies (Joshi and Dhawan, 2007), but also an explanation of the limited distribution of the narrow-ranged species, such as wild plants with extremely small populations (PSESP) (Yang et al., 2017). We assessed the wild populations of T. sutchuenensis using chloroplast simple sequence repeats (cpSSR) and nuclear restriction site-associated DNA sequencing (nRAD-seq) and found relatively high chloroplast genetic diversity but low nuclear genetic diversity, which provides molecular explanation in its distribution in the Daba mountains as a “refugia”.
T. sutchuenensis was widely distributed in China in the late Pliocene based on its fossil in the Shanxi province (Figure 1C) (Cui et al., 2015), and its distributed area and population size decreased from the Last Interglacial (LIG, ∼120–140 ka) to the Last Glacial Maximum (LGM, ∼22 ka) (Qin et al., 2017). Our results found that T. sutchuenensis wild populations were from a common ancestral population at 203 ka, which means that T. sutchuenensis southward-migrated and reached the Daba mountains before the population concentration of T. sutchuenensis. MaxEnt results from Qin et al. (2017) found that the current distribution size is smaller than the past (LIG, LGM) and the future. One related species, Cupressus chengiana, in the eastern Qinghai–Tibet plateau experienced a population contraction event during the late Quaternary (Xu et al., 2010). Rare and relict plant species such as Ginkgo biloba, Davidia involucrate, Liriodendron chinense, Camptotheca acuminate, and Taxus chinensis are found in the distribution area of T. sutchuenensis (Ma et al., 2017). Cool and dry climate in the Late Pliocene (Li et al., 2004), acting as a vicariance driver, forced a southward migration of Thuja in Shanxi province and survived in the Daba mountains (Cui et al., 2015). The relatively stable habitats within the Qinling–Daba Mountains have served as a refugia for more than 1,620 endemic plants in China (Ying, 1994), which have been supported by other plant species such as Cathaya argyrophylla (Wang and Ge, 2006), Saruma henryi (Zhou et al., 2010), and tree peony species (Xu et al., 2019).
Narrow-ranged species generally exhibited lower genetic diversity than widespread congeners (Cole, 2003). In this study, however, we found the cpSSR genetic diversity of T. sutchuenensis was not lower (He = 0.23–0.53) than other narrow-ranged species such as the Camellia huana (He = 0.339–0.605 for cpSSR, Ho = 0.353–0.605 for cpSSR) (Li et al., 2020), or widespread species, for example, Pinus kesiya (He = 0.654 for cpSSR) (Rai and Ginwal, 2018) and Ambrosia trifida (He = 0.203–0.645 for cpSSR) (Sharma et al., 2020). T. sutchuenensis have comparable genetic diversity with the other two congener species, T. koraiensis (SSR He = 0.548) and T. occidentalis (SSR He = 0.574–0.624) (Xu et al., 2013; Hou et al., 2018). This relatively high cpSSR genetic diversity of T. sutchuenensis might be due to current individuals from large populations that existed before in northern China, which indicates the likely effects of evolutionary history (from Northern China to the Daba Mountains) on genetic diversity. An alternative explanation is that chloroplast genomic inheritance is conservative and maintain genetic material, and chloroplast genome is mainly paternally inherited in Cupressaceae (Sakaguchi et al., 2014).
However, the nucleotide genetic diversity (pi = 0.082–0.091 for variant positions, pi = 0.0053–0.0062 for all positions) of T. sutchuenensis in this study is relatively low compared with other narrow-ranged species, such as Tetraena mongolica (He = 0.337–0.341; pi = 0.349–0.354) based on nRAD-seq (Cheng et al., 2020), or widespread species, such as Amorphophallus paeoniifolius (pi = 0.072–0.285) (Gao et al., 2017). Nuclear genetic diversity of T. sutchuenensis in our study was higher than six single-copy nuclear loci (pi = 0.00219) (Qin et al., 2020). Qin et al. (2020) reviewed 12 species of Cupressaceae and the nuclear genetic diversity was 0.00156–0.00989 with 6–13 loci (see Table 4 of Qin et al., 2020). Thus, methods (e.g., SSR, ISSR, and RAD) and the number of loci used to calculate genetic diversity influence the values of genetic diversity. For example, endemic species showed lower genetic diversity (Ho, He) within populations than regional or widespread species using STMS (sequence-tagged microsatellite sites)-based data, while there is no difference between them using RAPD (random amplified polymorphic DNA)-based data (Nybom, 2004). In addition, the genetic diversity of different species depends on their life form, breeding system, and geographic range (Nybom, 2004).
Generally, low genetic diversity result from high genetic differentiation and low gene flow, but the genetic differentiation was low between T. sutchuenensis populations (FST < 0.2), with high gene flow. This might result from the low population size of T. sutchuenensis (∼10,000 wild individuals in the nature). It is consistent with other narrow-ranged species, for example, the Cupressaceae species Platycladus orientalis (FST = 0.024–0.386) (Chang et al., 2019), endangered Camellia huana (FST = 0.2159 for cpSSR) (Li et al., 2020), threatened Rhododendron cyanocarpum (FST = 0.0314–0.0452) (Liu et al., 2020), and Tetraena mongolica (FST = 0.010–0.026) (Cheng et al., 2020). The genetic differentiation between populations was strongly influenced by genetic drift, gene flow, long-term evolution, mating systems, selection, and mutation (Hamrick and Godt, 1996). Due to gymnosperm species having long generation times, they usually show low genetic differentiation from wind pollination-mediated gene flow and breeding systems (Hamrick et al., 1992). T. sutchuenensis is a wind-pollinated hermaphrodite plant, and the low genetic differentiation resulted from the long-distance gene dispersal by pollen grains or seeds, because the extant T. sutchuenensis populations were geographically close. TreeMix results of nRAD-seq showed migration between close populations (QS and XY), and genetic and geographic distance indicating a weak effect of isolation by distance.
In fact, population structure results based on nRAD-seq of three cultivated and six wild populations failed to reveal any meaningful cluster patterns, which consisted with a previous study (Qin et al., 2020). However, population structure from six wild populations showed three clusters but with extensive gene flow and introgression between clusters. TreeMix results of cpSSR and nRAD-seq showed historic migrations between these wild populations, which is likely with the construction of roads in the valley between mountains (Figure 1D). In Cupressaceae, chloroplast genome is mainly paternally inherited, and compared with microsatellite regions, the low mutation rate of nRAD-seq loci can reveal the ancestral pattern of genetic structure (Andrews et al., 2016). This indicates the importance of constructing corridors in the conservation of genetic diversity (Yu et al., 2020).
In addition, the reintroduction of endangered species is an effective pathway to conserve genetic diversity. It has been demonstrated that reintroduction success can be enhanced by using plant material from large and stable source populations (Godefroid et al., 2011). Structure analysis can find the origination of reintroduced populations, we found that the reintroduced GMR and BQR populations cultivated by a cottage method were from the wild QS population. Interestingly, the reintroduced XYR population developed from the seeds of the XY population 15 years ago did not cluster with other any wild populations. Assessing the genetic diversity and genetic structure of endemic and endangered species is a fundamental criterion in developing scientific and rational conservation strategies. The reintroduction of T. sutchuenensis should collect plant material from populations with high genetic diversity such as the QS and BQ populations. In these populations, plant material from different age and size categories should be collected wherever possible in order to sample representative genetic diversity (Lauterbach, 2013).
In conclusion, we found the wild populations of T. sutchuenensis in the Daba mountains, i.e., its entire distributed region, were from a common ancestral population at 203 ka, which indicates T. sutchuenensis have colonized the Daba mountains before contraction at the Last Interglacial (LIG, ∼120–140 ka). The relatively high chloroplast and low nuclear genetic diversity of T. sutchuenensis populations might explain why it can and only survive in the Daba mountains as a “refuge.” After long-term (∼203 ka) adaptation, T. sutchuenensis prefers to grow on limestone cliffs, crest ridges, or steep slopes with special climate, soil (pH = 7.53), and sunshine (Ma et al., 2017), which leads to a narrow ecological niche of T. sutchuenensis and further to a low nuclear genetic diversity. These restrict its distribution only in the Daba mountains, although the seeds of T. sutchuenensis are wind-disseminated with membranous wings that facilitate population spread. Thus, enforcing the reintroduction of plant material from populations with high genetic diversity facilitates the conservation of T. sutchuenensis.
Data Availability Statement
The original contributions presented in the study are publicly available. These data can be found here: PRJNA777040.
Author Contributions
YL designed the research. ZY, XW, KW, WY, PD, JD, YL, KC, and YL performed the experiments. ZY, KW, and YL analyzed the data. ZY and YL wrote the paper.
Funding
This work was financially supported by the Biodiversity Survey, Observation, and Assessment of the Ministry of Ecology and Environment, China (2019HJ2096001006) and the Open Foundation of Key Laboratory of Forest Plant Ecology, Ministry of Education, Northeast Forestry University (K2020B01).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank Mr. Shiqiang Zhang and local farmers for their assistance in field sampling.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.801229/full#supplementary-material
References
Allendorf, F. W., and Luikart, G. (2009). Conservation and the Genetics of Populations. John Wiley & Sons.
Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G., and Hohenlohe, P. A. (2016). Harnessing the Power of RADseq for Ecological and Evolutionary Genomics. Nat. Rev. Genet. 17, 81–92. doi:10.1038/nrg.2015.28
Bai, W.-N., Yan, P.-C., Zhang, B.-W., Woeste, K. E., Lin, K., and Zhang, D.-Y. (2018). Demographically Idiosyncratic Responses to Climate Change and Rapid Pleistocene Diversification of the walnut Genus Juglans (Juglandaceae) Revealed by Whole-Genome Sequences. New Phytol. 217, 1726–1736. doi:10.1111/nph.14917
Catchen, J. M., Amores, A., Hohenlohe, P., Cresko, W., and Postlethwait, J. H. (2011). Stacks: Building and Genotyping Loci De Novo from Short-Read Sequences. G3 (Bethesda) 1, 171–182. doi:10.1534/g3.111.000240
Chang, E., Tian, Y., Wang, C., Deng, N., Jiang, Z., and Liu, C. (2019). Exploring the Phylogeography of Ancient Platycladus Orientalis in China by Specific-Locus Amplified Fragment Sequencing. Ijms 20, 3871. doi:10.3390/ijms20163871
Chau, M. M., Reyes, W. R., and Ranker, T. A. (2013). Ecological Factors Influencing Growth of the Endangered Hawaiian Fern Marsilea Villosa (Marsileaceae) and Implications for Conservation Management. Am. J. Bot. 100, 1532–1543. doi:10.3732/ajb.1200625
Cheng, J., Kao, H., and Dong, S. (2020). Population Genetic Structure and Gene Flow of Rare and Endangered Tetraena Mongolica Maxim. Revealed by Reduced Representation Sequencing. BMC Plant Biol. 20, 391. doi:10.1186/s12870-020-02594-y
Cohen, J. I., Williams, J. T., Plucknett, D. L., and Shands, H. (1991). Ex Situ Conservation of Plant Genetic Resources: Global Development and Environmental Concerns. Science 253, 866–872. doi:10.1126/science.253.5022.866
Cole, C. T. (2003). Genetic Variation in Rare and Common Plants. Annu. Rev. Ecol. Evol. Syst. 34, 213–237. doi:10.1146/annurev.ecolsys.34.030102.151717
Cornuet, J. M., and Luikart, G. (1996). Description and Power Analysis of Two Tests for Detecting Recent Population Bottlenecks from Allele Frequency Data. Genetics 144, 2001–2014. doi:10.3892/ijo_0000055110.1093/genetics/144.4.2001
Cui, Y.-M., Sun, B., Wang, H.-F., Ferguson, D. K., Wang, Y.-F., Li, C.-S., et al. (2015). Exploring the Formation of a Disjunctive Pattern between Eastern Asia and North America Based on Fossil Evidence from Thuja (Cupressaceae). Plos One 10, e0138544. doi:10.1371/journal.pone.0138544
Di Rienzo, A., Peterson, A. C., Garza, J. C., Valdes, A. M., Slatkin, M., and Freimer, N. B. (1994). Mutational Processes of Simple-Sequence Repeat Loci in Human Populations. Proc. Natl. Acad. Sci. 91, 3166–3170. doi:10.1073/pnas.91.8.3166
Earl, D. A., and Vonholdt, B. M. (2012). Structure Harvester: a Website and Program for Visualizing Structure Output and Implementing the Evanno Method. Conservation Genet. Resour. 4, 359–361. doi:10.1007/s12686-011-9548-7
Ellstrand, N. C., and Elam, D. R. (1993). “Population Genetic Consequences of Small Population Size: Implications for Plant Conservation,” in Annu. Rev. Ecol. Evol. Syst. 24, 217–242. doi:10.1146/annurev.es.24.110193.001245
Ewens, W. J. (1970). A Second Installment: Evolution and the Genetics of Populations . Vol. 2, the Theory of Gene Frequencies. Sewall Wright. University of Chicago Press, Chicago, 1969. Viii + 512 pp., Illus. $15. Science 2168, 722–723. The theory of gene frequencies. doi:10.1126/science.168.3932.722
Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V. C., and Foll, M. (2013). Robust Demographic Inference from Genomic and SNP Data. Plos Genet. 9, e1003905. doi:10.1371/journal.pgen.1003905
Excoffier, L., and Foll, M. (2011). Fastsimcoal: a Continuous-Time Coalescent Simulator of Genomic Diversity under Arbitrarily Complex Evolutionary Scenarios. Bioinformatics 27, 1332–1334. doi:10.1093/bioinformatics/btr124
Excoffier, L., and Lischer, H. E. L. (2010). Arlequin Suite Ver 3.5: a New Series of Programs to Perform Population Genetics Analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi:10.1111/j.1755-0998.2010.02847.x
Excoffier, L., Smouse, P. E., and Quattro, J. M. (1992). Analysis of Molecular Variance Inferred from Metric Distances Among DNA Haplotypes: Application to Human Mitochondrial DNA Restriction Data. Genetics 131, 479–491. doi:10.1016/1050-3862(92)90005-P10.1093/genetics/131.2.479
Fox, C. W., and Reed, D. H. (2011). Inbreeding Depression Increases with Environmental Stress: an Experimental Study and Meta-Analysis. Evolution 65, 246–258. doi:10.1111/j.1558-5646.2010.01108.x
Frankham, R., and Briscoe, D. A. (2010). Introduction to Conservation Genetics. Cambridge, UK: Cambridge University Press.
Frichot, E., Schoville, S. D., Bouchard, G., and François, O. (2013). Testing for Associations between Loci and Environmental Gradients Using Latent Factor Mixed Models. Mol. Biol. Evol. 30, 1687–1699. doi:10.1093/molbev/mst063
Gao, Y., Yin, S., Wu, L., Dai, D., Wang, H., Liu, C., et al. (2017). Genetic Diversity and Structure of Wild and Cultivated Amorphophallus Paeoniifolius Populations in Southwestern China as Revealed by RAD-Seq. Sci. Rep. 7, 14183. doi:10.1038/s41598-017-14738-6
Godefroid, S., Piazza, C., Rossi, G., Buord, S., Stevens, A.-D., Aguraiuja, R., et al. (2011). How Successful Are Plant Species Reintroductions. Biol. Conservation 144, 672–682. doi:10.1016/j.biocon.2010.10.003
Hamrick, J. L., and Godt, M. J. W. (1996). Effects of Life History Traits on Genetic Diversity in Plant Species. Phil. Trans. R. Soc. Lond. B 351, 1291–1298. doi:10.1098/rstb.1996.0112
Hamrick, J. L., Godt, M. J. W., and Sherman-Broyles, S. L. (1992). “Factors Influencing Levels of Genetic Diversity in Woody Plant Species,” in Population Genetics of forest Trees (Springer), 95–124. doi:10.1007/978-94-011-2815-5_7
Hou, L., Cui, Y., Li, X., Chen, W., Zhang, Z., Pang, X., et al. (2018). Genetic Evaluation of Natural Populations of the Endangered conifer Thuja Koraiensis Using Microsatellite Markers by Restriction-Associated DNA Sequencing. Genes 9, 218. doi:10.3390/genes9040218
Jombart, T. (2008). Adegenet: a R Package for the Multivariate Analysis of Genetic Markers. Bioinformatics 24, 1403–1405. doi:10.1093/bioinformatics/btn129
Jordan, R., Breed, M. F., Prober, S. M., Miller, A. D., and Hoffmann, A. A. (2019). How Well Do Revegetation Plantings Capture Genetic Diversity. Biol. Lett. 15, 20190460. doi:10.1098/rsbl.2019.0460
Joshi, P., and Dhawan, V. (2007). Analysis of Genetic Diversity Among Swertia Chirayita Genotypes. Biol. Plant. 51, 764–768. doi:10.1007/s10535-007-0156-z
Lauterbach, D. (2013). Ex situ-kulturen gefährdeter wildpflanzen - populationsgenetische aspekte und empfehlungen für besammlung, kultivierung und wiederausbringung. ANLiegen Natur 35, 32–39.
Lee, T.-H., Guo, H., Wang, X., Kim, C., and Paterson, A. H. (2014). SNPhylo: a Pipeline to Construct a Phylogenetic Tree from Huge SNP Data. BMC Genomics 15. doi:10.1186/1471-2164-15-162
Li, S., Liu, S.-L., Pei, S.-Y., Ning, M.-M., and Tang, S.-Q. (2020). Genetic Diversity and Population Structure of Camellia Huana (Theaceae), a limestone Species with Narrow Geographic Range, Based on Chloroplast DNA Sequence and Microsatellite Markers. Plant Divers. 42, 343–350. doi:10.1016/j.pld.2020.06.003
Li, X. Q., Li, C. S., Lu, H. Y., Dodson, J. R., and Wang, Y. F. (2004). Paleovegetation and Paleoclimate in Middle-Late Pliocene, Shanxi, central China. Palaeogeogr. Palaeoclimatol. Palaeoecol. 210, 57–66. doi:10.1016/j.palaeo.2004.03.007
Liu, D., Zhang, L., Wang, J., and Ma, Y. (2020). Conservation Genomics of a Threatened Rhododendron: Contrasting Patterns of Population Structure Revealed from Neutral and Selected SNPs. Front. Genet. 11, 757. doi:10.3389/fgene.2020.00757
Liu, J., Shi, S., Chang, E., Yang, W., and Jiang, Z. (2013). Genetic Diversity of the Critically Endangered Thuja Sutchuenensis Revealed by ISSR Markers and the Implications for Conservation. Ijms 14, 14860–14871. doi:10.3390/ijms140714860
Luikart, G., Allendorf, F., Cornuet, J., and Sherwin, W. (1998). Distortion of Allele Frequency Distributions Provides a Test for Recent Population Bottlenecks. J. Hered. 89, 238–247. doi:10.1093/jhered/89.3.238
Ma, F., Qin, A., Guo, Q., Jian, Z., Pei, S., Wang, X., et al. (2017). Geographical Distribution and Habitat Characteristics of Critically Endangered Species, Thuja Sutchuenensis. Chin. J. Ecol. 36, 1777–1784. doi:10.13292/j.1000-4890.201707.034
Nybom, H. (2004). Comparison of Different Nuclear DNA Markers for Estimating Intraspecific Genetic Diversity in Plants. Mol. Ecol. 13, 1143–1155. doi:10.1111/j.1365-294X.2004.02141.x
Ouborg, N. J. (2010). Integrating Population Genetics and Conservation Biology in the Era of Genomics. Biol. Lett. 6, 3–6. doi:10.1098/rsbl.2009.0590
Peakall, R., and Smouse, P. E. (2006). Genalex 6: Genetic Analysis in Excel. Population Genetic Software for Teaching and Research. Mol. Ecol. Notes 6, 288–295. doi:10.1111/j.1471-8286.2005.01155.x
Pickrell, J. K., and Pritchard, J. K. (2012). Inference of Population Splits and Mixtures from Genome-wide Allele Frequency Data. Plos Genet. 8, e1002967. doi:10.1371/journal.pgen.1002967
Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of Population Structure Using Multilocus Genotype Data. Genetics 155, 945–959. doi:10.1093/genetics/155.2.945
Qin, A., Ding, Y., Jian, Z., Ma, F., Worth, J. R. P., Pei, S., et al. (2021). Low Genetic Diversity and Population Differentiation in Thuja Sutchuenensis Franch., an Extremely Endangered Rediscovered conifer Species in Southwestern China. Glob. Ecol. Conservation 25, e01430. doi:10.1016/j.gecco.2020.e01430
Qin, A., Liu, B., Guo, Q., Bussmann, R. W., Ma, F., Jian, Z., et al. (2017). Maxent Modeling for Predicting Impacts of Climate Change on the Potential Distribution of Thuja Sutchuenensis Franch., an Extremely Endangered conifer from Southwestern China. Glob. Ecol. Conservation 10, 139–146. doi:10.1016/j.gecco.2017.02.004
Rai, K. C., and Ginwal, H. S. (2018). Microsatellite Analysis to Study Genetic Diversity in Khasi Pine (Pinus Kesiya Royle Ex. Gordon) Using Chloroplast SSR Markers. Silvae Genetica 67, 99–105. doi:10.2478/sg-2018-0014
Rosenberg, N. A. (2004). Distruct: a Program for the Graphical Display of Population Structure. Mol. Ecol. Notes 4, 137–138. doi:10.1046/j.1471-8286.2003.00566.x
Saccheri, I., Kuussaari, M., Kankare, M., Vikman, P., Fortelius, W., and Hanski, I. (1998). Inbreeding and Extinction in a Butterfly Metapopulation. Nature 392, 491–494. doi:10.1038/33136
Sakaguchi, S., Tsumura, Y., Crisp, M. D., Bowman, D. M. J. S., and Isagi, Y. (2014). Genetic Evidence for Paternal Inheritance of the Chloroplast in Four Australian Callitris Species (Cupressaceae). J. For. Res. 19, 244–248. doi:10.1007/s10310-012-0384-8
Sharma, H., Hyvönen, J., and Poczai, P. (2020). Development of Chloroplast Microsatellite Markers for Giant Ragweed ( Ambrosia Trifida ). Appl. Plant Sci. 8, e11313. doi:10.1002/aps3.11313
Smouse, P. E., Long, J. C., and Sokal, R. R. (1986). Multiple Regression and Correlation Extensions of the Mantel Test of Matrix Correspondence. Syst. Zoolog. 35, 627. doi:10.2307/2413122
Sun, Y., Li, L., Li, L., Zou, J., and Liu, J. (2015). Distributional Dynamics and Interspecific Gene Flow in Picea Likiangensis and P. Wilsonii Triggered by Climate Change on the Qinghai‐Tibet Plateau. J. Biogeogr. 42, 475–484. doi:10.1111/jbi.12434
Tamura, K., Stecher, G., Peterson, D., Filipski, A., and Kumar, S. (2013). MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Mol. Biol. Evol. 30, 2725–2729. doi:10.1093/molbev/mst197
Tang, C. Q., Yang, Y., Ohsawa, M., Momohara, A., Yi, S.-R., Robertson, K., et al. (2015). Community Structure and Survival of Tertiary Relict Thuja Sutchuenensis (Cupressaceae) in the Subtropical Daba Mountains, Southwestern China. Plos One 10, e0125307. doi:10.1371/journal.pone.0125307
Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M., and Shipley, P. (2004). Micro-Checker: Software for Identifying and Correcting Genotyping Errors in Microsatellite Data. Mol. Ecol. Notes 4, 535–538. doi:10.1111/j.1471-8286.2004.00684.x
Wang, H.-W., and Ge, S. (2006). Phylogeography of the endangeredCathaya argyrophylla(Pinaceae) Inferred from Sequence Variation of Mitochondrial and Nuclear DNA. Mol. Ecol. 15, 4109–4122. doi:10.1111/j.1365-294X.2006.03086.x
Wang, S., Zhu, W., Gao, X., Li, X., Yan, S., Liu, X., et al. (2014). Population Size and Time since Island Isolation Determine Genetic Diversity Loss in Insular Frog Populations. Mol. Ecol. 23, 637–648. doi:10.1111/mec.12634
Whitlock, M. C., and Lotterhos, K. E. (2015). Reliable Detection of Loci Responsible for Local Adaptation: Inference of a Null Model through Trimming the Distribution of FST. The Am. Naturalist 186 (Suppl. 1), S24–S36. doi:10.1086/682949
Xiang, Q., Fajon, A., Zhenyu, L., Likuo, F., and Zhengyu, L. (2002). Thuja Sutchuenensis: a Rediscovered Species of the Cupressaceae. Bot. J. Linn. Soc. 139, 305–310. doi:10.1046/j.1095-8339.2002.00055.x
Xu, H., Tremblay, F., and Bergeron, Y. (2013). Development and Multiplexed Amplification of SSR Markers for Thuja Occidentalis (Cupressaceae) Using Shotgun Pyrosequencing. Appl. Plant Sci. 1, 1200427. doi:10.3732/apps.1200427
Xu, T., Abbott, R. J., Milne, R. I., Mao, K., Du, F. K., Wu, G., et al. (2010). Phylogeography and Allopatric Divergence of cypress Species (Cupressus L.) in the Qinghai-Tibetan Plateau and Adjacent Regions. BMC Evol. Biol. 10, 194. doi:10.1186/1471-2148-10-194
Xu, X. X., Cheng, F. Y., Peng, L. P., Sun, Y. Q., Hu, X. G., Li, S. Y., et al. (2019). Late Pleistocene Speciation of Three Closely Related Tree Peonies Endemic to the Qinling-Daba Mountains, a Major Glacial Refugium in Central China. Ecol. Evol. 9, 7528–7548. doi:10.1002/ece3.5284
Yang, J., Cai, L., Liu, D., Chen, G., Gratzfeld, J., and Sun, W. (2020). China's Conservation Program on Plant Species with Extremely Small Populations (PSESP): Progress and Perspectives. Biol. Conservation 244, 108535. doi:10.1016/j.biocon.2020.108535
Yang, R., Feng, X., and Gong, X. (2017). Genetic Structure and Demographic History of Cycas Chenii (Cycadaceae), an Endangered Species with Extremely Small Populations. Plant Divers. 39, 44–51. doi:10.1016/j.pld.2016.11.003
Ying, J. (1994). An Analysis of the flora of Qinling Mountain Range: Its Nature, Characteristics, and Origins. J. Syst. Evol. 32, 389–410. (in Chinese with English abstract).
Yu, W., Wu, B., Wang, X., Yao, Z., Li, Y., and Liu, Y. (2020). Scale-dependent Effects of Habitat Fragmentation on the Genetic Diversity of Actinidia Chinensis Populations in China. Hortic. Res. 7, 172. doi:10.1038/s41438-020-00401-1
Yuan, N., Comes, H. P., Mao, Y.-R., Qi, X.-S., and Qiu, Y.-X. (2012). Genetic Effects of Recent Habitat Fragmentation in the Thousand-Island Lake Region of Southeast China on the Distylous Herb Hedyotis Chrysotricha (Rubiaceae). Am. J. Bot. 99, 1715–1725. doi:10.3732/ajb.1200054
Zhou, T.-H., Li, S., Qian, Z.-Q., Su, H.-L., Huang, Z.-H., Guo, Z.-G., et al. (2010). Strong Phylogeographic Pattern of cpDNA Variation Reveals Multiple Glacial Refugia for Saruma Henryi Oliv. (Aristolochiaceae), an Endangered Herb Endemic to China. Mol. Phylogenet. Evol. 57, 176–188. doi:10.1016/j.ympev.2010.07.001
Keywords: chloroplast SSR, genetic diversity, limited distribution, narrow-ranged species, nuclear RAD-seq, Thuja sutchuenensis Franch
Citation: Yao Z, Wang X, Wang K, Yu W, Deng P, Dong J, Li Y, Cui K and Liu Y (2021) Chloroplast and Nuclear Genetic Diversity Explain the Limited Distribution of Endangered and Endemic Thuja sutchuenensis in China. Front. Genet. 12:801229. doi: 10.3389/fgene.2021.801229
Received: 25 October 2021; Accepted: 30 November 2021;
Published: 23 December 2021.
Edited by:
Renchao Zhou, Sun Yat-sen University, ChinaReviewed by:
Yue Hong Yan, Orchid Conservation & Research Center, ChinaHengchang Wang, Wuhan Botanical Garden (CAS), China
Copyright © 2021 Yao, Wang, Wang, Yu, Deng, Dong, Li, Cui and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yongbo Liu, bGl1eWJAY3JhZXMub3JnLmNu