Study site and sampling design

Teak is an insect-pollinated, major outcrossing species [20]. The fruits are naturally dispersed by wind [21] or gravity, but dispersal is very limited due to the heavy weight and lack of wings on the seeds. This species is widely distributed throughout Myanmar, which lies between the latitudes of 9° 33'N and 28° 31'N and longitudes of 92° 10'E and 101° 10'E, although the northern and southern boundaries of the teak-growing area lie at the respective latitudes of 28° 30'N and 10°N approximately. In Myanmar, teak normally grows at altitudes ranging from sea level to 1,000 m and thrives best in warm, moist, tropical climates with an annual rainfall level between 1,300 and 3,800 mm. However, teak can grow at precipitation extremes of <800 mm and >5,000 mm, corresponding to the outer fringes of the central dry zone and the southern region of Myanmar, respectively.

Myanmar can be divided into five physiographic regions based on geology, topology, and terrain: the northern mountains, comprising the upper and northern sub-regions; western ranges; eastern plateau; central basin; and lowland Myanmar, which includes the lower region and southern sub-region [22]. These geographic regions comprise seven states and seven divisions. For example, the northern mountain region includes Kachin state; the western ranges include Yakhaing and Chin states; the central region comprises Mandalay, Magway, and a small part of the Sagaing division; the lower region includes the Yangon, Ayeyarwaddy, and Bago divisions; the southern sub-region includes Mon state, the southern part of Kayin state, and Tanintharyi division; and the eastern plateau includes Shan state, Kayar state, and the northern part of Kayin state.

Fresh leaves were collected from 480 individuals representing 20 natural populations that covered nearly the entire teak distribution area throughout Myanmar. Some samples were collected from fragmented and remnant natural forests (e.g., TZ, MLM). The sample size of each population ranged from 16 to 32 (Table 1), and sampled trees were spatially separated by at least 50 m to avoid sampling closely related individuals. The collected samples were dried with silica gel in plastic bags and kept at room temperature prior to DNA extraction.

DNA extraction

DNA was isolated from dried leaf tissues using the hexadecyltrimethylammonium bromide (CTAB) method [23]. Approximately 100 mg of dry leaves were pulverized and then mixed with 1 mL of CTAB buffer (100 mM Tris-HCl, pH 8.0, 20 mM EDTA, 2% CTAB) with 0.1% β-mercaptoethanol added immediately prior to use. The mixture was incubated at 65°C for 60 min and centrifuged for 2 min at 13,000 ×g, after which 600 µL of the supernatant were then transferred to a 1.5-mL microcentrifuge tube. The supernatant was mixed twice with phenol/chloroform/isoamyl alcohol (25:24:1) and shaken for 10 min, followed by centrifugation for 10 min at 12,000 ×g. DNA was precipitated from the aqueous phase by adding a 0.1 volume of 3 M sodium acetate and 2.5 volumes of ethanol. The DNA precipitate was washed twice with 70% ethanol and dissolved in water. The quality of randomly selected DNA samples was determined using a NanoVue Plus Spectrophotometer (GE Healthcare and Life Sciences, Little Chalfont, UK).

Table 1.

Locations and sample size of 20 natural teak populations in Myanmar. Letters in parentheses indicate abbreviated population names. The small m in parentheses indicates the altitude in meters.

DNA sequencing for the development of chloroplast markers

The complete teak chloroplast genome nucleotide sequence was previously determined (GenBank Accession No. NC020098). To detect intra-specific genetic variation, partial genome walking was conducted using eight individuals and 58 PCR genome walking primers (Appendix 1). One individual was selected from each of eight populations, HMB, SP, TDG, NKO, GGW, KTA, HLB and MLM, which nearly cover the entire geographical and altitudinal range of the teak-growing area of Myanmar (Table 1). DNA amplification was performed in a GeneAmp^® PCR System 9700 (Applied Biosystems, Carlsbad, CA, USA) using the following program: 95°C for 1 min; 30 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 60 s; and a final extension at 72°C for 10 min. PCR was performed in a volume of 10 µL, which contained 10 mM Tris-HCl, pH 9.0, 50 mM KCl, 1.5 mM MgCl₂, 0.2 mM of each dNTP, 0.25 µM of each primer, 20 ng of genomic DNA, and 0.1 U of Taq DNA polymerase (Life Technologies, Carlsbad, CA, USA). To remove excess dNTPs and primers, aliquots of each PCR product (5 µL) were purified with 2 µL of ExoSAP-IT reagent (Affymetrix, San Diego, CA, USA) during a 1 h incubation at 37°C, followed by exposure at 80°C for 15 min. Cycle sequencing was performed using the BigDye^® 3.1 Terminator Cycle Sequencing kit (Applied Biosystems) as follows: 96°C for 1 min; 25 cycles of 96°C for 10 s, 50°C for 5 s, and 60°C for 4 min; and a final hold at 4°C according to the manufacturer's instructions. The sequencing primer (3.3 pmol; same sequence as the PCR primer), 1.0 µL of ABI Dye Terminator Ready-Reaction sequencing premix, and 1.5 µL of 5x sequence buffer were added to the template. Excess dye terminators were removed via ethanol precipitation. The extension products were evaporated to dryness under vacuum, resuspended in Hi-Di^™ formamide (Applied Biosystems), incubated for 2 min at 95°C, and loaded onto an ABI PRISM^® model 3130xl DNA sequencer (Applied Biosystems) according to the manufacturer's instructions. For sequence analysis and alignment, we used the software package Sequencher^® 4.10.1 (Gene Codes Corporation, Ann Arbor, MI, USA).

