Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Phylogeography of the reticulated python (Malayopython reticulatus ssp.): Conservation implications for the worlds’ most traded snake species

  • Gillian Murray-Dickson ,

    Roles Data curation, Formal analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing

    gmdickson@rzss.org.uk

    Affiliation Royal Zoological Society of Scotland (RZSS) WildGenes Laboratory, Edinburgh, United Kingdom

  • Muhammad Ghazali,

    Roles Methodology, Writing – review & editing

    Affiliation Royal Zoological Society of Scotland (RZSS) WildGenes Laboratory, Edinburgh, United Kingdom

  • Rob Ogden,

    Roles Conceptualization, Writing – review & editing

    Affiliation Trace Wildlife Forensics Network, Edinburgh, United Kingdom

  • Rafe Brown,

    Roles Writing – review & editing

    Affiliation KU Biodiversity Institute, 1345 Jayhawk Blvd, Dyche, Lawrence, KS, United States of America

  • Mark Auliya

    Roles Conceptualization, Funding acquisition, Project administration, Writing – original draft, Writing – review & editing

    Affiliation Helmholtz Centre for Environmental Research – UFZ, Leipzig, Germany

Abstract

As an important economic natural resource in Southeast Asia, reticulated pythons (Malayopython reticulatus ssp.) are primarily harvested from the wild for their skins—which are prized in the luxury leather goods industry. Trade dynamics of this CITES Appendix II listed species are complex and management approaches on the country or regional level appear obscure. Little is known about the actual geographic point-of-harvest of snakes, how genetic diversity is partitioned across the species range, how current harvest levels may affect the genetic viability of populations, and whether genetic structure could (or should) be accounted for when managing harvest quotas. As an initial survey, we use mitochondrial sequence data to define the broad-scale geographic structure of genetic diversity across a significant portion of the reticulated python’s native range. Preliminary results reveal: (1) prominent phylogenetic structure across populations east and west of Huxley’s modification of Wallace’s line. Thirty-four haplotypes were apportioned across two geographically distinct groups, estimated to be moderately (5.2%); (2) Philippine, Bornean and Sulawesian populations appear to cluster distinctly; (3) individuals from Ambon Island suggest recent human introduction. Malayopython reticulatus is currently managed as a single taxonomic unit across Southeast Asia yet these initial results may justify special management considerations of the Philippine populations as a phylogenetically distinct unit, that warrants further examination. In Indonesia, genetic structure does not conform tightly to political boundaries and therefore we advocate the precautionary designation and use of Evolutionary Significant Units within Malayopython reticulatus, to inform and guide regional adaptive management plans.

Introduction

Habitat loss and degradation as a result of unsustainable per capita consumption of natural resources and rising human population levels [1] are having detrimental impacts on ecosystems and biodiversity, particularly in developing countries [24]. Overexploitations of natural resources and agricultural activity have been identified as some of the most prominent threats to biodiversity [5]. Unsustainable exploitation of wildlife [6], particularly reptiles [7,8], has been increasingly reported in Southeast Asia [913].

Human and pythons have a lengthy historical (and potentially evolutionary) association [14]. The international and commercial trade in python skins can be traced back to between 1910 [15] and the 1920s [1618] with contemporary uses for python-derived products spanning the fashion, food and traditional medicine industries [19,20]. Trade in reptile skin and leather products is valued at $339 million, approximately 5% of the legal, global wildlife trade [21], with five Southeast Asian python species (Malayopython reticulatus ssp., Python bivittatus ssp., P. curtus, P. brongersmai and P. breitensteini) being heavily exploited for this purpose. Among these the reticulated python (M. reticulatus ssp.) is the most economically important species [18] with approximately 350,000 skins legally exported annually for the high-end fashion market alone [20].

Formerly recognized in the genus Python, the reticulated python was genetically and morphologically allocated as a distinct clade, together with the Lesser Sunda python (Python timoriensis), and thereinafter included in the genus Broghammerus [22,23]. However, this genus [24] is considered invalid as it lacked accompanying data and analysis [25]. Reynolds et al. [26] therefore ascribed the new genus Malayopython, an action that has since been supported [27,28] and includes the two species M. reticulatus and M. timoriensis.

The present global range of Malayopython reticulatus is explained by the combination of the complex, regional geological history (particularly that of insular Southeast Asia) [29], the species’ excellent dispersal ability [30], and human introductions [31]. The species is distributed extensively across continental and insular Southeast Asia (Fig 1) with distinct insular ‘morphs’ (colour patterns) from Indonesia recognised in the pet and skin trade [32,33]. Although taxonomic subdivisions have not been extensively applied across the broader species range, morphologically and genetically differentiated populations (M. r. saputrai and M. r. jampeanus) have been identified on Selayar and Tanahjampea Island, respectively [34]. Phenotypic distinctiveness of the Philippines population has also been noted (Auliya and Brown, pers. obs.).

thumbnail
Fig 1. Sample geographic origins.

Distribution of the reticulated python across Southeast Asia, as indicated by the dotted line (after [3]). Pie chart size indicates the number of samples sequenced from each location and segments represent the frequency of each mitochondrial haplotype resolved.

https://doi.org/10.1371/journal.pone.0182049.g001

