SSR individual identification system construction and population genetics analysis for Chamaecyparis formosensis

Development of new SSR markers for C. formosensis

For a higher accuracy of court’s judgment on illegal felling, it is necessary to establish a complete forensic system. Although some SSR markers of cypress have been published^17,25,26,27, the reported detection rates for dried timber were only 20–40%^8,9. Therefore, it is necessary to develop more SSR markers as a contingency plan. When the sample is in a poor condition, more markers can be applied in order to achieve the threshold of combined power of discrimination (CPD) required for successful comparison between seized timbers and victim trees. In order to maximize potential loci, Next generation sequence (NGS) methods were used. In this study 3 DNA library were constructed. We used the Illumina MiSeq platform (2 × 301 bp; Illumina, San Diego, California, USA) to sequence the DNA libraries (Fig. 1. and Supplementary 1.).

A total of 70,325,072 raw reads were produced. The raw reads were deposited in the NCBI BioProject (PRJNA454510). After quality-trimming to the raw reads with CLC Genomics Workbench version 10 (QIAGENE, Aarhus, Denmark), 70,319,509 contigs were generated with the length between 133 and 146 bp on average. De novo assembly was conducted with the following parameters: contig number 208,467, minimum length of contigs 18 bp, maximum length contigs 108,928 bp, and average length contigs 491 bp. The sequence was assembled with software CLC Genomics Workbench version 10, and the length of the assembly sequence was 102,281,642 bp.

A sum of 78,250 SSR containing sequences was screened by MISA (v 1.0, MicroSAtellite)²⁸. We newly designed 100 candidate SSR primer pairs for testing in C. formosensis by BatchPrimer3²⁹.

There are 9 validated SSR markers that are polymorphic (success rate 9.00%) were registered in GenBank in NCBI (Table 1) and passed for cross-species tests (Supplementary 2 and 3).

Unlike the traditional SSR cloning method, with next-generation sequencing technology, it is easy to obtain a significant amount of SSR containing sequences from sequenced genomes³⁰. However, transforming candidate SSR primer pairs into validated SSR markers is still a time-consuming and expensive step. Qualified SSR markers need to succeed in PCR amplification, have good peak pattern quality with minor stuttering, and be free of non-amplifying (invalid) alleles. The turnover rate from candidate SSR primer pairs to validated SSR markers varies from species to species^30,31. The success rate in Chamaecyparis plants is between 5.24% and 9.27%^{8,17,25,26,32}.

Developing C. formosensis individual identification system

In this study, newly developed 9 validated SSR markers and other 27 validated SSR markers¹⁷ polymorphic SSR markers were analyzed against 92 individuals from 4 geographic areas (MM, HV, GW, SY, Fig. 2 and Supplementary 1). The results of developed 36 SSR markers are summarized in Table 2. Among the 92 individuals in this study, each number of alleles of SSR is between 2 and 27, with an average of 7.916. The levels of observed heterozygosity (Ho)³³ are from 0.000 to 0.891, with an average of 0.414. The levels of expected heterozygosity (He)³³ range from 0.103 to 0.906, with an average of 0.565. Significant (P < 0.001) deviations of Hardy–Weinberg equilibrium (HWE)³⁴ were detected in 23 SSR loci: Cred47, 225, 231, 236, 242, 248, 249, 250, 253, 260, 262, 276, 277, 280, 603, 610, 628, 640, 641, 674, 678, 682, 683. Ho is the actual proportion of heterozygous individuals in each locus within the population, whereas the He is the expected value estimated per HWE. Ho and He are among the most widely used parameters in estimating genetic diversity in a population. The population structural and even historical information can be obtained from Ho and He. When Ho = He, it means that the population is random mating. When Ho < He, it means that the population is inbreeding. When Ho > He, it means that the population is outcrossing³³. Most of these loci (36 tested) are Ho < He (except Cred211, Cred220, Cred225, Cred248, Cred276, Cred281, Cred297, Cred298), suggesting the population of C. formosensis has a low genetic divergence and is an inbred strain. HWE describes that under ideal conditions, there exists no mutations, no natural selection, no individuals moving in or out, the population is infinitely large, and random mating within the population. Therefore, gene frequency does not change over time or generation. However, there will always be one or more interfering factors (e.g. genetic drift, natural selection, mutation, gene flow, population bottleneck, founder effect, and inbreeding.) affecting gene frequency in nature³⁵. Therefore, HWE is difficult to achieve in nature. In this study, 23 loci out of 36 markers deviated from HWE (63.89% deviation rate). The reason for this deviation could be artificial selection, non-panmixia or genetic drift.

