High Genetic Diversity and Low Differentiation of Michelia coriacea (Magnoliaceae), a Critically Endangered Endemic in Southeast Yunnan, China

Michelia coriacea, a critically endangered tree, has a restricted and fragmented distribution in Southeast Yunnan Province, China. The genetic diversity, genetic structure and gene flow in the three extant populations of this species were detected by 10 inter-simple sequence repeat (ISSR) markers and 11 simple sequence repeat (SSR) markers. Examination of genetic diversity revealed that the species maintained a relatively high level of genetic diversity at the species level (percentage of polymorphic bands) PPB = 96.36% from ISSRs; PPL (percentage of polymorphic loci) = 95.56% from SSRs, despite several fragmental populations. Low levels of genetic differentiation among the populations of M. coriacea were detected by Nei’s Gst = 0.187 for ISSR and Wright’s Fst = 0.090 for SSR markers, which is further confirmed by Bayesian model-based STRUCTURE and PCoA analysis that could not reveal a clear separation between populations, although YKP was differentiated to other two populations by ISSR markers. Meanwhile, AMOVA analysis also indicated that 22.84% and 13.90% of genetic variation existed among populations for ISSRs and SSRs, respectively. The high level of genetic diversity, low genetic differentiation, and the population, structure imply that the fragmented habitat and the isolated population of M. coriacea may be due to recent over-exploitation. Conservation and management of M. coriacea should concentrate on maintaining the high level of genetic variability through both in and ex-situ conservation actions.


Introduction
The ultimate goal of conservation biology is to maintain the evolutionary potential of species by maintaining natural levels of genetic diversity [1][2][3][4]. To achieve this goal, understanding of the species' genetic diversity and genetic structure is necessary [5]. Within a species, genetic diversity of plant populations is largely determined by factors such as mating system, gene flow, genetic drift, evolutionary and life history [6][7][8]. Compared with widespread taxa, many rare and endangered species may become genetically depressed because of their small population size [9]. Habitat fragmentation exacerbates these problems, and has therefore been recognized as one of the greatest threats to the survival of many species in small and/or isolated populations [10]. Knowledge of genetic diversity and population age structure in rare plants not only enhances our understanding of population dynamics, adaptation and evolution, but also provides useful information for biological conservation [11,12].
The recently described Michelia coriacea (Magnoliaceae) is a critically endangered species, restricted to karst areas of Southeast Yunnan Province, China, bordering northern Vietnam [13][14][15]. This area has exceptionally high species diversity, and has been recognized as one of the most important biodiversity hotspots in China [16,17], but has been severely affected by habitat fragmentation due to human activities, causing populations of rare plants to become isolated from each other.
Currently, M. coriacea occurs as scattered individuals in evergreen woods on limestone slopes or roadsides at altitude of 1300-1600 m, and has a potential distribution range of approximate 4190 km 2 [18]. These individuals are distributed across only three extant populations: Daping and Shipen, Dongma and Tiechang, and Yang-Kai-Ping (defined as populations of DS, DT and YKP in this paper) (Figure 1). In DS population, some individuals of M. coriacea are over 100 years old, having been protected as symbols of -good fortune‖ by local villagers; however, they are surrounded and isolated by farmland and villages. In DT population, the maximum age of sampled trees was less than 30 years, and in YKP population the average tree age was 10-50 years [19]. All individuals in these populations are in habitats that have been degraded and fragmented due to heavy logging and vegetation destruction. Undoubtedly, an understanding of the genetic diversity of the critically endangered M. coriacea is urgently needed for planning its conservation action. The use of molecular analysis as an integral component in the conservation of rare and endangered species is becoming more widely used [5,11,20]. The molecular markers best suited for detecting genetic diversity should be relatively easy and inexpensive to detect and these markers should have evolved rapidly enough to be variable within populations. Inter-simple sequence repeats (ISSRs) and simple sequence repeats (SSRs) have been successfully employed and recognized as useful molecular markers in the analysis of a population's genetic variation, and for other purposes in various species [21,22]. In order to achieve valuable and comprehensive information on the genetic diversity and population differentiation of M. coriacea, we employed ISSR and SSR markers to address the key issue: what are the levels of genetic diversity and differentiation within, and among, populations of M. coriacea? In addition, implications for conserving this critically endangered species are discussed.

