High Genetic Diversity and Low Di ﬀ erentiation in Michelia shiluensis , an Endangered Magnolia Species in South China

: Research Highlights : This study is the ﬁrst to examine the genetic diversity of Michelia shiluensis (Magnoliaceae). High genetic diversity and low di ﬀ erentiation were detected in this species. Based on these results, we discuss feasible protection measures to provide a basis for the conservation and utilization of M. shiluensis . Background and Objectives : Michelia shiluensis is distributed in Hainan and Guangdong province, China. Due to human disturbance, the population has decreased sharply, and there is thus an urgent need to evaluate genetic variation within this species in order to identify an optimal conservation strategy. Materials and Methods : In this study, we used eight nuclear single sequence repeat (nSSR) markers and two chloroplast DNA (cpDNA) markers to assess the genetic diversity, population structure, and dynamics of 78 samples collected from six populations. Results : The results showed that the average observed heterozygosity (Ho), expected heterozygosity (He), and percentage of polymorphic loci (PPL) from nSSR markers in each population of M. shiluensis were 0.686, 0.718, and 97.92%, respectively. For cpDNA markers, the overall haplotype diversity (Hd) was 0.674, and the nucleotide diversity was 0.220. Analysis of markers showed that the genetic variation between populations was much lower based on nSSR than on cpDNA (10.18% and 77.56%, respectively, based on an analysis of molecular variance (AMOVA)). Analysis of the population structure based on the two markers shows that one of the populations (DL) is very di ﬀ erent from the other ﬁve. Conclusions : High genetic diversity and low population di ﬀ erentiation of M. shiluensis might be the result of rich ancestral genetic variation. The current decline in population may therefore be due to human disturbance rather than to inbreeding or genetic drift. Management and conservation strategies should focus on maintaining the genetic diversity in situ, and on the cultivation of seedlings ex-situ for transplanting back to their original habitat.


Introduction
Habitat destruction and environmental pollution caused by human disturbance [1] have limited the distribution and thus promoted the spatial isolation of many wildlife species, threatening their Forests 2020, 11, 469 3 of 15 between each individual). Leaves from each individual were placed into Ziploc bags sealed with silica gel, then brought back to our laboratory and stored at 4°C until DNA extraction.
Total genomic DNA was extracted from dried leaves following the supplied protocol in the Plant Genomic DNA Kit (DP305, TianGen Biotech Co., LTD, Beijing, China). The NanoDropTM 2000 Spectrophotometers (Thermo Scientific, Waltham, MA, USA) were used to test the concentration of extracted DNA solution, and the DNA work solution was prepared at an approximate concentration of 10 ng/µL.

Primers and Fragment Amplification
In total, 117 nSSR primer pairs were chosen from previous studies on Magnoliaceae; all were screened using eight samples from six populations. Screening revealed 16 pairs of primers that were able to amplify specific bands from selected samples. Among these, eight pairs showed high peak fluorescence signals and rich polymorphism in genotyping results. A total of 10 pairs of universal cpDNA primers were selected for screening; only two displayed rich polymorphisms in sequencing results. All primers (Table S2) were manufactured by Tsingke (Tsingke Biotech Co., Beijing, China).
Two PCR procedures were conducted for nSSR amplification, following the method of Schuelke (2000) [34]. The final volume used in the first PCR cycle was 10 µL; this contained 1 µL of genomic DNA (about 10 ng), 5 µL of 2X Es Taq Master Mix (Cwbiotech, Beijing, China), 1 µL each of forward and reverse primers, and ddH 2 O filled in the remainder. In the second PCR cycle, 5 µL of PCR product from the first cycle was used as a template in the final 30 µL volume, which contained 15 µL of 2X Es Taq Master Mix (Cwbiotech, Beijing, China), 3 µL fluorescent dye primer, and ddH 2 O filling in the remainder.
The DNA amplification protocol for nSSR was as follows: 5 min at 94°C for denaturation, followed by cycles of 30 s at 94°C, 30 s at specific Ta (Table S2), and 30 s at 72°C, with a final extension step of 72°C for 7 min. Twenty cycles were conducted for the first PCR procedure, and 35 for the second. The SSR-PCR products were run on the ABI 3730 DNA Analyzer (Biosystems, Foster City, CA, USA), using GS-500 as an internal size standard. Allele size was assessed using GeneMarker 2.2.0 [35]. For cpDNA amplification, cycling parameters were set as follows: initial denaturation at 94 • C for 3 min, followed by 35 cycles of denaturation at 94 • C for 1 min, annealing for 35 s at 56°C, extension at 72 • C for 1 min, and a final extension at 72 • C for 5 min. Products of cpDNA amplification were sequenced on the ABI 3730 xl DNA Analyzer (Biosystems, Foster City, CA, USA).

