Road killed specimens or tissue of Duberria l. lutrix were obtained from the South African National Biodiversity Institute tissue bank (SANBI – Cape Town, South Africa) and from several private collections and road kills. Road killed specimens or tissue samples were preserved in absolute ethanol and refrigerated at 4 °C. A total of 87 D. l. lutrix specimens were collected from 38 localities across South Africa, covering most of the subspecies distribution range (Fig. 1, Table 1). An additional Duberria specimen from Uganda was donated by the Californian Academy of Sciences (CAS), USA. Specimens collected during the present study were deposited in the Port Elizabeth Museum Reptile Collection (PEM R), Eastern Cape Province, South Africa (Table 1).

DNA was extracted from ethanol preserved muscle or liver tissue biopsies. A MacheryNagel DNA extraction kit was used for the DNA extraction, following the manufacturer’s protocol. Extracted DNA was stored at -20 °C until needed for the polymerase chain reaction (PCR). Three gene regions were targeted using the PCR, these included two mitochondrial (mtDNA) loci: nicotinamide adenine dinucleotide dehydrogenase subunit 4 (ND4) using the primer pairs listed in Arévalo et al. (1994) and Barlow et al. (2013) respectively; while for cytochrome b (cyt b) the primer pairs listed in Burbrink et al. (2000) and Ruane et al. (2015) were used; for the nuclear locus, β-spectrin nonerythrocytic intron 1 (SPTBN1) the primers pairs listed in Ruane et al. (2014) were used (see Table 2 for details of the primer pair combinations). The ND4 and cyt b loci have been extensively used in snake phylogeographic studies, while the nuDNA locus has been demonstrated to be a variable nuclear marker in other snake species (Burbrink et al. 2000; Barlow et al. 2013; Ruane et al. 2014, 2015). All specimens were sequenced for the two mtDNA loci, while a single sample per locality was sequenced for the nuDNA locus.

All PCR amplification was performed using standard protocols. PCR conditions were as followed: 94 °C for 4 min., 94 °C for 30 sec., the annealing temperature of the primers varied between 48 °C to 50 °C for 35 sec., 72 °C for 40 sec., for 35 cycles and a final extension at 72 °C for 10 min. A list of the six primer pairs used are provided in Table 2. PCR products were visualized using a 1% agarose gel that contained a 1% ethidium bromide solution. Following successful amplification of a locus, a BioFLUX gel purification kit was used to purify the amplicons, following the manufacturer’s protocol. The gel purified PCR amplicons were sequenced at the Central Analytical Facility (CAF), at Stellenbosch University, using an ABI 3700 automated DNA sequencer.

The mtDNA and nuDNA sequences were aligned in CLUSTAL W (Larkin et al. 2007) using the default settings. The two protein-coding mtDNA (ND4 and cyt b) loci were examined for the presence of pseudogenes by converting the DNA sequence to amino acids to detect the presence of stop codons. Since all loci on the mtDNA were linked, the two mtDNA loci were combined for the phylogenetic analysis. For ND4 and cyt b, 740 bp and 610 bp fragment was respectively sequenced for the 87 D. l. lutrix specimens. Sequences were deposited in GenBank (Table 1). The combined mtDNA data set yielded a total of 1350 bp. The DNA substitution models obtained in jModelTest 2.0 (Darriba et al. 2012) are provided in Suppl. material 1. In addition, ND4 and cyt b sequences for D. variegata, as well as an additional specimen of D. lutrix. sp from Kenya were downloaded from GenBank and included in the phylogenetic analysis (Kelly et al. 2008; Vidal et al. 2008). For the SPTBN1, the older samples failed to amplify, hence these were coded as missing during the combined analyses. For the nuDNA SPTBN1, a 760 bp fragment was amplified for 30 D. l. lutrix specimens and sequences were deposited in GenBank (accession numbers in Table 1). For the nuclear SPTBN1 locus, allelic heterozygotes were inferred using PHASE (Stephens et al. 2001; Stephens and Scheet 2005). PHASE implements a Bayesian method for reconstructing haplotypes from nuclear sequences that include multiple heterozygous base sites within individuals. To estimate allele frequencies, PHASE was run five times. The run with the best goodness-of-fit to an approximate coalescent model was retained, resulting in two nuclear haplotype sequences of alleles per individual.