All python species are listed on Appendix II of the Convention on International Trade of Endangered Species of Fauna and Flora (CITES) [35], except for the Indian python (Python molurus) which is listed on CITES Appendix I. Appendix II listing allows commercial trade that is regulated through a permitting process [36]. Data and trends in the trade of Malayopython reticulatus skins, compiled during a CITES commissioned study to examine commercial trade and review the current taxonomic status, distribution, biology and ecology in Asian pythons’, identified Indonesia as the primary country supplying skins of M. reticulatus to the global market [18]. With export permits documenting approximately 700,000 skins in 1987/1988, the harvest quota of M. reticulatus in Indonesia was limited to 445,000 specimens in 1989. However, the number of skins documented on export permits from Indonesia totalled 555,882 for the same year and therefore traceability issues (origin of skins and trade routes) were already being called to interest [18]. Trade of wild pythons was banned in Thailand in 1992 [Wild Animal Reservation and Protection Act, B.E. 2535 (A.D. 1992); www.ThaiLaws.com], the Philippine governments banned the export of M. reticulatus in 1986, and Viet Nam banned wild harvest of M. reticulatus in 1998 [37]. Today, Malaysia and Indonesia are the major countries of origin and export for skins of wild M. reticulatus [20] while the only country that claims to regularly export captive-bred specimens for the skin trade is Viet Nam. During the period 2000 and 2013, the country exported 776,916 skins and 380,870 metres of skins of M. reticulatus that were declared as captive-bred (CITES Trade Database—http://trade.cites.org).

Although harvest levels vary across regions, significant declines of the species have been reported in Bangladesh, Indonesia, Lao PDR, Malaysia, Myanmar and the Philippines–all of which has brought in to question the sustainability and legality of intense, large-scale commercial trade in M. reticulatus skins [18]. Whilst traders in Indonesia may not perceive a depletion of local populations, this may be explained by the increase in harvest areas over years, to fetch the annually allocated quotas [3]. Hunting areas have increased in West Kalimantan (province of Indonesian Borneo) as a result of a growing trade system and the improved mobility of hunters, but may also be caused by a local decline in populations due to the expansion of agricultural land [3]. Sustainability of annual harvest levels are also debatable, due to the pressures of continuous forest loss in that geographic realm [38].

When combined with the high rate of population growth in these areas (with the exception of Thailand, it is estimated at 1.17% above the world mean average in 2010) [39] and associated use of natural resources, it is unknown how current levels of harvesting will affect the viability of M. reticulatus populations throughout the many political and geographic regions of its distribution. Furthermore, little is known about how levels of genetic diversity are partitioned across the species range, nor how variable harvest rates will affect local population viability. The question of python population sustainability has prompted industrial stakeholders to identify a need for traceability systems for proving the legal origin of python skins in trade, ideally to the point of harvest [40,41].

Examples of DNA-based tools to inform monitoring and enforcement are widespread. They have been used effectively during wildlife forensic investigations (see [42] for review)—to identify parts of animal derivatives in an otherwise unidentifiable sample. Forensic procedures have been developed to identify animal parts or derivatives in the traditional Chinese medicine trade (e.g., snakes: [43,44]), identify turtle species from samples of shell [45], to verify origin [46] either where illegal activity or fraudulent claims are suspected e.g., Testudo graeca [47] and Morelia boeleni [48], and to authenticate population of origin where regional differences in population viability, management or harvest quotas exist (e.g., ivory, [49]; shark fins, [50]; regulated fisheries, [51]). Taxonomic clarification prior to monitoring trade was suggested for short-tailed pythons, the polytypic species Python curtus [52]. High levels of divergence were found between the three, recognised subspecies, all of which were being exploited as a single taxonomic unit across Borneo and Sumatra. The authors suggested elevating each taxon to species level, with recommendations that they no longer be managed as a single biological taxon. As with Python curtus ssp., consideration of genetic structure in establishing management plans has been recommended for populations of the intensely traded Nile monitor lizard (Varanus niloticus) [53,54], to avoid unnecessary erosion of individual genetic partitions [55]. Despite recognising Malayopython reticulatus as an important commercial resource across much of its native range, little consideration has been given to examining the distribution of genetic diversity, much less to incorporating the data into regional management plans.

The legality of wild-harvest to meet the demands of the python skin trade, and the existence of registered python farms, varies among countries [37,56,57]. Where wild-harvest is not permitted (i.e. Cambodia, Thailand, Viet Nam) but laundering is suspected, cross-border trade activities and confiscations have been reported [40,5862], suggesting that a portion of ‘captive-bred’ specimens may be of wild origin. In such instances, it would be beneficial to be able to identify either the farmed or wild source of a sample. However, to use genetic data as evidence for enforcement would require farmed individuals to be genetically distinguishable from their wild counterparts, or the establishment of a verifiable breeding registration scheme. As official python farms only exist across a small portion of the species range and farming practices vary considerably across the industry [56], it is unlikely that the level of genetic differentiation would be uniform for all farmed populations. Development of a DNA-based traceability tool would therefore require scrutiny on a farm-by-farm basis to assess its feasibility in each case. Alternatively, if wild populations are genetically structured across distinct geographic areas, genetic data can be used to trace the geographic origin of an individual and from this, determine whether its geographic origin was the same as reported on the associated permit [11,12]. Where wild-harvest is permitted, extensive testing of python skin shipments would facilitate assessment of whether quota levels were being adequately complied with, monitored and adapted. This information could also be used to help identify potential situations where wild caught (e.g. live reptile) individuals are fraudulently labelled as captive bred [7,63].

To enable development of DNA-based traceability tools it is first necessary to understand whether population genetic structure exists across the species range, and assess whether the genetic structure reflects geographic boundaries or partitions, whilst quantifying differences that exist across known localities. This approach provides reference data with which to compare future samples of ‘unknown origin’, and an assessment of the distribution of genetic diversity will help direct future taxonomic investigation across different regions [55,64,65]. This preliminary approach may also be useful for defining and devising an adaptive management strategy for Malayopython reticulatus. It should be noted however that access to comprehensive reference sample sets can be challenging. When sub-optimal, comparative observations such as morphological comparison (see [34]) are not possible and results remain provisional. To this end we present mitochondrial sequence data derived from a collection of historical and contemporary samples of M. reticulatus to investigate broad-scale geographic structure of genetic diversity across a significant portion of the species native range. Provisional genetic structure is examined with respect to biogeographical variability across the Indo-Australian archipelago, and the feasibility of using these data for identifying the geographic origin of individual samples is assessed. We identify distinct genetic lineages as candidate targets for conservation (such as Evolutionary Significant Units [66]), to facilitate and guide further fine-scale analysis to improve the resolution of potential traceability tools for conservation management and enforcement of laws relevant to illegal trade.

Materials and methods

Handling of live animals for samples collected in the Philippines was approved under the University of Kansas’ IACUC authorisation 158–04.

Sampling and DNA extraction

Eighty-nine Malayopython reticulatus samples were used in this study, representing populations from Thailand and Viet Nam, West Malaysia, Singapore, Sumatra, Java, Lesser Sundas, Borneo, Ambon, Halmahera and the Philippines (see Fig 1). Samples of tissue were collated from museum voucher specimens that had been collected across the Southeast Asian distribution of M. reticulatus between 1862 and 2014 (see S1 Table for details). Blood samples were also provided from the Singapore Zoological Gardens, collected from individual live pythons encountered in and around residential areas. Ventral scale clip samples and shed skins were taken from captive held specimens of known geographic origin. It should be noted that one sample obtained from a museum collection was labelled as ‘New Guinea’, however, the species is not known to naturally occur in New Guinea and so this locality is considered in error. Tissue samples preserved either in ethanol, or in formalin prior to transfer to ethanol, were also frozen at -20°C. Blood samples were stored immediately upon collection in Ethylenediaminetetraacetic acid (EDTA) tubes and subsequently stored at -20°C. Shed skin was dried and then stored at -20°C. Total genomic DNA was extracted from individual blood samples, dried shed skin and museum samples preserved in ethanol using the DNeasy blood and tissue kit (Qiagen Ltd) according to the manufacturer’s protocol. Samples previously in contact with formalin were first soaked overnight in phosphate-buffered saline (PBS) prior to DNA extraction using the QIAamp DNA investigator kit (Qiagen Ltd), including 1 ng carrier RNA per extraction as per the manufacturers protocol.

Mitochondrial DNA sequencing and analysis

A suite of primers was used to PCR amplify approximately 1200 base pairs (bps) of the mitochondrial cytochrome b (Cyt b) and adjacent, partial control region. Although regularly used for phylogeographic inference, a duplicated control region is present in the snake mitochondrial genome due to gene rearrangements [67]. This makes it difficult to validate amplified fragments and so primers which target the cytochrome b gene were selected instead [67,68]. Details of primers are given in Table 1. Fragment 1 and Fragment 2 were amplified from fresh tissue, blood and dried skin samples, whereas Fragment 3, Fragment 4 and Fragment 5 were amplified from museum-derived samples. Primers sets Mretic1a, Mretic2 and Mretic3 were designed to produce smaller, overlapping amplicons that could be successfully amplified from fragmented DNA. The PCR amplifications were performed in a total volume of 20 μl using a MJ Research PTC-100 thermal cycler (Waterton, MA). The final reaction mix contained 2 μL (at 10–50 ng/μL) of template DNA, 14 μL Maxima Hot Start PCR mastermix (which includes Maxima Hot Start Taq DNA polymerase, 2X hot start PCR buffer, 0.4 mM of each dNTP, 4 mM Mg2+) and 2ul of each forward and reverse primer at 10uM. The PCR profile consisted of an initial denaturation step of 5 min at 95°C, 40 cycles of PCR consisting of 30 seconds denaturation at 95°C, 30 seconds of annealing at 50°C and 60 seconds extension at 72°C, ending with 10 mins at 72°C. Products were cleaned using 1 μl of a 1:1 Exo1/FastAP solution (ThermoFisher Scientific) and incubated for 15 mins at 37°C and 15 mins at 85°C. Two microliters of the purified product were prepared for sequencing with BigDye® Terminator v3.1 Cycle Sequencing Kits, as per manufacturers instruction (ThermoFisher Scientific), and sequenced on an ABI3730 capillary sequencer (Edinburgh Genomics GenePool facility, Edinburgh).

Fragments 1 and 4 were sequenced in both the forward and reverse directions using PCR primers, whereas Fragments 2, 3 and 5 were sequenced in the forward direction only. Samples exposed to formalin were sequenced in both directions for all fragments to check for consistency. Sequence chromatograms were checked by eye to ensure unambiguous base-pair determination, and consensus sequences were edited using the trace editor in MEGA version 6 [69]. Fragment sequences were concatenated for each individual and subsequently aligned using CLUSTAL W [70]. The optimal nucleotide substitution model for aligned sequences was identified as TrN+I (proportion of invariable sites = 0.834, gamma = estimated from the data) using BIC as implemented in jModeltest v2.1.10 [71,72].

Hierarchical phylogenetic connectivity between individuals was inferred using both maximum likelihood (ML) and Bayesian approaches. Maximum likelihood inference was parameterised by the defined model and branch support was assessed over 1000 bootstrap replicates using PhyML version 3 [73], implemented within Geneious v 8.0.4 [74]. Bayesian inference was conducted using MrBayes [75] as implemented in Geneious, and analyses were conducted using HKY+I as the closest available substitution model. The model was run for 1,500,000 MCMC iterations after discarding 500,000 burn-in generation, and subsampled every 100 trees across 4 heated chains. The posterior distribution contained a total of 3601 samples, which were summarized by consensus. The resultant phylogenies were rooted with outgroups Python bivittatus (NCBI accession number JX401133.1) and Malayopython timoriensis (NCBI accession number EF545106.1). The probability of reciprocal monophyly under the null model of random coalescence [76] was calculated using a ‘Species Delimitation’ plugin, available for implementation in the Geneious software. Haplotype genealogy and geographic distribution were examined using a median joining network, constructed using PopArt software [77], under default parameters. The networks were constructed using both the full sequence data set and for individual amplified fragments (e.g. Fragment 2 only) to check for congruence between the clusters resolved.

Haplotype and nucleotide diversities (including the number of haplotypes, Nei’s haplotype (gene) diversity—Hd, number of segregating sites, and nucleotide diversity—π) were calculated for each putative island population and each haplotype group using DnaSP ver. 5.10.01 [78]. Haplotype richness was calculated using Contrib software [79] for individual populations only. A hierarchical analysis of molecular variance (AMOVA) was used to estimate genetic structure within and between geographically distinct groups, with deviation from null distribution tested with 1000 permutations. The pairwise fixation index Φ, (analogous to Wrights Fst) under the Tamura-Nei substitution model was calculated in Arlequin v 3.5 [80]. Population genetic analyses were completed using only sampling locations represented by > 6 individuals whereas diversity estimates were made for all island populations. Haplotype sequences were also compared to those available for M. reticulatus on the NCBI database and permitted comparison across the entire amplified sequence length (NCBI accession numbers U69860.1 and U69859.1).

Results

Analysis of 830 bps of concatenated Cyt b and partial control region sequence revealed 34 unique haplotypes throughout 81 samples of Malaopython reticulatus. Excluding sites with missing data, haplotypes were defined by 62 variable positions identified along the remaining 570 base pairs, of which 49 were parsimony informative (Table 2). Sequences are deposited on Genbank (accession numbers: MF576180 –MF576213).

thumbnail
Table 2. Variable base pair positions across the cytochrome b and control region fragment.

https://doi.org/10.1371/journal.pone.0182049.t002

Both Bayesian and ML inference produced similarly unresolved tree topologies (Fig 2) with strong ML bootstrap support and posterior Bayesian support for five localised haplotype clades, but no support for branching order among these clades. These haplotype groups included moderate to strongly supported clades of haplotypes from (1) Sulawesi, (2) Singapore, (3) Java/Ambon/Sumatra/Singapore (4) Borneo/Sumatra/Lesser Sundas, and (5) the Philippines. Although there was moderate to strong support for the monophyly of M. reticulatus and five well-supported but unresolved clades, 12 unique haplotypes exhibited no statistically supported affinities to each other or the five aforementioned haplotype clades. The addition of more distantly related outgroup species (P. bivitattus and P. regius) did not resolve the branching order among clades. Except for a single haplotype (Haplotype 10), haplotypes in the Philippines and Halmahera formed a distinct monophyletic clade which fell into a basal polytomy with the remaining haplotypes. Haplotypes from Borneo (haplotypes 2, 7, 16 and 20) also clustered independently although one sequence (Haplotype 2) was shared with Java, Sumatra and the Lesser Sundas. Two haplotypes found only on Sulawesi (Haplotype 18 and 34) demonstrate a high level of sequence identity to each other (99.3%) and a lower sequence identity to the remaining haplotypes (95.5% and 94.6%, respectively), including Haplotype 14 which was also found on Sulawesi. It is possible these divergent haplotypes represent samples of the subspecies M. reticulatus saputrai, a distinct population of only known from Selayar Island and South Sulawesi [34]. Pairwise estimates of evolutionary divergence found within either the Philippines (0.7%–1.2%) or Singapore/Sumatra/Java/Borneo (1.2%–1.6%) were significantly lower than comparisons made between these two groups (4.4%–6.2%). The exception was Haplotype 10 which was more divergent from haplotypes in the Philippines (4.8%–5.4%). Tamura-Nei estimates of divergence between groups indicate a higher level between the Philippines and Sulawesi (5.5%), of the Philippines and Singapore (5.2%) than between Singapore and Sulawesi (4.4%).

thumbnail
Fig 2. Phylogenetic analysis.

Majority rule consensus tree for Bayesian inference of phylogeographic relationships among mitochondrial haplotypes inferred across the range of M. reticulatus (the bootstrap ML consensus tree resolved the same topology with respect to all well supported nodes). The tree was rooted with two more distantly related members of the genus Python and one more closely related member of Malayopython [22] and nodal support is provided as Bayesian posterior probabilities. The number of individuals per haplotype is given in parenthesis and coloured shapes indicate sampling location.

https://doi.org/10.1371/journal.pone.0182049.g002

The median joining network identifies two geographically distinct groups of haplotypes (Fig 3). Fifteen percent of haplotypes are shared across multiple island locations (the Philippines being counted as a single location); most of the remaining haplotypes were found at only one island location (Fig 1). A total of 20 haplotypes were present across Thailand, Viet Nam, West Malaysia, Singapore, Sumatra, Java, Lesser Sundas, Borneo, Sulawesi and Ambon (in this study termed the ‘Western’ haplotype group). Bornean haplotypes clustered hierarchically within the Western haplotype group. Only one haplotype was shared with other locations and the most frequent haplotype (Haplotype 16) was only found on Borneo, where it represents 69% of individuals. Although not forming a distinct cluster in the network (they did form a strongly supported phylogenetic couplet in the phylogeographic analysis), haplotypes found on Sulawesi exhibit higher levels of differentiation than other Western haplotypes. Twelve haplotypes were present among individuals from the Philippines, Halmahera and a single sample erroneously labelled as ‘New Guinea’ (termed the Eastern haplotype group). Despite its geographic proximity to Ambon and Sulawesi, samples from Halmahera were more similar to those from the Philippines, with one haplotype shared between both locations (Haplotype19). Haplotype16 occurred at the highest frequency (n = 9) within the Eastern group and was only present across the Philippines. Within both Western and Eastern haplotype groups, haplotypes typically differed by 1–5 substitutions from their nearest neighbour. The mean within-group distance ranged from 0.6% (Eastern group) to 1.0% (Western group), with 5.5% divergence between the two. The probability of the clades being chosen via a random coalescent process was p < 0.05, suggesting that the division may represent natural, geographically endemic genetic variation, possible warranting formal taxonomic recognition. Median joining networks compiled using only the short sequence fragments (i.e. before concatenation) retained the two main haplotype groups (Western and Eastern) but resolution within each cluster was largely lost as the number of base pairs analysed was reduced.

thumbnail
Fig 3. Median-joining network.

Network of mitochondrial haplotypes. Each circle represents a unique DNA sequence and its frequency and geographic identity are denoted by the circle size and colour. Hatch marks represent 1 nucleotide substitution along the 570 bp sequence and dotted lines delineate the phylogenetically-defined groups (see text for discussion).

https://doi.org/10.1371/journal.pone.0182049.g003

When haplotypes were grouped as either ‘Western’ or ‘Eastern’, between-group differences accounted for 83.5% (p < 0.001) of the observed variation with a global AMOVA ΦST of 0.86 (p < 0.000001). All pairwise ΦST estimates between Singapore, Borneo, Sumatra, Halmahera and the Philippines, were significantly different, rejecting the null hypothesis of a single, panmictic population (or possibly, taxonomic entity). Overall haplotype diversity was Hd = 0.962 (± 0.008), with both the Western and Eastern groups showing comparable levels of this genetic diversity when estimated independently (Table 3). Haplotype and nucleotide diversities were highest for the Philippines (Hd = 0.915, π = 0.0104), then for Singapore (Hd = 0.876, π = 0.0032), Borneo (Hd = 0.526, π = 0.0015), and finally Sulawesi (Hd = 0.500, π = 0.0008). Haplotype and nucleotide diversity for the remaining insular populations could not be estimated due to the small sample sizes (n < 4).

Discussion

Resolving the geographic distribution of the naturally occurring genetic variation of Malayopython reticulatus is an ongoing issue consisting of several distinct challenges. In addition to being shaped by the evolutionary diversification of the lineage, in response to the dynamic geographical template representing major landmasses in this realm [81,82], the contemporary geographical range of M. reticulatus has been influenced by several factors including its exceptional dispersal behaviour and ability to cross large stretches of water [30], environmental conditions [83], and anthropogenic translocation [31]. The latter issue has been highlighted by [84]who notes ‘the reticulated python … is carried around as a food animal and rat-catcher’ and ‘some of its occurrences on Wallacean islands’ are ‘possibly … due to human agency’ (see comments [14,85]). Trade dynamics of M. reticulatus seem likely to have contributed to the introduction of non-native populations, and to have likewise facilitated the potential for genetic homogenisation between what would otherwise have been isolated populations.

Despite the potential for confounding factors, we revealed significant structure among mitochondrial haplotypes across widely distributed Southeast Asian populations of M. reticulatus. The 34 haplotypes resolved across the species expansive geographical range were structured between two broad geographic realms. The first comprised landmasses on the Sunda Shelf (continental Southeast Asia, Sumatra, Java, Borneo), the Lesser Sundas, and the oceanic islands of Sulawesi, Ambon and Palawan of the Philippines; the second included only the Philippines and Halmahera. In the context of this study, hese two geographic regions have been termed the Western and Eastern haplotype groups, respectively; the distinction reflects the dominant geographic location of the haplotypes found, although the ‘Western’ group is nested within the paraphyletic ‘Eastern’ group. Although mitochondrial haplotypes within the Western haplotype group showed a high degree of similarity to previously published M. reticulatus sequences (98%–99%), those found in the Eastern group appear more divergent (94%–96% sequence similarity). This level of divergence has previously been reported for three subspecies of Python curtus that were elevated to distinct species [52]; each subspecies contained less than 1% sequence divergence within the taxon, between 3% and 8.9% divergence between the sub-species, and 10%–12.4% sequence divergence from P. curtus sister species, M. reticulatus.

The shared ancestry of the land-bridge islands of the Sunda Shelf (Sumatra, Java and Borneo) with continental Southeast Asia is a likely result of the region’s paleogeographic setting [29,86], which was characterised by dramatically lowered sea levels and expanded continental land area during the Pleistocene glacial maxima [4,81,8789]). Of the three islands, Sumatra and Java remain geographically proximal to both each other and to continental Southeast Asia. Despite the rising sea level which separated Sumatra and Java from the mainland, surface currents between the land masses are presumed not to have been strong enough to act as a barrier to the dispersal of M. reticulatus as demonstrated by the colonisation of Krakatau Island [30]. Active marine dispersal from Java to Sumatra and vice versa could be facilitated by the circulation of sea surface currents triggered by monsoon winds, i.e. December–February south-eastern currents from Sumatra to Java, then June–August north-western currents from Java to Sumatra [90,91]. The species may therefore have crossed marine barriers as has been reported for the estuarine crocodile (Crocodylus porosus) [92], however, it is also noted that the ‘human agency on boats is highly significant’ with regard to the dispersal of terrestrial vertebrates in the Sunda Strait [93]. Conversely, individuals from Borneo form a nested cluster within the Western haplotype group and share only a single haplotype with other locations in this group. This difference could be evidence of an effective barrier to dispersal imposed by the islands greater distance from continental Southeast Asia, Sumatra or Java but the lack of resolution of branching order within this clade prevents confident interpretation of these possibilities.