nSSR Analysis
Genetic diversity parameters were generated using GeneAlEx 6.5.2 [36]. These were: the number of observed alleles (Na), number of effective alleles (Ne), the observed (Ho) and expected (He) heterozygosity, the Shannon diversity index (I), the percentage of polymorphic loci (PPL), and the private alleles (Np). A Mantel test was also carried out in GenAlEx 6.5.2 [36]. PowerMarker v3.25 [37] was used to calculate polymorphism information content (PIC) and linkage disequilibrium (LD). Population genetic structure analysis using mixed models and allele-dependent frequency models was conducted using STRUCTURE 2.3.4 [38]. The parameters in STRUCTURE were set as follows: K = 1-6, Number of Interaction = 10, Burn-in period = 1,000,000, Markov chain Monte Carlo (MCMC) = 200,000. The results from STRUCTURE were compressed and uploaded to Structure Harvester online service; from this, we obtained the most likely K value (∆K) and L value (K) using Evanno's method [39]. Optional clustering was summarized using CLUMPP v1.1.2 [40]. The neighbor-joining (NJ) tree clustering algorithm was conducted using Mega 7 [41] and edited on iTOL [42]. Additionally, gene flow between populations (Nm), F-statistics for inbreeding coefficient (F IS ), global inbreeding coefficient (F IT ), and the coefficient of genetic differentiation (F ST ) were generated using Popgene 32 [43]. Arelquin v3.5 [44] was used to calculate the pairwise F ST for each population and for analysis of molecular variance (AMOVA).

