Mitogenome analysis reveals a complex phylogeographic relationship within the wild tiger population of Thailand

Fig. S2 The alignment of the control region of the WEFCOM and Malayan tigers reported in this study with previous reports. These sequences contain different lengths of the control region resulting in different sizes of the complete mitochondrial genome. The different lengths of the control region are a result of the variable number of tandem repeats in the repetitive sequence box 3 (RS-3). HVS = hypervariable site, RS = repetitive sequence, CCR = central conserved region, CSB = conserved sequence box.


INTRODUCTION
Morphological and genetic studies of wild tigers Panthera tigris have resulted in the designation of 6 extant subspecies: Amur P.  (Maźak 1981, Luo et al. 2004).Thailand is one of 7 countries situated within the geographic range of tigers and contains sufficient habitat to support 2 source sites and up to 250 tigers (DN P 2010, Walston et al. 2010).The majority of tigers in Thailand are of the subspecies P. t. corbetti, and P. t. jacksoni occurs in the south of the country separated by the Isthmus of Kra (Luo et al. 2004).Most of the remaining tiger habitat in Thailand is within the Western Forest Complex (WEFCOM), a World Heritage Site which covers an area of approximately OPEN PEN ACCESS CCESS 18 000 km 2 (Simcharoen et al. 2007).However, despite government protection under the Wild Animal Reservation andProtection Act B.E.2535 (1992) of Thailand, the tiger population here is still declining due to habitat loss and human persecution.
Due to the critical importance of regions such as WEFCOM to tiger conservation in Thailand, National Park staff have been monitoring the population through regular camera trap studies, and they have also begun to collect non-invasive samples (faeces and hair) for genetic analysis.Recent genetic studies of mitochondrial DN A variation in the WEFCOM population and in southern Thailand (unpubl.data) have revealed the presence of 2 new haplotypes that have not been previously reported.One haplotype (GenBank Accession No. KC879297) was discovered in tiger samples from WEFCOM, and the other (KC879296), which was detected in an individual from southern Thailand, represents a new haplotype of P. t. jacksoni.
Here we describe a new study as an extension of this work, to characterize the complete mitogenome of tigers from WEFCOM and southern Thailand, and to demonstrate the phylogenetic relationship between these haplotypes and tigers in other populations.The results of this and subsequent genetic studies will act as the basis for a national genetic database, which could be used for both research and forensic applications in Thailand.

Primer design
Most of the primers used in this study were newly designed from the alignment of cymt (GenBank DQ151150) and F2 numt sequences (GenBank DQ151151) (Kim et al. 2006).Primers originally developed for Ursidae (Delisle & Strobeck 2002) were used to amplify the region between the 2 ends of F2 numt (ATP6 − ND5).Full details of the primers used are given in Table S1 in the Supplement, at www. intres.com/ articles/ suppl/n030 p125_ supp.pdf.

Samples and DNA extraction
Hair samples (n = 3) were collected from WEFCOM during field surveys by N ational Park staff in 2009.These 3 samples shared the same haplotype based on aligned sequences from ND6, Cyt b and the control region; we therefore chose only 1 sample for further analysis (ID = WEFCOM).Three muscle samples were also obtained from 3 tiger carcasses confiscated by police from poachers in N arathiwat province, southern Thailand, in 2009.One sample was obtained from each individual, and 1 sample (ID = Ti2) was identified as a new haplotype using the aligned sequences from ND6, Cyt b and the control region.This sample was therefore chosen for subsequent analysis.A map of Thailand showing the sampling sites is shown in Fig. S1 in the Supplement at www. int-res.com/ articles/ suppl/ n030 p125_ supp.pdf.
DNA extractions were performed following a modified method of Carter & Milton (1993).Briefly, the hair or muscle was cut into small pieces, digested in lysis buffer, and incubated overnight at 56°C.The supernatant was then collected and mixed with silica diatom.After shaking and centrifugation, the resultant pellet was washed with washing buffer.The pellet was then dried and DN A eluted with elution buffer.