Interestingly, [85] surmised that M. reticulatus may not be native to the Philippines, and that human introductions might be responsible and explain the existence of the species in this archipelago. Divergence between individuals from the Philippines and members of the Western haplotype group (continental Southeast Asia, Sumatra, Java, Borneo, Lesser Sundas, and Sulawesi) is demarcated by Huxley’s modification of Wallace’s line, a theoretical, biogeographic boundary summarizing the eastern edge of the Sunda Shelf (Fig 4). Huxley’s modification of Wallace’s line separates the Philippines archipelago from the Sunda Shelf [94], with the ‘exception’ of Palawan which has been variably classified or associated with the Sunda Shelf [81,9598]). This theoretical barrier corresponds to the distinction between the Eastern and Western haplotype groups described here, and the presence of the Palawan haplotype (Haplotype 10) in the Western (Sunda Shelf) group is an interesting example of a large vertebrate, distributed according to the view that Palawn is a final extension of Sundaland [97]. The Eastern versus Western pattern of divergence elucidated here has previously been documented in a variety of taxa but is far from universal in terrestrial vertebrates [81,97,99,100].

thumbnail
Fig 4. Biogeographic realms.

Chinese trade routes to the North Moluccas. Dotted line indicates the eastern route used during the Yuan period (1279–1368). The dashed line indicates the route used by Chinese Merchants connecting Chinese overseas settlements to northern Java likely used during the Ming period (1368–1644). Faunal boundaries are indicated by the Wallace line [101] and Huxley's line [94]. Map revised after Ptak 1992 [102].

