Journal of Threatened Taxa | www.threatenedtaxa.org | 26 June 2022 | 14(6): 21127–21139

 

 

ISSN 0974-7907 (Online) | ISSN 0974-7893 (Print) 

https://doi.org/10.11609/jott.7572.14.6.21127-21139

#7572 | Received 15 July 2021 | Final received 22 March 2022 | Finally accepted 15 May 2022

 

 

Identification of confiscated pangolin for conservation purposes through molecular approach

 

Wirdateti 1, R. Taufiq P. Nugraha 2, Yulianto 3   & Gono Semiadi 4

 

1,4  Research Centre for Ecology and Ethnobiology, 2,3 Research Centre for Applied Zoology,

National Research and Innovation Agency, Jl. Raya Jakarta-Bogor Km. 46, Cibinong 16911, Indonesia.

1 teti_mzb@yahoo.com (corresponding author), 2 tragulus@gmail.com, 3 yulianto.mzb@gmail.com, 4 semiadi@gmail.com

 

 

 

Abstract: Over the past decade, the pangolin has emerged as one of the most prominent illegally traded mammals, and high extraction rates of Manis javanica from Indonesia have become a world concern. With the rise of the illegal trade, tools for uncovering the origins of pangolins for law enforcement are needed. Use of genetic markers for species and population identification has become a versatile tool in law enforcement efforts related to illegal wildlife trade and the management of endangered species. This study aims to uncover the origin of confiscated pangolins via a molecular approach using COI mtDNA markers. Forty-eight samples came from confiscated pangolins in Jakarta, Surabaya, Jember, Pangkalan Bun, Medan, Lampung, Riau, and Palembang, as well as four samples from the wild population in Riau, Pangkalan Bun, and East Java. Grouping using phylogenetic trees showed two groups with a bootstrap value of 90% based on wild samples. The first group consists of Sumatra and Kalimantan populations, while the second group consists of a Javan population. From a total of 44 confiscated samples, 12 were identified as Javan, nine from Kalimantan, and 23 from Sumatra. Genetic distance value (d) among individuals was d= 0.012 ± 0.002, with haplotype diversity (Hd) 0.864 ± 0.0444. The analysis of molecular variance (AMOVA) shows a clear genetic difference among populations (75%) and within populations (25%). The results showed that animals confiscated in one location may come from several different populations. These results can be used to track the flow of the pangolin trade in Indonesia, and support conservation management for the release of confiscated animals.

 

Keywords: COI, confiscated, illegal wildlife trade, Manis javanica, Pangolin, population.

 

 

 

 

Editor: G. Umapathy, CSIR-Centre for Cellular and Molecular Biology (CCMB), Hyderabad, India. Date of publication: 26 June 2022 (online & print)

 

Citation: Wirdateti, R.T.P. Nugraha, Yulianto & G. Semiadi (2022). Identification of confiscated pangolin for conservation purposes through molecular approach. Journal of Threatened Taxa 14(6): 21127–21139. https://doi.org/10.11609/jott.7572.14.6.21127-21139

 

Copyright: © Wirdateti et al. 2022. Creative Commons Attribution 4.0 International License.  JoTT allows unrestricted use, reproduction, and distribution of this article in any medium by providing adequate credit to the author(s) and the source of publication.

 

Funding: This study was funded by Southeast Asian Regional Centre for Tropical Biology (SEAMEO BIOTROP) research grant 2014, contract number: 060.12/PSRP/SPK-PNLT/III/2014.

 

Competing interests: The authors declare no competing interests.

 

Author details: Wirdateti is a researcher focusing on molecular genetics and ecology of wildlife species, she has working closely with multiple law enforcement agency for molecular identification of wildlife on illegal wildlife crime case. R. Taufiq P. Nugraha is a researcher working on wildlife physiology and conservation, he also has a great interest on wildlife management, illegal wildlife crime and zoonotic diseases. Yulianto is a laboratory technician specializing on molecular analysis and comparative wildlife reproduction.  Gono Semiadi is a researcher who has been working on wildlife management, zoos management and biopolicy.

 

Author contributions: W, RTPN, Y and GS have an equal contribution to this work.

 

Acknowledgements: Our gratitude goes to the Research Center for Biology-LIPI, Southeast Asian Regional Centre for Tropical Biology (SEAMEO BIOTROP), The Jakarta Natural Resources Conservation Center (BKSDA), The Pangkalan Bun-Central Kalimantan Natural Resources Conservation Center (BKSDA), The Jember-East Java Natural Resources Conservation Center (BKSDA), The Lampung-South Sumatra Natural Resources Conservation Center (BKSDA), The Zamrud National Park Siak Regency, Riau Province, and to Dewi Citra Murniati, Pramesa Narakusomo for assisting data analysis.

 

 

 

INTRODUCTION

 

There are eight extant pangolin species (Manis sp.) distributed in Asia and Africa. Four species are known In Asia: Manis pendactyla in China, M. crassicaudata in India, and two in southeastern Asia, the Sunda pangolin (M. javanica) also occurring in Indonesia apart from other southeastern countries, and M. culionensis in the Philippines (Feiler 1998; Gaubert & Antunes 2005; Gaubert et al. 2018; Kumar et al. 2018). In Indonesia, the Sunda Pangolin is one of several species listed as protected under the Minister of Environment and Forestry Regulation Number P.106 of 2018 concerning protected plant and animal. Under the International Union for Conservation of Nature (IUCN), this species is ‘Critically Endangered’, while CITES (Convention on International Trade in Endangered Species of Wild Fauna and Flora) list Sunda Pangolin in Appendix I since 2016. In Indonesia, pangolins can be found in Sumatra, Java, Kalimantan, and other surrounding islands. Over the past decade, pangolins have emerged as one of the world’s highest illegally traded mammal species surpassing other iconic species such as tigers, rhinos, and elephants (Kumar et al. 2018a).

The illegal trade in the eastern Asian and southeastern Asian markets was primarily driven by the demand for pangolin scales that were allegedly used by Traditional Chinese Medicine (TCM) and as accessories/ornaments, for spiritual and ritualistic purposes (Boakye et al. 2004; Challender 2011; Mahmood et al. 2012; Kumar et al. 2018 Xing et al. 2020). Scales of pangolins are the most valuable part, followed by meat (Li & Wang 1999; Pantel & Chin 2009; Challender 2011). The decline of the pangolin population in mainland Indo-China region due to excessive utilization caused traders to expand the range of pangolin search to all types in southeastern Asia, such as Indonesia, Malaysia, and India, as well as Africa. Factors responsible for pangolin population vulnerability are a low reproduction rate, predation, habitat loss, and poaching.

