A non-invasive buccal swab method was applied to collect samples from 62 ball python specimens in Togo, including 21 samples from five of the nine known python farms currently located in Lomé and 41 samples from 12 wild populations located outside Lomé (Fig. 2). The number of samples differed between python farms and between locations of wild populations (ranging from 1–11) (Fig. 3, Suppl. material 1: Tables S1, S3).

Figure 2.

Southern Togo with sampling locations for wild Python regius (green) and housed in farms (red). Five specimens were sampled at each of three python farms, and at two farms each an additional three specimens were sampled. Number of wild specimens sampled per location: Hangoume, n = 3; Dagbati, n = 2, Kpove, n = 3, Assahoun, n = 2; Tado, n = 3; Tsevié, n = 3; Agbave, n = 1; Nyidove, n = 1, Amoussoukope, n = 4; Ountivou, n = 6; Towouganou (Zio), n = 11; Aveta, n = 2.

Figure 3.

Map of the spatial distribution of haplotypes of Python regius in Togo. Circle sizes correspond to haplotype frequency. Source of the map: https://www.esri.com.

Total genomic DNA was extracted from swab samples using the Blood and Tissue kit by Qiagen (Hilden, Germany) following the manufacturer’s protocol. DNA extracts are available from the ZFMK Biobank, Bonn. The COI segment (674 bp) with primers HCO2198-JJ and LCO1490-JJ (Astrin and Stüben 2008) was amplified via the polymerase chain reaction (PCR) in reaction volumes of 20 μl, including 2.5 μl of undiluted DNA template, and using the Multiplex PCR Master Mix (Qiagen). Thermal cycling was performed on GeneAmp PCR System 2700 instruments (Life Technologies, Carlsbad, USA) as follows: hot start Taq activation: 15 min at 95 °C; first cycle set, touch down (15 repeats): 35 s denaturation at 94 °C, 90 s annealing at 55 °C (−1 °C/cycle) and 90 s extension at 72 °C. Second cycle set (25 repeats): 35 s denaturation at 94 °C, 90 s annealing at 40 °C and 90 s extension at 72 °C; final elongation 10 min at 72 °C. Amplicons were purified with the ExoSAP-IT Purification Kit (USB Corporation, Cleveland, Ohio) and sequenced in both directions using the PCR primers at Macrogen Europe’s commercial Sanger sequencing service (Amsterdam, NL). All newly found haplotypes were deposited in GenBank (accession numbers: Suppl. material 1: Table S1).

Molecular data were first phylogenetically analysed based on mtDNA to gain insights into spatial pattern of genetic variation and the level of genetic divergence among populations of the species in Togo.

The number of haplotypes was calculated using DnaSP 6 (Rozas et al. 2017). A median-joining (Bandelt et al. 1999) and a median (Bandelt et al. 1995) haplotype network was constructed in PopART v. 1.7 (Leigh and Bryant 2015) and SplitsTree v. 4.14.8 (Huson and Bryant 2006), respectively, to visualise relationships among haplotypes of the 12 localities (hereafter referred to as “populations”). Haplotype frequency distribution was mapped using ArcGIS 10.3.1 (ESRI, Redlands, CA, U.S.A.).

We also tested for significant differences in nucleotide and haplotype diversity between sample locations using permutation tests implemented in the R script genetic_diversity_diffs v1.0.6 (Alexander et al. 2016; R Core Team 2019) with 10,000 iterations. The script resamples from the combined haplotype frequencies over all populations, in order to test whether the observed differences in haplotype diversity between specific populations, are greater than expected by chance. For this analysis, we grouped all farms together due to their small geographic distances, as we did with the neighbouring localities Agbave and Nyidove from which only one sample was available each.

Phylogenies were reconstructed by Maximum Likelihood (ML) and Bayesian Inference (BI) methods, using RAxML v. 8.2.10 (Stamatakis 2014) as well as MrBayes v. 3.2.7a (Ronquist et al. 2012). We used Xenopeltis unicolor (NCBI accession no. LC075326) as outgroup and included additional COI-sequences of P. regius available from GenBank (AB177878; KX012740; KX012757; KX012789–KX012791; NC007399). Optimal partitioning schemes and substitution models for each partition were selected using the Bayesian Information Criterion in a greedy search algorithm and linked branch lengths in PartitionFinder v. 1.1.1 (Lanfear et al. 2012, 2014). BI analyses were then run using models and codon partitions as selected by PartitionFinder. ML tree construction was performed under the GTR-GAMMA model applying the codon partition scheme preferred by PartitionFinder. Nodal support of the tree was assessed by rapid bootstrapping algorithm with 1,000 bootstrap replicates.