https://doi.org/10.1371/journal.pone.0182049.g004

Northern Halmahera is part of an arc system, which is shared by the Philippines and northern New Guinea [82,103]. The islands now reveal similar distribution patterns of e.g., insects and rodents [104], with Halmahera located at the southern Philippine trench [105107]. According to [108] the ‘east Halmahera-east Mindanao terrane’ is considered as an ‘eastward continuation of a Papuan arc complex’, termed as the Melanesian arc. This biogeographic affinity might explain the haplotype lineage found on Halmahera, that appears closely related to the Philippines; similar southern Mindanao biogeographic relationships have previously been inferred to involve Sulawesi [82,100,109]. This however, does not account for the divergence seen between Halmahera and Ambon which are both geographically close and located on the same biogeographic arc.

Geographical variants of mitochondrial DNA found on Ambon were more closely related to those present in Singapore and on Sumatra than nearby Halmahera. As part of the Seram Island Group (including the islands of the Outer Banda Arc, i.e. Seram, Ambon, Boano, Kelang, Maniapa, Harukuku and Saparua, see [110]), Ambon Island does have endemic species such as the python Morelia clastolepis [111] and monitor lizard Varanus cerambonensis [112], whereas other species are thought to have been introduced to Ambon e.g., marsupials, such as the northern common cuscus (Phalanger orientalis) and the common spotted cuscus (Spilocuscus maculatus) [113]. Historic anthropogenic movement of M. reticulatus individuals may account for the discordance between haplotype and geographic proximity across the Wallacea region as according to [31], M. reticulatus was also introduced to Ambon. Chinese trade routes that passed the Moluccas and Chinese overseas settlements in northern Java were in use during the Ming period and as the only Asian producers of cloves during this time, vessels would likely have visited the Moluccas for this commodity ([102], see Fig 4), and in particular, Ambon due to its accessibility along the shipping route. This may well have led to the introduction of M. reticulatus to Ambon following use aboard their ships as rat-catchers and food [102].

Although less divergent than the Philippines, and tentatively considered part of the Western haplotype group, Sulawesi shares no haplotypes with other group members. Representative of the Wallacean biogeographic faunal region and separated from the Sunda Shelf by deep underwater channels, divergence between faunal assemblages either side of Wallace’s line is well documented and a number of endemic species have been identified on Sulawesi, often further restricted to within one of seven areas of endemism [96,109]. A recent study also suggests that the population Malayopython r. jampeanus is genetically more closely related to the nominate form M. r. reticulatus, than to the geographical closer located M. reticulatus saputrai [34], suggesting genetic exchange between the Lesser Sundas and some of the Selayar Islands in the Flores Sea [114]. During the late Cretaceous, Borneo was still connected to western Sulawesi [87,115], thus enabling the exchange of terrestrial organisms; however, the formation of the Makassar Strait between both islands isolated terrestrial biota, explaining the high level of endemism on Sulawesi today (e.g., [116]). The Makassar Strait flow is considered a strong (sub-) surface current [117] that might prevent dispersal of M. reticulatus between both islands and explain the genetic distinction of Sulawesi haplotypes.

Despite the major tectonic collision events that occurred ca. 25 Million years ago, the current distribution patterns of terrestrial organisms in Southeast Asia were likely only formed during the last one Million years [29,86]. Formation of the world’s two largest archipelagos, Indonesia and the Philippines have, in part, resulted from very different tectonic processes [107]. Although major landmasses of the Greater Sundas and Palawan Island have been placed on the margin of the Eurasian plate, on the Sunda Shelf, most islands of the Philippines were paleotransported from the southeast [4,106]. Although the distribution of mitochondrial haplotypes for M. reticulatus spatially coincides with biogeographic expectations derived from the distribution of land and sea [86] and the geological mechanisms responsible for the formation of the archipelago [81], there are insufficient data to infer the exact number of colonisation events or the order in which they might have occurred [12]. Furthermore, comparison with haplotype sequences previously resolved for M. reticulatus is prevented by a paucity of equivalent data with associated locality information in Genbank. Nevertheless, the current dataset represents a step forward in that it is by far the most geographically comprehensive for M. reticulatus, Importantly, the phylogeographic structure inferred here suggests the existence of naturally occurring, properly documented genetic variation across the species distributional range. Although the sample size here is not conducive to fine-scale population genetic analysis, structuring within each of the main clades would be anticipated with the use of higher resolution markers, sampled from throughout the genome.

Malayopython reticulatus is currently managed as a single taxonomic unit across its range, and there is no consideration for genetic differences between harvested populations. Even though major countries of origin and export of M. reticulatus skins (Malaysia and Indonesia) have national/provincial harvest quotas [20], this management system is difficult to verify in terms of transparency and traceability [118]. There is no tool available to authenticate the likely origin of a skin except from the information included on the CITES exports permits, and listing the country of export does not automatically infer the country of origin. Illegal cross-border trade is known to exist, undoubtedly allowing unscrupulous traders to ‘circumvent national quotas’ [118], and loopholes in export reporting mechanisms have been repeatedly outlined in earlier studies [7,20,41,119,120]. In 2016, the provincial harvest quota was highest for South Sulawesi, followed by North Sumatra and West Kalimantan (Borneo) (see Table 4). The need to investigate the impact of regional harvest has been highlighted for the populations M. reticulatus saputrai and M. r. jampeanus [34]; however, trade in both is on-going with the former, to date, outnumbering provincial quotas set for Sumatra and Kalimantan (Anon. in litt. to Auliya; 22 Nov. 2016; see Table 4). Although there is no information available to address the impact that regional trade may have on the sustainability of different, taxonomically distinct, range restricted, populations of M. reticulatus, specific management tools are warranted to permit the identification of python skin origin, to ensure harvest quotas are adhered to, in order to maintain the viability of genetically distinct populations.

thumbnail
Table 4. Regional harvest quotas for the reticulated python.

https://doi.org/10.1371/journal.pone.0182049.t004