The level of poaching and illegal trade of pangolin in Indonesia is in stark contrast to the biological data, information, and studies on pangolins. Until now, the population, reproduction, most of the biological information of this species in nature are unknown. In contrast, the rapid decline in the population will undoubtedly continue every year, mainly due to hunting and habitat loss. Pangolins are particularly vulnerable to over exploitation because they are easy to hunt and have a slow reproductive rate (Yang et al. 2007; Challender 2011). Large-scale commercial harvesting and international trade have been going on since the early 20th century. Dammerman (1929) in Vincent (2015) reported the export of several tonnes of Sunda Pangolin scales from Indonesia on the island of Java to China in the period 1925–1929 involving at least 4,000–10,000 pangolins per year, even though the species is legally protected. Likewise, for the period 1958–1964, Harrisson & Loh (1965) in Vincent (2015) documented export licenses of more than 60,000 kg of pangolin scales which most likely came from Indonesian Borneo (Kalimantan) through Malaysia from Sarawak to Singapore and Hong Kong. Furthermore, data obtained from press and law enforcement authorities have shown that around 30,000 pangolins were caught in southeastern Asia between 2000 and 2007 (Chin & Pantel 2009 in Mahmood et al. 2012), indicating that M. javanica was mainly from Indonesia. The high hunting rate for Sunda Pangolins can be seen from the description of the results of confiscations from 2012 to 2015, where there were 45 confiscations: 12 in 2012, 10 in 2013, 17 in 2014, and seven in 2015. Sumatra is the location where most seizures occurred, with 21 confiscations totaling 4,046 individuals; Java had 14 confiscations with 6,736 individuals, and Kalimantan region seven with 793 pangolins destined for China (Vincent 2015). Data from tirto.id states that between 1999 and 2017, at least 192,567 pangolins were involved in illegal trade. Moreover, it estimated that the actual number is much higher due to many confiscation data not adequately recorded.

One of the main problems for law enforcement in the illegal wildlife trade of pangolins is the lack of information regarding the origin of confiscated pangolins (national and transboundary), since they can only be visually identified as Sunda Pangolin. This data is crucial for surveillance and conservation management to protect this species from extinction, e.g., choosing the right location to release confiscated animals. A DNA-based approach to species and population identification may prove to be a powerful tool for wildlife law enforcement agencies (Ogden et al. 2009; Zhang et al. 2015; Rajpoot et al. 2016). Genetic profiling of Indonesian Pangolin using mitochondria (mtDNA) reveals the genetic structure of the Sunda Pangolin population based on cytochrome b gene and control region (D-loop) (Kumar et al. 2018a; Wirdateti & Semiadi 2013, 2017). Nevertheless, this study is only conducted on a small part of the mtDNA gene, while mtDNA using a single marker is prone to bias (Ballard & Whitclok 2004). Recently, a whole-genome sequence of Sunda Pangolin originating from Malaysia provides a genome infrastructure for genetic research related to conservation and management (Cho et al. 2016), providing broader insight into genome conservation to reveal possible illegal trade routes and mixing of pangolin lineages in southeastern Asia (Nash et al. 2017). In the present study, we conducted identification of confiscated Sunda pangolins using COI genes to provide information for their management and conservation.

 

 

MATERIALS AND METHODS

 

Sample Collection

A total of 48 samples were taken from confiscated (44) and wild pangolin (4) in several places (Java, Sumatra, and Kalimantan; Table 1). DNA materials were collected as tissue from meats, and scales, and were preserved in absolute ethanol. Wild samples are taken from scales of dead pangolins found in their natural habitat. Confiscated samples were collected in 2008 from Medan, from Kalimantan in 2013 (Pangkalan Bun), from Java in 2010 and 2014 (Jember, Jakarta, and Surabaya), from Sumatra in 2014 and 2018 (Lampung, Riau, and Palembang; Figure 1).  Wild samples were acquired from Central Kalimantan, Riau, and East Java. Some of these samples (26 samples originating from 2013 and 2014 confiscation) had been analyzed using Cytochrome b (Wirdateti et al. 2013).

 

DNA Amplification

Total genome DNA was extracted using Qiagen Dneasy Blood and Tissue Kit Mini Stool (Qiagen) for tissue samples. For scale samples, and tissue with low yields, we extracted DNA using conventional phenol-chloroform (Kocher et al. 1989). This study used the COI gene mtDNA as a marker to determine the population origin of the confiscated pangolins by using a specific primer on Sunda Pangolin as long as 870 bp. The primer was designed as follows COI Treng F: TGGAAACTGACTAGTGCCCC; COI Treng R: GCTCCCATGGAGAGAACGTA. Primers were designed using a sequence template from COI Pangolin. The primers were designed using Primers3 (v.0.4.0) and Pick primers tools.

The amplification uses 30 µl polymerase chain reactions (PCR) containing 1 µl DNA template, 17 µl PCR mix reaction (FirstBase, Singapore), 2.5 µl primer F and R respectively, and distilled water (MQ) up to 30 µl. PCR reaction started with a 3-min denaturation at 95°C, followed by 35 cycles of denaturation at 94°C for 30 seconds, annealing at 56°C for 45 seconds and extension at 72°C for 30 seconds. The final incubation was at 72°C for 10 min.

 

Sequencing

PCR products were sequenced using the same forward and reverse primer as in amplification at FirstBase, Singapore using the Sanger method. PCR products were purified using the kit SureClean Plus (Bioline USA Inc.) according to the manufacturer’s manual and sequenced using BigDye Terminator v3.1 Cycle Sequencing Kit DNA Analyzer (Applied Biosystems) following Vendor’s protocol.

 

Data Analysis