The combined mtDNA data set was subjected to a Bayesian inference (BI). The BI analysis was performed in MrBayes v3.2.6 (Ronquist et al. 2012; Huelsenbeck et al. 2016), sampling every 10,000 generations for a total of 10 million generations. Convergence and burn-in were determined using Tracer v1.6 (Rambaut et al. 2014). Nodes were considered well supported when they had a posterior probability (pP) of > 0.95. All trees were visualised using TreeGraph 2 (Stöver and Müller 2010). For the combined DNA analysis, a single specimen per locality was used of each of the three loci and the same phylogenetic approaches listed above were followed to reconstruct evolutionary relationships. In addition, we also performed a maximum parsimony (MP) analyses on the combined mtDNA data and the total evidence data sets. MP analyses were executed in PAUP*4 v. beta 10 (Swofford 2002). For the MP analyses, trees were generated using the heuristic search option with tree bisection and reconnection (TBR branch swapping using 100 random taxon additions). Phylogenetic confidence in the nodes recovered from MP was estimated by bootstrapping (Felsenstein 1985) analysing 1000 replicates of the data set. Only bootstrap values >75% were regarded as statistically well supported (Felsenstein 1985).

Both Vidal et al. (2008) and Kelly et al. (2009) demonstrated that the two snake species Amplorhinus multimaculatus and Ditypophis vivax are sister to Duberria. Hence the former two species were used as outgroups. Mitochondrial DNA sequences for the latter two species were downloaded from GenBank whereas the nuDNA locus was coded absent. Uncorrected sequence divergence values for both the ND4 and the cyt b loci were calculated in PAUP*4 version beta 10 (Swofford 2002). We did not combine the two mtDNA to calculate the sequence divergence value since we would not be able to compare this value with other phylogeographic studies since most authors perform this analysis on loci individually.

Haplotype networks were constructed for the combined mtDNA data, as well as for the nuclear SPTBN1 using TCS version 1.3 using a 95% connection limit (Clement et al. 2000). An analysis of molecular variance (AMOVA) was conducted in ARLEQUIN v3.5.2.2 (Excoffier and Lischer 2010) using the combined mtDNA data. The preliminary analyses of the combined mtDNA topology revealed the presence of five clades, hence a hierarchical AMOVA was performed; 1) across all sample localities and; 2) for haploclade one detected using the combined mtDNA data set. The remaining four clades were not analysed further in AMOVA since the low sample sizes limited any statistical inferences. Two neutrality test using Fu’s Fs (Fu 1997) and Tajima’s D (Tajima 1989) were calculated in ARLEQUIN using 10,000 permutations.

