Pronounced genetic differentiation in Fokienia hodginsii revealed by simple sequence repeat markers

Abstract Fokienia hodginsii is a Tertiary relict conifer of the monotypic genus Fokienia (Cupressaceae s.l.). Currently, the species is distributed in southern China, northern Vietnam, and northern Laos and listed as a “near threatened” species by the IUCN. In this study, a total of 427 individuals of F. hodginsii were sampled from China and Vietnam to characterize its genetic diversity and population differentiation. Based on the profiles of 12 simple sequence repeat (SSR) markers, we observed a high level of genetic diversity in F. hodginsii at the species level (H e =0.635), albeit slightly lower than that of its sister species Chamaecyparis obtusa. Signals of bottleneck events were detected in the populations GXDMS, GXHJ, V‐PXB, and V‐HB, probably due to Pleistocene glaciations or overexploitation in recent years. Pronounced genetic differentiation (F st = 0.157) was found in this species. The inbreeding index (F is = 0.176 ± 0.024) indicated that F. hodginsii has a mixed mating system. Significant correlation was found between the pairwise genetic differentiation and geographic distance (r = 0.882, p = 0.01), suggesting that genetic differentiation among the populations follows the model of isolation by distance (IBD). STRUCTURE analysis and principal coordinate analysis revealed that these populations were divided into four groups: the western China group located mainly in the Yunnan–Guizhou Plateau, the central China group located mostly in the Luoxiao Mountains and Nanling Mountains, the eastern China group located in the Wuyi Mountains and the Vietnam group containing two populations in Vietnam. The different terrains and elevations of populations may be the most likely factors leading to the differentiation between the western China group and the central China group, while the geographic isolation caused by the lack of appropriate habitats may greatly contribute to the differentiation between the central China group and the eastern China group. Based on the results, some conservation suggestions for this species are provided, such as establishing seed orchards and multiple nature reserves.


| INTRODUC TI ON
Under current rapid global climate change, many endemic species are facing a high risk of extinction due to limited natural ranges resulting from genetic stochasticity or demographic, environmental, or other factors (Caughley, 1994;Gitzendanner & Soltis, 2000;Lande, 1993). It is vital to understand the genetic characteristics of these species, such as genetic diversity and population structure, for their management and the development of effective conservation strategies (Eckert, Samis, & Lougheed, 2008;Lesica & Allendorf, 2010).
The gymnosperm family Cupressaceae Bartling comprises approximately 22 genera and 150 species. Most of these species are Tertiary relict species that arose in the Jurassic (possibly as early as the Triassic), thrived in the Jurassic, and decreased in members continuously up to the present. It is also the only family of gymnosperms that is present on all continents except Antarctica (Yang, Ran, & Wang, 2012). However, except for Juniperus, Sabina, and Cupressus, most species in this family are locally endemic, and ensuring their survival under future climate change will require public and scientific attention.
The genus Fokienia Henry et Thomas (Cupressaceae s.l.) contains only one extant species, Fokienia hodginsii (Dunn) Henry et Thomas (Farjon, 2005; Figure 1). Fossil records show that Fokienia was widely distributed in the Northern Hemisphere in ancient periods: fossils in forms with foliage and attached seed cones of Fokienia were reported from the Paleocene in Saskatchewan, central Canada (McIver & Basinger, 1990); the Oligocene in Jilin, northeastern China (Guo & Zhang, 2002); and the Miocene in Zhejiang, eastern China (He, Sun, & Liu, 2012). However, this genus is currently distributed in only southern China, northern Vietnam, and northern Laos (Zheng & Fu, 1978).
In China, it occurs at elevations between approximately 1,000 and 1,800 m as a minor constituent of the subtropical evergreen (mixed) forest (Zheng & Fu, 1978). This conifer is a good landscape tree species with a beautiful shape and straight trunk  and is commonly cut down for building materials because of its light texture and material stability (Huang, Huang, Guo, & Zheng, 2015).
Currently, this conifer is listed as "near threatened (NT)" as part of the International Union for Conservation of Nature Red List (IUCN 2004) (Vuong, 2009).
Most recent studies on F. hodginsii mainly focused on seed breeding, nursery technology, plantation cultivation, essential oil extraction and development and utilization of other resources Zhao, 2005). Only one paper mentioned the progress in genetics of F. hodginsii, according to Tam, Trang, and Hoa (2011), who investigated the genetic diversity and population structure of F. hodginsii in Vietnam by applying ISSR markers and showed that F. hodginsii maintained a low level of genetic variability and a high level of genetic differentiation. They supposed that human disturbance may play a key role in the present status of F. hodginsii by leading to the degradation and fragmentation of its habitats.
Simple sequence repeat (SSR; microsatellite) markers, codominant markers with good reproducibility and high variability, are one of the best tools to understand species genetic diversity and population structure (Wang, Huang, & Long, 2013). Based on transcriptome sequencing, we synthesized 108 SSR primers that were successfully amplified in F. hodginsii (Ding et al., 2017). Applying these SSR markers, we aimed to investigate the levels of genetic diversity and population structure of this species, which could provide some reliable information for the protection of this endangered species.