All nucleotide sequence results were stored in a database using BioEdit software. The complement sequence between the forward primer and the reverse was edited with Chromas Pro software. All sequences were compared with the NCBI Genbank BLAST Database (www.ncbi.nlm.nih.gov/BLAST). DNA alignment was done using Clustal X (Thompson et al. 1997), and data analysis was conducted using MEGA 6.0 software (Tamura et al. 2013) and DNaSP ver. 5.0 (Librado & Rozas 2009). The MEGA 6.0 calculates the genetic distance and site variations among 48 samples, and the phylogenetic trees were used to determine each confiscated pangolins’ position based on the wild samples data. The analysis of DNA polymorphism includes the calculation of haplotype (h), haplotype diversity (Hd), and diversity of the nucleotides (π) using the DNaSP ver 5.0 software. Identification of Sunda pangolin was conducted using comparisons of Asian pangolin species, M. pendactyla (China), M. crassicaudata (India), and M. culionensis (Philippines) in GeneBank NCBI (NCBI Reference Sequence: NC_016008.1; NC_036434.1; and NC_036433.1, respectively). The phylogenetic tree formed was constructed using ML (Maximum Likehood) methods with bootstrap precision of 5,000.

For the selection of the best-fit model of nucleotide substitution using Bayesian inference (BAY) was conducted with the software IQ-TREE 1.6.12 (Nguyen L et al. 2015). The best-fitting of nucleotide substitution model for gene was determined with jModelTest v.2.1.6 (Kalyaanamoorth et al. 2017) selected by the Bayesian Information Criterion (BIC). The nucleotide frequencies for COI: A = 0.2509, C = 0.2915, G = 0.185, T = 0.2726; proportion of invariable sites I = 0.7542. The result was shown in FigTree v1.4.4 (Rambaut 2018). Bootstrap percentages (BP) were computed using 5,000 replicates.

The analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was conducted to investigate the hierarchical structure of mitochondrial marker variation to test the significance of the three pangolin populations using Arlequin v.3.5. 2.2 (Exocoffier & Lisher 2010). The significance of this structure was tested with 20,000 random permutations to test the significance of the three pangolin populations (Weir & Cockerham 1984).  AMOVA was performed by grouping samples according to their geographical location according to our result   from the previous analysis (MEGA). We calculated genetic differentiation among pangolin populations as pairwise fixation indices (Fst) in Arlequin. We used the pairwise FST values distances as the input data and 200 permutations were performed to determine the level of significance.

 

 

RESULTS

 

A.  Genetic Variations

The COI fragment from all samples was 866 bp long obtained using COI primer Treg F and COI Treg R designed from a sequence available on the GenBank NCBI. Only four samples had known origins: the wild samples obtained from Central Kalimantan (Pangkalan Bun), East Java (Jember), and Sumatra (Riau), while the other 44 samples came from the confiscated, market, and private collection with unknown origin. The use of wild samples is essential as a comparison to provide information of the unknown sample’s origin. Nucleotide blast in GeneBank NCBI revealed similarities (homology) sequence of 98.75 % to 99.75% for all samples with Sunda Pangolin (M. javanica). Furthermore, the genetic variation analysis of several parameters for the identification of confiscated samples can be seen in Tables 2 and 3.

The results of polymorphic sites based on variations in nucleotide sites (V), singleton base (difference of one base) (S), informative sites (P), and genetic distance (d) show differences from each population of Sumatra, Java, and Kalimantan (Table 2, Table 3). Overall sequence alignment along 866 nucleotides from 48 samples contained of 54 variation sites (V), 16 singleton variation sites (S), and 38 informative sites (P). Genetic distance between individuals d = 0.012 ± 0.002 (1.2% ± 0.2%), which is formed from base mutations or site variations in the 866 bp nucleotide sequence.  The use of Cyt b on M. javanica (Sunda Pangolin) showed higher variations of 83 site variations, 20 singletone variation sites (S), and 226 conserved sites with a sequence length of 331 bp nucleotide. (Kumar et al. 2018a).  Results of analysis of 20 sequences along 373 nt Cyt-b mtDNA showed 32 site variations, 21 sites informative, and 11 singleton sites from confiscated samples (Wirdateti et al. 2013). While the identification of confiscated pangolin species in Africa using the COI gene showed, the genetic distance (d) was from 0.001 to 0.055 (0.1% to 5%) among all species with M. javanica and P. tricuspis (Mwale et al. 2017).

DNA polymorphism analysis based on site variations showed 21 haplotypes (h) from the entire study sample with haplotype diversity Hd = 0.864 ± 0.0444. Nucleotide diversity (Pi) was Π = 0.01138 ± 0.00140, with the average nucleotide difference between individuals (k) = 9,801. This value gives the genetic distance between confiscated pangolin individuals of about 1.1%, indicating pangolins are in the same species but different populations.

To strengthen the quality of this study, we calculated the analysis data using AMOVA and Statistic test (Fst) to get strengthening the quality of the study, we calculated the analysis data using AMOVA and the statistic test (Fst). 

The results of AMOVA for total populations are shown in Table 3. The AMOVA for three populations shows a significant Fixation Index Fst value of 0.7525, indicating that at least  the pair-wise populations reveals significant heterogeneity (p < 0.05) (Table 4).

We found significant genetic differentiation (pairwise Fst) among population and within population base on localities pangolin samples. And statistic test with distance method are show values of pairwise Fst calculated between populations are genetically distinct (Table 5). The values of Fst between Java and Sumatra (0.81201); Java and Kalimantan (0.75713); Kalimantan and Sumatra (0.33619). Comparisons of pairs of these populations are significant (p = 0.000).

 

B. Phylogeny

The phylogeny tree was formed using the MEGA 6.0 program (Kumar et al. 2015) with the ML (Maximum Likelihood) method with a bootstrap of 5000, as shown in Figure 2. As a comparison, other species from Asia, M. pendactyla, M. culionensis, and M. crassicaudata from the NCBI GeneBank sequence (NCBI Reference Sequence: NC_016008.1; NC_036434.1; and NC_0364333.1) were used. Phylogenetic analysis was used to identify confiscated pangolins based on samples of pangolins from the wild (Figure 2). From 44 confiscated pangolins and four wild samples, two main groups with a bootstrap value of 90% were formed, representing the Sumatra-Kalimantan population and the Java population.

