SSR marker-based genetic diversity and structure analyses of Camellia nitidissima var. phaeopubisperma from different populations

Plant materials and DNA extraction

In March 2023, we gathered newly-grown healthy leaves of C. nitidissima var. phaeopubisperma from 14 golden camellia planting bases in Guangxi and Guangdong, China. All the 14 planting bases, except two located in Hepu County in Guangxi (with a distance of 19 km), have a distance greater than 50 km from each other. All golden camellia individuals from each planting base formed a natural population, and 14 natural populations were included in this study. At least five representative golden camellia individuals were randomly chosen from each natural population at each planting base, and a total of 84 individuals were sampled (see Table S1 in the Supplemental Materials for sample codes and collection sites). The geographical and climatic parameters of the sampling sites as well as the sample size is shown in Table 1. Soil composition at the 14 sampling sites is presented in Table S2. The climatic and soil conditions vary across these planting bases. All individuals of C. nitidissima var. phaeopubisperma at the 14 plant bases originated from Fangchenggang, Guanggxi, growing from cuttings taken from wild plants in Fangchenggang for more than 15 years.

The collected leaves were quickly dried with silica gel and stored at 4 °C for DNA extraction and PCR amplification. Total DNA in leaves was extracted (see Appendix 1 in the Supplemental Materials for DNA extraction procedure) using the magnetic bead-based kit for plant tissue DNA extraction from Wuhan Tianyi Huayu Gene Technology Co., Ltd. (Wuhan, China) and stored at −20 °C for PCR test. DNA concentration and quality were determined by NanoDrop 8000 and 1% agarose gel electrophoresis, respectively.

Primer screening and PCR amplification

The amplification effect and polymorphism of 32 pairs of primers (see Table S3 in the Supplemental Materials) screened in our previous study (Gao, 2022) as well as the efficiency with which they evaluated golden camellia diversity were analyzed. The result showed that the evaluation effect on the diversity of golden camellia obtained by the use of only eight pairs of primers (see Table 2) was the same as by using the above 32 pairs of primers. Therefore, the eight pairs of screened primers were used for amplification of DNA isolated from the 84 plant samples. The SSR fluorescent primer system used for PCR amplification had a volume of 10 μL, consisting of 1.0 μL of fluorescent PCR products, 0.5 μL of GeneScan™ 500 LIZ, and 8.5 μl of Hi-Di™ Formamide. The amplification procedure is shown in Table S4 in the Supplemental Materials. Following the PCR test, the amplification product was detected by fluorescence-labeled capillary electrophoresis. See Appendix 2 for specifics about sample preparation for capillary electrophoresis.

SSR detection and data arrangement

Formamide and relative molecular weight internal standard were mixed thoroughly at a ratio of 250:1 (v/v). A total of 9 μL of the resultant mixture was loaded into the microplate, and 1 μL of PCR product diluted tenfold was added. The ABI 3730xL DNA sequencer was used for capillary electrophoresis, and the original data obtained by it were analyzed by the Fragment (Plant) fragment analysis software in GeneMarker. The position of the internal standard in each lane and the peak position of each sample were compared and analyzed to obtain the size of fragments. The internal standard score for each of test data must be at least 90. During allele calling, the height of a single peak without any accompanying peak(s) within the Panel range should be no less than 250, and the height of the main peak in case of multiple peaks should be at a minimum of 500; the peak(s) must be sharp and bilaterally symmetrical and no inverted peak(s) or nested peaks (peaks of more than two kinds of fluorescence appearing at the same position) were allowed. Raw data in .fsa format were exported from the ABI 3730xL DNA sequencer and imported into the GeneMarker analysis software after sorting and filing by loci detected for genotype data reading, and the genotype raw data in Excel format and the genotyping peak chart in PDF format were exported by the name of loci.

Analyses of genetic diversity and differentiation