| Sample collection and DNA extraction
A total of 427 individuals of F. hodginsii were sampled from 24 locations across twelve provinces of China and Vietnam (Table 1; Figure 2). A Garmin GPS unit (GPSMAP 62sc, Taiwan) was used to record the sample geographic locations with a margin of 10 m. For each population, fresh leaves were collected from 5 to 23 randomly selected fully grown individuals, which were at least 30 m apart from each other. Then, the leaf tissues were dried by silica gel and stored in zip-lock plastic bags for DNA extraction. Voucher specimens for each population were all deposited in the Herbarium of Sun Yat-sen University (SYS).
Total DNA was extracted from dried leaf tissue using the modified CTAB method (Doyle & Doyle, 1987). For each population, two individuals were randomly selected for PCR amplifications with all 108 primers designed by Ding et al. (2017). Fluorescence was added to the 3′ end of the 12 SSR markers (

| Data analyses
Linkage disequilibrium (LD) between pairs of loci and deviation from Hardy-Weinberg equilibrium (HWE) for each locus/population combination were tested using ARLEQUIN version 3.1 (Schneider, Roessli, & Excoffier, 2000). Parameters of genetic variation were calculated using GenAlEx v6.41 (Peakall & Smouse, 2006), including the total number of alleles (N a ), the effective number of alleles (N e ), the expected and observed heterozygosities (H e and H o , respectively), the Shannon information index (I) and the fixation (inbreeding) index (F is ). Additionally, FSTAT version 2.9.3.2 (Goudet, 2002) was used to calculate the allelic richness (A R ), the unbiased estimate of Wright's F-statistic (including total-population inbreeding coefficients (F it ), the overall intrapopulation inbreeding coefficient (F is ) and the interpopulation genetic differentiation coefficient (F st ), Weir & Cockerham, 1984), and pairwise F st between paired populations. Based on pairwise F st , gene flow between populations (N m ) was further estimated with the following formula: N m = (1 − F st )/4F st (Wright, 1969). Four abiotic-climate variables, namely, minimum temperature, maximum temperature, average temperature, and precipitation, from the sampled locations were obtained from the WorldClim database (Version 1.4; https://www.worldclim.org/) and used to calculate the differentiation matrix. Mantel tests (Mantel, 1967) between the matrix of the pairwise population differentiation in terms of F st /(1 − F st ) and the differentiation matrix of geographic distances or abiotic-climate variables were performed with GenAlEx with 1,000 random permutations (Rousset, 1997).
Taking into account the geographic location of each population and the genetic differentiation within and among populations, Spatial Analysis of Molecular Variance (SAMOVA) software (Dupanloup, Schneider, & Excoffier, 2002) was used to define the best number of groups; then, ARLEQUIN version 3.11 was used for the analysis of molecular variance (AMOVA; Excoffier, Smouse, & Quattro, 1992), in which three levels of genetic differentiation were calculated: genetic differentiation within populations, genetic differentiation among populations within groups, and genetic differentiation among groups.  (Piry, Luikart, & Cornuet, 1999)  In addition, a Bayesian clustering approach implemented in STRUCTURE v2.3.4 (Evanno, Regnaut, & Goudet, 2005)  TA B L E 4 Genetic diversity at the 12 microsatellite loci the Jaccard distance between populations using MVSP software (Kovach, 1999).

| Genetic diversity
According to the LD analysis for these 12 polymorphic loci, no pairs of loci showed linkage disequilibrium after a sequential Bonferroni correction for multiple tests, indicating that the 12 markers can be considered independent markers for population genetics studies. The genetic variation across the 24 natural populations is summarized in  Table S1).