Polymorphism information content, or power of information content (PIC), is an index of the relative ability of the SSR marker’s genetic variability. The higher the polymorphism of marker’s genotype, the higher the PIC value³⁶. Polymorphic markers were highly informative (PIC > 0.50), reasonably informative (0.50 > PIC > 0.25), and slightly informative (PIC < 0.25). Power of discrimination (PD)³⁷ refers to the ability of genetic markers to distinguish individuals within a population. Obviously, in a population with more allele types and evenly distributed genotypes, the low probability of two random individuals having the same genotype, and the system can identify the greater probability of two random individuals. Probability of identity (P_I)³⁸ is the probability of two individuals with the same genotype. PD = 1-P_I . The value of PIC, PD and P_I of individual markers reflects its identification ability in the individual identification system. The greater PIC and PD, the lower P_I in value, suggesting the higher identification ability of the marker, and vice versa. The levels of PIC range from 0.097 to 0.876, with an average 0.528. The levels of PD range from 0.102 to 0.885, with an average 0.567. The levels of P_I range from 0.114 to 0.897, with an average 0.431. There were 19 out of 36 markers with PIC greater than 0.5, and the mean of these 36 markers PIC values was greater than 0.5, suggesting the markers have a high identification ability. The results of PD and P_I correspond to those of PIC. The highly informative markers presented in PIC also show higher identification ability in PD and P_I.

Significant linkages (P < 0.001) were detected among Cred35/229/277 (Group 1), Cred47/298 (Group 2), Cred231/249/253/262 (Group 3), Cred281/297 (Group 4), Cred603/683 (Group 5) and Cred640/678/682 (Group 6) with GENEPOP 4.2³⁹, suggesting the abovementioned group located in the same linkage group (Table 2). When identifying several independent polymorphic genetic markers simultaneously (polymorphic markers located in different linkage groups), the combined probability of identity (CP_I) is the product of the P_I of each genetic marker. At this time, CP_I will be greatly reduced, and the combined power of discrimination (CPD) will become very high. As defined above, CP_I + CPD = 1. The credibility of the individual identification system is calculated based on “Random match probability in population size and confidence levels’” published by Budowle et al.⁴⁰. Confidence levels (CL) = (1 – CP_I)^N, where N is number of individuals.

The individual identification system was applied to illegal felling cases⁸. When the seized timber and the victim tree are identified as the same particular plant, under the considerations of fairness and objective, the court usually adopts 99.99%, 99% or 95% confidence level as the credibility standard⁴¹, ISO ISO/IEC 17,025). In this study, the locus with the lowest P_I within a linkage group was used to calculate the CP_I (Table 3). The CP_I decreased along with accumulation of loci and the P_I of each locus were sorted in ascending order. The system can accumulate up to 28 loci without linkage. When reaching its maximum capability, even under most strict standard (confidence level 99.99%) dictated by court, this system can be used to identify 60 million C. formosensis, with CP_I as low as 1.652 × 10^–12, CPD as high as 0.999999999998348 (almost equal to 1), beyond the known population size of 32.06 ± 3.20 million C. formosensis⁴². Under ideal conditions, a minimum of 6 loci can be applied to the system, with an identifiable C. formosensis population of 2,900 under 95% confidence level. The CP_I is as low as 1.728 × 10^–5, and CPD is as high as 0.999982712603209 (Table 3).

