Genetic diversity and population structure of Cinnamomum balansae Lecomte inferred by microsatellites

2.1 Plant materials

With Illumina sequencing for the development of microsatellites, we collected three plant tissues, including roots, leaves, and stems of C. balansae (one individual) from a wild population (Ba Vi National Park, Hanoi city). All sample tissues were immediately frozen in liquid nitrogen and stored at –80°C until RNA extraction.

To assess the genetic variation in and population structure of C. balansae, the study area for this species was carried out in nine sites, one each in Hanoi, Phu Tho, Tuyen Quang, Yen Bai, Son La, Hoa Binh, Ninh Binh, Thanh Hoa, and Nghe An, representing its natural range in Vietnam (Table 1, Figure 1). Sampling locations were recorded using a global positioning system (GPS). The inner bark from 161 adult trees was randomly sampled for 9 known populations in central and northern Vietnam. Samples were placed into plastic bags containing silica gel in the field and then transferred to the Laboratory of Molecular Biology of Vietnam-Russia Tropical Centre and stored at −30°C until DNA extraction.

Figure 1

Map showing the collection sites of nine populations of C. balansae in Vietnam.

2.2 DNA/RNA extraction, cDNA sequencing, and microsatellite loci identification

For Illumina sequencing, total RNA from each sample was isolated using the OmniPlant RNA Kit (DNase I) (www.cwbiotech.com). Both quantity and quality of total RNA were checked through Nanodrop ND-2000 spectrophotometer (Thermo Electron Corporation, USA) and Agilent 2100 Bioanalyzer. An equal amount of total RNA from each sample was pooled together and sent to Breeding Biotechnologies Co. Ltd for cDNA library construction and transcriptome sequencing using Illumina HiSeq™4000. Cleaned mRNA was used for cDNA library construction extracted from 200 µg total RNA using Oligo (dT). The cDNA first strand was prepared from random hexamers using mRNA as a template, and the other strands were from the buffer, dNTPs, RNase H, and DNA polymerase I, and then cleaned with AMPure XP beads. The cleaned double-stranded cDNA was subjected to terminal repair, the sequencing linker was ligated, and then the fragment size was selected with AMPure XP beads. The cDNA library was acquired by PCR enrichment. After validating the library on an Agilent 2100 Bioanalyzer, Breeding Biotechnologies Co., Ltd sequenced the cDNA libraries on an Illumina HiSeq™ 4000. The Trimmomatic v3.0 [30] was used for raw reads filtration. The reads showing adaptor contamination, length < 36 bp, and low-quality value (quality <20) higher than 15% were eliminated. Trinity [31] with default parameters were used for de novo assembly of the cleaned reads. Then, the TIGR Gene Indices clustering tool (TGICL) v2.1 [32] was used to cluster and eradicate redundant transcripts and identified unigenes for further analysis.