We initially tested the cross-amplification of microsatellite loci previously described by Jordan et al. (2002) and Shaney et al. (2016). Out of these microsatellites, 11 loci cross-amplified in our samples and were genotyped for each sample using fluorescently labelled primers (for accession numbers see Table 1). Some of the primers were modified in their length to adjust the annealing temperature for a multiplex reaction. Before multiplexing, all primer pairs were tested in singleplex PCRs to check for correct amplification of the desired fragments. PCRs were performed in a 12,5µl reaction mix containing 5-40 ng DNA, 200 µM dNTP´s, 0,5 µM each primer pair (reverse primer of each pair was labelled, see Table 1), 5× Q5 reaction buffer and 1 U Q5 High-Fidelity polymerase (New England Biolabs). Amplification conditions were as follows: 98 °C for 30 s followed by 33 cycles at 98 °C for 10 s, 64 °C annealing temperature for 30 s and 72 °C for 30 s, with a final extension at 72 °C for 2 min. Amplification products were separated on an ABI 3730 automatic sequencer. Raw data were analysed and evaluated using GeneMapper ID v3.2.1 (Applied Biosystems).

All markers were checked for scoring errors, large allele dropout and the possible presence of null alleles using Micro-checker v2.2.3 (van Oosterhout et al. 2004). Departure from Hardy–Weinberg equilibrium (HWE) was tested for each locus by Fisher’s exact tests (Guo and Thompson 1992) using a Markov chain Monte Carlo (MCMC) approach with 10,000 steps and 1,000 iterations as implemented in GenePop v. 4.2 (Raymond and Rousset 1995). GenePop was also used to test linkage disequilibrium (LD) for all pairs of loci. P-values of LD were corrected using the FDR-method (Benjamini and Hochberg 1995). The number of alleles, polymorphic information content (PIC) (Botstein et al. 1980), expected (He) and observed (Ho) heterozygosity were calculated per locus using Cervus v. 3.0.7 (Marshall et al. 1998; Kalinowski et al. 2007). FSTAT (Goudet 1995) was used to estimate and test (with 10,000 permutations) F_ST values between wild populations as well as between wild populations and farm individuals. GenAlEx (Peakall and Smouse 2006, 2012) was used to compute the probability of identity (PI) with increasing locus combinations, i.e., the probability that two randomly chosen individuals in a population have identical genotypes.

Isolation by distance (IBD) was examined for the “wild” populations with reduced major axis (RMA) regression and Mantel test on matrices of genetic and geographic distances using IBD v.1.52 (Bohonak 2002) and 10,000 randomisations. Pairwise geographic distances between localities were calculated with der Geographic Distance Matrix Generator (Ersts 2019) and log-transformed. Genetic distance matrix (FST/[1-FST]; Rousset 1997) was obtained with GenePop, based on the eight polymorphic loci.

Individual assignment tests were performed to assign farm individuals to the population they have the highest probability of belonging to using a Bayesian approach according to Rannala and Mountain (1997) implemented in the program GeneClass2 (Piry et al. 2004). This method has been shown to outperform other assignment techniques, and has been specifically demonstrated to obtain high assignment accuracy for smaller sample sizes (Cornuet et al. 1999). Moreover, this assignment method has the advantage that it does not assume that the true population of origin is among the sampled populations (Manel et al. 2002). We included the eight polymorphic loci only, assumed a threshold of 0.05 and applied a Monte Carlo resampling of 100,000 simulated individuals (Paetkau et al. 2004). For broad comparison, we also performed self-assignment tests for the individuals from wild populations using the same parameter setting. We considered individuals with an assignment score > 0.7 to be assigned with confidence (Dellicour et al. 2011).

The extraction of DNA from buccal swabs of 60 individual specimens from Togo proved successful; for two samples DNA extraction failed. The aligned sequence data set of these samples contained 674 bp with 13 variable characters of which 10 were parsimony informative. Translation of the gene segment revealed no frameshift mutations or premature stop codons. Both BI and ML trees show that samples are split into three major clades, which, however, show no geographic structure and are only weakly differentiated (Suppl. material 1: Fig. S1). These clades are consistently recovered in the median network (Suppl. material 1: Fig. S2). Sequence comparison revealed a total of nine haplotypes, including a haplotype found by Dong and Kumazawa (2005) (AB177878); the nine haplotypes were deposited in GenBank (MN295674–82). One haplotype (hap1) was shared by all localities except Nyidove (close to the Ghanaian border); three haplotypes (hap7, hap8 and hap9) were unique to their geographic population (Ountivou, near the Benin border; Towouganou [Zio district] and Aveta), and not represented among the samples from any farm (Fig. 3). Remaining haplotypes (which differed from the most common haplotype by only 1–4 mutations) were found at moderate to low frequency (Figs 3, 4). Five among nine haplotypes (based on 21 individuals) were represented among farm animals.