Genetic Diversity and Population Structure from nSSR Markers
A total of 115 alleles were detected at the eight loci (Table S3). The Na at each locus ranged from 7 to 21, with an average of 14.375 per loci. The average I was 1.506; the lowest was 0.820 (WS18) and the highest was 1.767 (LT116). The PIC values ranged from 0.531 to 0.909, with an average of 0.809 per locus. The Ho values ranged from 0.482 (WS18) to 0.780 (MA3-7), with an average value of 0.686. Four of the eight loci displayed higher He than Ho, with the highest at LT116 (0.797) and the lowest at WS18 (0.445), and an average of 0.718. There was no evidence of significant LD. The average estimated F IS , F IT , F ST , and Nm were 0.042, 0.176, 0.139, and 1.546, respectively (Table S4).
At the population level, the average Na, Ne, I, Ho, and He were 6.021 (4.625-9.500), 4.236  Table 1). The PPL in each population ranged from 87.50% to 100%, with an average of 97.92%. The PPL of the WZ accession was 87.50, whereas those of all other populations were 100. All populations contained private alleles, with 24 found in DL, followed by seven in YC, four in YG, three each in BS and WZ, and two in CJ. The inbreeding coefficient (F IS ) for CJ was 0.207; those for all of the other populations were less than 0.1. This result indicates that intense breeding had occurred within the CJ population. An examination of the correlation between population genetic diversity parameters and population size showed that Na, Ne, I, and Np were significantly related to the population size; however, Ho and He were not related (Figure 1). At the species level, Na, Ne, I, Ho, He were 14.375, 7.019, 2.120, 0.689, and 0.825, respectively. These parameters show high genetic diversity in M. shiluensis at both the species and population levels.
Analysis of nSSR genetic structure shows a K of 167.38 for K = 2. The results show that all 78 individuals comprised two clusters (Figure 2), which separate individuals from DL from the other populations; there was evidence of minor genetic mixing between DL and other populations.  Analysis of nSSR genetic structure shows a △ K of 167.38 for K = 2. The results show that all 78 individuals comprised two clusters ( Figure 2), which separate individuals from DL from the other populations; there was evidence of minor genetic mixing between DL and other populations.
The AMOVA analysis of the M. shiluensis populations revealed a high proportion of variance within the population (89.92%, P < 0.001), and a lower proportion of variance across populations (10.18%, P < 0.001; Table 2).   Analysis of nSSR genetic structure shows a △ K of 167.38 for K = 2. The results show that all 78 individuals comprised two clusters ( Figure 2), which separate individuals from DL from the other populations; there was evidence of minor genetic mixing between DL and other populations.
The AMOVA analysis of the M. shiluensis populations revealed a high proportion of variance within the population (89.92%, P < 0.001), and a lower proportion of variance across populations (10.18%, P < 0.001; Table 2).    The AMOVA analysis of the M. shiluensis populations revealed a high proportion of variance within the population (89.92%, p < 0.001), and a lower proportion of variance across populations (10.18%, p < 0.001; Table 2). The tree diagram obtained by the NJ clustering method shows that the 78 individuals from six accessions were divided into two clusters ( Figure 3). Forty-six individuals from DL and one from BS comprise the green cluster, and three individuals from DL and the remaining individuals comprise the red cluster. This result is very similar to that obtained through STRUCTURE analysis. The tree diagram obtained by the NJ clustering method shows that the 78 individuals from six accessions were divided into two clusters ( Figure 3). Forty-six individuals from DL and one from BS comprise the green cluster, and three individuals from DL and the remaining individuals comprise the red cluster. This result is very similar to that obtained through STRUCTURE analysis.
The Mantel test ( Figure 4) confirmed that genetic distance and geographic distance were significantly correlated (P = 0.04, R = 0.66), indicating that an increase in geographic distance leads to increased genetic differentiation among populations. However, this correlation becomes insignificant upon removal of the YC population, which is located far from the main distribution area (P = 0.20, R = 0.33).   The Mantel test (Figure 4) confirmed that genetic distance and geographic distance were significantly correlated (p = 0.04, R = 0.66), indicating that an increase in geographic distance leads to increased genetic differentiation among populations. However, this correlation becomes insignificant upon removal of the YC population, which is located far from the main distribution area (p = 0.20, R = 0.33).