One of the problems with SSR marker is the appearance of null alleles. One possible cause of SSR null alleles is poor primer annealing caused by the nucleotide sequence divergence of the flanking primer on one or both sides (for example, point mutation or indel in the primer sequence)⁴³. In addition, due to the competitive nature of PCR, smaller alleles usually have a higher amplification efficiency than larger alleles. Therefore, only the smaller of the two alleles can be detected from heterozygous individuals. The null alleles caused by differential amplification can usually be seen by loading more samples or adjusting the contrast⁴⁴. The third cause of null alleles may be due to inconsistent quality or the low quantity of DNA templates. Some loci are relatively easy to amplify, yet others cannot be amplified within the same DNA preparation⁴⁵. When a null allele is present, the observed genotype represents one of the several possible true genotypes⁴⁶. SSR markers inevitably produce null alleles, and each SSR marker has a different background for null alleles.

Dakin et al. (2004)⁴⁴ reviewed 233 publications by examining how authors detect and deal with null alleles and the methods used to estimate the frequency of null alleles across articles. The authors demonstrate that the frequency of simulated null alleles is usually overestimated, which will lead to underestimating the usability of this marker. It was misunderstood that the existence of null alleles will reduce the availability of paternity testing, individual identification, and population genetic research. However, it has been demonstrated that null alleles do not change the overall result on assignment testing^43,44. Compared with the presence of null alleles, increasing the number of loci and the degree of genetic differentiation has a more significant impact on the accuracy of assignment testing. This argument is valuable for studying SSR markers and populations prone to invalid alleles, as it allows researchers to use loci affected by invalid alleles^43,44.

In Huang et al.⁸ where C. taiwanensis individual identification system was applied to an illegal felling conviction case, CPD calculations exclude any markers that show homozygous PCR results per ISO/IEC 17,025. The CPDs are calculated only from the possibilities of the markers found in timber and tree samples simultaneously. Null alleles and PCR fail will only reduce the identification rate but will not cause seized timber and victim tree from different individuals to be identified as the same source. However, this is not to say that efforts should not be made to use loci that display low-frequency null alleles. On the contrary, markers that are less prone to invalid alleles should always be preferred because they are less ambiguous and are more potent in assignment testing. However, before many individual identification markers are developed and optimised, the impact of null alleles should not be overemphasized, as it reduces the usability of markers^43,44.

Population genetics analysis

F_is, by definition F_is = 1-Ho/He, is the inbreeding coefficient of an individual concerning the local subpopulation. When Ho < He, then F_is > 0, indicating that the population is an inbreeding. The 36 polymorphic SSRs were used to evaluate the genetic diversity parameters of the four groups (MM, HV, GW, SY) (Table 4). The number of alleles (A) for each locus is 4.417 and 5.444. Ho and He are ranged from 0.376 to 0.506 and from 0.474 to 0.583, respectively. All the groups show positive inbreeding coefficients, suggesting these four groups are inbreeding lines.

The fixation index (F_st) estimates population differentiation due to genetic structure⁴⁷. A higher F_st value means a higher degree of difference between populations. When F_st is less than 0.05, there is no differentiation among populations. When F_st is between 0.05 and 0.15, there is low differentiation among populations. On the other hand, the estimation of the number of migrants (Nm) is gene flow value⁴⁷. If Nm is more than one, genes frequently exchange, which counteracts the genetic drift and prevents the population differentiation⁴⁸. If Nm is greater than four, the population is a random mating⁴⁹. The analyses of F_st and Nm of the four geographic areas were conducted by GeneAlex 6.503⁵⁰ (Table 5). The F_st value between HV and GW was 0.035, suggesting no population differentiation in these two populations. The highest F_st value (0.074) was found between HV and MM. The F_st values ranged from 0.056 to 0.065 were found between the rest geographic areas, indicating a low differentiation in these geographic areas. The highest Nm value (6.832) was found between HV and GW, whereas the lowest value (3.141) was between HV and MM. The Nm values of four geographic areas were greater than 1 (between 3.141 and 6.832), suggesting a frequent gene exchange between the four geographic regions, which offsets genetic drift and prevents population differentiation. For GW/MM (Nm = 4.022), GW/HV (Nm = 6.382), the Nm values of the population are more significant than four, suggesting that these populations are random mating.