Based on PCR amplification results, GenAlEx 6.501 was used to estimate various genetic parameters, including observed number of alleles per locus (Na), effective number of alleles per locus (Ne), Shannon information index (I), polymorphism information content (PIC), observed heterozygosity (Ho), expected heterozygosity (He), fixation index (F), intra-population inbreeding coefficient (Fis = (Hs − Hi)/Hs, where Hi is the number of expected heterozygotes for individuals and Hs is the mean of expected heterozygotes for each population), inter-population inbreeding coefficient (Fit = (Ht − Hi)/Ht, where Ht is the overall number of expected heterozygotes), genetic differentiation coefficient (Fst = (Ht − Hs)/Ht), and gene flow value (Nm = 0.25(1 − Fst)/Fst) (Wright, 1931).

Genetic structure analysis

By referring to the method described in Wu et al. (2024b), STRUCTURE 2.3.4 was used to analyze the genetic structure of 84 plant samples and 14 populations with K = 1~20, 10,000 burn-in cycles (the model reaching a steady state), and 100,000 Markov chain Monte Carlo (MCMC) iterations resulting in a progressive convergence toward more reliable statistical results (Porras-Hurtado et al., 2013). Each K value was run 20 times and the optimal ΔK value (which denotes the best clustering result) was obtained. The genetic structure graph was generated according to the optimal K value. Using find.cluster in R 4.4.1 (adegenet), the K-means algorithm was run to generate a graph, in which different K values were the horizontal coordinate and the obtained BIC values were the vertical coordinate; the K value at the first inflection point was selected to determine the number of clusters of the 84 plant samples. The powermarker software was used to estimate the genetic distances between the populations (samples). The Unweighted Pair Group Method using Arithmetic Averages (UPGMA) clustering based on Nei’s genetic distances was performed on all samples and clustering charts were generated. GenAlEx 6.501 was used to run principal coordinates analysis (PCoA) on all samples based on Nei’s genetic distances. PCoA can reveal the difference between two samples by intuitively comparing the linear distance between the samples on coordinate axes. The shorter the linear distance between two samples, the smaller the difference between them, or vice versa. Here we used this dimensionality reduction technique to map high-dimensional golden camellia data to the low-dimensional coordinate space to clearly show the difference in genetic relationships among golden camellia from different growing areas, so that we could identify the influence of geographic distances on golden camellia genetic relationships.

Genetic diversity

The statistics of genetic diversity and heterozygosity for the eight SSR loci are shown in Table 3. By using the eight pairs of primers, 72 alleles (Na) were observed in the 84 plant samples, with a range of 5–13 alleles and a mean of nine alleles for all loci. The number of effective alleles (Ne) ranged from 1.576 to 6.045 (3.206 on average). For all loci the Ne was lower than the Na, indicating the interaction between alleles. The I, Ho, and He ranged from 0.732 to 2.028 (1.345), from 0.274 to 0.679 (0.482), and from 0.366 to 0.835 (0.611), respectively. The Ho was lower than the He and the F was bigger than 0 in all the eight loci, indicating a departure from the Hardy-Weinberg equilibrium due to homozygous excess and heterozygote deficiency. The PIC ranged from 0.334 to 0.816 (0.579), 0.25 < PIC < 0.5 for three loci (moderately polymorphic) and PIC > 0.5 for five loci (highly polymorphic).

The Fis ranged from −0.012 to 0.233 (0.095), its value being positive for seven loci, also indicating deviation from the Hardy-Weinberg equilibrium due to homozygous excess and heterozygote deficiency. The Fit had a range of 0.144–0.315 (0.207).

The parameter values representing the genetic diversity of the 14 populations are listed in Table 4. The average value of the Na, Ne, I, Ho, and He was 3.493, 2.400, 0.910, 0.445, and 0.489, respectively. For populations LP, CF, HZ, WM, CT, SK, SL, EL, WZ, and GG, the He was higher than the Ho, indicating the presence of a certain degree of inbreeding within the populations. The population HZ had the highest value of Na, I, and He, suggesting that its genetic diversity was higher than that of other populations. The F ranged from −0.035 to 0.234 (0.071) for all the populations; populations XA, CZ, HP, and BB had a negative F value, indicating the presence of excessive random mating within each of them (they were subjected to inbreeding repression or under the pressure of natural selection).

