Species, sex and geo-location identification of seized tiger (Panthera tigris tigris) parts in Nepal—A molecular forensic approach

Tiger (Panthera tigris) populations are in danger across their entire range due to habitat loss, poaching and the demand for tiger parts. The Bengal tiger (Panthera tigris tigris) is an endangered apex predator with a population size estimated to be less than 200 in Nepal. In spite of strict wildlife protection laws, illegal trade of tiger parts is increasing; and Nepal has become one of the major sources and transit routes for poached wildlife parts. Identification of wildlife parts is often challenging for law enforcement officials due to inadequate training and lack of available tools. Here, we describe a molecular forensic approach to gain insight into illegally trafficked tiger parts seized across Nepal. We created Nepal’s first comprehensive reference genetic database of wild tigers through the Nepal Tiger Genome Project (2011–2013). This database has nuclear DNA microsatellite genotype and sex profiles, including geo-spatial information, of over 60% (n = 120) of the wild tigers of Nepal. We analyzed 15 putative cases of confiscated poached tiger parts and all were confirmed to be of tiger. Ten samples were identified as male and five were female. We determined probable geo-source location for 9 of the 14 samples with 6–8 nuclear DNA microsatellite loci using inferences from four different statistical assignment methods. Six samples were assigned to Bardia National Park and one of these was an exact match to a female tiger previously profiled in our fecal DNA reference database. Two tiger samples were assigned to Shuklaphanta Wildlife Reserve and one to Chitwan National Park. We are unable to definitively assign five tiger samples which could be offspring dispersers or might have come from tiger population outside of Nepal. Our study revealed that the western region, particularly Bardia National Park, is a poaching hotspot for illegal tiger trade in Nepal. We present feasibility of using molecular forensic based evidence to incriminate criminals in a court of law in the fight against wildlife crime.