PCR conditions
Each DN A template was amplified with a singledirection PCR prior to specific fragment amplification.PCRs were performed in a 20 µl reaction volume containing 10 µl of 2× Phusion buffer (Finnzymes), 200 µM dNTP, 0.25 µl of each 1 µM reverse primer, 0.2 µl of Phusion High-fidelity DN A polymerase (Finnzymes), 2 µl of DN A template (diluted to 5 ng µl −1 ), and the relevant amount of purified water.PCR cycling conditions were as follows: a pre-denaturation step at 98°C for 3 min, then 35 cycles of 98°C for 30 s, annealing at 58°C for 30 s, extension at 72°C for 2 min, then 72°C for 5 min.The higher denaturation temperature of 98°C is recommended for this highfidelity Taq enzyme.
The second round of PCR was performed in 50 µl volumes containing 25 µl of 2× buffer, 200 µM dNTP, 0.1 µl of each of 10 µM forward and reverse primer, 0.5 µl of Phusion High-fidelity DN A polymerase (Finnzymes), 1 µl of the PCR product from the first round of single-direction PCR and purified water.PCR was performed for 40 cycles under similar conditions to those described above, apart from reducing the extension time to 1 min.PCR products were purified using the silica diatom method as described above (Carter & Milton 1993) and sequenced in either forward or reverse directions.This 2-step PCR protocol was repeated twice for each mitochondrial region except the control region, which was repeated 3 times.

Data analysis
Amplified mtDN A fragments were assembled into complete mtDN A sequences using BioEdit (Hall 1999).Haplotypes identified in this study were then analysed jointly with haplotypes of 4078 bp of concatenated sequence from extant tiger subspecies (Luo et al. 2004); see Fig. 1.Three felid species related to tigers, viz.snow leopard Panthera uncia (GenBank N C_010638), leopard P. pardus (GenBank EF551002) and clouded leopard Neofelis nebulosa (GenBank DQ257669), were used as outgroups in the analysis (Table S2 in the Supplement).Phylogenetic relationships between the 2 tiger samples in this study and published P. tigris haplotypes were assessed using neighbourjoining (N J), maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference.The Bayesian information criterion implemented in jModelTest v2.1.4(Posada 2008) was used to select the best-fit substitution model for ML and Bayesian analysis.The Hasegawa-Kishino-Yano model + gamma distribution (base = [0.32500.2911 0.1338], nst = 2, tratio = 20.6989,rates = gamma, shape = 0.1520, ncat = 4, pinvar = 0) was selected as the best-fit model.Using these model parameters, Bayesian inference was executed in MrBayes v3.2.5 (Ronquist & Huelsenbeck 2003, Ronquist et al. 2012) for 2 000 000 generations with trees sampled every 100 generations and the first 25% generations discarded as burn-in.N J, MP and ML analyses were performed in PAUP* v4a142 (Swofford 2002) with 2000 repetitions for MP and N J, and 100 repetitions for ML.Heuristic search with random additions of taxa and tree-bisection-reconnection branch swapping was employed.Trees were then viewed and manually annotated with FigTree v1.4.2.The median-joining phylogenetic network was computed using N etwork 4.6.1.2(Bandelt et al. 1999) to analyse the relationship between tiger haplotypes from this study and those from Luo et al. (2004).