Gene flow and genetic differentiation

As shown in Table 3, the Fst ranged from 0.090 to 0.166 (0.123), meaning that all loci had a moderate level of genetic differentiation(0.05 < Fst < 0.25). The Nm ranged from 1.254 to 2.529 (1.865); its value was bigger than 1 for the eight loci, suggesting frequent gene exchange between populations.

Our analysis of molecular variance (AMOVA) result (Table 5) shows that, genetic variations among populations, among individuals, and within individuals accounted for 3%, 19%, and 78%, respectively of the total variation. F statistics for 14 populations were given in Table 6, which shows that the Fis, Fit, Fst, and F’st were 0.197, 0.223, 0.033, and 0.073, respectively, suggesting that the populations were at a moderate level of genetic differentiation. The Nm was 7.367, greater than 1, suggesting that inter-population genetic differentiation was not a result of genetic drift.

The gene flow (Nm) values and genetic differentiation coefficients (Fst) for the 14 populations are presented in Table 7. The Nm ranged from 1.759 (HP-BB) to 10.817 (CT-SK), 3.931 on average and bigger than 1 in all cases, indicating frequent gene exchange between those populations, from which it can be deduced that genetic differentiation caused by genetic drift may be prevented among the populations. The Fst ranged from 0.023 (SK-CT) to 0.130 (XA-SK), 0.0720 on average, suggesting the presence of moderate genetic differentiation among the populations.

Genetic structure and relationships

The PCoA result of the 84 plant samples is shown in Fig. 1, from which we know that samples from the same population were distributed in a scattered manner, and genetic distances between most individuals were independent of the population they belonged to, suggesting that the genetic difference among the 84 plant samples was mainly attributed to individuals.

In order to have a deeper understanding of the genetic relationship between the 84 plant samples and between the 14 populations of C. nitidissima var. phaeopubisperma, Bayesian clustering was performed in the STRUCTURE software to analyze their genetic structure. The change of Delta K with the number of groups (K) is shown in Fig. 2, from which it can be found that Delta K reached the maximum level when K = 5, suggesting that the 84 plant samples from the 14 populations could be categorized into five ancestral groups. Genetic maps of group structure were generated with K = 5 (see Figs. 3A and 3B). In Fig. 3A, all the 84 plant samples from the 14 growing areas inherited genetic properties of the five ancestral groups, and their hereditary materials were contributed by each and all groups. However, the proportion of genes inherited by each sample from the ancestral groups differed, genes from ancestral groups 1–5 mainly inherited by 19, 12, 17, 17, and 19 samples, respectively. In Fig. 3B, the 14 natural populations inherited genes from the five ancestral groups to varying degrees, but the proportion of genes inherited by each of these 14 natural populations from the five ancestral groups also differed. Populations BB, CZ, HP and WM inherited more genes from group 1 (48%, 65%, 55%, and 63%, respectively) than from other groups; population CT inherited more genes from group 3 (47%); population XA inherited more genes from group 1 and group 3 (42% and 44%); populations CF and LP inherited more gene from group 4 (50% and 61%); and population RL inherited more genes from group 5 (52%) (refer to Table S5 in the Supplemental Materials for more information).

K-means clustering was conducted on amplified fragment length data of SSR primers for 84 plant samples. Bayesian information criterion (BIC) values were obtained via choosing different K values (the number of clusters). The value of K at the first inflection point was chosen as the optimal number of clusters, and it seemed most appropriate to divide the 84 samples into three clusters (see Fig. 4). The result of UPGMA clustering is shown in Fig. 5, which exhibits that cluster 1 was comprised of 37 samples, cluster 2 consisted of 30 samples, and cluster 3 contained 17 samples (see Table S6 in the Supplemental Materials for specific samples contained in each cluster). Samples from populations CF, HF, HP, and XA were grouped as a cluster, samples from populations BB, CT, LP, RL, SL, and WM fell into two different clusters, and samples from the remaining four populations scattered in the three clusters.