Genetic Diversity and Genetic Structure Investigated by ISSR Markers
A total of 110 bands were presented from the 10 selected primers among 58 individuals of the three populations. Of these, 106 (96.36%) were polymorphic (Table 1). At the population level, the percentage of polymorphic bands (PPB) within each population ranged from 60.00% (DS) to 81.82% (DT) with an average of 72.73% (Table 2). Populations YKP, DT and DS contained nine, eight and zero private bands respectively. Assuming Hardy-Weinberg equilibrium, the mean Nei's gene diversity (H) was estimated to be 0.226 within populations and 0.283 at the species level, and the average Shannon information index (H pop ) within populations and species levels was 0.345 and 0.436, respectively (Table 2).  In the STRUCTURE analysis of ISSR data, the value of ΔK was 903.1 for K = 2, 412.9 and 233.5 for K = 3 and K = 4, respectively, and <155 for all values of K higher than 4. Therefore K = 2 best represents the data. Structure analysis run with K = 2 showed a clear separation between YKP and the other two populations ( Figure 2).
The two-dimensional PCoA plot ( Figure 3) shows that the first principal coordinate accounts for 39.59% of total variation and separates the YKP population from the DS and DT populations. The second principal coordinate (18.88% of total variation) separated most individuals from DS from those of DT, but there was some overlap, with certain individuals from DS grouping with those from DT. AMOVA analysis among populations showed that 77.16% variation occurred within populations (Table 3). Deduced from the G st value, the level of gene flow (N m ) was estimated at 1.086 (N m > 1), indicating no differentiation among populations.

Genetic Diversity and Genetic Divergence Based on SSR Markers
A total of 45 alleles at 11 SSR loci were revealed across 58 individuals of M. coriacea. The number of alleles per locus ranged from 2 at locus MC64 to 7 at locus MC49, with an average of 4.091 per locus (Table 4). At species level, the percentage of polymorphic loci was 95%. Values for observed heterozygosity (H o ) ranged from 0.172 at locus MC45 to 0.759 locus MC8 (average 0.412), and for expected heterozygosity (H e ) from 0.247 at locus MC45 to 0.782 at locus MC8 (average 0.505). PIC values ranged from 0.417 to 0.750, with an average value of 0.502 per locus (Table 4). Table 4. SSR primers used for DNA amplifications, locus, repeat motif, primer size (bp), PCR annealing temperature (T a ), number of alleles (A), observed heterozygosity (H o ), excepted heterozygosity (H e ) and polymorphism information content (PIC). At the population level, the percentage of polymorphic loci per population ranged from 72.73% (DS) to 91.90% (DT) with an average of 82.15%. The average value for H o was 0.425 (ranging from 0.335 to 0.533) and that for H e was 0.470 (ranging from 0.429 to 0.498). Unlike private bands detected from ISSRs, private SSR alleles were present in all three populations: DS, DT and YKP contained 2, 3, and 2 private alleles respectively ( Table 5). The estimated genetic differentiation among populations (F st ) was 0.090 and the estimated level of gene flow (N m ) among populations was 2.543. AMOVA analysis showed that 13.90% of the variance existed among populations (Table 3). In the STRUCTURE analysis of SSR data, the value of ΔK was 61.2 for K = 3, 59.3 for K = 4, and <51 for all values of other K. The scores for K = 3 and K = 4 were very similar; therefore we performed runs with each K value. However, results either from K = 3 or K = 4 failed to differentiate between populations (Figure 4a,b). Similarly, PCoA analysis did not clearly separate the populations. Individuals from DS tended to score lower for coordinate 1 (accounting for 25.57% of total variation), whereas individuals from DT and YKP tended to have higher and lower scores respectively for coordinate 2 (accounting for 19.84% of total variation). However, all three populations overlapped for both coordinates, meaning that discrimination based on this data and PCoA analysis of this data was not possible ( Figure 5).
High genetic diversity can be maintained in rare plants that are adapted to an existence comprising small isolated populations (i.e., are naturally rare; [26]). However, high genetic diversity within very small populations can also be observed if very recent population size reduction has occurred, especially where that reduction has occurred within a generation or two for the species concerned; in such cases the surviving individuals are effectively samples from the larger population that existed before [11,20,27,28]. In such cases, unlike naturally rare plants, a dramatic loss of diversity in future generations is to be expected via genetic drift.
Habitat destruction appears to be recent for M. coriacea, as evidenced both by witness accounts of local people and the literature record. The species was first discovered in 1986 and then described as a new species in 1987 [29]. Initially, it was known from four localities in SE Yunnan, i.e., DS, DT, YKP and a fourth, Guangnan [29]. Despite repeated searches of the Guangnan area between 2004 and 2006, we could not find any plants of M. coriacea, indicating its extinction from Guangnan within the last 30 years. However, there is no evidence for range contraction during the same period for the three extant populations of M. coriacea in this study. In addition, our molecular data clearly indicated apparently high levels of gene flow for M. coriacea (N m = 1.086 for ISSR data, and 2.543 for SSRs), which would account for low differentiation between populations (F st = 0.090 for SSRs and G st = 0.187 for ISSRs). It should be noted that indirect estimates of Nm values must be interpreted with caution [30,31], and this data therefore should be viewed as general indicators of the magnitude of genetic exchange.
Low genetic differentiation and high gene flow between populations can result from long-distance gene dispersal either by pollen or by seed [10]. Long range seed dispersal has been implicated in maintaining links between populations in the rare Michelia formosana [32]. M. coriacea also has animal-dispersed seeds, and insect pollination may also have contributed to gene flow over distance in this species. However, seed or pollen flow appears unlikely between the modern, isolated populations. Therefore, gene flow between plants at the three extant sites via seeds and/or pollen was probably extensive and relatively unhindered in the past, before the populations became isolated.