SSRs were detected in the Simple Sequence Repeat Identification Tool program (SSRIT, http://www.gramene.org/db/markers/ssrtool). Only unigenes longer than 1 kb were included in the EST-SSR detection. The parameters were set to detect perfect, di-, tri-, tetra-, penta-, and hexa-nucleotide motifs with a minimum of 6, 5, 5, 5, and 5 repeats, respectively. Identification of microsatellite loci showed simple sequence repeats and was deposited in GenBank.

Total genomic DNA was extracted from the samples using the Plant Genome DNA Extraction Kit (BioTeke, Beijing, China). The quality of total DNA was verified by 1% agarose gel electrophoresis and NanoDrop2000 (BioTeke Instruments, Winooski, VT, USA). The concentration was then diluted to 20 ng/µL and stored at −20°C.

2.3 Microsatellite amplification

Polymerase chain reaction (PCR) was performed in a 25 μL of reaction mixture containing 2.5 µL of template DNA, 12.5 µL of 2X Taq Master mix, 1 µL of each primer, and 8 µL of deionized water. All SSR primers were designed using Primer 5.0 software [33]. The main stated selection criteria for the screening of primers were as follows GC content of 45–55 with 55% as the optimum; primer length of 18–24 bp with 20 bp as the optimum length; annealing temperature from 55 to 65°C, with 60°C as the optimum temperature; and PCR product size ranging from 100 to 300 bp. Fifteen SSR polymorphic markers were used in the present study (Table 2). The PCR amplification was carried out on a Mastercycler X50s Eppendorf, under the following thermal cycler: initial denaturation at 94°C for 2 min, followed by 40 cycles consisting of 1 min at 94°C for denaturation, 30 s alignment at the annealing temperature for each primer pair at 55°C and 1 min at 72°C for extension, and 10 min at 72°C for final cycle to complete the extension of any remaining products before holding the samples at 4°C until they were analyzed. The amplification products were separated on 8% polyacrylamide gels in 1 × TAE buffer using the Cleaver Scientific electrophoresis system (UK) and visualized with GelRed^TM 10000X. The sizes of the PCR products were detected and analyzed using GelAnalyzer software with a 20 bp DNA ladder (Invitrogen).

2.4 Data analysis

Genetic diversity was estimated using GenAlEx [34] based on SSR allele frequencies, including the number of alleles per locus (N _a), effective alleles (N _e), private alleles (N _p), the proportion of polymorphic loci (P%), the coefficient of total inbreeding (F _IT), the genetic differentiation (F _ST), Nei’s genetic distances among the populations, the observed heterozygosities (H _o), and the expected heterozygosities (H _e). Gene flow between populations (N _m) was calculated using the F _ST value: N _m = (1/F _ST − 1)/4. The inbreeding coefficients (F _IS) value was also calculated for each population by Fstat [35]. Tests for deviation from the Hardy-Weinberg equilibrium at each locus and the linkage disequilibrium for each locus pairwise combination in each population were also performed by Genpop v.4.6 [36]. Arlequin v.3.0 [37] was used to perform an analysis of molecular variance (AMOVA) for the estimation of genetic divergence as PhiPT, and its associated significance based on 10,000 permutations. Testing for recent bottleneck events for each population using Bottleneck v.1.2 [38] was applied with the stepwise mutation model (SMM) and two-phase model (TPM), and their significance was tested. The proportion of the SMM was set to 70% under default settings. The significance of these tests was evaluated by the one-tailed Wilcoxon signed-rank test.

Principal coordinate analysis (PCoA) was implemented based on G′_ST values [39] using GenAlEx. Neighbor-joining (NJ) tree was constructed genetic relationships among populations based on a matrix of F _ST values using Poptree2 [40]. Bayesian analysis of population structure was performed with Structure v.2.3.4 [41] to determine the optimal value of the genetic cluster (K). Setting the admixture model with correlated allele frequency and the ancestry models, 10 separate runs of the number of clusters (K) in the dataset were performed from 1 to 10 for each K value at 500,000 Markov Chain Monte Carlo (MCMC) repetitions and a 100,000 burn-in period. The optimal value of K was detected using Structure Harvester [42] based on the ∆K value by Evanno et al. [43].

A Mantel test [34] was applied to determine the relationships between genetic and geographic distances among the C. balansae populations based on 999 permutations.

3.1 Genetic diversity

A total of 92 different alleles were detected from 15 SSR markers in C. balansae in the tropical forests of Vietnam. Allele numbers ranged from 4 (VD51 and VD81) to 9 (VD08), with an average of 6.1 (Table 2). Among these, 35 alleles were private. Loci VD08 and VD49 had the greatest number of private alleles (4). One private allele at 2 loci (VD31 and VD51), 3 private alleles at 3 loci (VD07, VD75, and VD108), and 2 private alleles at 8 remaining loci were recorded.

Genetic diversities for C. balansae populations are presented in Table 3. The average percentage of polymorphic loci was 68.15%, ranging from 60% in the 4 populations of BV, SL, TQ, and NA to 86.67% in the TH population. The average value of alleles (N _a) per population was 1.9, ranging from 1.6 (SL) to 2.2 (YB and TH).

The highest number of private alleles (N _p = 8) was detected in the CP population, followed by three populations of YB, BE, and QH with N _p = 6, two populations of MC and DH with N _p = 4, the BV population (N _p = 2), and two populations (TC and SD) with N _p = 1. The N _e value varied from 1.4 in 5 different populations (MC, DH, SD, BE, and QH) to 1.7 (YB population), with an average of 1.5 alleles. The observed heterozygosity (H _o) and expected heterozygosity (H _e) averaged 0.246 (0.143–0.339) and 0.262 (0.198–0.338), respectively. The QH population had the lowest genetic diversity (H _o = 0.143 and H _e = 0.198) and the YB population had the highest values (H _o = 0.339 and H _e = 0.338). Five populations of BV, SD, TC, CP, and QH had positive F _IS values, indicating an excess of homozygosity and inbreeding in these populations. The highest F _IS value (0.212) was found in the CP population and was significant (p < 0.05). A significant deficit of heterozygosity was detected in two populations (TC and CP), based on the Bottleneck analysis under the TPM and SMM (p < 0.05), and suggested evidence of a bottleneck.

3.2 Population structure and Mantel test

The average F _IT value was 0.624, ranging from 0.06 (VD108) to 0.992 (VD25), indicating a deficiency of heterozygosity in the populations (Table 2). The genetic differentiation among populations (F _ST) varied from 0.285 (VD108) to 0.992 (VD25) with an average of 0.601, indicating high differentiation among all the populations and reflecting low gene flow (N _m = 0.211). Our results of the AMOVA showed that the proportion of genetic variation among the populations was higher (74%) than within populations (Table 4). The PhiPT value and gene flow among populations were 0.744 and 0.086, respectively. The PhiPT value for the variation among populations was significant (p < 0.001).

The Mantel test revealed a significantly high correlation between genetic distance and geographical distance (R ² = 0.7331) (Figure 2), indicating high genetic differentiation among populations.

Figure 2

Mantel test of genetic distance and geographical distance of nine C. balansae populations.

Genetic clusters were distinguished through the NJ analysis by using Poptree2 and the PCoA based on G′_ST values (Figure 3). Two different clusters were generated. The first cluster included three populations (CP, BE, and QH) in the central region, while the second included the remaining six populations (DH, SD, BV, MC, YB, and TC) in the northern region. An admixture model was performed to evaluate the population structure of the 161 individuals of C. balansae based on the Bayesian analysis using the Structure program. The results showed that the optimum number of genetic groups (K) was two, with the highest value of delta K (1,145.4) obtained from Structure Harvester might an optimum number of genetic clusters (Figure 4). Each individual in each population was assigned to one of the two genetic clusters. The ancestry coefficients of each individual and each population are shown in Figure 4. The red cluster represented all individuals from the populations of BV, MC, DH, TC, SD, and YB, while the green cluster represented CP and QH populations. All the individuals in the BE population were separated and closer to the red cluster than the green cluster.

Figure 3

Principal Coordinates (PCoA) (a) and NJ tree for nine populations of C. balansae with bootstrap values above branches (b).

Figure 4

The Delta K distribution graph (a) and bar plot of admixture assignment for nine C. balansae populations to cluster (K = 2 and highest ∆K value = 1145.4) based on Bayesian analysis (b).