| Genetic structure
The results from F-statistics showed that the overall intrapopulation inbreeding coefficient (F is ) was 0.176 ± 0.024, the total-population inbreeding coefficient (F it ) was 0.308 ± 0.019, the interpopulation genetic differentiation coefficient (F st ) was 0.157 ± 0.019, and the gene flow (N m ) was estimated to be 2.013 ± 0.698 (Table 4). All pairwise F st values were highly significant (p < 0.001), ranging from 0.009 (between FJDYS and FJFHS) to 0.234 (between V-HB and ZJJD; Table 5). Correlation analyses showed that the genetic differentiation was most correlated with geographic distance (r = 0.882,

| Genetic bottleneck assessments
The Wilcoxon test and sign test indicated that bottleneck events may have occurred in the populations GXDMS, GXHJ, V-PXB, and V-HB via the infinite allele model and the two-phased mutation model (Table 8).

| Genetic diversity
Genetic diversity is crucial for species, as it may influence the ability of species to cope with environmental change (Frankham, Ballou, & Briscoe, 2002;Frankham, 1995aFrankham, , 1995b. In this study, microsatellite markers were used to estimate population genetic diversity and to investigate the genetic structure of F. hodginsii. Slightly lower genetic diversity was found in F. hodginsii (H e = 0.635 ± 0.005) than in Chamaecyparis obtusa (H e = 0.780), the sister species of F. hodginsii (Matsumoto, Uchida, Taguchi, Tani, & Tsumura, 2010). Compared to other species (Nybom, 2004), the expected heterozygosities (H e ) of F. hodginsii are similar to those of regional species (H e = 0.65) and long-lived woody perennial species (H e = 0.68). Allelic diversity (N a ) and expected heterozygosity (H e ) are also commonly used to estimate the genetic diversity in natural populations (Freeland, Kirk, & Petersen, 2011;Hamilton, 2009 In China, the populations GXDMS and GXHJ, where only 5-7 individuals were collected, had the lowest genetic diversity (H e = 0.590 and 0.606, respectively), and signals of bottleneck events were also detected in these two populations (Table 8). These phenomena may be explained by insufficient sampling. However, as a Tertiary relict species, this conifer was strongly influenced by the Pleistocene glaciations, resulting in the populations contracting sharply. In China, it has been more than 2,600 years since this conifer was used to build boats and houses, and due to extensive deforestation, the lower distribution limit of this conifer has moved up by 500 m since the 1980s (Hou, Cheng, Lin, & Yu, 2004). During our field investigations, we also observed substantial evidence of deforestation near the F. hodginsii populations, and in many places where ample specimens were recorded, few or no individual were found, especially in the populations of GXDMS and GXHJ. Further, the geographic locations of these two populations were near Vietnam, indicating that the low genetic diversity observed in GXDMS and GXHJ may be caused by the same factors that account for the low genetic diversity observed in Vietnam.

| Genetic differentiation
Most conifers have high levels of genetic diversity within populations and low levels of differentiation among populations (Hamrick, Godt, & Sherman-Broyles, 1992). According to the AMOVA results in this study, the genetic diversity of F. hodginsii is primarily maintained within populations (84.66%, p < 0.01), while the genetic differentiation among populations of F. hodginsii (F st = 0.157 ± 0.019) is weak; however, the value of F is was 0.176 ± 0.024, indicating a mixed mating system in which inbreeding occurred frequently. The genetic differentiation among populations of F. hodginsii (F st = 0.157 ± 0.019) is also in accordance with that of other mixed-breeding species of seed plants (79.2%, Nybom & Bartish, 2000), slightly higher than that of wind-dispersed species (F st = 0.13), and much lower than that of entomophilous species (F st = 0.21) (Nybom, 2004). This pattern is also in accordance with previous observations that the dispersal of Fokienia is mainly through the wind, though sometimes also through insects (Jin et al., 2012;Lu et al., 2011;Wang & Ran, 2014). Such patterns were also observed in Cupressus funebris, for which the genetic diversity within populations is 88.15%, F st = 0.1580 and F is = 0.1579 (Lu et al., 2014). For the species C. obtusa, much higher genetic diversity was maintained within populations (91.7%), and genetic differentiation among populations was lower (F st = 0.039). The F is value estimated for C. obtusa was only 0.034, indicating a random mating system. Therefore, the different levels of genetic differentiation among the three species may be caused primarily by the differentiation of mating systems.
In this study, a significant correlation was found between genetic differentiation (F st /(1 − F st )) and geographic distance (r = 0.882, p = 0.01), suggesting that the genetic differentiation among populations follows the model of isolation by distance (IBD), that is, the differentiation among populations is strongly associated with geographic distance. Such a phenomenon was also observed in C. obtusa (r 2 = 0.3997 and p = 0.001, Matsumoto et al., 2010). It is also known that the dispersal of Fokienia is mainly through the wind (Jin et al., 2012;Lu et al., 2011;Wang & Ran, 2014); thus, its capability for long-distance dispersal could be limited as the geographic distance increases.
Although significant correlations were also found between genetic differentiation and climatic variables in the sampled locations, such as average temperature (r = 0.178, p = 0.04) and precipitation (r = 0.256, p = 0.01), their correlations were rather weak compared to those with geographic distance (r = 0.882, p = 0.01). It was observed that the flowering period of F. hodginsii is delayed with a decrease in temperature and precipitation (Hou et al., 2006); therefore, climatic factors may also actively increase the genetic differentiation among populations to a lesser extent.

| Population structure
The STRUCTURE model based on 12 loci identified three as the most likely number of genetic clusters, as the highest ΔK value was at K = 3.
The assignment results for K = 3 showed that the two populations in Vietnam were clustered with the eastern China group. In contrast, the results for K = 4 showed that the Vietnam populations were separated from the eastern China group and clustered as a fourth group.
However, the populations located in Vietnam are located far away from those in eastern China, and the climatic conditions are much different between the two regions. It is surprising that the two populations in Vietnam were clustered with the eastern China group and not the western China group, which is much closer to Vietnam in terms of geographic distance. More molecular data need to be analyzed to understand this pattern.
In this study, the assignment results for K = 4 were the same as the results from SAMOVA and PCoA. Therefore, it is reasonable to divide all populations into four groups: the eastern China group,

| Conservation implications
Genetic diversity plays an important role in determining the survival and adaptability of a species (Liao et al., 2015). The high genetic diversity maintained within F. hodginsii and the initial significant genetic differentiation among its populations found in this study are encouraging. However, we found recent bottleneck events in the populations GXDMS, GXHJ, V-PXB, and V-HB, suggesting that individual populations may suffer from a dramatic decline in population size. As a Tertiary relict species, the range of this conifer contracted sharply during the Pleistocene glaciations, and our field investigations also showed that the F. hodginsii populations have been overexploited since the 1980s, especially in the last ten years. For the conservation of this species, measures should be taken to increase the number of individuals and avoid the destruction caused by human activities. Ex situ conservation and breeding can also be considered to maintain the greatest within-species genetic variation, especially for the populations GXHJ and GXDMS, with higher inbreeding coefficients. Establishing seed orchards is also a good method, which could preserve favorable genes and prepare for breeding in the future. According to the results from STRUCTURE, the optimum number of groups is 4; thus, we also should establish seed orchards for these four groups to preserve their genotypes. In addition, establishing multiple F. hodginsii nature reserves, such as the Daiyunshan National Nature Reserve and Nanling National Nature Reserve, is needed, and the communities containing F. hodginsii should be classified as absolute protection areas to avoid human destruction.

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S' CO NTR I B UTI O N S
Liao, W.B. and Fan, Q. designed the research. Guo, W. and Huang, Y.SH. collected the samples. Yin, Q.Y., Huang, Y.L. and Zhou, R.CH. generated the data. Yin, Q.Y., Chen, S.F. and Zhou, R.CH. analyzed and interpreted the data. Yin, Q.Y. wrote the manuscript, and Chen, S.F. and Zhou, R.CH. edited the manuscript.

DATA ACCE SS I B I LIT Y
The primers used in this study are shown in Table 2, and all other data supporting the findings are available within the article and supplementary information file.