The phylogeny tree shows that the first group being represented by, Sumatra and Kalimantan, came from four populations: Population 1 representing the Sumatra, and populations 2, 3, and 4 representing the Kalimantan.  This grouping is based on the wild sample of each population. Wild sample from Zamrud National Park in Riau (MZBR 1423; 1424) representing Sumatra, and the Pangkalan Bun wild sample (MZBR 1163) representing the population of Kalimantan. This first group shows no clear differences between the populations of Sumatra and Kalimantan with a low genetic distance d = 0.004 ± 0.001 (Table 3), and low bootstrap value (30%). In contrast, the Javan population is clearly separated from the Sumatran and Kalimantan groups based on the wild samples from Jember (MZBR 1179), and it was showed high genetic distance from Sumatra and Java (d = 0.024 ± 0.005 with Kalimantan; d = 0.023 ± 0.004 with Sumatra) than the Sumatra and Kalimantan (Table 3). Based on this grouping, 44 confiscated samples were identified as 12 samples from the Java population, 23 samples from Sumatra, and nine samples estimated to be from Kalimantan. Each population variation site (polymorphic sites) can be seen in Table 2. The population of Java has a fairly high diversity compared to Sumatra and Kalimantan.  As many as 23 sites varied (different nucleotides) from 13 samples in the Java population; in the Sumatran population, from 25 samples, only 14 sites were varied, while in the Kalimantan population from 10 samples, only seven site variations were found. The result indicates that populations with high nucleotide differences (site variations) in the sequence range provides high haplotype diversity on the confiscated pangolins population being tested in this study. The analysis results show that the genetic diversity in the Java population is quite high, namely nine haplotypes from 13 individuals. In comparison, the Sumatran population has eight haplotypes from 25 individuals, and Kalimantan has four haplotypes from 10 individuals or 21 haplotypes were formed in this confiscated sample (Figure 3). Besides high haplotypes, higher genetic diversity in the Javan population was also indicated by a higher genetic distance (d = 0.006 ± 0.001) than the Sumatran and Kalimantan populations (d = 0.003 ± 0.001; d = 0.002 ± 0.001). However, the haplotype diversity of the confiscated samples was still quite high (Hd = 0.864 ± 0.0444). Haplotypes can identify the origin of the population from confiscated samples based on on-site variations and groupings in the phylogeny tree. Individuals who share the same haplotype can provide information on the origin of confiscated pangolins, trade routes, assist in controlling and monitoring hunting sites, and policymaking on conservation directions.

 

The phylogeny using Bayesian (BAY):

These three populations are clustered into two distinct groups; the first group includes the Java population; the second one includes the other two populations. The second group seemingly came from either Sumatra population or Kalimantan (Borneo) population.

We included the posterior probabilities obtained by BAY in the tree obtained by ML and supported by bootstrap values (Figure 4). The TN+F+I (Tamura Nei, parameter F: Nucleotide Frequencies; I: Invariance Sites) model nucleotide substitution was selected as the best-fit model according to BIC (Bayesian Information Criterion scores and weights) of evolution for all gene fragments using JModeltest (Kalyaanamoorthy et al. 2017).  The nucleotide frequencies for COI: A = 0.2509, C = 0.2915, G = 0.185, T = 0.2726; proportion of invariable sites I = 0.7542. This topology is almost similar to the NJ tree in the MEGA Program, where Kalimantan and Sumatra are in one group (Figure 5).

 

 

DISCUSSION

 

To see the position of the Sunda Pangolin species, other Asian pangolin species M. pendactyla, M. culionensis and M. crassicaudata were used as a comparison. The results showed that all samples used in this research belonged to the M. javanica species group. The results of the analysis showed that M. javanica was separated from M. culionensis (Philippines) by a genetic distance (d) of 3.6%, and from M. crassicaudata (India) by a genetic distance (d) of 14.4%, the separation between these two species had a high bootstrap value (93%). The genetic distance (d = 3.6%) between M. javanica and M. culionensis indicates that the two species are closely related. This result is also supported by the results of previous studies, which stated that the Palawan Pangolin M. culionensis (Philippines) is often considered a subspecies of M. javanica; the species was later raised to the species level based on morphological differences with M. javanica (Feiler 1998; dan Antunes 2005). Likewise, Gaudin et al. 2009 (Gaubert et al. 2018) stated that the thick-tailed pangolin, the Sunda and Palawan pangolins (M. javanica and M. culionensis) are sister species. The Chinese Pangolin species are located in population 4 or one clade with the Sunda Pangolin in the phylogenetic tree (Figure 2). Meanwhile, another sample in population 4 came from confiscated in Medan, Jakarta, Lampung, and Palembang. The other studies on both species using the same COI marker showed the separation between M. pentadactyla and M. javanica in the phylogenetic tree (dan Antunes 2015; Hassanin et al. 2015). So, the presence of M. pendactyla in subgroup 4 (Kalimantan), possibly indicates that the sample from NCBI is not M. pendactyla, or the confiscated sample is of unknown origin.

The confiscated pangolins are identified as Sunda Pangolin with a genetic distance of d = 0.012 ± 0.002 (1.2%), and nucleotide diversity Π = 0.01138 ± 0.00140, which indicated the value of differences within one species. However, the use of the mtDNA COI gene in this species has not shown a clear separation between the Sumatran and Kalimantan populations, as shown in the phylogenetic tree, while the Javan population is clearly genetically separated (Figure 2). Based on the group formed, the location of the confiscation does not always indicate the origin place of the pangolin.  It can be seen that several confiscated samples from central Kalimantan (Pangkalan Bun) were in the same group as Sumatra (1165, 1164, 1157, etc. Table 1.), while confiscated samples from Sumatra (Medan 270, Lampung 057) and Jakarta (1034, 1030) were in the same group as Kalimantan (wild).  The same result can also be seen in the one Kalimantan confiscated sample (1166) clustered in the Javan group. Previous research using mtDNA (mitochondria) levels also showed that some samples from Medan, Kalimantan were clustered with the Javan population, then, confiscated samples from Sumatra and Java were clustered in the Kalimantan population (Nash et al. 2017). The grouping of each individual also gave the same results as the haplotype phylogeny, which gave a clear difference in the Java population, with nine haplotypes from 13 individuals with a bootstrap value of 99% (Figure 3, Table 2). Zhang et al. (2015) revealed that the analysis of confiscated scales using multiple levels of mitochondrial DNA also gave an unclear separation in the population of M. javanica species. The results above can illustrate that the illegal trade of pangolin in Indonesia is run through several routes, namely Sumatra, Java, and Kalimantan.