Introduction
The Bengal tiger (Panthera tigris tigris) is the most prevalent and endangered tiger subspecies found in the Indian subcontinent [1,2]. In just over a century, wild tiger populations have dramatically fallen (>97%) and currently about 3,200 tigers are left in the wild [1,2]. Widespread deforestation, habitat fragmentation and loss, prey depletion, and poaching are the major causes of tiger population decline [3][4][5][6]. Moreover, diseases such as canine distemper virus (CDV) could pose further threat to already declining tiger populations [7]. Poaching and illegal wildlife trade present the greatest threats to the survival of tigers across their range [8][9][10]. With increasing demand for wildlife parts and illicit trade of wildlife, Asia has become one of the major hubs for wildlife crimes [4,11,12]. Nepal is one of the major sources of illegal wildlife parts [13,14]. Each year lucrative, illicit wildlife commodities including tiger parts such as skin, bone, claws, teeth, blood, genitals, and meat are illegally trafficked to Chinese markets [14,15] and used in traditional Chinese medicine [15][16][17]. Wildlife crime in Nepal is driven by a multitude of factors, including poverty, lucrative financial rewards from illicit wildlife products in black markets, porous borders, lack of conservation awareness, and poor law enforcement for the prevention of wildlife crimes [18,19].
Often the only evidence recovered from wildlife crime scene is traces of blood, pieces of meat, skin, fur, bones and other forms of biological materials. Many of the seized specimens are virtually impossible to identify based on morphological analysis alone. Lack of precise identification of confiscated specimens has often resulted in no convictions in court of law in Nepal. Hence, there is an urgent need to use advanced forensic techniques to identify seized wildlife specimens to determine their species and region of origin.
Molecular forensics can provide information on species, sex and individuals of poached animals with relatively high accuracy [20][21][22][23][24][25][26][27][28][29]. Using a geo-spatial reference genetic database, it is possible to determine the source of the seized parts using various statistical and computational analyses.
Applications of DNA profiling in wildlife forensics were successfully applied for whales and dolphins [30,31]; mule deer, white-tailed deer, moose, caribou and American black bear [32]; wild boar and wolves [33,34]; African elephants [20,23,35] and leopards [24]. For tigers, molecular forensic techniques have been, for example, used to investigate a tiger bone smuggling case in China [36] and illegal sale of tiger meat from a circus tiger [37]. In Nepal, utility of genetic tools in both population studies and forensics have recently been explored for tigers, through the Nepal Tiger Genome Project (NTGP) (2011-2013) [38].
Genetic tools combined with geo-location baseline data on species can discern geographic origin of confiscated wildlife parts and their derivatives [20,23,24,29,32,35,39]. This kind of information is highly valuable for conservation managers to prioritize their efforts in effective anti-poaching activities. Best practices in anti-poaching efforts require scientifically sound wildlife monitoring efforts and the effective sharing of information amongst conservation managers and relevant stakeholders [40]. The Government of Nepal (GoN) has asserted that the incidences of seized tiger parts have dramatically increased over the last years. Nepal is both a source and a conduit for wildlife parts trafficking into China/Tibet and other areas. Hence there is a need to utilize effective tools such as molecular forensics to track poaching and illegal trade in tiger parts within the country and beyond. The GoN has proposed establishing a system to better track records of poaching incidences, confiscation and enforcement of the Convention of International Trade in Endangered Species (CITES) of Wild Fauna and Flora by increasing the in-country capacity for identification of animal parts and their derivatives using effective and modern DNA-based forensic tools [41].
The primary goal of this study is to apply a molecular forensic approach to identify species, sex and the population of origin of 15 putative tiger parts confiscated by the Central Investigation Bureau (CIB) of Nepal Police. To identify the geo-location source of these samples, we used inferences from four different assignment methods and a previously developed geo-spatial genetic reference database of wild tigers of Nepal [38]. The database includes genetic and geo-location information on 120 wild tigers of Nepal, covering 60% of the estimated 200 tigers across the Terai Arc Landscape (TAL) sampled across the main tiger habitats ( [42,43].

Source of forensic samples
The laboratory of the Center for Molecular Dynamics Nepal (CMDN) received a total of fifteen seized putative tiger parts, which included thirteen skin pieces, and two blood smeared knives (S1 Fig) seized by the CIB during investigative operations (2014-16) in southern Nepal (S1 Table). These samples were labeled, photographed and stored with desiccant at room temperature. This study has been authorized through study permit provided by the Department of National Parks and Wildlife Conservation, Ministry of Forests and Soil Conservation-Government of Nepal.

Geo-spatial baseline genetic database of wild tiger in Nepal
NTGP created Nepal's first comprehensive tiger reference genetic database by collecting and analyzing fecal samples (n = 770;CNP = 420, BNP = 116, SWR = 79, PWR and wildlife corridors = 155) across the TAL region of Nepal (Fig 1 and S2 Table). This database includes species, sex and individual DNA microsatellite profiles of 120 wild tigers [44].
Genetic analyses of forensic samples DNA extraction. DNA was extracted using DNeasy Blood and Tissue Kit following the manufacturer's instruction (Qiagen Inc., Germany) [46]. A 1 cm 2 piece of skin (forensic) was digested overnight at 56˚C in 180 μL ATL and 20 μL Proteinase K solution. Extracted DNA was precipitated with ethanol and purification was conducted in a spin column following the kit's protocol. For the blood smeared knives, samples were taken by swabbing the smeared surface with a sterile cotton swab saturated with phosphate buffer saline. DNA from swab samples were extracted following the same protocol used in skin samples. Each batch of DNA extractions included a negative extraction control. Precaution was taken to avoid contamination. Extracted DNA was stored at -20˚C.
Species identification. Tigers were identified using a PCR assay that used tiger specific mtDNA Cytochrome-b (CYT-B) primers [47]. A total of 7 μL PCR reaction was prepared containing 3.5 μL of 2X Qiagen multiplex mastermix (Qiagen Inc., Germany), 0.7 μL Q-solution, 200 nM each CYT-B primer sets and 1.5 μL of purified DNA template. The thermocycling condition was 95˚C for 15 min followed by 35 cycles of 94˚C for 30 sec, 59˚C for 90 sec and 72˚C for 90 sec with the final extension at 72˚C for 10 min. Amplified 162 bp target PCR product was visualized under 1.5% agarose gel electrophoresis (S2 Fig). Sex identification. Sex of identified tiger samples was determined by amplifying the Amelogenin (AMEL) gene of sex chromosomes [48]. A 7 μL PCR reaction was prepared containing 3.5 μL of 2X Qiagen multiplex mastermix (Qiagen Inc., Germany), 0.7 μL Q-solution, 200 nM each AMEL forward and reverse primers and 1.5 μL of purified DNA template. The thermocycling condition was 95˚C for 15 min followed by 45 total cycles of 94˚C for 30 sec, 53˚C for 60 sec and 72˚C for 60 sec with final extension at 72˚C for 10 min. Amplified PCR products were visualized in 3% agarose gel electrophoresis. Female samples yielded a single (214 bp) PCR band, while two bands (194 bp and 214 bp) identified male samples. Each sample was run in triplicate. Samples that had amplification of X and Y alleles on at least two out of three replicates were identified as male. Samples that had only X allele amplification in all three replicates were assigned as female [49]. Individual identification. Individual tigers were identified using a panel of polymorphic microsatellite markers (n = 10) developed from the domestic cat (Feliscatus) and tiger genomes [50][51][52]. The primers are arranged in a single multiplex panel. PCR amplification was carried out in a 7 μL reaction volume containing 3.5 μL of Qiagen multiplex master-mix (Qiagen Inc., Germany), 0.7 μL of Q solution, each ten microsatellite primers (FCA391: 0.2 μM, PttD5: 0.07 μM, FCA232: 0.14 μM, FCA304: 0.07 μM, F85: 0.30 μM, FCA441: 0.14 μM, FCA043: 0.09 μM, F53: 0.49 μM, FCA205: 0.21 μM and PttA2: 0.07 μM) and 2.5 μL of purified DNA template. The PCR thermocycling condition was 95˚C for 15 min with a 10 cycles of touchdown step (94˚C for 30 sec, initial annealing at 62˚C reduced by 0.5˚C in each cycle for 90 sec and extension at 72˚C for 60 sec), followed by 25 cycles of 94˚C for 30 sec, 57˚C for 90 sec and 72˚C for 60 sec, and a final extension at 72˚C for 10 min. Amplified product (0.7 μL) was genotyped by adding 0.3 μL of LIZ-500 size standard in ABI 310 genetic analyzer (Applied Biosystems, USA). Microsatellite alleles were scored using GENEMAPPER, version 4.1 (Applied Biosystems, USA). To finalize the consensus genotypes, a multi-tube approach was used where at least three identical homozygote PCR results were required for a homozygote genotype call; each allele was observed in two independent PCRs to record a heterozygous genotype [53]. By utilizing the matching tool in GenAlEx version 6.503 [54], individual tigers were identified from genotype data.
Genetic structure analysis. The 8 loci microsatellite genotypes from previously identified tiger individuals were used as the reference baseline data [38]. A Bayesian clustering approach was applied to determine the genetic structure in tiger population of Nepal using STRUC-TURE, version 2.3.4 [55,56]. The Bayesian clustering was implemented using 2,000,000 Markov chain Monte Carlo (MCMC) replications after initial 500,000 burn-in with repetitions of the analysis for 10 independent times per each assumed population (K = 1 to 6). This analysis was run using the admixture model with correlated allele frequencies and was tested under supervised (with LOCPRIOR) and unsupervised (without LOCPRIOR) learning algorithms. In the supervised method, sampling location information consisting of three geographical locations (CNP, BNP, SWR) were provided as a priori information, whereas it was not provided in unsupervised method. The rate of change of likelihood (delta K) value was estimated by Evanno method [57] to determine the optimal numbers of genetic clusters or populations present in our reference genotype data for both supervised and unsupervised analysis using STRUCTURE HARVESTER, version 0.6.94 [58]. The percentage of inferred ancestry (Qscores) for each sample from 10 independent runs were aggregated into an average Q-scores using CLUMPP, version 1.1.2 [59].
Geo-source assignment of unknown tiger samples. Confiscated tiger samples were assigned to potential source populations (reference baseline data obtained through NTGP) using four different approaches, including Bayesian clustering analysis with STRUCTURE, version 2.3.4 [55,56], likelihood and Bayesian-based assignment in Geneclass, version 2 [60], and multivariate discriminant analysis of principal components (DAPC) using the R package adegenet, version 2.1.1 [61]. First, STRUCTURE analysis using the same parameters as described above was performed by combining the genotype data from the forensic samples (unknown origin; PopInfo = 0) with the tiger samples from NTGP's reference database (known origin; PopInfo = 1). The assignment of unknown samples to sub-populations was determined based on Q-score membership values. Tigers were identified as residents of an area when their average Q value score was ! 70%. We classified individuals with Q scores < 70% as admixed and unassigned [62]. Next, we estimated the probability of assignment for each forensic tiger sample to different source populations with assignment threshold score of 0.05 using Paetkau's frequency based method [63] and with a missing allele frequency of 0.01 using Rannala and Mountain's Bayesian method [64] in Geneclass, version 2 [60]. We assigned a tiger sample to a population of origin when the probability of assignment was >0.99 and classified samples as unassigned when the probability was below this threshold. To complement these analyses, we also conducted a discriminant analysis of principal components (DAPC) using the R package adegenet, version 2.1.1 [61], to cluster genetically similar individuals based on their multi-locus genotype. To establish conservative criteria for forensic assignment, we required a sample to meet the assignment criteria defined above for three of the four methods. These four methods were chosen because they have been shown to perform well in accurately assigning individuals in multiple studies [29,65,66].

Genetic analyses
A total of 120 individual tigers (CNP = 64; BNP = 36; SWR = 20) were identified by matching tool implemented via GenAlEx version 6.503 [54]. We calculated probability of identification (P ID ) to be 1.5E-06 and probability of identity among siblings (P ID(sibs) ) to be 3.2E-03 (<0.01) for 8 loci and achieved an overall genotyping success of 40%. Of all the evaluated microsatellite markers, eight were selected for genetic analysis of all samples with fairly high overall PCR amplification success (84%) and genotyping accuracy (82%) (S3 Table). Two microsatellite markers (FCA205 and PttA2) were included in the PCR multiplex but were not considered in the individual identification as one (Locus FCA205) did not amplify in the majority of samples and the other (PttA2) was not polymorphic. Similarly, the mean allelic dropout rate was 2.46% and false alleles were around 15.85% [44].
All forensic samples (n = 15) were of tiger (S2 Fig). Sex identification revealed ten males and five females (S3 Fig). Genotyping success was 100% for all forensic samples (S4 Table), except for one of the blood-smeared knife samples (F-NP-0012) in which only half the number of loci successfully amplified. A multiple genotype (microsatellite loci) pattern was observed in F-NP-0012, where three different alleles were detected at three loci FCA391, FCA232 and FCA043 hence, this sample was excluded from further analyses. By using matching tool in GenAlEx, version 6.503 [54], the genotype of one of the forensic samples (F-NP-0004), determined to be originating from a female tiger [ Table 1]), matched 100% to previously profiled female tiger at the BNP site in the NTGP baseline tiger reference database.