Figure 4.

Haplotype network based on 674 bp of COI gene from 60 specimens of Python regius from Togo. Each circle represents a haplotype with its size proportional to the frequency of the haplotype. Ticks on branches connecting the haplotypes indicate nucleotide mutations. Localities are indicated by different colour.

Consistent with the phylogenetic trees, the haplotype network does not inform on a specific geographic structure and pairwise comparison of genetic diversity between the localities revealed no significant differences (Fig. 4, Suppl. material 1: Table S2).

Seven out of 60 individual DNA samples could not be amplified for any of the microsatellite loci and were excluded from subsequent analyses. At least eight microsatellites were polymorphic and four of them showed higher polymorphism with PIC values > 0.7 (Table 2). None of the loci showed indications for the presence of null alleles, large allele dropout or stuttering. The sample set from the “Adaptation farm” showed a higher number of homozygotes at locus 5-MS16. However, because these farm samples may originate from different populations, we considered this excess as uncritical for subsequent analyses. Tests for LD between all pairs of loci across populations yielded no significant genotypic disequilibrium; none of the loci deviated significantly from HWE. For polymorphic loci, the mean number of alleles was seven (2–25 alleles per locus), the observed heterozygosity ranged from 0.02 to 0.91, and the expected heterozygosity from 0.01 to 0.93. The probability of two random animals having identical genotypes was estimated at 2.0 × 10^-8 (Table 2).

Information on population differentiation can ideally serve the management of genetic populations or, in our case, monitor national management regimes. The level of genetic differentiation (F_ST) between localities where wild individuals were sampled and specimens sampled from farms varied from 0 to 0.09 (Suppl. material 1: Table S3); none of the pairwise F_ST-values was significant, and there was no evidence of isolation by distance across populations (Z = -11588.18, r = 0.24, p ≤ 0.08). Notably, F_ST-values among pairs that included a farm were remarkably low, indicating a potentially higher level of individuals from different populations among the farm samples that minimise the degree of genetic differentiation between farms and wild populations.

Self-assignment of individuals from wild populations was successful and correct in 94% of the samples, with assignment probabilities ranging between 0.66 and 0.96 (mean 0.86); only two of these correctly allocated samples were assigned with lower confidence (0.66 − 0.70); (Suppl. material 1: Table S4). Sample MA26 (from Tsevié) was assigned with a slightly higher confidence to a different population than to the one where it was sampled from; sample MA50 (from Zio) couldn’t be assigned to any of the wild population. In “real life” assignment, two of the farm individuals (MA15, MA19) could be assigned with high confidence to one of the wild reference populations, while a further five samples (MA3, MA11, MA12, MA15, MA20) were assigned with a probability > 0.6 (Suppl. material 1: Table S5).

There are no significant geographical or climatic barriers in the sample area and gene flow between populations of ball pythons has likely occurred over the last 10,000 years. However, unlike our study, genetic studies focussed on wild populations of other savannah inhabitant reptile species in West Africa have revealed phylogenetically distinct clades [e.g., the African helmeted turtle (Pelomedusa subrufa) (Vargas-Ramírez et al. 2010; Wong et al. 2010), the Puff adder (Bitis arietans) (Barlow et al. 2013), and the Egyptian cobra (Naja haje) (Trape et al. 2009)]. Likewise, phylogenetically distinct clades have also been previously reported for ball python prey species [e.g., Misonne’s soft-furred mouse (Praomys misonnei)] (Nicolas et al. 2011).

Furthermore, stronger population subdivisions than those observed would be expected, at least across a wider geographic scale, given the low assumed dispersal capacity of ball pythons. Although an important research priority, to date field research provides some information on habitat use (Luiselli and Akani 1999), while studies on movements and home range sizes of ball pythons in West Africa are almost lacking. Little published information is currently available on daily and seasonal activities (Luiselli et al. 2007), however based on field experience and discussions with collectors in Benin, the species likely exhibits narrow home range (Toudonou, pers. comms.). Overall, it is suggested that the ball python represents a rather sedentary species (Luiselli, pers. comm.) that likely reflects other congeneric sedentary species, such as the Short-tail python (Python brongersmai), of which specimens have been shown to cover a mean distance of less than 200 metres in ca. 50 days (Erdelen et al. 1997, in: Auliya 2006).