AMOVA analysis with genetic structure testing showed significant genetic differentiation (pairwise Fst) among pangolin populations. Although the phylogenetic tree shows several genealogical branches or geographic clusters, the results of cluster analysis, sequence statistics, and AMOVA indicate a significant division between these three populations. The cluster analysis suggests that these three populations can be clustered into two groups, one includes Java populations, and the second population includes Kalimantan and Sumatra. Fst values between the Java population and Sumatra, between the Java and Kalimantan populations, and between the Sumatra and Kalimantan populations show significant genetic differences (Table 4), indicating that at least two populations exist of pangolins in Indonesia. The AMOVA results show that among the population percentage of variance is 75.25%, and within a population is 24.75%, and the Fixation index (Fst) value is 0.7525 which indicates a significance (p< 0.005) (Table 4). A considerable Fst value indicates a genetic structure with a high degree of variation between populations, and each population is geographically separated where the allele frequencies are different. While within the population shows a small diversity value in the genetic structure, the possibility of mating or breeding is high among the population due to the low effective population size (Ne). If Fst is small, it means that allele frequencies in each population are the same; if it is large, the allele frequencies are different (Hosinger & Weir 2019).

The sample size that is not large enough or irregular or small will affect the genetic structure. The FSt statistic test showed significant results both between populations and within populations. Based on the comparison of pairs of population sample test, the Java and Sumatra populations gave a higher value (Fst = 0.812, p <0.001) than Java and Kalimantan (Fst = 0.757, p <0.001), and the lowest values were Sumatra and Kalimantan (Fst = 0.336, p <0.001) (Table 5). The FSt values above indicate a robust genetic structure for the Javan population, with high differentiation with Sumatra and Kalimantan. The amount of genetic differentiation among populations has a predictable relationship with the rate of evolutionary processes (migration, mutation, and drift). Large populations with a lot of migration tend to show little differentiation, whereas small populations with little migration tend to be highly differentiated (Wright 1931). The results of other studies also showed that the intraspecific p-distance in M. javanica and P. tricupis was higher (COI: 0.037 to 0.030) than African pangolins, which averaged between 0.001–0.055. It has a higher maximum intraspecific divergence indicating a geographic sub-structure (Mwale et al. 2016).

Like the previous analysis, the results of the statistical distance test through Alerquin, showed that the populations of Sumatra and Kalimantan were closer and also shown in the BI phylogram tree. Bootstrapping does not support the separation of the two populations, namely node 82, following the Kalimantan population to be a sub-population, and this result is also shown from different haplotypes or no shared or nested haplotypes. In contrast, the genetic structure of the Javan population shows that the population is separated from the other two populations through Bayesian analysis with nodes 100 bootstrap, AMOVA and Fst statistical tests, with high distance values. The sample size that is not large enough also affects gene flow (Nm) in genetic structure (not shown), these results indicate a significant difference in the Java population (7,025; p <0.001) than in Sumatra (2,220, p <0.001) and Kalimantan (2,911. P <0.001). Test Differentiation Based on Haplotype Frequencies (Raymond & Rousset 1995) was significant between populations (p <0.05). The strong and significant genetic structure indicates substantial limitations on genetic and demographic connectivity (Hedgecock et al. 2007) among pangolin populations in Indonesia.

The Bayesian inference phylogenetic analysis results can be seen from the phylogram (Figure 4) using the IQ Tree program. The value at each branch point node is the result of the bootstrap support value in supporting topological credibility. Some results of bootstrap on several nodes/branch points have unsupported values with indistinguishable branches (polytomy) so that the position of external nodes or individuals may be incorrect. The Bayesian Inference (BI) phylogenetic results are not much different from the previous analysis, namely MEGA in terms of population divergence on valid bootstrap support (Hoang et al. 2017). The Java population still represents a separate group from Sumatra and Kalimantan with valid bootstrap support of 100 and the position of the Kalimantan population from Sumatra. However, the sample numbers MZBR 1417 and 1424 were separated from Sumatra and Kalimantan with a bootstrap of 100, while the Kalimantan population was separated from Sumatra by a bootstrap of 82. The bootstrap node value of 82 did not support the phylogenetic tree in BI. A phylogenetic tree has supporting nodes with a bootstrap value of 95 for Bayesian values ​​(Huelsenback & Hilis 1993). The branching or divergence of each individual in the population seems to show better resolution and description, although a very valid bootstrap value has not supported it for several nodes. Although it doesn’t produce a valid bootstrap support value, the topology with a better resolution may be due to the Effective population size (Ne), which is analyzed heuristically to minimize polytomy. The advantages of Bayesian Inference (BI) resolution over MEGA can be caused by complex parameters in BI, the use of the MCMC (Monte Carlo Markov Chain) numerical algorithm, and prior and posterior distributions.

The Java population represents a monophyletic group with the same common ancestor and lineages and forms a natural group with a valid bootstrap support value of 100. Although the AMOVA data clearly shows the population structure, this result cannot be clearly explained by the separation of the Sumatran and Kalimantan populations.

The results above show that mitochondrial COI markers in this study have not provided sensitive information for each population or intra-species. But a DNA-based approach to species and population identification may prove to be a powerful tool for wildlife law enforcement agencies (Ogden et al. 2009; Zhang et al. 2015; Rajpoot et al. 2016). Several experts have used mitochondrial markers as validation for species identification, including cytochrome b (Cyt b), 12S ribosomal RNA (12S rRNA), 16S ribosomal RNA (16S rRNA), and Cytochrome oxidase subunit I (COI) genes which are routinely used for species identification in wildlife forensics (DeSalle et al. 1993; Hsieh et al. 2001; Guha & Kashyap 2006; Alacs et al.2009; Kumar et al. 2016, 2018). Likewise, for the identification of confiscated pangolins, the use of several mitochondrial COI, cyt b genes, and D-loop can distinguish several confiscated species, namely P. tricuspis, P. tetradactyla, S. gigantea, S. temminckii, M. javanica, and M. pentadactyla with high bootstrap values ​​>70%, and the distance between all species was around 0.100–0.188 for COI and 0.10–0.20 for Cyt b, and 0.048–0.125 for the D-loop (Mwale et al. 2017). Thus, COI, Cyt b, and D-loop markers were more effectively used for identification or as inter-species markers.