Genetic Diversity and Population Structure from cpDNA Marker
The aligned sequences consist of 1142 bp from two chloroplast DNA regions: trnH-psbA (400 bp) and trnK59-matK (742 bp). A total of 12 nucleic acid substitutions were observed within the region, and a total of six haplotypes were identified (Table S5). Haplotype diversity is highest in YG, followed by YC and DL; the lowest was found in BS, WZ, and CJ. Haplotype diversity ranges from 0 to 0.800, with an average of 0.674, reflecting the degree of difference among the haplotypes in each population. The YG accession shows moderate nucleotide diversity (0.00333), while DL (0.00074) and YC (0.00117) both show lower diversity. Haplotype and nucleotide diversity are both zero in three populations (BS, WZ, and CJ), which reflects the absence of variation between individuals within each of these populations (Table 3). The haplotype H3 was shared among all populations ( Figure 5), and the remaining haplotypes were private in populations. The haplotype H3 was central in the haplotype network, and can mutate into any of the other five haplotypes through base substitution. Haplotypes H1 and H2 were exclusive to DL; H4 and H5 were exclusive to YG; and H6 was exclusive to YC. Based on its widespread distribution in all populations, we surmise that H3 may represent the ancestral haplotype of M. shiluensis.
The analysis of the ML tree ( Figure 6) shows that the haplotypes H3, H4, and H5 comprised one cluster; H1 and H2 comprised another cluster; and H6 comprised a third cluster.
Total genetic diversity (H T = 0.856) was higher than that within each population (H S = 0.137). The permutation test revealed a higher value for N ST (0.728) than G ST (0.534), with P < 0.05, indicating a clear phylogenetic structure among the populations. 27 2 1 2 78 0.674 0.00220 The haplotype H3 was shared among all populations ( Figure 5), and the remaining haplotypes were private in populations. The haplotype H3 was central in the haplotype network, and can mutate in   The analysi AMOVA results ( Table 2) show that chloroplast DNA diversity of all groups was 77.56% among populations and 22.44% within populations.
The neutrality test (Table 4) shows that both Tajima's D and Fu's FS were positive, but not significant (P = 0.525 and 0.904, respectively). This result indicates that chloroplast DNA diversity was less affected by selection than nuclear DNA diversity. The analysis of the mismatch distribution using multimodal plots (Figure 7) shows that a recent expansion in population was not supported, indicating that the population may be in dynamic equilibrium. The values of SSR and HRag were 0.059 (P = 0.25) and 0.190 (P = 0.41), respectively.
The results of the correlation between geographic distance and genetic distance (Mantel test; Figure 8) show no significant correlation between the genetic distance and geographic distance of the M. shiluensis population at the cpDNA level (R=0.01, P=0.96), regardless of whether YC is removed (R=0.63, P=0.13). Table 4. Neutrality test and mismatch distribution from cpDNA fragments of Michelia shiluensis. AMOVA results ( Table 2) show that chloroplast DNA diversity of all groups was 77.56% among populations and 22.44% within populations.
The neutrality test (Table 4) shows that both Tajima's D and Fu's FS were positive, but not significant (p = 0.525 and 0.904, respectively). This result indicates that chloroplast DNA diversity was less affected by selection than nuclear DNA diversity. The analysis of the mismatch distribution using multimodal plots (Figure 7) shows that a recent expansion in population was not supported, indicating that the population may be in dynamic equilibrium. The values of SSR and HRag were 0.059 (p = 0.25) and 0.190 (p = 0.41), respectively.

High Level of Genetic Diversity
The level of genetic diversity plays an important role in the adaptive evolution and long-term survival of species and populations [13][14][15]. Generally, higher genetic diversity is found in widely distributed species than in those with restricted ranges [49,50]. However, many studies have shown that high genetic diversity can be maintained in endangered species despite dispersed populations [17,51]. In this study, a relatively high degree of genetic variation was found in the sampled populations of M. shiluensis. The nSSR data showed the population-level He to be 0.718, which is higher than that found in Magnolia tomentosa (He = 0.675) [52], Michelia coriacea (He = 0.47) [53], or Magnolia wufengensis (He = 0.184) [54], but lower than that of Magnolia stellata (He = 0.773) [9]. It is worth noting that the He found in this study is significantly higher than that found in some narrowranging species (He = 0.420) and in some widely distributed species (He = 0.620) [55]. Studies have shown that genetic diversity is positively correlated with population size [56][57][58]. Surprisingly, in this study, neither Ho (P = 0.953) nor He (P = 0.166) was significantly related to the population size ( Figure 1). This outcome has been found in some special cases [59,60]. However, because the population size in this study has remained relatively uniform, this result might require further The results of the correlation between geographic distance and genetic distance (Mantel test; Figure 8) show no significant correlation between the genetic distance and geographic distance of the M. shiluensis population at the cpDNA level (R = 0.01, p = 0.96), regardless of whether YC is removed (R = 0.63, p = 0.13).

High Level of Genetic Diversity
The level of genetic diversity plays an important role in the adaptive evolution and long-term survival of species and populations [13][14][15]. Generally, higher genetic diversity is found in widely distributed species than in those with restricted ranges [49,50]. However, many studies have shown that high genetic diversity can be maintained in endangered species despite dispersed populations [17,51]. In this study, a relatively high degree of genetic variation was found in the sampled populations of M. shiluensis. The nSSR data showed the population-level He to be 0.718, which is higher than that found in Magnolia tomentosa (He = 0.675) [52], Michelia coriacea (He = 0.47) [53], or Magnolia wufengensis (He = 0.184) [54], but lower than that of Magnolia stellata (He = 0.773) [9]. It is worth noting that the He found in this study is significantly higher than that found in some narrowranging species (He = 0.420) and in some widely distributed species (He = 0.620) [55]. Studies have