A heat map of genetic distances between the 14 populations is shown in Fig. 6 (see Table S7 in the Supplemental Materials for specific values). The maximum and minimum genetic distances were 0.3601 (WM-SL) and 0.0739 (LP-HZ), respectively. The UPGMA clustering analysis based on Nei’s genetic distances shows that the 14 populations could be divided into two clusters (see Fig. 7): one cluster consisted of populations CF, HZ, LP, RL, WZ, SK, and SL, and the other cluster comprised populations GG, BB, CT, HP, XA, CZ, and WM. Populations with a short geographic distance did not necessarily fall into the same cluster, suggesting that the grouping was not dependent of geographic distances.

Analysis of genetic diversity of SSR loci in Camellia nitidissima var. phaeopubisperma

Constituting an important part of biodiversity, genetic diversity is the prerequisite for the survival, adaptation, development, and evolution of species, and understanding the genetic diversity of a species not only provides the scientific basis for its conservation and utilization, but is of great importance to genetic improvement and innovation of its germplasm resources (Li et al., 2024). C. nitidissima var. phaeopubisperma is the most widely used variety in the section Chrysantha and has been introduced to different areas with different environmental conditions. However, no previous report has been made on the status of genetic differentiation of this variety after introduction and cultivation. In this study, we performed SSR genetic diversity analyses on 84 samples of C. nitidissima var. phaeopubisperma from 14 locations in Guangxi and Guangdong, China, and found that differences in the Na, Ne, PIC, I, Ho, and He were present among the eight screened SSR loci, suggesting that there was a large difference in genetic variations among SSR loci in C. nitidissima var. phaeopubisperma. This discrepancy is of great significance to germplasm resource evaluation of C. nitidissima var. Phaeopubisperma since it holds great potential to obtain SSR loci for efficient appraisal of genetic diversity. Five out of eight SSR loci were highly polymorphic (PIC > 0.5), and three loci were moderately polymorphic (0.2 < PIC < 0.5). According to Purvis & Franklin (2005), loci with PIC > 0.7 are best genetic markers. In our study, the PIC was higher than 0.7 for three out of eight SSR loci, meaning that it was possible to choose SSR loci for efficient investigation on genetic diversity of C. nitidissima var. phaeopubisperma from the tested loci. The He ranged from 0.366 to 0.835, indicating substantial genetic diversity among the tested samples (Wu et al., 2024a). We also found that the Ho was lower than the He for all the eight loci, indicating deviation from Hardy-Weingerg equilibrium in most loci of the progeny populations, possibly due to the influence of non-random factors on parental populations in the process of gene transfer from generation to generation (Wu et al., 2024b). Moreover, the genetic differentiation coefficients were higher than 0.05 and lower than 0.25, suggesting a moderate level of differentiation (Stebbins & Wright, 1977); the gene flow values were greater than 1, suggesting frequent gene exchange among populations and maintenance of current genetic structure. The results of both of these parameters indicated the presence of a moderate level of subdivision of populations (Wright, 1965). The populations we studied were characterized by low heterozygosity. Low heterozygosity across a species increases its vulnerability to extinction (Kramer & Havens, 2009; Oakley & Winn, 2012; Poudel et al., 2014) and goes against its sustainable development; however, a low level of heterozygosity may favor the selection and breeding of high-quality germplasms (Wu et al., 2024b).

Analysis of genetic diversity of populations of Camellia nitidissima var. phaeopubisperma