RESULTS
The size of the complete mitogenome in the WEFCOM and Ti2 individuals (GenBank KJ508412 and KJ508413) was 16 954 and 17 008 bp, respectively (Table S3 in the Supplement).The difference in mitogenome size between this and other studies was mainly due to variation in the tandem repeats TACACACG and TATACACG in the repetitive sequence box 3 of the control region (Fig. S2 in the Supplement).Thirty variable sites were identified in the whole mitogenome excluding the control region (21 C/T and 9 A/G nucleotide substitutions), and there were 22 new variable sites (11 in WEFCOM and 11 in Ti2) when comparing our results to other studies (Luo et al. 2004, Driscoll et al. 2009, Zhang et al. 2011, Kitpipit et al. 2012, Mondol et al. 2013, Sun et al. 2015, Xue et al. 2015) (see Table S4 in the Supplement).The WEFCOM haplotype carried 1 Panthera tigris altaica-specific single nucleotide polymorphism (SNP; 14711A) as described by Luo et al. (2004), whereas Ti2 contained no subspecies-specific SNPs.
In each sample, we found 2 distinct sequences with numt contamination in the region of COX2-ATP6.No heteroplasmy was found in either sequence.The first sequence (amplicon 9.1) generated from primers described for Ursidae (Delisle & Strobeck 2002) contained an F2 numt sequence from the COX1 to COX2 gene as found in the report by Zhang et al. (2011) and could not be assembled with the flanking fragments.The second sequence (amplicon 9.2), amplified by the newly designed primers in this study, contained F2 numt contamination in the COX2 to ATP8 genes, similar to the results of Kitpipit et al. (2012) (see Table S1 and Figs.S3 & S4 in the Supplement), and perfectly assembled to the adjacent fragments.The last domain of sequences covering part of the ATP6 gene also differed (Figs.S3 & S4).We therefore excluded the mitochondrial fragment between positions 7399 and 12 707 (based on the WEFCOM mitogenome) from further analysis as it may have contained additional numt sequences.
Results from the MP, N J and ML trees and the Bayesian inference analyses on 4078 bp of concatenated mtDN A indicate that the WEFCOM haplotype is nearly equidistant to the fixed haplotype of P. t. altaica (8 different nucleotides) and a haplogroup consisting of all reported P. t. corbetti mtDNA haplotypes (9 different nucleotides; Fig. 1A).The high Bayesian posterior probability and bootstrap support of this monophyletic cluster suggests a recent common ancestor of tigers now distributed in  2004).(A) Results from maximum parsimony (MP), neighbour-joining (NJ), maximum likelihood (ML) and Bayesian analyses were similar, so only the MP tree is shown here.The numbers above each branch represent the statistical support in percentages for the MP, NJ, ML, and Bayesian analyses, respectively.Dashes represent statistical support values < 50.The numbers below each branch represent the number of changes.(B) The median-joining network demonstrates the evolutionary relationship between our 2 samples and previous reports (Luo et al. 2004).Each circle represents a haplotype, and the marks on each branch indicate mutation points.Circle size is not proportional to the haplotype frequency.Both the MP phylogenetic tree and the median-joining network support an additional Indochinese haplotype (WEFCOM), which evolved from a common ancestor between P. t. altaica and P. t. corbetti.The new Malayan haplotype (Ti2) descended from the same historic ancestor as P. t. jacksoni Indochina and Far Eastern regions of Asia, and substantial population genetic substructure within Indochinese tigers.The southern haplotype, Ti2, clustered with JAX1 and JAX2 haplotypes and formed a monophyletic group within the Malayan tiger subspecies P. t. jacksoni.These results were supported by the median-joining network analysis and suggested a complex phylogeographic pattern within P. t. corbetti in that the WEFCOM haplotype is different from other P. t. corbetti haplotypes, and that the Ti2 haplotype belongs to P. t. jacksoni (Fig. 1B).
The analysis of the longer mtDN A sequence showed that the WEFCOM haplotype appears to be more closely related to P. t. altaica than to other P. t. corbetti haplotypes with high statistical support in all algorithms.However, the number of changes from the WEFCOM haplotype to the haplogroup of P. t. altaica (15 nucleotides) is larger than that to the haplo group of P. t. corbetti (13 nucleotides; Fig. 2).This result further indicates a close and complex relationship between the 2 subspecies.All analyses using the MP, ML and Bayesian approaches consistently supported the finding that the lineage of P. t. jacksoni diverged earlier than the splitting of P. t. corbetti and P. t. altaica.