Reference population genetic structure
Bayesian clustering analysis performed in STRUCTURE (with LOCPRIOR) predicted two as the most probable number of genetic clusters (K = 2) based on Evanno methodology, where one group consisted of samples from CNP, while the other constituted of samples from BNP and SWR (Fig 2 and S4 Fig). However, K = 3 was the most supported structure from mean likelihood (mean LnP(K) = -2097.47 and SD = 9.6), and Q value plots where the three clusters corresponded to tiger samples from CNP, BNP and SWR. The latter structure pattern with three genetic clusters (K = 3) was also corroborated based on Q-score plots across CNP, BNP and SWR, making an effective reference structure that would be essential for forensic assignment of unknown samples.
The unsupervised approach in STRUCTURE also predicted two as the most probable number of genetic clusters (K = 2) based on Evanno method dividing the data into one cluster consisting of CNP samples and another with samples originating from BNP and SWR (Fig 3 and  S5 Fig). However, this model seemed to support K = 4 from mean likelihood (mean LnP(K) = -1999.7 and SD = 0.34). The CNP population showed two sub-structures at K = 4, but no significant spatial patterns were found, rather the samples were panmictic across the CNP site. Therefore, we usedthe supervised model (with LOCPRIOR) with K = 3 subpopulation clusters for the population assignment of unknown forensic samples.