Chloroplast SNP genotyping

We successfully designed three new teak chloroplast SNP markers based on polymorphisms identified in our analysis of intraspecific variation in the teak chloroplast genome. SNP genotyping was conducted via a primer extension-based method with the ABI PRISM SNaPshot Multiplex Kit (Applied Biosystems). Based on the flanking sequences of the target SNPs, a locus-specific primer pair (LSP) and extension primer (EP) were designed for each SNP (Appendix 2). The LSPs were designed to anneal at 58–60°C. Multiplex PCR was performed in a 10 µL total volume with the QIAGEN multiplex PCR Kit (QIAGEN, USA); each reaction contained 2x QIAGEN Multiplex PCR Master Mix, 0.25 µM of each primer, and 20 ng of template DNA. PCR amplification and purification were performed using the same protocols as described above. A single-base extension reaction was conducted using the ABI PRISM SNaPshot Multiplex Kit (Applied Biosystems). The reaction was conducted in a total volume of 10 µL, which contained 2 µL of ExoSAP-treated PCR product, 0.2 µM of each EPs, 1x SNP reaction buffer, and 0.50 µL of ABI PRISM SNaPshot Multiplex Kit. The thermal cycler program was run according to the manufacturer's recommendations. To dissolve the unincorporated fluorescent ddNTP, the SNaPshot product cleanup was performed with shrimp alkaline phosphatase (SAP; USB Corp., USA): 1 µL (1U) was added to 5 µL of SNaPshot product; the mixture was briefly centrifuged and incubated at 37°C for 1 h for digestion, followed by 80°C for 15 min to terminate the reaction. For analysis, 1 µL of SAP-treated product was mixed with 9.8 µL of HiDi^™ formamide and 0.2 µL of GeneScan 120 LIZ size standard (Applied Biosystems). Data were generated after performing capillary electrophoresis on an automated sequencer (ABI 3130 Genetic Analyzer; Applied Biosystems) and analyzed using GeneMapper analysis software, version 4.0 (Applied Biosystems).

Nuclear microsatellite genotyping

For nrSSR analysis, eleven nuclear microsatellite primer pairs [14] that yielded high levels of polymorphism and clear peaks (Teak-A06, Teak-B07, Teak-B03, Teak-H10, Teak-Da12, Teak-F01, Teak-F02, Teak-A11, Teak-F05, Teak-B02, and Teak-C03) were used in this study. Multiplex PCR was performed according to the manufacturer's protocol (QIAGEN Multiplex PCR kit; QIAGEN, Venlo, The Netherlands) in a final volume of 10 µL. Four fluorescently labeled (labels: FAM, NED, PET, and VIC) primers were multiplexed. Amplifications were performed in a thermal cycler (Bio-Rad T100^™; Bio-Rad, Tokyo, Japan) with an initial denaturation step at 95°C for 15 min; 30 cycles of denaturation at 94°C for 30 s, annealing at 57°C for 90 s, and extension at 72°C for 1 min; and a final extension at 60°C for 30 min for all primers. Multiplexed PCR products were diluted 10-fold with sterilized water, and 1 µL of the diluted product was mixed with Hi-Di^™ formamide and GeneScan 600 LIZ size standard (Applied Biosystems). PCR product electrophoresis was performed on an ABI 3130xl sequencer (Applied Biosystems), and alleles were scored using GeneMapper software. If many artifacts and peaks with abnormally high fluorescence intensities were detected, independent PCR reactions were conducted to confirm the amplification.

Data analysis