DISCUSSION
Tiger populations with a sufficient number of breeding individuals to support the wider landscape are currently thought to occur in only 8 countries, including Thailand (Goodrich et al. 2015).Therefore, it is vitally important to understand the dynamics of the wild tiger population in Thailand.In the present study, we provide the first report of the complete mitogenome of the Malayan tiger and of novel haplotypes in the wild Indochinese tiger population in Thailand.We also demonstrate the phylogenetic relationship of these 2 subspecies to other tiger populations in Asia.
The results of our phylogenetic analyses are in agreement with the subspecies classification previously described (Luo et al. 2004), but indicate phylogenetic differences between the WEFCOM haplotype and other Panthera tigris corbetti haplotypes found across Asia.We found that the WEFCOM haplotype contained the 14711A SNP of P. t. altaica, a common SN P shared among P. t. altaica and P. t. corbetti (15756C), and no subspecies-specific SN Ps described for P. t. corbetti (Luo et al. 2004).Our results suggest that the WEFCOM haplotype is   (Sun et al. 2015) and the WEFCOM individual in this study provides further support for a recent common ancestor for these 2 haplotypes, and hints at further local Indochinese substructure.Thus, tiger populations in northern Thailand may contain some of the haplotypes present in the ancestral population when tigers expanded out of China and into Russia and mainland Asia.Thailand is now the only country in the region with a reasonable population of P. t. corbetti tigers, as they have largely been poached from other neighbouring countries such as China, Laos and Cambodia.Thus, the existence of new haplotypes in Thailand also highlights its critical geographic position between the populations of P. t. altaica in the Russian Far East and P. t. jacksoni in Malaysia.
The discovery of 2 distinct sequences containing numt in this study is notable, as it suggests the presence of additional numt in multiple fragments, as has been reported in species such as the domestic cat, human and chimpanzee (Lopez et al. 1994, 1996, Tourmen et al. 2002, Antunes et al. 2007, Hazkani-Covo & Graur 2007).Therefore, future studies involving the complete mitochondrial genome of tigers must take care to check for confounding sequences between the 2 ends of F2 numt, which should be excluded from any downstream phylogenetic analysis.

CONCLUSIONS
Here, we have provided an account of the complete mitogenome of Panthera tigris corbetti and P. t. jacksoni, and have provided further information on the phylogenetic relationship between these 2 tiger subspecies in Thailand.Our results also provide evidence for the close relationship between P. t. altaica and P. t. corbetti, and support the importance of WEFCOM to the Asian tiger population.

Fig. 1 .
Fig. 1.Phylogeny of mtDNA haplotypes of extant tiger Panthera tigris subspecies based on 4078 bp of concatenated sequences following Luo et al. (2004).(A) Results from maximum parsimony (MP), neighbour-joining (NJ), maximum likelihood (ML) and Bayesian analyses were similar, so only the MP tree is shown here.The numbers above each branch represent the statistical support in percentages for the MP, NJ, ML, and Bayesian analyses, respectively.Dashes represent statistical support values < 50.The numbers below each branch represent the number of changes.(B) The median-joining network demonstrates the evolutionary relationship between our 2 samples and previous reports(Luo et al. 2004).Each circle represents a haplotype, and the marks on each branch indicate mutation points.Circle size is not proportional to the haplotype frequency.Both the MP phylogenetic tree and the median-joining network support an additional Indochinese haplotype (WEFCOM), which evolved from a common ancestor between P. t. altaica and P. t. corbetti.The new Malayan haplotype (Ti2) descended from the same historic ancestor as P. t. jacksoni

Fig. 2 .
Fig. 2. Maximum parsimony (MP) phylogenetic tree based on 10 140 bp of mitochondrial sequence (Zhang et al. 2011, Kitpipit et al. 2012, Sun et al. 2014) excluding the control region (CR) and the numt-contaminated fragments.The numbers above each branch represent the bootstrap values for the MP, neighbour-joining, maximum likelihood and Bayesian analyses, respectively.Dashes represent bootstrap values < 50.The number below each branch represents the number of changes.The placement of the WEFCOM haplotype was consistent in all methods with high statistical support (Xue et al. 2015)o P. t. altaica and that there may be more complex genetic structure in Indochinese tigers than previously thought.As expected, the Ti2 haplotype sampled in southern Thailand formed a monophyletic group with JAX haplotypes and indicates that P. t. jacksoni diverged from a common Indochinese ancestor earlier than the divergence between P. t. corbetti and P. t. altaica.A recent genetic study of tigers in the Sunda region suggests that the tiger subspecies now found in much of southeast Asia and the Russian Far East (P.t.corbetti, P. t. jacksoni and P. t. altaica) most likely evolved from a common ancestor in a second wave of expansion and divergence from P. t. amoyensis in China(Xue et al. 2015).Amur tigers (P.t. altaica) and Indochinese tigers (P.t. corbetti) have also been shown to be sister taxa separated by only a few mitochondrial steps.The close phylogenetic relationship between 1 wild P. t. altaica (GenBank KF297576) from China