High Level of Genetic Diversity
The level of genetic diversity plays an important role in the adaptive evolution and long-term survival of species and populations [13][14][15]. Generally, higher genetic diversity is found in widely distributed species than in those with restricted ranges [49,50]. However, many studies have shown that high genetic diversity can be maintained in endangered species despite dispersed populations [17,51]. In this study, a relatively high degree of genetic variation was found in the sampled populations of M. shiluensis. The nSSR data showed the population-level He to be 0.718, which is higher than that found in Magnolia tomentosa (He = 0.675) [52], Michelia coriacea (He = 0.47) [53], or Magnolia wufengensis (He = 0.184) [54], but lower than that of Magnolia stellata (He = 0.773) [9]. It is worth noting that the He found in this study is significantly higher than that found in some narrow-ranging species (He = 0.420) and in some widely distributed species (He = 0.620) [55]. Studies have shown that genetic diversity is positively correlated with population size [56][57][58]. Surprisingly, in this study, neither Ho (p = 0.953) nor He (p = 0.166) was significantly related to the population size (Figure 1). This outcome has been found in some special cases [59,60]. However, because the population size in this study has remained relatively uniform, this result might require further consideration. The cpDNA data show that the overall nucleotide diversity was relatively low (0.00220), but the haplotype diversity (Hd) was 0.674 (Table 3); this is higher than that of both M. stellata (Hd = 0.527) [9] and Michelia maudiae (Hd = 0.44) [27], but lower than that of Michelia formosana (Hd = 0.953) [61]. High haplotype diversity further indicates that M. shiluensis is highly genetically diverse and is not greatly affected by genetic drift [9]. The lower nucleotide diversity may be due to highly conserved sequences and low substitution rates within Magnoliaceae [62].
Although the six sampled M. shiluensis groups formed a single large population comprising several smaller subpopulations [63], the genetic diversity in each subpopulation was relatively high. The reasons for high genetic diversity may be as follows: first, we speculate that gene exchange between M. shiluensis populations has historically been frequent, and thus, the existing population inherited rich ancestral genetic variation [64]; second, M. shiluensis is a perennial plant, and therefore, the recent sharp decline of individuals has not increased the chance of inbreeding, which reduces genetic diversity [65]; third, plants of the Michelia genus tend to cross-pollinate, which can reduce the loss of genetic diversity through large gene flows [55,66,67].

Lack of Genetic Differentiation between Populations
The results of molecular analysis of variance (AMOVA; Table 2) showed that genetic variation among the M. shiluensis populations was 10.18% based on nSSR markers, while that based on cpDNA markers was 77.56%. This inconsistency may be related to the differential focus of the markers or high levels of gene flow [68]. Comparative analysis between bi-parental markers and maternally inherited markers can provide comprehensive insights into population dynamics because cpDNA mutations reflect past changes, while nSSR mutations can provide inferences for recent population events [69]. Thus, differences in genetic patterns and rates of evolution often produce large contrasts between nuclear and organelle genetic diversities under conditions of genetic differentiation [28].
According to Wright [70], if Nm is less than 1, genetic drift and differentiation may occur. In this study, Nm reached 1.546 (Table S3), thus, it is very likely that the M. shiluensis population has not yet undergone either genetic drift or any apparent degree of genetic differentiation. From the perspective of M. shiluensis haplotype distribution ( Figure 5), haplotype H3 is distributed in all populations, which greatly reduces the level of differentiation between populations. In addition, the characteristic red seeds of Magnoliaceae are easily seen and eaten by birds [71]. As a result, the seeds can be widely dispersed via long-distance bird flight, thereby increasing the range of gene flow. Based on the above reasons, we can infer that a large amount of gene flow may be the reason for the current low level of genetic differentiation and indistinct geographic structure in the M. shiluensis population [72].
Generally, the level of genetic differentiation between populations increases with increasing geographic distance [73]. However, the two markers showed a different pattern when applied to geographically isolated M. shiluensis accessions (Figure 4, Figure 8). The results of nSSR showed that the genetic distance increased significantly with the increase in geographic distance (p = 0.04); however, removal of the YC population eliminated significance (p = 0.20). The YC population is far from the main distribution area and isolated by a strait (Figure 5). We speculate that the latter was the source of the significant difference we found. Results based on cpDNA showed no correlation between genetic and geographic distances regardless of whether YC is removed (P=0.96, 0.13). The result of the Mantel test may be due to the fact that, in the primary distribution range of M. shiluensis on Hainan Island [23], the maximum distance between groups of less than 80 km is insufficient to limit the range of bird activity; the resulting pattern is consistent with a lack of geographical isolation [9]. In addition to bird-mediated seed dispersal, gravity is another agent of dispersal, which results in a spatially-restricted genetic structure [9].