Forensic assignment of unknown samples
Using K = 3 as the reference population structure, Q membership scores of tiger samples of unknown origin were used to assign them to a geographical location. Eight forensic samples had Q-scores above 70% and were assigned to a specific reference population in Nepal, while the remaining six samples were determined to be of admixed genetic makeup (Table 1). Among the eight assigned samples, one sample (F-NP-0009) was assigned to CNP, five samples (F-NP-0004, -0006, -0010, -0011 and -0013) were assigned to BNP, and two samples (F-NP-0002 and -0014) were assigned to SWR. For the admixed individuals, 3 showed greater than 50% ancestry with SWR and most of the remaining ancestry with BNP while the other 3 showed greater than 50% ancestry with BNP and most of remaining ancestry with SWR (2) or CNP (1). Results from Geneclass2 analyses confirmed our findings in STRUCTURE. Using the frequency-based and Bayesian assignment methods, six samples were assigned to BNP, two were assigned to SWR and one was assigned to CNP with a probability >99% (Table 1).   As observed in the STRUCTURE analyses, five samples did not meet our assignment probability threshold (>99%) and showed admixed ancestry (Table 1). DAPC analysis suggested that one forensic sample originated from CNP while most others were assigned to the BNP cluster confirming the assignment results obtained using the other methods (Fig 4). Based on our criteria of requiring three or more assignment methods to confirm the results (see methods), we concluded that 6 forensic samples originate from BNP, 2 from SWR and 1 from CNP. 5 samples were classified as unassigned (Table 1).

