1. Introduction

Conservation of rare and endangered species and their sustainable utilization relies on research into their genetic diversity, population differentiation, and gene flow (Yang et al., 2010, Zaya et al., 2017). Genetic diversity is an important part of conservation biodiversity (Ramanatha Rao and Hodgkin, 2002) and essential for maintaining both the ability to cope with environmental changes and the evolutionary potential of species (Bhattacharyya and Kumaria, 2015, Ellstrand and Elam, 1993). Genetic variation may be correlated with the geographical distribution range of species (Solórzano et al., 2016). Species with a wide geographical range might maintain higher levels of genetic diversity than species with limited range (Chung et al., 2018, Levy et al., 2016). In particular, rare or endemic species with restricted distribution range and small population size are characterized by poor genetics (Solórzano et al., 2016, Zaya et al., 2017).

Camellia huana T. L. Ming et W. J. Zhang is a small shrub with yellow flowers belonging to Camellia sect. Archecamellia Sealy (Theaceae) (Fig. 1). The distribution of C. huana is very narrow, being restricted to evergreen, broad-leaved forests in the limestone hills of Tian'e County in northern Guangxi and Luodian County, Ceheng County and Xingren County in southern Guizhou, China. C. huana is the most northerly species of sect. Archecamellia and shows a discontinuous distribution with other sect. Archecamellia species. It may be better adapted to northern climates (Logan et al., 2019) and is of great significance in the introduction of yellow-flowered Camellia to the north. Although the total number of C. huana plants was once more than 310, 000, by 2013, only about 30, 000 individuals existed in natural populations (An, 2005, Xie et al., 2014). A great number of C. huana plants have been transplanted because C. huana is an ornamental plant and its flowers and leaves can be used as tea. C. huana is listed as an Endangered (EN) species on the Threatened Species List of China's Higher Plants (Qin et al., 2017).

In the present study, we combined biparentally inherited microsatellites (simple sequence repeats; SSRs) with uniparentally inherited chloroplast DNA (cpDNA) to investigate the genetic diversity and structure of C. huana. Microsatellites are effective molecular tools for assessing genetic variation, characterized by their high variability, low cost of detection, large distribution and co-dominant inheritance (Abbasov et al., 2018, Spooner et al., 2007). cpDNA displays conserved gene order, low mutation rate and a lower effective population size than nuclear genes, and is thus suitable for providing supplementary and more accurate information on population variation and differentiation (Huang et al., 2011). In C. huana, cpDNA is maternally inherited and dispersed by seeds (Birky, 2008). This study aimed to (1) estimate the level of genetic diversity within and between C. huana populations, (2) infer the genetic structure and differentiation of C. huana, and (3) develop effective conservation strategies and management measures.

2. Materials and methods 2.1. Plant materials

The entire distribution area of C. huana was surveyed in 2018 and 2019, and only five natural populations were found, in Tian'e County, Guangxi. Four populations were from the Longtan Nature Reserve (populations DLT, NRT, HLT and JSP). A total of 358 individual C. huana samples were collected from these five natural populations and seven other populations directly transplanted from wild habitats: three in the south of Guizhou province, and the remaining four in the north of Guangxi. Twenty-nine to thirty individuals were collected randomly from each population. One plant was sampled approximately every 10–50 m in the natural populations, depending on the size of each population. The location and number of sampled individuals for each population are provided in Table 1 and Fig. 2a.

2.2. Molecular procedures

Healthy, fresh leaves were harvested and stored with silica gel in zip-lock plastic bags until DNA isolation. Genomic DNA was extracted according to a modified CTAB method (Doyle and Doyle, 1990). Extracted DNA was dissolved in 50 μL 1 × TE buffer for use as PCR templates. Primer pairs described by Wei et al. (2017) and Xi et al. (2012) were adopted for amplifying a cpDNA fragment from the small single-copy (SSC) region comprising the following three fragments: SSC1, CP30 and SSC3. PCR products were sequenced using the same primers plus four internal sequencing primers: SSC1-1, CP30-1, CP30-2 and SSC3-1 (Supplemental file 1: Table S1).