Previous studies found that no significant correlation exists between the genetic diversity and population size of C. nitidissima (Chen et al., 2022; Wei et al., 2008). There is a high degree of differentiation among populations of this species. Compared to natural populations, ex-situ populations derived from transplant of wild seedlings show more genetic variations (Wei et al., 2008; Zhu et al., 2023). However, the genetic diversity of populations originating from seeds or cuttings from a few high-quality trees is low (Zhu et al., 2023). The 14 populations in our study had a moderate level of genetic differentiation, consistent with the result of Chen et al. (2022), but lower than the result of Wei et al. (2008), likely due to the fact that the populations in our study had a higher level of gene flow than those in Wei et al. (2008). Our result validated the conclusion of Zhu et al. (2023). Both Li et al. (2020) and Tang et al. (2006) stated that genetic distances of golden camellia are significantly correlated with their geographic distances. In the present study, Nei’s genetic distances-based UPGMA clustering analysis performed on 14 natural populations of C. nitidissima var. phaeopubisperma revealed that these populations could be divided into two clusters. One cluster was composed of populations CF, HZ, LP, RL, WZ, SK, and SL, and the other cluster consisted of populations GG, BB, CT, HP, XA, CZ, and WM (see Fig. 7). It is clear that populations with a short geographic distance did not necessarily fall into the same cluster, suggesting that population clustering was independent of geographic distances. This result was in discord with those of Li et al. (2020) and Tang et al. (2006). Li et al. (2020) sampled from four sites in Fangchenggang, Guangxi, China and Tang et al. (2006) sampled from six sites in Guangxi, China (four in Fangchenggang, and two in Fusui, Nanning, Guangxi), whereas our study had a higher number of sampling sites (14 in total) and covered more environmental conditions, being more representative of the current genetic status of Camellia nitidissima var. Phaeopubisperma. Although gene exchanges existed among different populations, there are also genetic subdivision in these populations, this together with the difference in sample size may be the cause of the discrepancy regarding whether there is any correlation between genetic distances and geographic distances. The inter-population gene flow value was 7.367 in our study, bigger than 1, suggesting that frequent gene exchange was present among the populations, inconsistent with the result of C. tunghinensis chang in Tang et al. (2020), in which outcrossing mainly occurs between individuals in sub-populations, and there are low levels of gene flow and differentiation between sub-populations. As shown by our AMOVA result (see Table 5), the total genetic variation of C. nitidissima var. phaeopubisperma was mostly attributed to within-individuals (i.e., the individuals themselves). The low level of inter-population genetic variation in our study may be due to frequent gene exchange among populations. We also conducted STRUCTURE analysis to further explore the genetic relationships among the 14 populations, and the result shows that the 84 samples from the 14 populations were derived from five ancestor groups. Due to the impact from gene flow, the 14 populations inherited genes from the five groups to varying degrees.

Influences of species and material resources of golden camellia on results of genetic diversity

C. niidissima is rare, and the specimens of C. niidissima and C. chrysantha (Hu) Tuyama used in many documents may be C. nitidissima var. phaeopubisperma; for example, the Latin name of golden camellia novel food was C. chrysantha (Hu) Tuyama as provided in the Announcement of the Ministry of Health of the People’s Republic of China (The Ministry of Health of the People’s Republic of China, 2010), however, after tracing the origin of the samples submitted for inspection to the species level, it was confirmed that these samples were actually C. nitidissima var. phaeopubisperma. In the present study, the populations growing in different areas originated from seedlings growing from branch cuttings of wild plants identified to be C. nitidissima var. phaeopubisperma. By contrast, natural populations of C. nitidissima were studied in Wei et al. (2008), Li et al. (2020), and Tang et al. (2006), and C. tunghinensis chang was studied in Tang et al. (2020). The research materials are different between our study and these previous studies. This difference may lead to different results of genetic diversity, but in essence the discrepancy is due to different levels of gene flow as described above.

Our study has its limitations. For example, the use of transcriptomics alone faces challenges as biological organism might develop adaptive responses which finally can not be explained only from transcriptomics analysis. The integration of transcriptomics and metabolomics might be helpful for further explanation (Bueno & Lopes, 2020; Cadena-Zamudio et al., 2022; Fu et al., 2022).