Despite uncertainty of the impact current harvest levels will have, particularly on ecological and long term on genetic variability, no precautionary measures have been established to prevent local declines and extinctions. [121] provided substantial biological and ecological information on the species harvested on the island of Sumatra. Although their study indicated apparent sustainability, the authors also note that ‘far more information will be needed before we can confidently assess sustainable levels of offtake of Indonesian pythons’, and the same conclusions were drawn within a study published 17 years later [122]. Both studies provide regional-based results (Sumatra) and the latter authors conclude that ‘the harvest appears to be sustainable’, ignoring the entire issue of trade dynamics (legal and illegal within a country and between countries), and thus genetic variation as a result of local adaptation. Therefore, the identification of sustainability indices will be required (see [123]). The phylogeographic inference presented here suggests that geographically distinct clades (in particular, the Philippines and continental Southeast Asia, but also Borneo) warrant separate consideration as Evolutionary Significant Units. Furthermore, the level of divergence estimated for populations in the Philippines would also seem to justify a taxonomic review, and it is conceivable that phylogenetically distinct species may be present. Divergent, evolutionary distinct lineages should be taken into consideration when establishing harvest quotas to avoid over-exploitation and the erosion of intra-specific genetic variation [11]. This requires a thorough understanding of the fine scale species structure and the source/skin dynamics of populations at different scales [124].

Because genetic divergence, lineage isolation, and even speciation are not necessarily accompanied by morphological divergence (see [125,126]), the use of DNA-based tools is a logical development for verifying the origin of traded python skins. Skins that are destined for incorporation to the luxury goods industry are heavily treated with a cocktail of chemicals, and so it is necessary to ensure that the same genetic signals can be detected in processed skins as in untreated ones before developing forensic tests and databases. Chemically treated skins can yield low concentrations of heavily degraded and fragmented DNA [127] and therefore, it is necessary to ensure that any DNA-based test devised be feasible for use with poor quality samples, from which only degraded DNA may be extracted. Although phylogeographic analysis of shorter mitochondrial sequence fragments did not retain the same number of haplotypes that were evident when using concatenated sequences, it did continue to distinguish the Western and Eastern haplotype groups, as well as the Sulawesi population. This demonstrates that the sequence data may be used to assign individuals of ‘unknown origin’ [11], back to either Western or Eastern genotypes. This result suggests that heavily processed samples of ‘unknown origin’ within the python skin trade (where amplification is restricted to small fragments of DNA; [128]) could be localised to geographic origin by comparison to the data presented here.

Conclusions and recommendations

Although the results presented here offer a very provisional insight to the genetic structure of Malayopython reticulatus across the species contemporary range, they offer an encouraging conservation genetics baseline that justifies implementing a precautionary approach [129] to population management, and further consideration and investigation is necessary.

If genetically distinct populations of Malayopython reticulatus are to be continually exploited for the skin trade, management strategies should not solely be tailored to obtaining the maximum economic yield. Populations should instead be managed following an adaptive management scheme to ensure their long-term sustainability (see [130]). Among the complex geography of continental and insular Southeast Asia, trade dynamics currently outbalance sustainable measures [3,18], and thus effective management practises for preventing local over-exploitation are not established to address illegal trade activities and trace origins of sourced populations [20,41,131]

Genetic variation is fundamental for species conservation [96]. To safeguard genetically distinct populations, prevent genetic erosion and retain their ability to undergo evolutionary change, it is important to manage populations in such a way that they do not fall under a viable threshold [132]. Based on the results of this work, the following conclusions and recommendations are offered:

  1. Whilst caution is to be exercised in inferring taxonomic divisions within the current dataset due to the limited availability of quality reference samples, the genetic distinctiveness of the Philippine population suggests that the status of M. reticulatus within this archipelago warrants further investigation and possible review.
  2. Distinct genetic structure across the ecoregions and faunal regions suggests that variation in heritage should be incorporated into the implementation of regional management schemes (also see [55]), and may provide a basis for the development of regional skin traceability systems.
  3. To regulate trade dynamics, wild harvest rates and exploitation levels require regular monitoring of the status of resource population(s) to reduce uncertainties.
  4. Further fine scale genetic analysis is also warranted to delineate local genetic partitions, and the potential application of adaptive genetic markers to establish conservation units should be considered (see [133]).
  5. To encourage and foster scientific collaborative networks with the countries of origin to permit construction of comprehensive, reference sample databases on which to base the development of robust traceability protocols.

Supporting information

S1 Table. Details of samples included in this study.

The table contains information about samples included in this study along with the corresponding haplotypes inferred from sequence variation across a mitochondrial cytochrome b fragment.

https://doi.org/10.1371/journal.pone.0182049.s001

(PDF)

Acknowledgments

M. Auliya expresses sincere thanks to Guy Morgan (BSR, Business for Social Responsibility, Paris) for establishing communication with LVMH (Moët Hennessy Louis Vuitton). Thanks to LVMH for support from Daniel Peterlin, Alexandre Campelli, Emmanuel Mathieu, Alessandra Colombo, Isabelle Laronde, Sylvie Bénard, Jerome de Rinaldo and Delphine Raimbert and further thanks to the following scientific research institutions, zoological museums, tanneries and individuals for providing voucher specimens and tissue samples: Berlin, Museum für Naturkunde, Mark-Oliver Rödel, Frank Tillack; Bonn, Zoologisches Forschungsmuseum Alexander Koenig, Ursula Bott, Claudia Koch, Dennis Rödder, Wolfgang Böhme; Florence, Dolmen Spa, Luca Giananti; Frankfurt, Senckenberg Forschungsinstitut und Naturmuseum, Linda Acker; Hamburg, Biozentrum Grindel und Zoologisches Museum, Alexander Haas, Jakob Hallermann; Kansas, University of Kansas, Cameron D. Siler, Lori Bryn Schlenker, Matt Buehler; Kuching, Universiti Malaysia Sarawak (UNIMAS), Indraneil Das; Leiden, naturalis, Esther Dondorp, Caroline Perpermans, Pim Arntzen; London, Natural History Museum, Patrick Campbell; Osnabrück, professional breeder, Karsten Woellner; Paris, Muséum national d'Histoire naturelle, Annemarie Ohler; Singapore, Singapore Zoological Gardens & Wildlife Reserves Singapore, Sonja Luz, Timothy Tan; Singapore, Heng Long tannery; Vienna, Natural History Museum Vienna, Heinz Grillitsch, Silke Schweiger. Thanks to: Edinburgh, Royal Zoological Society of Scotland WildGenes laboratory, Helen Senn, Jo Howard-McCombe, Jennifer Kaden; German CITES Management Authorities, Franz Boehmer, Stefan Raths; CITES Management Authority of Sarawak, Engkamat Anak Lading for issuing CITES import/export permits. Samples were collected under the following permits: Philippine samples (supported by NSF DEB 0743491, 1418895 to RMB)—Foreign CITES Permit No. DE 208–02; U.S. CITES Permit No. 132US705634/9; Singapore samples—CITES Permit No. 14SG101157CE.