Reports of high extraction rate of pangolins from Indonesia have become a concern to the world. However, counter measures and origin of these pangolins is not clearly understood. One of the main problems of pangolin confiscation in Indonesia is identifying the source and distribution of these confiscated pangolins; there is no data on the genetic distribution map of Sunda Pangolin in Indonesia. A distribution map will help the conservation of pangolin by allowing stakeholders to monitor the population and prevent its illegal trade. The latest report states that about 30% of the proportion of pangolin confiscated in Sumatra came from Kalimantan (Nash et al. 2017). Identification using one or two genes certainly cannot reveal the origin of the pangolin in the same species (intraspecies). This study is only conducted on a small fraction of mtDNA genes, mtDNA as a single marker, is prone to bias (Ballard & Whitclok 2004). With this argument, it is necessary to reveal the whole genome mtDNA and approximately 15,000 bp nucleotides as genetic markers for identification at the population level, especially for Indonesian pangolin. The data can be used as a baseline for mapping Indonesian pangolin genetic diversity to assist the conservation and handling of confiscated animals. The main problem with confiscated pangolin is that life confiscated animals will be released back into the wild as soon as possible; in many cases, these animals have been released back to the nearby confiscation area or region the pangolin while the pangolin itself might not come from the same population. This will undoubtedly affect each population’s gene pool, as the results of this study show that there are pretty clear differences between Sumatra, Kalimantan, and Java populations. In this regard, the information provided by this research is essential for policymakers and stakeholders to better understand the management and conservation of pangolins.

           

 

CONCLUSION

 

The use of COI gene markers in this study has not been able to provide effective information on confiscated samples based on population origin, especially between Sumatran and Kalimantan, owing to low genetic distances. However, it can provide a clear separation between Sumatran and Kalimantan populations with Java populations based on phylogenetic trees and a higher genetic distance values, and the Javan population had a stronger genetic structure than the other two populations. Based on the distribution of haplotypes from confiscated samples can identify the origin of confiscated pangolin from Java, Sumatra, and Kalimantan populations. Even though genetic distances and nucleotide differences between Sumatra and Kalimantan are very low, they can be distinguished from the haplotype type. This study’s findings showed that the seized material came from several organized hunting locations from illegal traders in the range of pangolin distribution areas, as shown that samples from one confiscation location originated from more than one population. Further analysis is required with the addition of wild samples with known geographical origin as a comparison reference. Policymakers can apply this information to release live pangolin, manage and supervise wild pangolins, and carry out effective law enforcement.

 

 

Table 1. List of samples used in this study.

 

Catalog number

Year

Type

Sample location

1

MZBR. T01 (1352)

2018

Confiscated

KSDA Lampung

2

MZBR. T02 (1355)

2018

Confiscated

KSDA Lampung

3

MZBR. T03 (1353)

2018

Confiscated

KSDA Lampung

4

MZBR. T04 (1354)

2018

Confiscated

KSDA Lampung

5

MZBR. T05 (1359)

2018

Confiscated

KSDA Lampung

6

MZBR. T06 (1360)

2018

Confiscated

KSDA Lampung

7

MZBR. T07 (1361)

2018

Confiscated

KSDA Lampung

8

MZBR. T08 (1363)

2018

Confiscated

KSDA Lampung

9

MZBR. T09 (1356)

2018

Confiscated

KSDA Lampung

10

MZBR. T10 (1367)

2018

Confiscated

KSDA Lampung

11

MZBR. T11 (1322)

2018

Confiscated

KSDA Lampung

12

MZBR. T12 (1340)

2018

Confiscated

KSDA Lampung

13

MZBR. T13 (1421)

2018

Confiscated

Palembang Market

14

MZBR. T14 (1334)

2018

Confiscated

KSDA Lampung

15

MZBR.1038

2012

Confiscated

KSDA Bogor 1

16

MZBR. T15 (1341)

2018

Confiscated

KSDA Lampung

17

MZBR. T16 (1418)

2018

Confiscated

Palembang Market

18

MZBR.17  (1420)

2018

Confiscated

Palembang Market

19

MZBR.18  (1416)

2018

Confiscated

Palembang Market

20

MZBR.19  (1422)

2018

Confiscated

Palembang Market

21

MZBR.1034

2012

Confiscated

KSDA Bogor 1

22

MZBR.20  (1423)

2018

Wild

Zamrud National Park

23

MZBR.1036

2012

Confiscated

KSDA Bogor 1

24

MZBR.21  (1424)

2018

Wild

Zamrud National Park

25

MZBR.22   (1417)

2018

Confiscated

Palembang Market

26

MZBR.1040

2012

Confiscated

KSDA Bogor 1

27

MZBR.1165

2013

Confiscated

Pangkalanbun, Central Kalimantan

28

MZBR.0270

2008

Confiscated

Sukabumi, West Java

29

MZBR.1180

2014

Confiscated

Jember, East Java

30

MZBR.0273

2008

Confiscated

Medan, North Sumatra

31

MZBR.1181

2014

Confiscated

Jember, East Java

32

MZBR.1030

2012

Confiscated

Tegal Alur, Jakarta

33

MZBR.0272

2008

Confiscated

Medan, North Sumatra

34

MZBR.1182

2014

Confiscated

Jember, East Java

35

MZBR.1183

2014

Confiscated

Jember, East Java

36

MZBR.1179

2014

Wild

Jember, East Java wild

37

MZBR.0276

2008

Confiscated

Medan, North Sumatra

38

MZBR.1166

2013

Confiscated

Pangkalanbun, Central Kalimantan

39

MZBR.1057

2012

Confiscated

Tanggamus, Lampung

40

MZBR.1069

2012

Confiscated

Surabaya, East Java

41

MZBR.1070

2012

Confiscated

Surabaya, East Java

42

MZBR.1071

2012

Confiscated

Surabaya, East Java

43

MZBR.1072

2012

Confiscated

Surabaya, East Java

44

MZBR.1157

2013

Confiscated

Pangkalanbun, Central Kalimantan

45

MZBR.1163  

2013

wild

Pangkalanbun, Central Kalimantan

46

MZBR.1164

2013

Confiscated

Pangkalanbun, Central Kalimantan

47

MZBR.0275

2008

Confiscated

Medan, North Sumatra

48

MZBR.1162

2013

Confiscated

Pangkalan Bun, Central Kalimantan

 

Table 2. Genetic variation using mtDNA COI markers in 48 pangolin samples.

Parameter

Total samples

(Confiscated, Wild)

Java

Sumatra

Kalimantan

 

n = 48

N = 13

n = 25

n = 10

Polymorphic sites

 

 

 

 

Variable (polymorphic) sites  

Singleton variable sites

Parsimony informative sites

Genetic distance (d)

  54

  16

  38

d = 0.012 ± 0.002

23

8

15

0.006± 0.001

14

9

5

0.002 ± 0.001

7

7

-

0.003 ± 0.001

DNA Polymorphism

 

 

 

 