Fig 4. Discriminant analysis of principal components (DAPC) to infer genetic assignment of confiscated tiger parts to source populations in Nepal.
Scatter plots of tiger genotypes (known and unknown origin) in relation to discriminant functions were generated using R package adegenet, version 2.1.1 [61]. Each point represents one tiger individual and ellipses around each cluster represent 95% confidence. The barplot graphs eigen values of the first three principal components in relative magnitude. BNP, Bardia National Park; CNP, Chitwan National Park; SWR, Shuklaphanta Wildlife Reserve; FS, forensic tiger samples. https://doi.org/10.1371/journal.pone.0201639.g004

Discussion
Wildlife crime is a serious trans-national threat to biodiversity in South Asia/Indian subcontinent and worldwide. There is an immediate need to develop tools that can assist in tracking and ultimately preventing wildlife crime which is driven by organized criminal syndicates that control the burgeoning and highly lucrative illicit trade in wildlife parts [67]. CIB has been actively involved in the investigation of wildlife crime across all regions of Nepal and has had major successes in bringing the perpetrators to justice by closely working with other national and international stakeholders. A high number of seizure cases of tiger body parts from Mid/ Far western Nepal compared to other regions indicates that this area might be one of the major trafficking routes for illegal wildlife trade in Nepal. The recently established molecular forensic capability in Nepal provides powerful, new tools in the fight against wildlife crime in Nepal.
The use of molecular forensic tools plays an important role in strengthening law enforcement efforts to address wildlife crime and biodiversity conservation. It has already been successfully used in some countries to improve investigation of wildlife cases, identification of poaching hotspots and trafficking routes for African elephants and Indian leopards, as well as prosecution of wildlife criminals by connecting individuals to crime scenes [20-24, 34, 67].
TAL is a biologically diverse habitat. Once alluvial grasslands and subtropical deciduous forests of TAL has supported the highest recorded density of tigers in the world [68]. Poaching now represents a great threat to tiger survival in this landscape [1]. Dwindling tiger numbers in some protected areas calls for urgent conservation action and requires changes in the current protection strategy for tigers. Tiger populations in the western part of TAL may already have reached a dangerously low tipping-point, and further decline in local population numbers may push the sub-population towards local extinction, similar to what happened to the rhino population across Babai River floodplain in BNP [18,69].
Application of molecular forensics is relatively new and rapidly evolving, and is gaining acceptance in wildlife conservation efforts [70,71]. The information obtained through the molecular forensics is more accurate and informative than conventional methods [21,72]. To identify source population of poached tigers, we created a wild tiger geo-spatial, genetic baseline database through the NTGP [38], which serves as a reference dataset for the forensic analysis. The molecular tracking of tiger parts has provided valuable information on wildlife crime "hot spots" within Nepal and presented scientific evidence to incriminate criminals in a court of law. We identified one tiger skin samples (F-NP-0004; female) that had an exact DNA genotyping profile with one of our female reference NTGP tigers from BNP. We collected scats of that individual in the winter of 2012; this particular individual was probably poached as recently as one to two years ago. Of all the forensic evidence we analyzed, most were traced to BNP (6/15) ( Table 1), highlighting this area as a wildlife crime hotspot. This is crucial information for Nepal's law enforcement officials. After dissemination of these results (Table 1), the concerned stakeholders, including the Nepal Police and the Nepal Army increased their patrols and intelligence activities particularly in the BNP region. As a result, a few of the criminal networks involved in wildlife crime in that area were discovered and necessary actions were taken to prevent further threat from these poachers [73,74]. Our analysis confirmed that all the seized parts were from tiger and majority of them (10 of 15, >65%) were male.
The molecular assay we have used in our baseline tiger genetic database and forensic samples have high probability of identification (P ID ; 1.5E-06) and probability of identification among siblings (P ID (sibs) ; 3.2E-03) with eight microsatellite loci showing high polymorphism [38]. Similarly, our species and sex identification PCR assay was highly effective [44].Other studies have shown that assignment tests can be very effective and accurate when pairwise F st values between source populations exceed 0.05, which was true in our study [44,65]. Using our analysis framework, we were able to assign 9 of the 14 individuals to a source population in Nepal. Five individuals were unassigned and showed admixed genetic profiles between CNP, SWR and BNP. These individuals may have been offspring of migrants. Previous research detected 3 migrants from SWR to BNP and showed higher gene flow from SWR into BNP than BNP to SWR [44]. This is also reflected in the STRUCTURE Q value plots that show a larger proportion of admixed individuals in BNP (Fig 2). Thus, these admixed individuals are more likely to have been poached in BNP than SWR [44]. Forensic samples may also originate from source populations outside of Nepal, which are not represented in our database. To test this hypothesis, we generated Nei's pairwise genetic distance values for all individuals and built a UPGMA phylogenetic tree to search for outlier samples (S1 File). Two of the unassigned samples (F-NP-0001 and -0008) in this analysis did not cluster with the 3 main groups and may represent individuals from unsampled source populations in India. Overall, the results of our forensic study which used four different analytical methods (STRUCTURE Q, GeneClass 2-Assignment and Bayesian likelihood and DAPC) ( Table 1) were broadly congruent, indicating that this analysis approach is a highly valuable tool against wildlife crimes in Nepal and beyond. However, our tiger genetic database currently only includes genotypes from tiger populations in Nepal. Hence, we acknowledge that some of the samples assigned to regions in Nepal could have been poached in tiger habitats just across the border in India where tiger allele frequencies are likely similar. There is now an urgent demand to have trans-boundary collaboration between tiger habitat countries [75,76] with the main goal to build and share standardized tiger genetic data to develop a comprehensive and cross-referencing anti-poaching genetic database. We also recommend a molecular forensic data validation process to evaluate accuracy of the technique and integrate this in courts of law as admissible evidence in wildlife crime cases by conducting multi-year scat surveys and building a robust baseline genetic database of wild tigers not only in Nepal but throughout the tiger range countries in South Asia.
There is an ongoing effort to have court of law in Nepal accept DNA-based evidence for incriminating wildlife criminals. We are working with all the relevant stakeholders, including the United States Fish and Wildlife Service (USFWS) and the International Criminal Police Organization (INTERPOL) to bring awareness on the utility of molecular forensics in Nepal and beyond.

Conclusions
The world has experienced an unprecedented spike in illegal wildlife trade, threatening to overturn decades of conservation gains. Wildlife trade is the fourth most lucrative illegal trade in the world after drugs, human trafficking and arms trade [77]. To combat wildlife crime through providing scientifically accurate and reliable information, Nepal has now generated a genetic database of wild tiger populations in Nepal, which includes geo-location data that allows molecular forensic tools to determine the source location of seized tiger parts. With efforts like NTGP, developing countries such as Nepal have made significant progress in building local capacity to gather important biodiversity information on various wildlife species. This forensic study provides critical information for effective law enforcement in the fight against wildlife crime.