PCR amplification of cpDNA was carried out in 50 μL reactions containing 5 μL 10 × PCR buffer (Mg²⁺ plus), 4 μL dNTP mix (2.5 mM), 0.5 μL each primer (50 μM), 0.5 μL rTaq DNA Polymerase (all reagents from TaKaRa, China), 2.0 μL DNA (20 ng) and 37.5 μL double-distilled water. PCR amplification was performed in a thermocycler under the following conditions: initial 5 min denaturation at 94 ℃; followed by 35 cycles of 1 min at 94 ℃, 30 s annealing at specific annealing temperatures, 1 min elongation at 72 ℃; and final extension for 10 min at 72 ℃. All PCR products were directly purified and sequenced. All chloroplast haplotype sequences generated in this study were deposited in GenBank under accession numbers MN380425–MN380432.

Microsatellite markers were selected from previously developed microsatellites in Camellia pingguoensis D. Fang (Lu et al., 2014) and Camellia flavida H. T. Chang (Liufu et al., 2014). PCR amplification was carried out in 20 μL reactions containing 2 μL 10 × PCR buffer (Mg²⁺ plus), 0.4 μL dNTP mix (10 mM), 0.2 μL each primer (50 μM), 0.2 μL rTaq DNA Polymerase (all reagents from TaKaRa, China), 1.0 μL DNA (20 ng) and 16 μL double-distilled water. PCR amplifications were performed in a thermocycler under the same conditions as cpDNA fragment amplification. Preliminary screening identified 12 microsatellite loci (Supplemental file 1: Table S2) for further use in this study. The selected primer pairs were labeled with a fluorescent dye at the 5′ end, and their PCR products were separated and visualized using an ABI3730 automated sequencer (Applied Biosystems, USA). Microsatellite marker profiles were read using the Genemarker software (Softgenetics, USA).

2.3. cpDNA sequence analysis

All sequences were aligned using MUSCLE (Edgar, 2004). Population genetic parameters, including the number of haplotypes (n), haplotype diversity (Hd), nucleotide diversity (Pi) and number of polymorphic sites (S), were calculated for aligned cpDNA sequences using DNAsp version 5.0 (Librado and Rozas, 2009).