For chloroplast data analysis, the gene diversity of each population was calculated using CONTRIB ver. 1.02 [24]. Haplotypes occurring at a frequency of <5% of the total samples were defined as rare haplotypes [25]. Two measures of population differentiation, G_ST and N_ST, were compared using the programs PERMUTE and CPSSR version 2.0 [26]. G_ST incorporates haplotype frequencies, whereas N_ST accounts for differences between haplotypes. To infer the genealogical relationships among haplotypes, a network analysis based on parsimony methods was conducted using TCS [27]. The proportion of each haplotype per population was computed and mapped to visualize the haplotype distribution. The CONTRIB version 1.02 software package [24] was used to assess the conservation value of a particular population; this estimates the contribution of each population to the total diversity in terms of genetic diversity and differentiation.

For nrSSR analysis, FSTAT ver. 2.9.3 [28] was used for the calculation of genetic diversity parameters such as the expected heterozygosity (H_E), allelic richness (R) and fixation index (F_IS) of each population. Moreover, genetic diversity parameters (H_E and R), fixation indices (F_IS), genetic differentiation among populations (F_ST) of each group, and the significance of differences in each parameter among regions were measured using the same software. The Lositan software package [25] was used to detect outlier loci (non-neutral loci) in order to evaluate the neutral genetic diversity of teak. Genetic clusters of 20 natural populations and F_ST estimates between the clusters were generated using the program STRUCTURE ver. 2.3 [29], which was run for 15 iterations using the admixture model and a number of clusters (K) ranging from 1 to 10 with a burn-in length of 100,000 and a run length of 100,000 iterations. The most suitable number of clusters was inferred from the delta (K) statistic [30]. The average proportion of ancestry among individuals in each population was calibrated. Bottleneck version 1.2.02 [31] was used to examine the bottleneck effect on studied populations.

Putative genetic boundaries within the sampling regions were identified using Barrier software version 2.2 [32] according to the Euclidean distance in the cpSNP dataset and Nei's genetic distance in the nrSSR dataset. For an analysis of molecular variance (AMOVA), populations were grouped based on the combined cpSNP and nrSSR barrier results, and hierarchical AMOVAs were performed to determine the different genetic hierarchies among groups, among populations within a group, and among individuals within populations for both markers using Arlequin ver. 2.0 [33]. Significance testing was conducted with 1,000 permutations at each level. Spatial genetic structure was investigated using the logarithm of geographic distance in the range of 32–889 km and the pair-wise genetic distance F_ST/(1 - F_ST), which was based on Mantel's test as performed with 1,000 permutations using GeneAlex 6.5 [34].

Chloroplast polymorphism

A total of 43,734 base pairs, excluding missing data and gaps, were sequenced in eight individuals using 58 walking primers (Appendix 1). Three variable sites, all SNPs, were detected; two were in intergenic spaces (IGS), a G/T mutation in the psbK-psbI IGS and AAA/TTT in the trnK-rps16 IGS and one, C/A, occurred in the rpl16 intron. Fragments of these three variable sites were used to develop chloroplast markers with which to investigate geographic variations in Myanmar teak (Appendix 2).

Geographic variation revealed by cpSNP markers

A cpSNP variation analysis identified four haplotypes in 20 natural populations (Table 2). The frequency of H1 was highest (84%), and this haplotype was the most common in all populations except TDG and HMB, whereas H2 was the rarest haplotype (3%). Although the frequencies of H3 (6%) and H4 (7%) were relatively low, these haplotypes were detected in many populations. No private haplotype was found in any population. The gene diversity ranged from 0.000 to 0.565, with a mean value of 0.173 (Table 3, Fig. 1). The highest level of gene diversity was observed in the KTA population from the upper region, followed by the HMB and TDG populations from the lower and central regions, respectively. At the regional level, the lower region and eastern regions harbored the highest (0.249) and lowest gene diversity (0.075), respectively. Half of the 20 populations, mostly those from the central region, were monotypic and fixed for the common haplotype (H1), whereas the remaining populations were polytypic and contained two to four haplotypes. The overall F_ST value using cpDNA markers was 0.293 (Table 3). The measure of population differentiation, G_ST = 0.346, was higher than the N_ST value of 0.302. The geographic distribution of the chloroplast haplotype frequency in each population and the star-shaped parsimony haplotype network, in which each haplotype is distinguished by a single mutation, are shown in Fig. 1. The rare haplotype H2 occurred mainly in the southern region; it was also present in HMB, TDG, and POL from the lower, central, and upper regions, respectively, but not in the eastern region. The H3 haplotype was absent from the central populations, and H4 was not found in the southern and eastern regions. The highest frequency of H4 was observed in TDG, the only polytypic population in the central region. The assessment of populations to prioritize for conservation depicted stronger contributions of the HMB, TDG, KTA, and POL populations to the total diversity in terms of their high levels of genetic diversity (Fig. 2).