Variation Between Populations
Populations DT (PPB = 81.82%, PPL = 91.90%) and YKP (PPB = 76.36%, PPL = 81.82%) contained relatively higher genetic diversity than did population DS (PPB = 65%, PPL = 72.73%). Because trees at DS are usually >100 years old, whereas those at DT and YKP are < 50 years old [19], range contraction within the past 50 years can be eliminated as a cause for this difference. Instead, patterns of gene flow, migration and population structure prior to range contraction can impact on modern populations [10], and may explain differences in genetic diversity between these M. coriacea populations. The fact that heterozygosity for SSR alleles was slightly higher in DS than the other two populations indicates that lower diversity here had not yet brought about inbreeding when the current generation was founded. However, the absence of seedlings from all three populations has prevented us from assessing whether range contraction will lead to an immediate inbreeding effect in the next generation.
It should be noted that there are some discrepancies between the results (e.g., levels of examined population's heterozygosity, STRUCTURE and PCoA clusterings) from different markers, suggesting that the manner of polymorphism differs because of marker specificity. Contrasting patterns are commonly found when multiple markers are used to detect genetic structure, as for example in soybean [33], cashew [34], olive [35] and Ficus carica [36]. Such contrasts between markers could be due to any or all of the following causes: (1) these two marker systems sample genomic regions with different evolutionary modes; (2) ISSR is based on consensus markers, whereas the microsatellite markers used in this study were specifically developed for M. coriacea; and (3) ISSR and SSR have different modes of inheritance (i.e., dominant versus codominant). In addition, some marker loci could be under selection, or tightly linked to loci under selection, meaning that a small minority of markers for either system might not be truly neutral, although neutrality of both ISSR and SSR markers is often assumed. To test directly for these possibilities is beyond the scope of the current study, but using two types of marker, and interpreting results of each with a degree of caution, can give a reasonable assessment of the true genetic situation for the examined populations. In the current study, population YKP was clearly separable from the other populations based on ISSR, but not SSR, markers. Hence, YKP might be the most distinct population, but SSR data indicates that this conclusion must be treated with caution.

Conservation
The maintenance of a maximal amount of genetic diversity is one of the major objectives for conserving endangered and threatened species [11], and the loss of genetic variation may largely limit adaptability to changing environments [37,38]. In the case of M. coriacea, AMOVA results shows that each of the three examined populations currently maintains a high level of within-population genetic diversity. Therefore, the endangered status of this species currently reflects a small number of extant individuals and very poor natural regeneration. Conservation measures should therefore focus on establishing large numbers of seedlings, both in and ex situ, involving many different individuals as parents, to preserve as much as possible of the existing genetic diversity in subsequent generations. Given the current lack of recruitment, establishing seedlings in situ will require management and co-operation from local communities.
In addition, the observed high diversity among relatively few adult trees also places a conservation value on each individual adult tree as a mini-reservoir of genetic diversity, in a way that would not apply for genetically depauperate, relict populations. Hence, conservation efforts to preserve each existing tree would also have a value in maintaining genetic diversity.
Conservation strategies may differ between sites. In DS, the oldest trees are protected by villagers as good luck symbols [19], so encouraging the locals to grow seedlings from these trees to bring luck to future generations might be effective. This is, however, the least diverse population. Individuals from DT and YKP occur in secondary vegetation and enjoy no protection of any kind [19]. Of these two populations, YKP might be the most distinct and hence of higher conservation value deduced from ISSR results.
Within each population, however, SSR data also indicated clusters of individuals more similar to one another than they were to others. For example, individuals DS1 and DS2 were genetically highly distinct from others is DS, and likewise YKP3, YKP4 and YKP6 were genetically distinct from others in YKP. Therefore, seedlings with a higher degree of heterozygosity could be generated by crossing individuals that, according to our results, are genetically dissimilar.
In the short term, conservation priorities must be to preserve as many as possible of the surviving individuals. In the medium term, new generations must be ensured, and in the long term, genetic health must be maintained. Therefore, in-situ efforts to conserve remaining habitats need to be combined with ex-situ research on seed propagation, with a view to establishing a new generation of plants both in cultivation and the wild. Should in-situ conservation be more successful at one site than others, then seedlings from all three sites could be established at the protected site.

The Study Species and Sampling Procedures
Michelia coriacea is diploid (2n = 38) [39]. The mature trees can reach up to 20-30 m in height and 50-60 cm in diameter [18,19,40], and live for at least 100 years. The species bears scented whitish or creamy yellow flowers with 6-7 tepals which are in anthesis from February to April [13,14].
The plant sampling and material collection were carried out in March and July 2007, and in August 2008. Prior to this, comprehensive field surveys had located only 15, 33 and 53 individuals at DS, DT and YKP, respectively [19]. It should be noted that seedling regeneration is not good, because no seedlings were detected in DS, and only a few seedlings were found in the other two populations. Therefore, we sampled all 15 mature individuals in DS, whereas 21 and 22 mature individuals were sampled from DT and YKP respectively (Figure 1 and Table 6). Healthy, young leaves were collected from each sampled individual, and dried in silica gel for subsequent DNA extraction. Table 6. Known populations of Michelia coriacea examined for ISSR and SSR analyses with population code, locality, altitude, location coordinates, numbers of individuals sampled, sample code and voucher numbers.

DNA Extraction and PCR Amplification
Total DNA was extracted from the dried leaves following the modified CTAB method described by Doyle [41]. The purified total DNA was quantified by gel electrophoresis and its quality verified by spectrophotometry. DNA samples were stored at −20 °C.
One hundred ISSR primers from the University of British Columbia (UBC) were initially screened on eight randomly selected individuals for PCR. Of these, ten primers which consistently amplified well, and produced polymorphic bands, were selected for analyzing (Table 1). To ensure repeatability of amplification and scoring, all loci were amplified 2 times independently and run on separate gels. PCR reactions were performed in a reaction volume of 15 μL, containing 9.92 μL ddH 2 O, 0.4 μL formamide, 1.5 μL 10× PCR buffer (Mg 2+ Plus), 0.25 mM of each dNTP, 0.9 U Taq polymerase, 0.6 μM primers, 0.6 μL template DNA (approximate 50 ng). Cycling conditions consisted of an initial denaturation step at 97 °C for 4 min, followed by 35 cycles of denaturation at 94 °C for 1 min, annealing temperature for 1 min (Table 1), extension at 72 °C for 1.5 min, and finally 10 min at 72 °C for final extension. Amplification products were detected by using 1.6% agarose gel stained with 0.5 pg/mL ethidium bromide, and were electrophoresed in 1× Tris-borate-EDTA (TBE) buffer (pH 8.0) at 85 V for 2 h. The bands were visualized under UV light.
SSR genotyping was performed according to the methods described in Zhao et al. [40] (Table 2); then final extension at 72 °C for 7 min. The amplified product was then separated on 6% denaturing polyacrylamide sequencing gels using silver staining. A 20 bp DNA ladder standard (Fermentas) was used as standard for scoring.

ISSR Data
Amplified ISSR fragments showing consistent amplification were scored manually as present (1) or absent (0) for each sample to form a binary matrix. The POPGENE program Version 1.31 [42] was employed, assuming Hardy-Weinberg equilibrium, to obtain the genetic diversity parameters within each population: percentage of polymorphic bands (PPB), private bands (referring to bands found only within one population), observed allele number per locus (A o ), effective allele number per locus (A e ), Nei's gene diversity (H) and Shannon index of diversity (H pop ) [31]. Genetic diversity parameters were also calculated both at species and population levels.
Genetic differentiation among populations was analyzed using Nei's gene diversity statistics [27]. The gene flow was estimated by using the equation N m = (1 − G st )/4G st , where N m is the number of migrants per generation [43]. To infer population structure and assign individuals to populations, the program STRUCTURE version 2.3.1 was used [44], following the methods described by Ma et al. [45]. We adopted the admixture model with correlated allele frequencies. No prior knowledge of the species was included in the analyzed data set. To determine the optimal number of groups (K), we ran STRUCTURE with K varying from 1 to 10, with five runs for each K value. Previous studies have found that, in many cases, the posterior probability for a given K increases slightly, even after the real K is reached [46]. Therefore, we used an ad hoc statistic, ΔK, to determine the true value of K [47]. Our parameters were 100,000 burn-in periods and 10,000 MCMC repetitions after burn-in. Furthermore, principal co-ordinate analysis (PCoA) in GenAlEx 6.0 was employed to examine further the genetic relationships among detected populations on the basis of the same ISSR data [48].
The AMOVA (analysis of molecular variance) was performed through Arlequin 3.0 program [49] to describe variance components and their significance levels for variation among individuals within and among the populations.

SSR Data
To estimate overall genetic diversity, the following measures were calculated for each SSR primer locus and each population using POPGENE program Version 1.31 [42]: observed allele number per locus (A o ), effective allele number per locus (A e ), observed heterozygosity (H o ), expected heterozygosity (H e ), and the number of private loci, i.e., those found only in a single population. The same program was used to test for deviations from Hardy-Weinberg equilibrium (HWE) and pair-wise linkage disequilibrium (LD). To measure the marker polymorphism, the polymorphism information content (PIC) for each SSR Marker was calculated according to the formula PIC = 1 − ΣP i 2 , where P i is the frequency of the allele for each SSR marker locus [50].
Wright's F st was estimated by a weighted analysis of variance with the GENEPOP 4.0 program [51]. In addition, the amount of gene flow (N m ) among populations was estimated by an indirect genetic differentiation method based on F st value, N m = (1 − F st )/4F st . Partitioning of genetic variation within and among populations was further performed by analysis of molecular variance using the Arlequin 3.0 program [49]. For these SSR data, genetic relationships among three populations were also examined using STRUCTURE and GenAlEx 6.0 with the same methods described as ISSR markers above.

Conclusions
In summary, our results from both ISSR and SSR markers show that a relatively high level of genetic diversity and low levels of genetic differentiation among populations exist in the critically endangered plant Michelia coriacea. Furthermore, Bayesian model-based STRUCTURE and PCoA analysis could not reveal a clear separation between populations, although YKP was differentiated from the other two populations by ISSR markers. We therefore presume that the fragmented habitat and the isolated populations of M. coriacea may be due to very recent contraction of its range. Based on this data, conservation strategies for this critically endangered species were also proposed.