There are other factors that can influence snake movement and dispersal, such as the seasonal flooding of python habitat during the wet season (Auliya 2006). Activity levels can also vary between savannah and forest populations (Luiselli and Angelici 1998), and males of congeners have reported higher activity levels when searching for females during the breeding season (Slip and Shine 1988). Although such aspects may have certainly influenced ball python population genetics in Togo, it is important to note that this species is reported to select specific microhabitats that provide numerous shelter opportunities and a constant prey source (Barker and Barker 2006).

Similarly, with regards to the ball python specimens sampled within python farms, haplotype diversity indicates that python farm ranching activity has historically targeted several populations but does not provide a clear spatial pattern or “trade chain” in this regard. The extent to which this observed lack of well-defined haplogroups is the result of regional trade activity or continuous historical gene flow/long-distance dispersal (or indeed vice versa) cannot be determined based on the mtDNA sequence data alone. However, it does raise a number of important questions regarding the impact of ongoing commercial trade activity on remaining wild populations.

In particular, it is currently unclear whether any of the haplotypes identified in Togo during this study actually originated and / or have current core distributions in neighbouring range States. This is a distinct possibility given that researchers have reported illegal cross-border hunting of ball pythons (de Buffrénil 1995) and other reptile species, such as the Müller’s sand boa (Eryx muelleri) (Vignoli et al. 2015).

Our initial genetic findings line up with results obtained from a recent questionnaire-based study focused on python hunters in Togo, which reported that the collection and subsequent release of ball python ranched specimens was carried out in a relatively ad hoc and diffuse manner, without an effective monitoring process in place (D’Cruze et al., in prep.). Specifically, this study raised concerns that ball pythons may be released in insufficient numbers, in inappropriate habitats and geographic locations that may be at least partly responsible for reported decreases in local wild ball python populations over the last five years (D’Cruze et al., in prep.).

It is beyond the scope of this study to assess the full impact that trade and associated ranching activity has had on the conservation status of wild ball python populations in Togo and other neighboring range states in West Africa. However, the challenges associated with implementing proper wild release protocols, and the multiple risks posed to focal species, their associated communities and ecosystem functions in both source and destination areas (including disease introduction and genetic pollution) are well known and should be mitigated (IUCN/SSC 2013).

In particular, ball pythons are known to harbour a number of infectious diseases (e.g., cryptosporidium, Yimming et al. 2016; nidoviruses, Uccellini et al. 2014) and parasites [including ticks (e.g., Amblyomma latum, Corn et al. 2011)]. For a recent review of the scientific literature (published between 2009 and 2019), which identified a list of nearly 150 underlying pathogens associated with this species, please see Green et al. (2020). The fact that ball pythons are being ranched under intensive commercial captive conditions prior to their release in large numbers (as has been reported in the CITES trade database), together with concerns regarding a lack of proper biosecurity measures at these facilities (D’Cruze et al., in prep.), arguably makes a higher level of monitoring and evaluation a key priority.

Moreover, the translocation of individuals of non-local origin may lead to introgression that disrupts spatial genetic structure, alters local genetic diversity, and ultimately threatens local adaptations (e.g., Söderquist et al. 2017 and references therein). Hybridisation between translocated farmed and wild individuals might also be a driver of biodiversity loss (Randi 2008; Sutherland et al. 2006; Söderquist et al. 2017).

Our study provides only an initial insight into ball python genetic diversity in Togo. Yet, regardless of the low sample size per population, more than 90% of the samples from wild captures of ball python were correctly assigned to their population of origin, while 10% of farm samples could be assigned with high probability to one of the wild populations (and ca. 1/3 with a probability > 0.6). This could indicate that farms harvest populations, which have not been sampled in our study. However, the relatively low percentage of individuals assigned is similar to those reported in other studies and the limited assignment accuracy may result mainly from population pairs with a F_ST < 0.05 (Manel et al. 2002).

We also acknowledge that identifying the full population status and or exact geographic region of origin of ranched ball pythons in the wild requires a much denser sampling from a wider geographical area (including populations that have been subject to harvest in Ghana and Benin, in addition to reference samples of non-harvested populations across the species’ range), and most importantly, a higher sample size per (sub-)population (i.e., 25–30 samples; Kalinowski 2005; Hale et al. 2012). This would enable us to identify geographic regions where substantial genetic admixture of farmed genotypes with wild populations may or may not have occurred. Moreover, when a larger sample size is applied, polymorphism may be revealed at the three additional loci (12-MS27, 3-KE959105, 14-KE961431) that were monomorphic in the samples used in this present study, which consequently may contribute “decisively” to forensic assignment.

It is important to note that larger than expected dispersal rates for this species may be partly responsible for the genetic admixtures reported in this study and more ecological studies are needed in this regard. However, based on current knowledge (cf. above), the ball python is likely a relatively sedentary species, assumed low dispersal capacity.