Bothrops asper exhibits robust genetic partitioning that accounts for at least ten distinct mitochondrial phylogroups, included in two separate lineages. The groups occupy different geographic regions and show private haplotypes, indicating a clear division among them. We found evidence of transition zones only in Peten (Guatemala) and in central Panama, but in both cases, contact occurs between sibling groups within a single B. asper lineage (MY/NCA and PICA/CHOCO1, respectively). Conversely, we did no observe mitochondrial admixture between B. asper lineages, nor even in Isthmian Central America where they converge. Additional sampling in this region, especially in western Panama, could help to establish the extent to which these lineages are introgressing and the precise geographical boundaries in this apparent admixture zone.

Phylogeographic history

The observed molecular variation might, in part, result from the complicated geological and climatic history that shaped the distribution of the species during the Pliocene-Pleistocene, in addition to the species own ability to spread into the new empty niches that arose during that period.

As stated before, the sister relationship between B. asper and the B. atrox group has been established previously [27, 34, 65]; and a monophyletic B. asper was recovered in our BEAST estimations of divergence time and the putative ancestral distributions of our clades. Our estimates of divergence time suggest that the B. asper-B. atrox group ancestor was distributed in South America no older than 4.48 Mya (Fig 4). Also, BPA analysis resulted in Colombia as that the most probable 'state of locality' for the deep branches in the phylogeny of B. asper (Fig 3). We interpret these combined results as evidence that the diversification of B. asper occurred in the northwestern region of South America by the early Pliocene.

The allopatric distribution of B. atrox group and B. asper lineages on each side of the Eastern Andes Cordillera, suggests that the final uplift of this mountain range played a significant role in the cladogenesis of these lanceheads. Thus, this Andean Cordillera reached no more than 40% of its present elevation during much of the Neogene, but intense mountain building followed in the late Miocene and, especially in the Pliocene when elevations increased rapidly [70–71]. As a result, we hypothesize that a B. asper stock diverged to the west of the Eastern Andes, likely within the foothills along the Pacific coast of northern South America, as supported by the relatively high number of haplotypes observed in that region.

During that time, this ancient stock was further isolated by several marine incursions in northern South America, especially at the Magdalena Valley and the Maracaibo basin to the west, and the presence of three sea corridors in lower Central America to the north: the Atrato seaway, the Panama Portal, and the San Carlos Basin [72]. These marine transgressions resulted from changes in sea level that occurred during the warm climate periods of the Pliocene [73–75], but also during the Quaternary interglacial periods when warmer and wetter climates dominated again. Isolation during the cycles of sea incursions allowed rapid species evolution in the trans-Andean region, as has been recognized for several organisms by Nores [76], and is feasible that similar forces managed the divergence within B. asper.

Considering its South American origin, and given the low posterior probabilities for the ancestral locality that we retrieved for the root of Clade A (Fig 3), we postulated that the early separation between lineages in B. asper was driven by a dispersal event north into Caribbean Mesoamerica. In addition, the divergence between lineages in mid-Pliocene coincides with estimations of the final closure of Isthmian Central America, a dynamic event that, as previously mentioned, has shaped the region’s biogeography.

There is still considerable debate regarding the dynamics and timing of the final closure of Isthmian Central America. Some authors argue that the Isthmus arose as a series of islands in a shallow sea with a concluding land connection established in the late Pliocene [9, 10]; whereas others suggest that it emerged as a continuous land bridge since the Miocene [11, 72, 77], only broken by the three previously mentioned sea gates that eventually followed a north-to-south closure at the end of the Neogene [78]. Regardless of the paleogeographic model of emergence, colonization of Caribbean Mesoamerica by B. asper may have entailed crossing water corridors, a challenge that lancehead pitvipers seem suitable to carry out [65]. In support of this view, present-day surface currents in the southern Caribbean are mostly directed toward the northwest [79], and passive transportation from the north Caribbean coast of South America to northern Mesoamerica has been demonstrated experimentally [80].