STRUCTURE analysis^51,52 was used to analyze the population genetic structure of C. formosensis (Fig. 3), and the Delta K value was calculated to obtain the optimal number of clusters. K and Delta K are shown in Fig. 3a. The individuals of C. formosensis were most likely to be three clusters (Fig. 3b): the SY located in Eastern Taiwan is an independent cluster, the MM located in Southwestern Taiwan is another cluster, whereas the two HV and GW geographic areas are in the same genetic cluster. The results of F_st (Table 5), Nm (Table 5) and STRUCTURE analyses (Fig. 3) show that C. formosensis of the four geographical areas belongs to the same genetic population. The F_st and STRUCTURE analyses suggest that the samples fall into three clusters. The hypothesis that Taiwan Island is one of the plant refuges during the Quaternary glaciation^53,54 may help to explain the results. The study of historical biogeography and phylogeny of cypress⁵⁴ suggested that C. formosensis in Taiwan diverged from Chamaecyparis in Japan 2.9 million years ago. The arrival of the Quaternary glaciation led to species extinction and the continued retreat of species to lower latitudes⁵⁵; thus Taiwan Island became a refuge for many ancient species, such as Juniperus morrisonicola (Cupressaceae)⁵⁶, Abies kawakamii (Pinaceae)⁵⁷, Castanopsis carlesii (Fagaceae)⁵⁸. After the glaciation, species spread from the refuge to the surrounding areas and created species diversity across the latitude gradient⁵⁹. In our results, the low polymorphism of C. formosensis probably indicates that they originally derived from the same large population during the glaciation. After the glacial retreat, these four C. formosensis clusters spread out from the refuge and formed the four populations due to geographic isolation.

Studies^55,60 also show that, based on molecular evidence, many plants (eg. Cunninghamia konishii, Cyclobalanopsis glauca, Trochodendron aralioides) in Taiwan island have high genetic diversity, higher than that of mainland China and Japanese archipelago. This remarkable high genetic diversity is associated with the Ice Age history in Taiwan^55,60. The low genetic diversity of C. formosensis differs from most Taiwanese plants but is similar to another endangered plant of the genus Cypress, C. taiwanensis, in Taiwan Island (A = 6.507, Ho = 0.392, He = 0.501)⁸. Compared to C. obtusa, an endangered cypress plant in the Japanese archipelago, C. formosensis is also inbreeding (F_is = 0.034), but the degree of genetic diversity (A = 23.9) is significantly lower. One possible explanation of the low genetic diversity is that a large population of C. formosensis was divided into several smaller populations after ancient glacial retreat in Taiwan, and then they were recently overexploited by humans (REF).

GENECLASS v. 2.0⁶¹ was applied to analyse the provenance of 92 individuals independently. The probability of samples returning to the correct provenance is 95.00% (MM), 88.00% (HV), 69.57% (GW), and 100.00% (SY), with an overall mean correct rate of 88.04% (Table 6). Three HV individuals were misassigned to GW and four GW individuals were misassigned to HV, corresponding to the observation that HV and GW are the same clusters. However, three GW individuals were misassigned to MM, possibly because the geographic location of GW is between HV and MM. Therefore, GW has characteristics of north and south at the same time. Likelihood, one MM was mis-assigned to GW, further supporting the inference that there is partial gene exchange between MM and GW. Our data show that the populations in eastern (SY) and western Taiwan (the rest populations) have distinct genotype differences. Within the western populations, the northern (HV) and the southern ones (MM) have obvious differences. Therefore, when seizing timbers in the future, the genotype can be served as a prefilter to infer the geographic area of the victim tree if the provenance is found to be MM, HV or SY. A further inspection is required if the provenance is GW because of the existence of gene exchange between nearby geographic areas.