Haplotype (h):

Haplotype diversity Hd

Variance of Hd

 

Nucleotide diversity, Pi:

 

Average of nucleotide differences, k

h = 21

Hd = 0.864 ± 0.0444

V  = 0.00195

 

Π = 0.01138 ± 0.00140

 

K = 9.801

9

 

 

 

0.00812±  +/-     0.00458

11.647268

8

 

 

 

0.00256± +/-     0.00163

2.869185

4

 

 

 

0.00336± +/-     0.00217

3.301497

 

 

Table 3. Genetic distance between and within population.

 

Population

 Between population (d ± SE)

Within population

( d ± SE )

     1               2                 3   

1.

Java

                    0.005        0.004

  0.006 ± 0.001

2.

Kalimantan

   0.024                         0.001

  0.003 ± 0.001

3.

Sumatra

   0.023       0.004

  0.002 ± 0.001

 

 

Table 4. Analysis of molecular variance (AMOVA) based on sequences of cytochrome c oxidase subunit I (COI) of pangolin populations from Java, Sumatra, and Kalimantan.

Source of variation

Sum of squares

d.f.

Variance components

Percentage of variation

Fixation index (FST)

Among populations

2

166.210

5.53432 Va

     75.25

0.75254*  (p <0.05)                     

Within populations

45

81.894

1.81986 Vb

24.75

 

Total

47

248.104

7.35418

 

 

 

 

Table 5. Pairwise Fst calculated among pairs of pangolin population from Java, Sumatra, and Kalimantan based on sequences of cytochrome c oxidase subunit I (COI).

 

1

 2

3

 1. Java

 0.00000

 

 

 2. Sumatra

0.81201

0.00000

 

 3. KALIMANTAN

0.75713

0.33619

0.00000

 

 

For figures - - click here for full PDF

 

 

REFERENCES

 

Alacs, E.A., A. Georges, N.N. FitzSimmons & J. Robertson (2009). DNA detective: a review of molecular approaches to wildlife forensics. Forensic Science, Medicine, and Pathology 6: 180–194. https://doi.org/10.1007/s12024-009-9131-7

Ballar, J.W.O. & M.C. Whitlock (2004). The incomplete natural history of mitochondria. Molecular Ecology 13(4): 729–744

Boakye, M., D. Pietersen, D.A. Kotze, D. Dalton & R. Jansen (2004). Ethno medicinal use of African pangolins by traditional medical practitioners in Sierra Leone. Journal of Ethnobiology and Ethnomedicine 10: 76.

Challender, D.W.S., S.R. Harrop & D.C. MacMillan (2015). Understanding markets to conserve trade-threatened species in CITES. Biological Conservation 187: 249–259. https://doi.org/10.1016/j.biocon.2015.04.015

Challender, D.W.S. (2011). Asian pangolins: increasing affluence driving hunting pressure. TRAFFIC Bulletin 23: 92–93.

Choo, S.W., M. Rayko, T.K. Tan, R. Hari, A. Komissarov, W.Y. Wee, A.A. Yurchenko, S. Kliver, G. Tamazian, A. Antunes, R.K. Wilson, WC. Warren, KP. Koepfli, P. Minx, K. Krasheninnikova, A. Kotze, DL. Dalton, E. Vermaak, I.C. Paterson, P. Dobrynin, F.T. Sitam, J.J. Rovie-Ryan, W.E. Johnson, A.M. Yusoff, S.J. Luo, K.V. Karuppannan, G. Fang, D. Zheng, M.B. Gerstein, L. Lipovich, S.J. O’Brien & G.J. Wong (2016). Pangolin genomes and the evolution of mammalian scales and immunity. Genome Research 26: 1312–1322.

DeSalle, R., A.L. Williams & M. George (1993). Isolation and characterization of animal mitochondrial DNA. Methods in Enzymology 224: 176–204.