According to our divergence time estimations, the separation within the South American clade B started in late Pliocene, when a group reached the Pacific coast of Isthmian Central America (PICA). The route involved in this expansion to Mesoamerica is unclear, but most likely it also involved passive transportation through the Atrato seaway when global climate was warm [81]. By then, the Talamanca mountain range in Isthmian Central America isolated this population from the Caribbean lowlands [82, 83]. This hypothesis explains why Caribbean and Pacific B. asper populations in Mesoamerica are not sister taxa. Close affinities between lineages along the Pacific slope of Isthmian Central America and those in the Choco region in Colombia have also been reported in other taxa: frogs [84, 85], birds [86] and snakes [23, 87], suggesting a comparable history of colonization.

Hence, we postulate that B. asper invaded Mesoamerica in at least two, chronologically separated independent events. This trend coincides with the two-step dispersal-pulse hypothesis proposed by Savage [88] to explain the unequal distribution of amphibians and reptiles of South American origin inhabiting Middle America. According to this author, the first and earlier episode took place some 3.4 MYA when the sea level lowered [89], and a handful of South American taxa were able to invade Mesoamerica before the final closure of the isthmus. Consequently, few South American taxa have distributions that extend as far north as Mexico. Savage’s [88] second dispersion episode presumably occurred almost a million years later, after the final closure of the Isthmian bridge. As an outcome of this late invasion, the majority of genera of South American origin have distributions that do not reach beyond southern Nicaragua or Costa Rica. Chronologically, Savage’s first pulse corresponds to the split of B. asper lineages, whereas his second pulse matches the invasion of B. asper from South America to PICA.

Alternatively, a single dispersal event into Middle America is also possible, especially if incomplete lineage sorting [90] or secondary (interspecific) gene flow is contemplated. Both are distinct phenomena but produce very similar patterns of shared genetic diversity that in turn could affect species integrity.

During range expansion, some alleles could increase their frequencies due to gene drift, thus, “surfing” into fixed spatial sectors at the expanding front. Lineages that arise as a consequence of demographic expansion are temporary but tend to last longer in species with limited dispersal abilities [91]. Genetic surfing has been shown to explain patterns of spatial assorting of gene variation in loci with weak effective dispersal, such as mitochondrial DNA. Thus, genetic surfing should be more frequent in species that exhibit male-biased dispersal, and that have relatively low vagility, as it has been recently reported in the coral snake Micrurus tener [92]. Although the dispersal patterns of B. asper are unknown, male-biased dispersal is expected in this species, as it has been described in several other species of vipers, comprising both crotalines [93] and viperines [94].

Our mismatch distribution analyses for demographic expansion provide evidence that MV and CCO groups represent recent populations that diverged in the early Pleistocene. Thus, it is possible that the sorting of these haplotype groups resulted from gene surfing and not due to divergence in allopatry. However, no further evidence for demographic expansion was recovered in other populations, not even in those located at the northern extreme of the species distribution. Furthermore, we did not observe divergent mitochondrial haplotypes in sympatry in the putative ancestral populations, which is a crucial component to evidence of gene surfing, therefore limiting our interpretation. Further sampling of these populations and the inclusion of genome-wide markers in the analyses could resolve the question of whether the spatial sorting of haplotype groups in B. asper resulted from allopatric divergence or genetic surfing during demographic expansion.

Another possible scenario is that mtDNA introgression from B. atrox into B. asper occurred in the past, as this could explain the paraphyletic relation of this last species retrieved in our analyses (Fig 2). Introgression between sibling species is more widespread than previously thought and is reported in a great variety of plant and animal taxa [95, 96] including snakes [97, 98]. Ancient introgression could occur even if no evidence of shared haplotypes is available, as has been recently reported by Ruane et al. [98] for milk snakes genus Lampropeltis. For species that exhibit male-biased dispersal, as expected in B. asper, rates of introgression in mtDNA markers are often higher compared with those from biparentally inhered nuclear DNA markers [95]. Unfortunately, we did not include nuclear markers in our analysis, and none of our mitochondrial haplotypes were shared among our phylogroups, thus precluding further comparisons to evaluate introgression.