References

  1. 1. Laurance WF, Sayer J, Cassman KG. Agricultural expansion and its impacts on tropical nature. Trends Ecol Evol. 2014;29: 107–116. pmid:24388286
  2. 2. Sodhi NS, Koh LP, Brook BW, Ng PKL. Southeast Asian biodiversity: An impending disaster. Trends Ecol Evol. 2004;19: 654–660. pmid:16701328
  3. 3. Auliya M. Taxonomy, Life History and Conservation of Giant Reptiles in West Kalimantan (Indonesian Borneo) [Internet]. 2006. Available: https://www.amazon.com/Taxonomy-History-Conservation-Reptiles-Kalimantan/dp/3937285520
  4. 4. Brown RM, Diesmos AC. Philippines Biology. In: Gillespie R, Clague D, editors. Encyclopedia of Islands. University of California Press, Berkley; 2009. pp. 723–732.
  5. 5. Maxwell SL, Fuller RA, Brooks TM, Watson JEM. The ravages of guns, nets and bulldozers. Comment Nat. 2016;536: 143–145. pmid:27510207
  6. 6. TRAFFIC. What’s Driving the Wildlife Trade? A Review of Expert Opinion on Economic and Trade Control Efforts in Cambodia, Indonesia, Lao PDR, and Vietnam. 2008; 120. Available: www.traffic.org/general-reports/traffic_pub_gen24.pdf
  7. 7. Lyons J a., Natusch DJD. Wildlife laundering through breeding farms: Illegal harvest, population declines and a means of regulating the trade of green pythons (Morelia viridis) from Indonesia. Biol Conserv. 2011;144: 3073–3081.
  8. 8. Nijman V, Shepherd CR, Mumpuni , Sanders KL. Over-exploitation and illegal trade of reptiles in Indonesia. Herpetol J. 2012;22: 83–89.
  9. 9. Scheffers BR, Corlett RT, Diesmos A, Laurance WF. Local demand drives a bushmeat industry in a Philippine forest preserve. Trop Conserv Sci. 2012;5: 133–141.
  10. 10. Böhm M, Collen B, Baillie JEM, Bowles P, Chanson J, Cox N, et al. The conservation status of the world’s reptiles. Biol Conserv. 2013;157: 372–385.
  11. 11. Welton LJ, Siler CD, Linkem CW, Diesmos AC, Diesmos ML, Sy E, et al. Dragons in our midst: Phyloforensics of illegally traded Southeast Asian monitor lizards. Biol Conserv. Elsevier Ltd; 2013;159: 7–15.
  12. 12. Siler CD, Lira-Noriega A, Brown RM. Conservation genetics of Australasian sailfin lizards: Flagship species threatened by coastal development and insufficient protected area coverage. Biol Conserv. Elsevier Ltd; 2014;169: 100–108.
  13. 13. Auliya M, Altherr S, Ariano-Sanchez D, Baard EH, Brown C, Brown RM, et al. Trade in live reptiles, its impact on wild populations, and the role of the European market. Biol Conserv. Elsevier Ltd; 2016;204: 103–119.
  14. 14. Headland TN, Greene HW. Hunter-gatherers and other primates as prey, predators, and competitors of snakes. Proc Natl Acad Sci. 2011;108: E1470–E1474. pmid:22160702
  15. 15. Jenkins M, Broad S. International Trade in Reptile Skins [Internet]. Cambridge: TRAFFIC International; 1994. Available: https://portals.iucn.org/library/sites/library/files/documents/Traf-013.pdf
  16. 16. Pope CH. The Giant Snakes. Routledge & Kegan Paul, London; 1962.
  17. 17. Taylor EH. Philippine adventures: An autobiographical memoir. Recollections of an Herpetologist. University of Kansas Museum of Natural History (Monograph No. 4), Lawrance, KS, USA;
  18. 18. Groombridge B, Luxmoore R. Pythons in South-East Asia: A review of distribution, status and trade in three selected species. A report to the CITES Secretariat. Lausanne, Switzerland; 1990.
  19. 19. Klemens MW, Thorbjarnarson JB. Reptiles as a food resource. Biodivers Conserv. 1995;298: 281–298.
  20. 20. Kasterine A, Arbeid R, Caillabet O, Natusch D. The trade in South-east Asian python skins. Geneva; 2012.
  21. 21. Martin DL. Identification of reptile skin products using scale morphology. In: Huffman JE, Wallace JR, editors. Wildlife Forensics: methods and applications. 1st ed. Wiley-Blackwell, Chichester; 2012. pp. 161–199.
  22. 22. Rawlings LH, Rabosky DL. Python phylogenetics: inference from morphology and mitochondrial DNA Biological Journal of the Linnean Society 2008;1: 603–619. Available: http://onlinelibrary.wiley.com/doi/10.1111/j.1095-8312.2007.00904.x/full%5Cnpapers2://publication/uuid/BDD46442-DCA1-4E96-88EA-5A7AF7AC1F7F
  23. 23. Pyron R, Burbrink FT, Wiens JJ. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol Biol. 2013;13: 93. pmid:23627680
  24. 24. Hoser RT. A reclassification of hte Pythoninae including the description of two new genera, two new species and nine new subspecies. Crocodilian J Vic Assoc Amatuer Herpetol. 2004;4: 21–40, 31–37.
  25. 25. Kaiser H, Crother BI, Kelly CMR, Luiselli L, O’Shea M, Ota H, et al. Best Practices: In the 21 st Century, Taxonomic Decisions in Herpetology are Acceptable Only When Supported by a Body of Evidence and Published via Peer-Review. Herpertological Rev. 2011;44: 8–23.
  26. 26. Graham Reynolds R, Niemiller ML, Revell LJ. Toward a tree-of-life for the boas and pythons: Multilocus species-level phylogeny with unprecedented taxon sampling. Mol Phylogenet Evol. Elsevier Inc.; 2014;71: 201–213. pmid:24315866
  27. 27. Barker DG, Barker TM, Davis MA, Schuett GW. A review of the systematics and taxonomy of Pythonidae: An ancient serpent lineage. Zool J Linn Soc. 2015;175: 1–19.
  28. 28. Uetz P, Freed P, Hošek J. The Reptile Database [Internet]. 2017 [cited 21 Mar 2017]. Available: http://www.reptile-database.org
  29. 29. Hall R. The plate tectonics of Cenozoic SE Asia and the distribution of land and sea. Biogeogr Geol Evol SE Asia. 1998; 99–131.
  30. 30. Thornton I. Krakatau: the destruction and reassembly of an island ecosystem. Cambridge, Massachusetts: Harvard University Press; 1996.
  31. 31. Brehm A. Band 22. Reptilien II: Giftlose Schlangen—Giftschlangen. In: Meyer A, editor. Brehms Tierleben. Gutenberg-Verlag, Hamburg; 1927.
  32. 32. Barker D, Barker T. Big & Beautiful: The Reticulated Python. Reptiles. 1997;6: 48–67.
  33. 33. Yuwono FB. The Trade of Live Reptiles in Indonesia. In: Erdelen W, editor. Conservation, Trade and Sustainable Use of Lizards and Snakes in Indonesia. Mertensiella 9, Rheinbach, Germany; 1998. p. 144.
  34. 34. Auliya M, Mausfeld P, Schmitz a., Böhme W. Review of the reticulated python (Python reticulatus Schneider, 1801) with the description of new subspecies from Indonesia. Naturwissenschaften. 2002;89: 201–213. pmid:12135085
  35. 35. The Checklist of CITES Species Website. In: UNEP-WCMC (Comps.) [Internet]. 2015 [cited 30 Sep 2015]. Available: http://checklist.cites.org
  36. 36. CITES. The CITES Appendicies [Internet]. 2015 [cited 30 Sep 2015]. Available: https://www.cites.org/eng/app/index.php
  37. 37. Thomson J. Captive breeding of selected taxa in Cambodia and Viet Nam: A reference manual for farm operators and CITES authorities. 2008.
  38. 38. Luiselli L, Bonnet X, Rocco M, Amori G. Conservation Implications of Rapid Shifts in the Trade of Wild African and Asian Pythons. Biotropica. 2012;44: 569–573.
  39. 39. Woodruff DS. Biogeography and conservation in Southeast Asia: How 2.7 million years of repeated environmental fluctuations affect today’s patterns and the future of the remaining refugial-phase biodiversity. Biodivers Conserv. 2010;19: 919–941.
  40. 40. Webb GJW, Manolis C, Jenkins RWG. Improving international trade in reptile skins based on sustainable use. United Nations Conf Trade Dev. 2012; 1–9. Available: http://unctad.org/en/PublicationsLibrary/ditcted2011d7_en.pdf
  41. 41. Ashley D. Traceability systems for a Sustainable International Trade in Reptile Skins based on Sustainable Use [Internet]. 2014. Available: http://unctad.org/en/PublicationsLibrary/ditcted2013d6_en.pdf
  42. 42. Alacs EA, Georges A, FitzSimmons NN, Robertson J. DNA detective: A review of molecular approaches to wildlife forensics. Forensic Sci Med Pathol. 2010;6: 180–194. pmid:20013321
  43. 43. Yau FCF, Wong KL, Shaw PC, But PPH, Wang J. Authentication of snakes used in Chinese medicine by sequence characterized amplified region (SCAR). Biodivers Conserv. 2002;11: 1653–1662.
  44. 44. Cao S, Guo L, Luo H, Yuan H, Chen S, Zheng J, et al. Application of COI barcode sequence for the identification of snake medicine (Zaocys). Mitochondrial DNA Part A. 2016;27: 483–489. pmid:24857374
  45. 45. Lo CF, Lin YR, Chang HC, Lin JH. Identification of turtle shell and its preparations by PCR-DNA sequencing method. J Food Drug Anal. 2006;14: 153–158.
  46. 46. Ogden R, Linacre A. Wildlife forensic science: A review of genetic geographic origin assignment. Forensic Sci Int Genet. Elsevier Ireland Ltd; 2015;18: 152–159. pmid:25795277
  47. 47. Mucci N, Mengoni C, Randi E. Wildlife DNA forensics against crime: Resolution of a case of tortoise theft. Forensic Sci Int Genet. Elsevier Ireland Ltd; 2014;8: 200–202. pmid:24315609
  48. 48. Austin CC, Spataro M, Peterson S, Jordan J, McVay JD. Conservation genetics of Boelen’s python (Morelia boeleni) from New Guinea: Reduced genetic diversity and divergence of captive and wild animals. Conserv Genet. 2010;11: 889–896.
  49. 49. Wasser SK, Joseph Clark W, Drori O, Stephen Kisamo E, Mailand C, Mutayoba B, et al. Combating the illegal trade in African elephant ivory with DNA forensics. Conserv Biol. 2008;22: 1065–1071. pmid:18786100
  50. 50. Benavides MT, Horn RL, Feldheim KA, Shivji MS, Clarke SC, Wintner S, et al. Global phylogeography of the dusky shark Carcharhinus obscurus: Implications for fisheries management and monitoring the shark fin trade. Endanger Species Res. 2011;14: 13–22.
  51. 51. Martinsohn JT, Geffen AJ, Maes GE, Nielsen EE, Waples RS, Carvalho GR. Tracing Fish and fish products from ocena to fork using advanced molecular technologies. In: Hoorfar J, Jordan K, Prugga R, editors. Food Chain Integrity: A holistic approach to food traceability, safety, quality and authenticity. Woodhead Publishing; 2011.
  52. 52. Koegh J, Barker D, Shine R. Exploited but poorly known: systematics and biogeography of commercially harvested python. Biol J Linn Soc. 2001; 113–129.
  53. 53. de Buffrénil V. Les varans Africains: Varanus niloticus et Varanus exanthematicus. Donnees de synthese sur leur biologie et leur exploitation. Geneva (Switzerland) PNUE; 1993.
  54. 54. de Buffrénil V, Castanet J. Age Estimation by Skeletochronology in the Nile Monitor (Varanus niloticus), a Highly Exploited Species. J Herpetol. 2016;34: 414–424.
  55. 55. Dowell SA, de Buffrénil V, Kolokotronis S- O, Hekkala ER. Fine-scale genetic analysis of the exploited Nile monitor (Varanus niloticus) in Sahelian Africa. BMC Genet. 2015;16: 32. pmid:25884730
  56. 56. Natusch DJD, Lyons J a. Assessment of Python Breeding Farms Supplying the International High-end Leather Industry. 2014.
  57. 57. Nossal K, Livingstone DG, Aust P, Kasterine A, Ngo Viet C, Nguyen V, et al. Trade in python skins: impact on livelihoods in viet nam. Geneva; 2016.
  58. 58. Baird IG. Field Report No. 3: Wildlife trade between the southern Lao PDR provinces of Champasak, Sekong and Attapeu and Thailand, Cambodia and Viet Nam. Kuala Lumpur; 1993.
  59. 59. Stuart BL. Traffic Bulletin July 2004. Traffic Bull. 2004;20: 53.
  60. 60. WCS. Commercial wildlife farms in Vietnam: A problem or solution for conservation? Hanoi, Viet Nam; 2008.
  61. 61. Information and analysis bulletin on animal poaching and smuggling. In: On the Trail [Internet]. 2016. Available: www.robindesbois.org/wp-content/uploads/ON_THE_TRAIL_11.pdf
  62. 62. Information and analysis bulletin on animal poaching and smuggling. In: In the Trail [Internet]. 2016 p. 12. Available: www.robindesbois.org/wp-content/uploads/ON_THE_TRAIL_12.pdf
  63. 63. Nijman V, Shepherd CR. Wildlife Trade From ASEAN To the EU: Issues with the trade in captive-bred reptiles from Indonesia. TRAFFIC Europe Report for the European Commission. Brussels, Belgium; 2007.
  64. 64. Spinks PQ, Thomson RC, Hughes B, Moxley B, Brown R, Diesmos A, et al. Cryptic variation and the tragedy of unrecognized taxa: The case of international trade in the spiny turtle Heosemys spinosa (Testudines: Geoemydidae). Zool J Linn Soc. 2012;164: 811–824.
  65. 65. Card DC, Schield DR, Adams RH, Corbin AB, Perry BW, Andrew AL, et al. Phylogeographic and population genetic analyses reveal multiple species of Boa and independent origins of insular dwarfism. Mol Phylogenet Evol. Elsevier Inc.; 2016;102: 104–116. pmid:27241629
  66. 66. Moritz C. Applications of mitochondrial DNA analysis in conservation: a critical review. Mol Ecol. 1994;3: 401–411.
  67. 67. Kumazawa Y, Ota H, Nishida M, Ozawa T. Gene rearrangements in snake mitochondrial genomes: highly concerted evolution of control-region-like sequences duplicated and inserted into a tRNA gene cluster. Mol Biol Evol. 1996;13: 1242–1254. pmid:8896377
  68. 68. Kocher TD, Thomas WK, Meyer a, Edwards S V, Pääbo S, Villablanca FX, et al. Dynamics of mitochondrial DNA evolution in animals: amplification and sequencing with conserved primers. Proc Natl Acad Sci U S A. 1989;86: 6196–200. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=297804&tool=pmcentrez&rendertype=abstract pmid:2762322
  69. 69. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30: 2725–2729. pmid:24132122
  70. 70. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994;22: 4673–80. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=308517&tool=pmcentrez&rendertype=abstract pmid:7984417
  71. 71. Guindon S, Gascuel O. A Simple, Fast, and Accurate Algorithm to Estimate Large Phylogenies by Maximum Likelihood. Syst Biol. 2003;52: 696–704. pmid:14530136
  72. 72. Darriba D, Posada D. jModelTest 2.0 Manual. Nat methods2. 2012;9: 772. Available: http://code.google.com/p/jmodeltest2
  73. 73. Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Syst Biol. 2010;59: 307–321. pmid:20525638
  74. 74. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28: 1647–1649. pmid:22543367
  75. 75. Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics. 2001;17: 754–755. pmid:11524383
  76. 76. Rodrigo A, Bertels F, Heled J, Noder R, Shearman H, Tsai P. The perils of plenty: what are we going to do with all these genes? Philos Trans R Soc B Biol Sci. 2008;363: 3893–3902. pmid:18852100
  77. 77. Leigh JW, Bryant D. POPART: full-feature software for haplotype network construction. 2015; 1110–1116. 10.1111/2041-210X.12410
  78. 78. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25: 1451–2. pmid:19346325
  79. 79. Petit RJ, Mousadik A El, Pons O. Identifying populations for consevation on the basis of genetic markers. Conserv Biol. 1998;12: 844–855.
  80. 80. Excoffier L, Lischer HEL. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10: 564–567. pmid:21565059
  81. 81. Brown RM, Siler CD, Oliveros CH, Esselstyn JA, Diesmos AC, Hosner PA, et al. Evolutionary Processes of Diversification in a Model Island Archipelago. Annu Rev Ecol Evol Syst. 2013;44: 411–435.
  82. 82. Lohman DJ, de Bruyn M, Page T, von Rintelen K, Hall R, Ng PKL, et al. Biogeography of the Indo-Australian Archipelago. Annu Rev Ecol Evol Syst. 2011;42: 205–226.
  83. 83. Ims RA, Hjermann DØ. Condition-dependent dispersal. In: Clobert J, Dachin E, Dhont AA, Nichols JD, editors. Dispersal. Oxford University Press, New York; 2001. pp. 203–216.
  84. 84. Heinsohn TE. The Realm of the Cuscus: Animal Translocation and Tiological Invasion East of Wallace’s Line. Australian National University. 1998.
  85. 85. Brown WC, Alcala AC. The zoogeography of the herpetofauna of the Philippine Islands, a fringing archipelago. Proc Calif Acad Sci. 1970;38: 105–130. Available: http://biostor.org/reference/3125
  86. 86. Hall R. The palaeogeography of Sundaland and Wallacea since the Late Jurassic. J Limnol. 2013;72: 1–17.
  87. 87. Moss SJ, Wilson MEJ. Biogeographic implications of the Tertiary palaeogeographic evolution of Sulawesi and Borneo. Biogeogr Geol Evol SE Asia. 1998; 133–163.
  88. 88. Voris HK. Maps of Pleistocene sea levels in Southeast Asia: Shorelines, river systems and time durations. J Biogeogr. 2000;27: 1153–1167.
  89. 89. Bird MI, Taylor D, Hunt C. Palaeoenvironments of insular Southeast Asia during the Last Glacial Period: A savanna corridor in Sundaland? Quat Sci Rev. 2005;24: 2228–2242.
  90. 90. Birowo S, Uktolseja H. Oceanographic features of Sunda Strait. Proceedings of Workshop on coastal resources management of Krakatau and the Sunda Strait region, Indonesia Indones Instit Sci United Nation Univ. 1983. pp. 54–75.
  91. 91. Daryabor F, Ooi SH, Samah AA, Akbari A. Dynamics of the water circulations in the Southern South China Sea and its seasonal transports. PLoS One. 2016;11: 158415. pmid:27410682
  92. 92. Campbell HA, Watts ME, Sullivan S, Read MA, Choukroun S, Irwin SR, et al. Estuarine crocodiles ride surface currents to facilitate long-distance travel. J Anim Ecol. 2010;79: 955–964. pmid:20546063
  93. 93. Rawlinson PA, Widjoya AHT, Hutchinson MN, Brown GW. The Terrestrial Vertebrate Fauna of the Krakatau Islands, Sunda Strait, 1883–1986. Philos Trans R Soc London B, Biol Sci. 1990;328: 3 LP–28. Available: http://rstb.royalsocietypublishing.org/content/328/1245/3.abstract
  94. 94. Huxley TH. On the classification and distribution of the Alectoromorphae and Heteromorphae. Proc Zool Soc London. 1868; 294–319.
  95. 95. Brown RM, Guttman SI. Phylogenetic systematics of the Rana signata complex of Philippine and Bornean stream frogs: reconsideration of Huxley’ s modification of Wallace’ s Line at the Oriental–Australian faunal zone interface. Biol J Linn Soc. 2002;76: 393–461.
  96. 96. Evans BJ, Brown RM, McGuire JA, Supriatna J, Andayani N, Diesmos A, et al. Phylogenetics of Fanged Frogs: Testing Biogeographical Hypotheses at the Interface of the Asian and Australian Faunal Zones. Syst Biol. 2003;52: 794–819. pmid:14668118
  97. 97. Esselstyn JA, Oliveros CH, Moyle RG, Peterson AT, McGuire JA, Brown RM. Integrating phylogenetic and taxonomic evidence illuminates complex biogeographic patterns along Huxley’s modification of Wallace’s Line. J Biogeogr. 2010;37: 2054–2066.
  98. 98. Piper PJ, Ochoa J, Robles EC, Lewis H, Paz V. Palaeozoology of Palawan Island, Philippines. Quat Int. Elsevier Ltd; 2011;233: 142–158.
  99. 99. De Bruyn M, Wilson JA, Mather PB. Huxley’s line demarcates extensive genetic divergence between eastern and western forms of the giant freshwater prawn, Macrobrachium rosenbergii. Mol Phylogenet Evol. 2004;30: 251–257. pmid:15022775
  100. 100. Linkem CW, Brown RM, Siler CD, Evans BJ, Austin CC, Iskandar DT, et al. Stochastic faunal exchanges drive diversification in widespread Wallacean and Pacific island lizards (Squamata: Scincidae: Lamprolepis smaragdina). J Biogeogr. 2013;40: 507–520.
  101. 101. Wallace AR. On the zoological geography of the Malay Archipelago. J Proc Linn Soc Zool. 1860;4: 172–184.
  102. 102. Ptak R. The Northern Trade Route to the Spice IslandsSouth China Sea—Sulu Zone—North Moluccas (14th to early 16th century). Archipel. 1992;43: 27–56.
  103. 103. Hall R. Cenozoic reconstructions of SE Asia and the SW Pacific: changing pattern of land and sea. In: Metcalfe I, Smith JMB, Morwood M, Davidson ID, editors. Faunal and Floral Migrations and Evolution in SE Asia-Australasia. Swets & Zeitlinger; 2001. pp. 35–56.
  104. 104. Holloway JD. Sulawesi biogeography—discussion and summing up. In: Knight WJ, Holloway JD, editors. Insects and the rainforests of Southeast Asia. Royal Zoological Society of London; 1990. pp. 95–102.
  105. 105. Michaux B. Terrestrial birds of the Indo-Pacific. In: Hall R, Holloway JD, editors. Biogeography and Geological Evolution of SE-Asia. Backhuys Publishing, Leiden; 1998. pp. 361–391.
  106. 106. Yumul GP, Dimalanta CB, Maglambayarw VB, Marquez EJ. Tectonic setting of a composite terrane: A review of the Philippine island arc system. Geosci J. 2008;12: 7–17.
  107. 107. Yumul GP, Dimalanta CB, Queaño K, E M. Philippines, geology. In: Gillespie RG, Clague DA, editors. Encyclopedia of Islands. Berkley, University of California; 2009. pp. 732–738.
  108. 108. Hall R. Plate boundary evolution in the Halmahera region, Indonesia. Tectonophysics. 1987;144: 337–352.
  109. 109. Evans BJ, Supriatna J, Andayani N, Setiadi MI, Cannatella DC, Melnick DJ. Monkeys and toads define areas of endemism on Sulawesi. Evolution (N Y). Blackwell Publishing Ltd; 2003;57: 1436–1443. https://doi.org/10.1111/j.0014-3820.2003.tb00350.x
  110. 110. Natus IR. Biodiversity and endemic centres of Indonesian terrestrial vertebrates. University of Trier. 2005.
  111. 111. Harvey MB, Barker DG, Ammerman LK, Chippindale PT. Systematics of Pythons of the Morelia amethistina Complex (Serpentes : Boidae) with the Description of Three New Species. Herpetol Monogr. 2017;14: 139–185.
  112. 112. Philipp K, Böhme W, Ziegler T. The identifiy of Varanus indicus: Redfinition and description of a sibling species coexisting at the type locality. Spixiana. München :Zoologische Staatssammlung München,; 1999;22: 273–287. Available: http://www.biodiversitylibrary.org/item/89816
  113. 113. Heinsohn TE. Marsupials as introduced species: Long-term anthropogenic expansion of the marsupial frontier and its implications for zoogeographic interpretation. Altered Ecol Fire, Clim Hum Influ Terr landscapes. 2010; 133–176. Available: internal-pdf://heinsohn-0312701185/Heinsohn.pdf
  114. 114. Hanifa BF, Nugraha AP, Nanda IF, Daryono BS. Phylogenetic analysis of Malayopython reticulatus (Schneider, 1801) from Southern Sulawesi based on morphological and molecular character. 2016;20008: 20008. 10.1063/1.4953482
  115. 115. Hall R. Cenozoic geological and plate tectonic evolution of SE Asia and the SW Pacific: Computer-based reconstructions, model and animations. J Asian Earth Sci. 2002;20: 353–431.
  116. 116. Riley EP. The endemic seven: Four decades of research on the Sulawesi macaques. Evol Anthropol Issues, News, Rev. Wiley Subscription Services, Inc., A Wiley Company; 2010;19: 22–36.
  117. 117. Mayer B, Damm PE. The Makassar Strait throughflow and its jet. J Geophys Res Ocean. 2012;117: 1–14.
  118. 118. Natusch D, Lyons JA, Mumpuni, Riyanto A, Khadiejah S, Mustapha N, et al. Sustainable Management of the Trade in Reticulated Python Skins in Indonesia and Malaysia. 2016. 10.2305/IUCN.CH.2016.SSC-OP.61.en
  119. 119. Holden J. By Hook or by Crook. A reference manual on illegal wildlife trade and prosecutions in the United Kingdom. The Royal Society for the Protection of Birds; 1998.
  120. 120. Phelps J, Webb EL, Bickford D, Nijman V, Sodhi NS. Boosting CITES. Science (80-). American Association for the Advancement of Science; 2010;330: 1752–1753. pmid:21205655
  121. 121. Shine R, Ambariyanto , Harlow PS, Mumpuni . Reticulated pythons in Sumatra: Biology, harvesting and sustainability. Biol Conserv. 1998;87: 349–357.
  122. 122. Natusch DJD, Lyons JA, Mumpuni , Riyanto A, Shine R. Jungle giants: Assessing sustainable harvesting in a difficult-to-survey species (Python reticulatus). PLoS One. 2016;11: 1–14. pmid:27391138
  123. 123. Milner-Gullaud EJ, Akçakaya HR. Sustainability indices for exploited populations under uncertainty. Tredns Ecol Evol. 2001;16: 686–692.
  124. 124. Weinbaum KZ, Brashares JS, Golden CD, Getz WM. Searching for sustainability: Are assessments of wildlife harvests behind the times? Ecol Lett. 2013;16: 99–111. pmid:23062121
  125. 125. Rawlings LH, Donnellan SC. Phylogeographic analysis of the green python, Morelia viridis, reveals cryptic diversity. Mol Phylogenet Evol. 2003;27: 36–44. pmid:12679069
  126. 126. Barley AJ, White J, Diesmos AC, Brown RM. The challenge of species delimitation at the extremes: Diversification without morphological change in philippine sun skinks. Evolution (N Y). 2013;67: 3556–3572. pmid:24299408
  127. 127. Ojeda GN, Amavet PS, Rueda EC, Siroski PA. DNA extraction from skins of wild (Hydrochoerus hydrochaeris and Pecari tajacu) and domestic (Sus scrofa domestica) species using a novel protocol. Genet Mol Res. 2012;11: 672–678. pmid:22535403
  128. 128. Dubey B, Meganathan PR, Haque I. DNA mini-barcoding: An approach for forensic identification of some endangered Indian snake species. Forensic Sci Int Genet. Elsevier Ireland Ltd; 2011;5: 181–184. pmid:20457097
  129. 129. Cooney R, Dickson B. Biodiversity and the Precautionary Principle: “Risk, Uncertainty and Practice in Conservation and Sustainable Use” [Internet]. Taylor & Francis; 2012. Available: https://books.google.co.uk/books?id=HoIDoJUJw7sC
  130. 130. Heino M, Enberg K. Sustainable Use of Populations and Overexploitation. Encyclopedia of Life Sciences. John Wiley & Sons, Chichester, U.K.; 2008.
  131. 131. Howes BJ, Pither R, Prior KA. Conservation implications should guide the application of conservation genetics research. Endanger Species Res. 2009;8: 193–199.
  132. 132. Frankham R, Ballou JD, Briscoe DA, McInnes KH. Introduction to Conservation Genetics. Cambridge University Press, Cambridge, U.K.; 2009.
  133. 133. Gebremedhin B, Ficetola GF, Naderi S, Rezaei H-R, Maudet C, Rioux D, et al. Frontiers in identifying conservation units: from neutral markers to adaptive genetic variation. Anim Conserv. 2009;12: 107–109.