Demographic History and Population Structure
Population genetic structure reflects the genetic relationships within and between populations. From the haplotype network diagram, we can infer that M. shiluensis is likely to have originated in central Hainan. Colonization in DL fostered the formation of the haplotypes H1 and H2, while haplotypes H4 and H5 were formed after colonization in YG. Based on the large number of individuals and the relatively complete age structure in the DL population, we conclude that H1 and H2 form dominant haplotypes; however, as the number of individuals remaining in YG is small, it is impossible to determine which one of these haplotypes is dominant. Results of nSSR analysis confirm the genetic distinctiveness between DL and other Hainan groups, but suggest that the YG population may be very similar to BS, WZ, and CJ. The large number of private alleles in the DL population may be an important adaptation to an atypical environment [74], or a response to climate change; this result may be related to the local adaptability this species developed in response to the low temperature conditions of the Cenozoic Era [75,76]. It is worth mentioning that M. shiluensis has also been found outside Hainan Province. Studies have shown that Hainan Island was connected to Guangxi Province 65 million years ago [77]. In this study, STRUCTURE results of the analysis of nSSR markers indicate that the YC population shares the same gene pool as BS, YG, WZ, and CJ, while the analysis of cpDNA markers indicate that all sampled populations share the H3 haplotype. Furthermore, we speculate that there may be residual M. shiluensis in the area around the border between Guangdong and Guangxi.

Implication for Conservation
One of the main goals in protecting endangered species is to maximize existing genetic variation [17], as genetic variation may largely promote adaptation to a changing environment [78,79]. As all of the sampled M. shiluensis groups showed similar levels of genetic diversity, all of these populations make important contributions to the viability and protection of the species. Strengthening in situ conservation is an ideal approach that enables existing genetic resources to continue developing in an existing suitable environment. According to our field investigation, natural reserves are well-established in the DL, WZ, YG, and YC populations. Although no natural reserves are established in BS and CJ, these populations have also received protective support from the local government. The protection of existing genetic resources is only a short-term goal; the medium-term goal is to establish a new generation to improve genetic diversity. Except for DL, the individuals in the remaining populations are mature trees with a high level of genetic diversity, which provides protective value to each individual. One strategy for propagating future generations is to collect seeds in order to cultivate young seedlings and establish breeding in an off-site location [11]. In addition, germplasm resources can be propagated by grafting, and genetically dissimilar individuals can be crossed by artificial pollination to produce highly heterozygous seedlings. Great care needs to be taken to avoid genetic contamination when transplanting seedlings back to their native land [80], as this could allow foreign genes to infiltrate the original population and potentially eliminate it due to a lack of competitiveness. Finally, future research can encompass building a niche model to explore possible habitats of M. shiluensis, further expand its distribution range, and generate more genetic variation through local adaptation.

Conclusions
In summary, our results from nSSR and cpDNA markers indicate that high levels of population genetic diversity and low levels of genetic differentiation still exist in the endangered plants we studied. Therefore, we infer that the fragmentation and isolation within the M. shiluensis population may be