Table 2.

Polymorphic sites and chloroplast DNA haplotypes according to a SNaPshot analysis.

Fig. 1.

Myanmar teak haplotype network and geographic variation patterns as revealed by chloroplast single nucleotide polymorphism markers. The population numbers and abbreviations correspond to those used in Table 1. Black lines show the administrative boundaries of states and divisions. Yellow lines show the regional boundaries and overlap some administrative boundary lines. The circle sizes are proportional to the sample size of each population.

Fig. 2.

Contribution (CT%) of each population to total diversity in terms of genetic diversity and differentiation from other populations. The population abbreviations correspond to those used in Table 1.

Genetic structure revealed by nrSSR

Of the 11 nuclear microsatellite loci, one outlier (Teak-B03) was excluded from the nuclear microsatellite dataset analysis. Regarding the average genetic diversity statistics, the expected heterozygosity (H_E) and allelic richness (R) were 0.652 and 5.76, respectively (Table 3). The range of genetic diversity indices in these populations among 10 loci varied with respect to both the expected heterozygosity (H_E), which ranged from 0.530 to 0.763, and the allelic richness (R), which ranged from 3.09 to 8.29. The minimum and maximum fixation indices (F_IS) were -0.148 for the HMB population and 0.076 for the MMK population. The genetic differentiation among 20 populations was 0.094. At the regional level, the genetic diversity estimates and F-statistics of each group are shown in Table 3. The lower region possessed the highest genetic diversity in terms of H_E (0.712) and R (7.28), whereas the southern region exhibited the lowest for both values (Table 3). The highest F_ST value (0.097) and lowest F_IS (-0.091) were found in the eastern populations, whereas the lowest F_ST value (0.017) and highest F_IS (0.040) were found among the lowland populations. Significant differences among regional groups were detected in terms of allelic richness (p = 0.012), expected heterozygosity (p = 0.011), and F_IS (p = 0.02).

When a Bayesian clustering approach was implemented via STRUCTURE analysis, a maximum delta K value of 4 was detected. Using this value, we grouped the 20 populations of Myanmar teak into four genetic clusters (Fig. 3). The greatest genetic differentiation between clusters (F_ST) was observed in cluster III (0.310), followed by clusters IV (0.264), II (0.128), and I (0.037). Cluster I predominated in the lower and central regions, cluster II was dominant among populations from the upper and eastern regions apart from the KTA and NKO populations, and cluster IV was present mainly in the southern region and TDG population. Cluster III was predominant in two populations (POL and NKO) that were from different regions but were geographically close. Of these 20 populations, 16 exhibited a significant bottleneck effect under the stepwise mutation model, but only two populations (HMB and YGN) showed significant bottleneck effect under infinite assumption model.

Table 3.

Statistical summary of the diversity revealed by the chloroplast single nucleotide polymorphism and nuclear simple sequence repeat markers in 20 teak populations. A, number of alleles; R, allelic richness; H_O, observed heterozygosity; H_E, expected heterozygosity; F_IS, fixation index; F_ST, genetic differentiation among populations. Populations with asterisks indicate those with bottleneck effects at different significance levels. ^* p < 0.05, ^** p < 0.001, ^*** p < 0.001.

Fig. 3.

Myanmar teak genetic structure and geographic variation as detected using nuclear simple sequence repeat markers. The population numbers and abbreviations correspond to those used in Table 1. Yellow lines show the regional boundaries and overlap some black lines that indicate the administrative boundaries of states and divisions. The sizes of the circles are proportional to the sample size of each population.

Geographic variation and genetic differentiation

The KTA, TDG, and HMB populations were not included in any cpSNP marker groups (Fig. 4a), and the KTA, GGW, and NKO populations also were separated from the nrSSR marker groups (Fig. 4b) according to the barrier results. The boundaries between the lower and southern regions as determined with nrSSR markers were used for grouping because of the high cluster IV F_ST value, although the genetic boundaries between those regions were incongruent. Therefore, an identification of genetic boundaries based on cpSNP markers and nrSSR generally yielded four groups (Fig. 4c) that were used for AMOVA. Those groups likely corresponded to geographic regions (e.g., upper and eastern regions group, central populations group, lower region group, and southern region group). The AMOVA results revealed that the same amount of genetic variation accounted for the different groups generated using both types of markers. Significant genetic variations were observed at all hierarchical levels (Table 4). No significant isolation by distance (IBD) pattern was detected using cpSNP markers (r = -0.049, p = 0.503), whereas the nrSSR markers revealed a significant pattern (r = 0.333, p = 0.002).