We used BEAST v1.8.3 (Drummond et al. 2012) on the mitochondrial data set to determine the age of divergent events within Duberria l. lutrix. No fossil calibrations points are available for Duberria. The genus Duberria belongs to the subfamily Pseudoxyrhophiinae which is part of the family Lamprophiidae. The Lamprophiidae is in turn sister to the Elapidae and both being sister to the superfamily Colubroidea. Hence published mutation rates from the superfamily were used in the divergence time estimation (Figueroa et al. 2016; Hsaing et al. 2015 For the two mtDNA loci (ND4 and cyt b) a strict substitution rate of 1.34% (SD=0.251) per million years was used (Daza et al. 2009) after having checked that the relaxed log-normal clock’s standard deviation approached zero. The mutation rate for the nuclear marker is unknown, hence this marker was excluded from the divergence time estimation. We ran the analyses under the substitution model as inferred above (TrN+Gamma), unlinked between both mitochondrial genes, and under a coalescent prior with constant population size, as the resulting tree described within species relationships. We ran the Markov chain for 50 million iterations and sampled every 10,000 iteration. Tracer v1.6 (Rambaut et al. 2014) was used to assess chain convergence to ensure minimal autocorrelation between iterations (effective sample size > 2000 for all sampled parameters) and to determine the burn in (10% of samples). TreeAnnotator in BEAST was used to determine a maximum clade credible tree.

The MP and BI analyses retrieved near identical tree topologies, hence only the BI topology is shown and discussed. For MP, of a total of 1350 characters, 278 characters were found to be parsimony informative, and recovered 167 trees with a tree length of 597 steps with a consistency index (CI) of 0.59 and a retention index (RI) of 0.88. The BI topology (Fig. 2) retrieved the two East African specimens (Kenya and Uganda) of Duberria as basal, while D. variegata appeared sister to a South African clade of D. l. lutrix, with low statistical support for the monophyly of the clade. Within D. l. lutrix, five geographically discrete, statistically well-supported clades were detected (>75%/>0.95 pP). Clade one consisted of specimens that occurred predominantly from above the Hottentots Holland Mountains, Agulhas plain and Overberg, the Cape Peninsula and the south-eastern Cape and adjacent interior, and was sister to clade two. Clade two comprised specimens from below the Hottentots Holland Mountains, including the Cape Peninsula and Boland region. Clade three was comprised of specimens from the Eastern Cape coast and samples from interior of KwaZulu-Natal Province, and was sister to clade four. Clade four comprised two specimens from Mpumalanga (Sabie), while clade five comprised two specimens exclusive to the Limpopo Province (Entabeni and Wolkberg).

Figure 2. 

A Bayesian inference topology derived from the combined mtDNA analyses (ND4 + cyt b) amongst the South African Duberria l. lutrix. The posterior probability value (pP) values are presented above each node. Only pP values >0.95 are shown. Values below each node are bootstrap values for MP. Only bootstrap values >75% are shown. An asterisk (*) below or above a node indicates the absence of statistical support. The insert shows the typical Duberria lutrix lutrix.

For the ND4 locus the maximum uncorrected sequence divergence between the first and second clades was 1.68%. Between the second and third clade the maximum uncorrected sequence divergence was 3.84%. Between the third and the fourth clades the maximum uncorrected sequence divergence was 5.95%. Finally, the maximum uncorrected sequence divergence between clades four and five was 6.75% (Suppl. material 1: Table S1). For cyt b locus the maximum uncorrected sequence divergence between the first and second clades was 0.98%. Between the second and third clade the maximum uncorrected sequence divergence was 2.62%. Between the third and the fourth clades the maximum uncorrected sequence divergence was 3.77%. The maximum uncorrected sequence divergence between clades four and five was 6.01% (Suppl. material 2: Table S2). The ND4 locus revealed slightly higher levels of uncorrected sequence divergence values during the present study in comparison to the cyt b locus.

A total of 35 haplotypes were retrieved for the 87 Duberria l. lutrix specimens using the combined mtDNA (Fig. 3). For details of the haplotype distribution consult the Suppl. material 1: Table S3. Five haplogroups were retrieved, revealing a pattern congruent with the combined mtDNA topology (Fig. 2). The AMOVA results among all 38 sample localities revealed that 94.26% (Va = 15.89; df = 37; SS = 1394.37) of the variation occurred among sample sites, while 5.74% (Vb = 0.96; df = 49; SS = 47.44) of the variation occurred within sample sites. These results are indicative of marked genetic differentiation, a result that is corroborated by the marked Φst (0.94) as well as the high F_ST values among sample localities that were statistically significant for 56 combinations, ranging from 0.05 to 0.99 (results not shown). Haploclade one (Fig. 3) consisted of 20 haplotypes comprising samples from above the Hottentots Holland Mountains, Agulhas plain and the Cape Peninsula separated by ten unsampled / missing mutations from samples from the south-eastern Western Cape, with 76.46% of the variation occurring among sample localities, (Va = 3.07; df = 22, SS = 167.63), 23.58% occurred within sample localities (Vb = 0.94; df = 29; SS = 27.53), with a high Φst (0.74) as well as high F_ST values among sample localities that were generally statistically significant. The remaining four haploclades, two, three, four and five respectively (corresponding to the identical clades observed in Fig. 2) contained fewer than five haplotypes, hence we did not undertake any further statistical analyses due to the low sample sizes.

Figure 3. 

A minimum spanning network derived from combined mtDNA (ND 4 + cyt b) analyses demonstrating the five haploclades for Duberria l. lutrix. The number inside the boxes correspond to the haplotypes in Table 3. The closed black dots represent unsampled or missing haplotypes.

Fu’s Fs values were positive for Agulhas, Kokstad, Pringle Bay and Swellendam and indicate an excess of intermediate polymorphisms due to recent population bottlenecks or balancing selection, however none of these were statistically significant. Two of the Fu’s Fs values were negative for Somerset West and Napier indicating an excess of low frequency polymorphisms consistent with population expansions or positive directional selection. However, only the Somerset West population was statistically significant (P < 0.02). The remaining sample localities had a Fu’s Fs of zero. For Tajima’s D, five sample localities, Swellendam, Pringle Bay, Napier, Somerset West and Kokstad were negative, while only Napier and Kokstad were statistically significant (P < 0.05). The remaining sample localities has a Tajima’s D of zero. Negative Tajima’s D indicates an excess of low frequency polymorphism, population expansion or purifying selection.

Ten haplotypes were retrieved for the 30 specimens using the TCS analyses (network not shown). For a list of the sample localities per haplotype consult Suppl. material 2. Three haploclades were retrieved. Haploclade one contained six haplotypes from all the remaining sample localities (from clades 1, 2, 3 and 4; Fig. 2). Haploclade two contained a single haplotype from Oudtshoorn (clade 1; Fig. 2). Haploclade three contained three haplotypes from localities in the KwaZulu-Natal (Kokstad) and the Limpopo provinces (Wolkberg and Entabeni) (from clades five and three respectively; Fig. 2).

The combined mtDNA and nuDNA sequence data yielded a total of 2110 bp. The MP and the BI analyses retrieved highly congruent tree topologies, hence only the BI topology is shown. For the MP analyses, 293 characters were found to be parsimony informative, tree length of 627 steps, 730 trees, with CI = 0.58 and RI = 0.77. The total evidence BI topology (Fig. 4) was congruent with the combined mtDNA topology (Fig. 2). A monophyletic Duberria was retrieved with the two East African Duberria specimens (D. l. atriventris) form a basal split, while D. variegata was sister to a clade containing all the D. l. lutrix samples from South Africa. These five clades were also evident in the combined mtDNA topology (Fig. 2).

Figure 4. 

A Bayesian inference topology for the total evidence data sets (ND4 + cyt b + SPTBN1) amongst the South African Duberria l. lutrix species complex. The posterior probability value (pP) values are presented above each node. Only pP values >0.95 are shown. Values below each node are bootstrap values for MP. Only bootstrap values >75% are shown. An asterisk (*) below or above a node indicates the absence of statistical support.

Clade five diverged from the remaining D. l. lutrix clades 3.42 Mya (4.13 Mya to 2.72 Mya; 95% HPD). These divergences dates fall into the Pliocene and Pleistocene epochs. Clade four diverged from clades one, two and three occurred 2.05 Mya (1.34 Mya to 2.95 Mya; 95% HPD). The divergence between clade three from clades one and two occurred 1.23 Mya (0.74 Mya and 1.81 Mya; 95% HPD). Cladogenesis between clades one and two occurred 0.46 Mya (0.27 Mya and 0.69 Mya; 95% HPD) (Suppl. material 3).

The phylogeographic results demonstrated the presence of five mtDNA clades across the sampled distribution range of Duberria l. lutrix in South Africa, implying that the snake represents a species complex. Furthermore, these five clades were characterised by the absence of shared mtDNA haplotypes and marked sequence divergences values for both the ND4 and cyt b loci, suggesting possible genetic isolation and limited dispersal. However, the nuclear DNA sequences data failed to retrieve patterns congruent with the mtDNA data, hence it is not possible to exclusively imply the presence of five possible taxa within the Duberria l. lutrix species complex, and requires further taxonomic delineation. In addition, the divergence time estimates suggest that cladogenesis in the D. l. lutrix species complex occurred during the Plio/Pleistocene epochs, a period that was characterised by increased aridification throughout South Africa, resulting in the contraction of mesic adapted species. Our results reflect the impact of Plio/Pleistocene climatic driven fragmentation on a snake species resulting in possible cladogenesis. The latter result is in line with what has been reported for the widely distributed puff adder in South Africa (Barlow et al. 2013).

During the late Miocene, climatic profiles changed dramatically, resulting in decreasing levels of precipitation and marked aridification, a trend that was enhanced in the Plio/Pleistocene (Siesser 1980; Marlow et al. 2000; Cowling et al. 2009; Engelbrecht et al. 2013). Furthermore, during the Plio/Pleistocene epochs, coastal regions experienced dramatic marine transgressions that have been estimated to vary between 150–200 meters in certain regions (Partridge and Maud 1987, 2000; Partridge 1997). These events likely resulted in the extinction of low-lying terrestrial taxa characterised by low vagility and habitat specificity and the contraction of animals to high-lying refugial areas along the coastal belt mountains and the adjacent interior. More recently, during periods of glacial maxima in the Holocene, the interior of South Africa is thought to have become more inhospitable to ectotherms, due to low winter temperatures and reduced precipitation levels. These factors possibly caused mesic adapted organisms such as D. l. lutrix to seek more favourable habitat along the high-lying mountainous coastal regions of the Cape Fold and Drakensberg Mountains (Barlow et al. 2013; Tolley et al. 2014). During the Last Glacial Maximum coastal areas would have had exposed areas of continental shelf, off the current south-west and western coasts of South Africa, due to the lowering sea levels. This would have provided favourable habitat, as well as, acting as dispersal corridors for many species (Schreiner et al. 2013). During the interglacial period these corridors would have been inaccessible due to the rising of sea levels. Rapid changes in elevation can provide significant biogeographic barriers to the dispersal of ectotherms; this is evident when one compares the geographic topology of clades one and three. Clade one extends from the western coast of the Cape Peninsula above the Hottentots Holland Mountains range until just before the south-eastern Cape coastline whilst clade three occurs below the Hottentots Holland mountain range extending into the Cape Flats. The rapid changes in elevation between the two clades limits the dispersal of D. l. lutrix. Similar phylogeographic breaks have been observed in other co-distributed reptile species (Daniels et al. 2007, 2009; Tolley et al. 2009, 2014; Barlow et al. 2013). Climatic fluctuations would have altered environmental conditions during the Plio/Pleistocene allowing for dispersal of D. l. lutrix around these mountain ranges due to the changes in sea levels. Further evidence for isolation induced by climatic conditions during the Plio/Pleistocene can be found between clades one and two. Clade one occurs predominately throughout the western half of the Greater Cape Floristic Region (GCFR), namely the Cape Peninsula, above the Hottentots Holland mountain ranges and throughout the Agulhas plains and Klein Karoo, whilst clade two occurs within the eastern half of the GCFR, along the south-eastern Cape coastline and interior. Potentially, the reason for the genetic isolation between the two clades may be due to changes in the climatic conditions across the GCFR. The western and eastern sections of the GCFR are characterised by distinct rainfall regimes with the western half being characterised by a winter rainfall, while the eastern half is characterised by aseasonal and/or summer rainfall regime (Siesser 1980; Cowling et al. 2009; Tolley et al. 2009). It has been hypothesised that the division between some clades of reptiles corresponds to the changes in the rainfall regimes (Tolley et al. 2009, 2014). This is further corroborated by clades two and four which are found in the Eastern Cape, KwaZulu-Natal and Mpumalanga provinces. The Bedford gap is situated between the two clades. However, it is uncertain whether the genetic isolation is due to a combination of the xeric conditions and changes in rainfall patterning or simply due to one of the two variables. As rainfall patterns change and environments become more xeric, the minimum annual temperature of the area decreases, which can potentially limit the dispersal capabilities, as well as, survival capabilities of ectotherms. Furthermore, this would explain why D. l. lutrix has not dispersed along the western coastline of South Africa where the average annual precipitation and minimum annual temperatures are lower. This can be further evaluated by observing the effects that the Breede River xeric corridor had on the phylogenetic patterning of D. l. lutrix. When examining the topology of the phylogenetic trees (Figs 2, 4) the Breede River xeric corridor did not display any pronounced impact on the phylogeographic patterning. This observation favours the change in rainfall patterning as a possible source for the genetic isolation observed within D. l. lutrix. However, climatic fluctuations during interglacial periods would have lessened the impact that this xeric valley had on the dispersal of this species as precipitation levels changed. Finally, the extinction of intermediary haplotypes, possibly due to the climatic fluctuations, in widely distributed species with low dispersal capabilities and gene flow may have resulted in these pronounced phylogenetic gaps. Widespread sampling of D. l. lutrix is required to affirm or reject these inferences.

Although the climatic fluctuation would have forced species to retreat into refugia (Barlow et al. 2013; Tolley et al. 2014), evidence for the relictual populations, which are proposed to be restricted to the interior in the Klein Karoo and Free State are not corroborated by our results. The haplotype network for the nuclear marker SPTBN1 showed the Klein Karoo locality to be isolated from the other localities in the western and south-eastern Cape coastline and interior. However, this can be potentially biased; as firstly, the SPTBN1 is a protein coding locus and secondly, the mutation rates for nuclear markers in ectotherms are much slower than the mitochondrial markers, limiting inferences derived from them.

The combined mitochondrial and nuclear data set retrieved five clades that show evidence for geographically distinct Duberria lutrix lutrix lineages. The mtDNA data shows, high levels of uncorrected sequence divergence values. Glaw et al. (2007) reported that within the Malagasy Geodipsas infralineata using the cyt b marker uncorrected sequence divergence between the two clades ranged between 4.7–4.8%. Similarly, Ruane et al. (2017) using morphology and DNA sequencing (of the cytochrome oxidase one locus) observed two clades within the widespread Malagasy snake Mimophis mahfalensis. It is noteworthy, that while the phylogenetic affinities within the Lamprophiidae has recently received attention (Vidal et al. 2008; Portillo et al. 2018), phylogeographic studies remain limited. The latter observation suggests that species diversity among widespread species in the family may have resulted in an underestimation of alpha taxonomic diversity.

These sequence divergence values are similar to values observed in other snake lineages within Colubroidea and the Malagasy Pseudoxyrhophiinae that are considered genetically different (Feldman and Spicer 2002; Gou et al. 2011; Kindler et al. 2013; Ruane et al. 2015). This indicates that the respective clades might potentially be composed of cryptic species; however, we are cautious of the pitfalls of using these divergence estimations as exclusive evidence for species boundaries. Our results partially support the observation by Branch (2002, 2003), that the clade to which specimens from Port St Johns belong may represent a cryptic lineage. However, considering the sparse sampling of D. l. lutrix during the present study this might be an underestimation of species diversity. We advocate larger sample sizes per sample locality and a more comprehensive geographic sampling of the species range throughout South Africa coupled with more sensitive nuclear DNA markers, such as microsatellites or single nucleotide polymorphisms (SNP’s) to examine patterns of biparental gene flow. However, it is frequently difficult to obtain large sample sizes for snakes, specifically for phylogeographic studies. Similar observation has been made in other studies on snakes (Martínez-Freiría et al. 2015).