The software PermutCpSSR (Pons and Petit, 1996) was employed for evaluating differentiation for unordered alleles (G_ST) and ordered alleles (N_ST) under 1000 stochastic permutations, to test for the presence of phylogeographic structures with significant differences between G_ST and N_ST. Analysis of molecular variance (AMOVA) (Excoffier et al., 1992) was implemented using Arlequin 3.0 (https://www.fousoft.com/arlequin.html), to study the genetic variation partitioned within and among populations. Meanwhile, pairwise fixation indices (F_ST) were also calculated using this method. Any correlation between the matrices of F_ST genetic differentiation and geographical distance were studied by performing a Mantel test (Mantel, 1967) using GenAIEx6.5 (http://biology-assets.anu.edu.au/GenAlEx/Welcome.html).

Network analysis (http://www.fluxus-engineering.com/sharenet.htm) was performed to construct a haplotype network using the median-joining method (Bandelt et al., 1999), analyzing the genealogical relationships among all the haplotypes. ArcMap GIS (http://desktop.arcgis.com) was then used to map the geographical distribution of populations and haplotypes.

2.4. SSR analysis

GenAIEx6.5 (http://biology-assets.anu.edu.au/GenAlEx/Welcome.html) was used for dataset editing and formatting. Genepop version 4.2.2 (http://genepop.curtin.edu.au/) was used for testing departure from Hardy–Weinberg equilibrium and estimating linkage disequilibrium, and the sequential Bonferroni method (Rice, 1989) was used to correct significance levels for each population and each locus. Differentiation between pairwise populations based on F_ST was calculated using Arlequin 3.0 (https://www.fousoft.com/arlequin.html), and gene flow between population pairs was estimated using Wright's method, Nm = (1-F_ST)/4F_ST (Wright, 1931). Genetic parameters, consisting of mean number of alleles (N_A), observed heterozygosity (Ho), expected heterozygosity (H_E), percentage of polymorphic loci (PPB), fixation index (F) and inbreeding coefficient (Fis) for each microsatellite locus, were calculated using GenAIEx6.5.

An admixture model based on the Bayesian clustering method was performed with SSR data in STRUCTURE 2.3 (http://web.stanford.edu/group/pritchardlab/structure.html) to determine whether species were genetically distinct. A total of 20 independent simulations were run for a K-value of 1-13 with 100, 000 burn-in steps followed by 500, 000 steps. Structure harvester (http://taylor0.biology.ucla.edu/structureHarvester/) was then utilized to estimate the most likely number (K) (Rosenberg et al., 2001) of genetic clusters based on both log-likelihood L (K) and delta K (ΔK) (Evanno et al., 2005) methods. Principal coordinate analysis (PCoA) was carried out using GenAIEx6.5. Mantel tests in GenAIEx6.5 were performed to estimate isolation by distance for SSR data, using a correlation of F_ST/(1-F_ST) (Mantel, 1967) with geographic distance for all pairs of populations.

Bottleneck 1.2.02 (Cornuet and Luikart, 1996) was used to test for recent bottleneck events in the species. The demographic history of populations was explored under the stepwise mutation model (SMM) and two-phased model (TPM), and by a heterozygosity excess test. Computation was performed using two methods (Sign and Wilcoxon tests) appropriate for a population number < 20.

3. Results 3.1. cpDNA sequences

cpDNA consensus sequences from 120 individuals (10 per population) were 5200 bp in length containing nine polymorphic sites with no indels (Table 2). We identified eight chloroplast haplotypes (Table 1). Haplotypes C1 and C2, and haplotypes C7 and C8 were unique to populations HLT and TBC, respectively; haplotype C4 was fixed for population NJT. Haplotype C5 was the most abundant haplotype with wide distribution in five populations (MGS, ECC, LJPT, DLT and NRT); haplotypes C3 and C6 were shared by four populations (ECC, XDZ, JSP and LBT) and two populations (MGS and WGT), respectively. Genetic diversity indices of total nucleotide (Pi) and haplotype (Hd) diversity in all populations were 0.00042 (individually ranging from 0 to 0.00021) and 0.759 (individually ranging from 0 to 0.533), respectively, as inferred from cpDNA (Table 1). We detected cpDNA sequence variation only within four populations (MGS, ECC, TBC and HLT).

To further study population genetic structures of C. huana, we examined genetic differentiation by detecting sequence variation in cpDNA sequences across the 12 natural and transplanted populations. Analysis of molecular variance (AMOVA) revealed that the majority of genetic variation (93.18%) existed between populations, while only 6.82% was partitioned within populations at the cpDNA level (Table 3). The F_ST value of 0.9318 also indicated relatively large differentiation between the sampled populations. In addition, no significant correlation was detected between genetic and geographic distance at the cpDNA level (R² = 0.0189, P = 0.180) (Supplemental file 1: Fig. S1a).

To further investigate the relationship between haplotypes, we constructed a network using the median-joining method for cpDNA (Fig. 2b). Haplotypes C5 and C8 took central positions while the other haplotypes were distributed on outside nodes in the network. This distribution pattern suggested that haplotypes C5 and C8 likely represent the ancient haplotypes for this distribution area. We found a significantly larger N_ST of 0.932 than G_ST of 0.856 (P < 0.05) for cpDNA data, suggesting pairs of different cpDNA haplotypes from the same populations have more similar sequences than pairs from different populations.

3.2. SSR data

Hardy–Weinberg equilibrium test results showed that four of the 144 locus–population combinations deviated from expectation after Bonferroni correction, and no significant linkage disequilibrium was found, suggesting that the 12 microsatellite loci represented independent information across all samples.

Diversity estimates were summarized for each locus (Supplemental file 1: Table S3). A total of 99 alleles at 12 SSR loci were revealed across 358 individuals (approximately 30 per population) from 12 C. huana populations. The total number of alleles (N_T) per locus ranged from three (at locus FLA14) to 12 (at loci FLA11 and TER14), with a mean of 8.25 alleles per locus. The number of effective alleles (N_E) per locus varied from 1.135 (at locus FLA14) to 3.430 (at locus TER2), with a mean of 2.272. The observed (Ho) and expected (H_E) heterozygosity ranged from 0.053 (at locus FLA14) to 0.857 (at locus FLA21) and from 0.078 (at locus FLA14) to 0.656 (at locus FLA11), respectively.

Intra-population microsatellite variation parameters are listed in Table 1. The number of mean alleles (N_A) per population varied from 2.833 (populations XDZ and LBT) to 5.417 (population MGS), with an average number of 3.931 across all populations. The average number of effective alleles (N_E) for all populations was 2.272, varying from 1.722 (population XDZ) to 2.908 (population MGS). The average values of Ho and H_E for all populations were 0.457 and 0.466, respectively; Ho values varied from 0.353 (population JSP) to 0.589 (population TBC) and H_E values from 0.339 (population JSP) to 0.605 (population MGS). The fixation index (F) varied from -0.130 (population XDZ) to 0.196 (population LJPT), with a mean value of 0.025. The above parameters show that population MGS possesses the highest genetic diversity.

By STRUCTURE analysis of SSR data, we estimated the optimum number of groups (K) according to the procedure described by Evanno et al. (2005). The highest delta K (ΔK) for 358 individuals was K = 2 (Fig. 3a), suggesting that the 12 populations can be split into two groups. Populations XDZ, JSP and LBT were assigned to Group I; the remaining populations (MGS, ECC, TBC, LJPT, WGT, NJT, DLT, NRT and HLT) composed Group II (Fig. 3b). The existence of two groups was also supported by principal coordinate analysis (PCoA) (Fig. 4). cpDNA-based haplotype geographical distribution revealed that there was only one haplotype (C3) in Group I (populations XDZ, JSP and LBT; Fig. 2). Although haplotype C3 also existed in a cultivated population (ECC) of Group II, it may have been introduced from Group I; therefore, haplotype C3 may be a unique haplotype of group I.

AMOVA analysis of SSR data showed that 21.59% genetic variation resided between populations, and 78.41% within populations. Estimates of gene flow for each pair between the 12 populations are presented in Table 4. The maximum gene flow was observed between populations DLT and ECC (9.184), whereas the minimum gene flow was observed between populations LJPT and JSP (0.347). No significant effect of isolation by distance was discovered (Supplemental file 1: Fig. S1b) as the correlation between genetic and geographic distance was insignificant (R² = 0.0234, P = 0.360), which was supported by the results of the Mantel test.

Mutation-drift equilibrium estimated under different models and methods was calculated using Bottleneck analysis (Supplemental file 1: Table S4). This revealed that C. huana did not experience a population contraction. Under the TPM, populations TBC and NRT showed a significant excess of heterozygosity using the Wilcoxon test (P < 0.05); under the SMM, four populations (TBC, NJT, DLT and JSP) had a significant excess of heterozygosity using the Sign test. This indicates five populations (TBC, NJT, DLT, NRT and JSP) deviated from mutation–drift equilibrium. However, all populations had normal L-shaped distributions under the Mode shift model, suggesting that C. huana has not undergone a recent severe bottleneck.

4. Discussion 4.1. Chloroplast and microsatellite diversity in Camellia huana

cpDNA and microsatellite markers in C. huana revealed a low level of genetic diversity (n = 8, Hd = 0.759, Pi = 0.00042 for cpDNA, N_A = 3.931, H_E = 0.466 for SSRs). Relationships between geographic range and level of genetic diversity, at both the population and species levels, tend to be strong (Solórzano et al., 2016). Theoretically, plant species with geographically restricted distribution often possess reduced genetic diversity compared with more widespread species (Coates et al., 2003, Yoichi and Tomaru, 2014, Levy et al., 2016, Chung et al., 2018). As expected, C. huana, with limited distribution, possessed relatively lower genetic diversity than many of its relatives with more widespread distribution based on similar markers, e.g., Camellia taliensis (W. W. Smith) Melchior (Liu et al., 2012, Zhao et al., 2014), C. flavida (Wei et al., 2017), C. sinensis (L.) O. Kuntze (Yao et al., 2012) and C. japonica L. (Ueno et al., 2000) (Table 5). Furthermore, the expected heterozygosity in C. huana measured by SSR markers (H_E = 0.466) was slightly higher than the average value based on SSRs in endemic species (H_E = 0.420) but considerably lower than that in regional or widespread species (H_E = 0.650 and H_E = 0.620) (Nybom, 2004). C. huana shows a discontinuous distribution with other yellow-flowered Camellia species, resulting in limited hybridization and introgression between C. huana and its relatives. This discontinuous distribution may also contribute to the low genetic diversity.

Similar to C. huana, many other yellow-flowered Camellia species, such as Camellia chrysanthoides H. T. Chang (Chen et al., 2019), Camellia micrantha S. Y. Liang et Y. C. Zhong (Chen et al., 2019), Camellia parvipetala J. Y. Liang et Su (Chen et al., 2019), Camellia nitidissima Chi (Li et al., 2019, Lu et al., 2020) and C. nitidissima Chi var. microcarpa Chang et Ye (Lu et al., 2019), exhibit a small distribution range and low level of genetic diversity (Table 5). Many natural populations have disappeared owing to excessive land use and transplantation. The number of individuals in the surviving natural populations has declined sharply (Chen et al., 2019, Li et al., 2019, Lu et al., 2019). The influence of human factors is not reflected in the results of the present study. In the future, the genetic diversity of these Camellia species may be further eroded by these adverse disturbances. It is therefore urgent to strengthen the protection of existing populations of these yellow-flowered Camellia species.

4.2. Genetic structure and differentiation of Camellia huana

Genetic differentiation between C. huana populations in this study was high (F_ST = 0.2159 for SSRs, F_ST = 0.9318 for cpDNA). The dispersal mechanism of a species can influence the level of gene flow and patterns of genetic structure among populations (Miao et al., 2017); in C. huana, this includes the migration and movement of both seeds and pollen. A previous study suggested that seeds of the congener C. japonica were mainly transmitted by gravity with limited movement (Ueno et al., 2000). Parentage analysis of seedlings of C. flavida in a 15 hm² plot revealed a seed dispersal distance of 0.65–51.05 m (mean 12.47 m), with 81.0% of dispersed seeds found within 20 m (Peng and Tang, 2017). Our fieldwork found that C. huana seeds are also large and heavy, so we speculate that they disperse by gravity without long-distance movement, contributing to the limited cytoplasmic seed-mediated gene flow. Low gene flow (Nm = 0.073 for cpDNA) between populations also confirmed this inference. Although the pollen transmission mechanism of C. huana is still unclear, the related species Camellia petelotii (Merrill) Sealy receives frequent sunbird and honeybee visitors (Sun et al., 2017). C. flavida has a pollen dispersal distance of 3.06–194.73 m (mean 29.03 m), with 72.2% of dispersed pollen found within 20 m (Peng and Tang, 2017). The natural populations of C. huana grow in evergreen, broad-leaved forests on limestone. Limestone landforms are natural "terrestrial islands" exhibiting spatial isolation on restricted land masses (Gao et al., 2015). Therefore, the limestone habitat may have led to the discontinuous distribution of C. huana and prevented the transmission of seeds and pollen, resulting in restricted gene flow among populations and increased genetic differentiation. Although we identified two genetic groups consistent with geographical distribution using SSR markers, there was only one mutational step between the chloroplast haplotype (C3) of group I and the chloroplast haplotype (C5) of group II. Whether there is a morphological differentiation between members of these groups remains to be determined.

4.3. Conservation implications

Genetic diversity is the driving force for adaptation to various changeable environments, and one of the main purposes of conservation genetics is to maintain the genetic diversity of threatened species (Lowe et al., 2005, Bhattacharyya and Kumaria, 2015). According to the results of our STRUCTURE analysis, the five natural C. huana populations should be split into two management units, and we should place more emphasis on in situ conservation. Since one management unit contains only the JSP population with about 200 individuals, we should give this population highest priority for conservation; the other management unit contains four populations (NJT, DLT, NRT and HLT). High levels of genetic differentiation among populations suggest that the transplanted populations originated from the neighboring limestone mountains, respectively. Therefore, some ex situ conservation measures are necessary for transplanted populations, with as many sample individuals as possible covering all populations.

Author contributions

Shuang Li, Shang-Li Liu, Si-Yu Pei and Man–Man Ning collected plant samples. Shuang Li, Shang-Li Liu and Si-Yu Pei performed the experiments. Shuang Li analyzed the data and wrote the manuscript. Shao-Qing Tang designed and supervised the study, and also wrote and revised the manuscript. All authors read and approved the final manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (grant number 31760088) and Biodiversity Survey, Observation and Assessment Program of the Ministry of Ecology and Environment of China (grant number 837208).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2020.06.003.