Excoffier, L. & H.E. Lischer (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Molecular Ecology Resources 10: 564–567.

Feiler, A. (1998). Das Philippinen-Schuppentier, Manis culionensis ELERA, 1915, eine fast vergessene Art (Mammalia: Pholidota: Manidae). Zoologische Abhandlungen (Staatliches Museum fürTierkunde Dresden) 50: 161–164.

Gaubert, P. & A. Antunes (2015). What’s behind these scales? Comments to “The complete mitochondrial genome of Temminck’s ground pangolin (Smutsiatemminckii; Smuts, 1832) and phylogenetic position of the Pholidota (Weber, 1904). Gene 563: 06-108.

Gaubert, P., A. Antunes, H. Meng, L. Miao, S. Peigne, F. Justy, F Njiokou, S. Dufour, E. Danquah, J. Alahakoon, E. Verheyen, W.T. Stanley, S.J. O’Brien, W.E. Johnnson & S.J. Luo. (2018). The complete phylogeny of pangolin: Scaling up resources for the molecular tracing of the most trafficked mammals on earth. Journal of Heredity 109(4): 347–359

Gaubert, P. & A. Antunes  (2005). Assessing the taxonomic status of the Palawan Pangolin Manis culionensis (Pholidota) using discrete morphological characters. Journal of Mammalogy 86: 1068–1074

Gerstein, L., S.J. Lipovich, O’Brien & G.J. Wong (2016). Pangolin genomes and the evolution of mammalian scales and immunity. Genome Research 26: 1312–1322.

Guha, S. & V.K. Kashyap (2006). Molecular identification of lizard by RAPD and FINS of mitochondrial 16S rRNA gene. Legal Medicine 8: 5–10.

Hassanin, A., J.P. Hugot & B.J. van Vuuren (2015). Comparison of mitochondrial genome sequences of pangolins (Mammalia, Pholidota). Comptes Rendus Biologies 338: 260–265.

Hedgecock, D., S. Launey, A. Pudovkin, Y. Naciri & S. Lapegue. (2007). Small effective number 667 of parents (nb) inferred for a naturally spawned cohort of juvenile European flat oysters Ostrea 668 edulis. Marine Biology 150: 1173–1182

Hoang, D.T., C. Olga, V.h. Arndt, Q.M. Bui & V. Le Sy (2018). UFBoot2: Improving the Ultrafast Bootstrap Approximation. Molecular Biology and Evolution 35: 518–522.

Hsieh, H.M., H.L. Chiang, L.C. Tsai, N.E. Huang, A. Linacre & L.C. Lee (2001). Cytochrome b gene for species identification of the conservation animals. Forensic Science International 122: 7–18.

Huelsenbeck, J.P. &  D.M. Hillis (1993). Success of Phylogenetic methods in the four-taxon case. Systematic Biology 42(3): 247-264

Kalyaanamoorthy, S., B.Q. Minh, T.K.F Wong, V.H. Arndt & S.J. Lars (2017). ModelFinder: fast model selection for accurate phylogenetic estimates. Nature Methods 14: 587–589.

Kocher, T.D., W.K. Thomas, A. Meyer, S.V. Edwards, S. Pääbo, F.X. Villablanca &  A.C. Wilson (1989). Dynamics of mitochondrial DNA evolution in animals: Amplification and Sequencing with Conserved Primers. Proceedings of the National Academy of Sciences USA 86: 6196– 6200. https://doi.org/10.1073/pnas.86.16.6196

Kumar, V.P., A. Rajpoot, M. Thakur, M. Shukla, D. Kumar & S.P. Goyal (2016).   Illegal trade of Indian pangolin (Manis crassicaudata): genetic study from scales based on mitochondrial genes. Egyptian Journal of Forensic Sciences 6: 524–534.

Kumar, V.P., A. Rajpoot, M. Shukla, P. Nigam & S.P. Goyal (2018a). Inferring the molecular affinity of Indian pangolin with extant Manidae species based on mitochondrial genes: a wildlife forensic perspective. Mitochondrial DNA Part B: Resources 3(2): 640–644.

Kumar, V.P., A. Rajpoot, A. Shrivastav, K. Kumar, P. Nigam, A. Madhanraj & S.P. Goyal (2018). Phylogenetic relationship and molecular dating of Indian Pangolins (Manis crassicaudata) with other extant Pangolins species based on complete cytochrome b mitochondrial gene. Mitochondria Part A. 29:1–6.

Li, H. & H. Wang (1999). Wildlife trade in Yunnan Province, China, at the border with Vietnam. TRAFFIC Bulletin 18: 21–30.

Librado, P.J.R. & J. Rozas (2009). DnaSP v5: A Software for Comprehensive Analysis of DNA Polymorphism Data. Bioinformatics 25(11): 1451-2.

Mahmood, T., R. Hussain, N. Irshad, F. Akrim & M.S. Nadeem (2012). Illegal mass killing of Indian pangolin (Manis crassicaudata) on Potohar Region, Pakistan. Pakistan Journal of Zoology 44: 1457–1461.

Mwale, M., D.L. Dalton, R. Jansen, M.D. Bruyn, D. Pietersen, P.S. Mokgokong & A. Kotze (2017). Forensic application of DNA barcoding for identification of illegally traded African pangolin scales. Genome 60: 272–284.

Nash, H.C., Wirdateti, G.W. Low, S.W. Choo,  J.L. Chong, G. Semiadi, R. Hari, M.H. Sulaiman, S.T. Turvey, T.A. Evans & F.E. Rheind (2018). Conservation genomics reveals possible illegal trade routes and admixture across pangolin lineages in Southeast Asia. Conservation Genetics 19: 1083–1095

Nijman, V. (2015). Pangolin seizure data reported in the Indonesian media. TRAFFIC Bulletin 27(2): 43–46.

Nguyen, L.T., H.A. Schmidt, A. von Haeseler & B.Q. Minh (2015). IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Molecular Biology and Evolution 32: 268–274.

Ogden, R., N. Dawnay & R. McEwing (2009). Wildlife DNA forensics - bridging the gap between conservation genetics and law enforcement. Endangered Species Research 9: 179–195.

Pantel. S. & S.Y. Chin (2009) In: Proceedings of the workshop on trade and conservation of pangolins native to South and Southeast Asia. TRAFFIC Southeast Asia. Petaling Jaya, Selangor

Posada, D. & T.R. Buckley (2004). Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests. Systematic Biology 53: 793–808.

Raymond, M and F. Rousset (1995). An Test Exact Test for Population differentiation. Evolution 49(6): 1280–1283.

Rajpoot, A., K.P. Kumar, A. Bahuguna & D. Kumar (2016). Forensically informative nucleotide sequencing (FINS) for the first time authentication of Indian Varanus species: implication in wildlife forensics and conservation. Mitochondria DNA Part A. 27: 1–9.

Tamura, K., G. Stecher, D. Peterson, A. Filipski & S. Kumar (2013). MEGA 6: Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology Evolution 30(12): 2725–2729.

Thomson, J.D., T.J. Gibson, F. Plewniak, F. Jeanmougin & D.G. Higgins (1997). The Clustal_X Windows Interface: Flexible Strategies for Mutiple Sequence Alignment Aided by Quality Analysis. Nucleid Acids Research 25(24): 4876–4882.

Wei, B. & C. Cockerham (1984). Estimating F-statistics for the analysis of population structure. Evolution 38: 1358–1370. http://www.jstor.org/stable/2408641

Wirdateti, G.S. (2017). Genetic variation of confiscated pangolins of Sumatra, Java, and Kalimantan based on control region mitochondrial DNA. Jurnal Veteriner 18(2): 181–191

Wirdateti, G.S. & Yulianto (2013). IdentifikasiTrenggiling (Manis javanica) MenggunakanPenanda Cytochrome B Mitokondria DNA (Identification Of Pangolin (Manis javanica Desmarest, 1822) Using Cytochrome B mtDNA Marker). Jurnal Veteriner 14(4): 467–474.

Wright, S. (1931).  Evolution in Mendelian populations. Genetics 16: 97–159.

Xing, S., T.C. Bonebrake, W. Cheng, M. Zhang, G. Ades, D. Shaw & Y. Zhou (2020). Chapter 14 - Meat and medicine: historic and contemporary use in Asia, pp. 227–239. In: Challender, DW.S., H.C. Nash & C. Waterman (eds.), Pangolins. Academic Press.

Yang, C.W., S. Chen, C. Chang, M.F. Lin, E. Block, R. Lorentsen, J.S.C. Chin & E.S. Dierenfeld (2007). History and husbandry of pangolins in captivity. Zoo Biology 26: 223–230.

Zhang, H., M.P. Miller, F. Yang, H.K. Chan, P. Gaubert, G. Ades & G.A. Fischer (2015). Molecular tracing of confiscated pangolin scales for conservation and illegal trade monitoring in Southeast Asia. Global Ecology and Conservation 4: 414–422.