Table 4.

Results of an analysis of molecular variance (AMOVA) involving the chloroplast single nucleotide polymorphism (cpSNP) and nuclear simple sequence repeat (nrSSR) markers. Groups are identified based on the genetic boundaries from the cpSNP and nrSSR marker barrier results.

Fig. 4.

a: Identification of genetic boundaries based on chloroplast single nucleotide polymorphism (cpSNP) markers. The segment values indicate the Euclidean distances between the separated populations. The segment thicknesses correspond to the Euclidean distance values. b: Identification of genetic boundaries based on nuclear simple sequence repeat (nrSSR) markers. The segment values represent Nei's genetic distances between the separated populations. c: Identification of genetic boundaries based on the combined cpSNP and nrSSR barrier results. Black lines indicate the congruent boundary segments for both the cpSNP and nrSSR barrier results. Red lines indicate the boundary segments demonstrated by the nrSSR barrier results.

Geographic variation revealed by chloroplast SNP markers

Due to teak's limited seed dispersal ability, we expected clear geographic patterns in Myanmar teak would be detected using cpSNP markers. However, the significantly higher N_ST relative to G_ST indicated the lack of a phylogeographic structure and randomized teak haplotype distribution at cpSNP markers. Using maternally inherited markers, tree species are typically shown to exhibit strong genetic differentiation among populations, with F_ST values exceeding 0.6 (average: 0.764) [18]. Therefore, the overall cpSNP-based F_ST of Myanmar teak populations (0.293) was quite low in term of tree species. The lack of phylogeographic structure, star-shaped haplotype network pattern, and low F_ST values indicated recent colonization from a small number of Myanmar teak founders and therefore insufficient time for the occurrence of mutation in long-lived species. The predominant distribution of a single haplotype (H1) throughout much of the range suggested a bottleneck effect at the time of postglacial recolonization. In fact, we detected significant bottleneck effects in 16 of the 20 populations under stepwise mutation model in this study. In contrast to the distribution of H1, the rare haplotype H2 occurred mainly in the southern region, but not the eastern region. The haplotype H3 was absent from the central region, whereas the haplotype H4 was not found in southern and eastern regions. Therefore, rare haplotypes exhibit geographic patterns, despite the wide distribution of the major haplotype throughout Myanmar.

Genetic diversity and genetic structure revealed by nrSSR markers

In the nrSSR analysis, the populations in the lower region exhibited the greatest genetic diversity and allelic richness (Table 3). The lower populations are located near the Bago Mountains, an area that is considered the home of Myanmar teak [35]. In temperate or boreal forests, populations that possessed high levels of genetic diversity, allelic richness, and geographic clustering of rare alleles acted as refuges during the LGM [36] and might have derived from an ancestral population or an admixture of different refuges [37]. In contrast, the southern populations exhibited the lowest diversity found in Myanmar (Table 3). A small population size and isolation generally result in low genetic diversity [37]. Forest fragmentation has occurred in the southern region and resulted in decreased population sizes and increased isolation, which might explain the low genetic diversity in southern populations.

Four genetic clusters that basically corresponded to geographic regions were detected in our Myanmar teak sample. The genetic structures of tree species are profoundly influenced by geographical isolation and environmental factors such as topography (including the effects of tectonic shifts), river course dynamics, and sea level fluctuations [38]. For instance, the southern populations were separated from the other groups by rivers and high mountain ranges. Those geographic features might have acted as genetic barriers between the southern and other populations.

Comparison of findings obtained using cpSNP and nrSSR markers

In this study, we observed four main genetic groups of Myanmar teak using both types of markers. Genetic variation, which accounted for the different groups (Table 4), was relatively low in terms of cpSNP markers when compared to the variation observed in other tree species such as Neobalanocarpus heimii [39] and Caryocar brasilense [40]. A genetic bottleneck in Myanmar teak might account for the low level of cytoplasmic variation. In the chloroplast dataset, the rare haplotype H2 was distributed throughout all southern populations; these populations possessed the lowest genetic diversity and represented the major component of cluster IV, with a high F_ST value according to nrSSR data. The haplotype H2 might have derived from the common ancestor H1, and cluster IV might similarly have derived from cluster I.