A total of 24 tissue samples were obtained from wild Kinosternon scorpioides turtles in six sampling locations in three Colombian regions: one in the Colombian Amazon Basin (Leticia), two in the Caribbean Coast or Caribe Basin (Cispatá Bay and Lorica swamp), three in the Orinoco basin (sierra de la Macarena-Guaviare Sub-basin, Puerto Carreño-Orinoco river floodplain and Tuparro National Park-Tomo Sub-basin) and one in San Andrés Island ( Figure 1 and Table 1 ).

Turtles were captured at night. Then, following the protocol designed by the Instituto de Investigación de Recursos Biológicos “Alexander von Humboldt” in Colombia (Vargas-Ramírez, 2017), approximately 0.3 cm³ of tissue was cut out from the back foot using a scalpel and stored in 90% ethanol. An iodine solution was applied after sampling to prevent infection. Sampling protocols were evaluated by the animal care committee (Cicual) from Universidad de los Andes.

We implemented a calibrated genealogy reconstruction analysis using the R35 Intron and the 35,507 SNPs matrix in a combined way, considering each dataset as one partition but with linked trees. For this purpose, we ran a Bayesian phylogenetic analysis in BEAST v1.8.2 (Drummond et al., 2012) using GTR+G as the best-fit evolutionary model inferred in jModelTest (Posada, 2008), an ‘log-normal relaxed clock model’ [substitution rate of 2.7x10-4 mutations per site per million year (Lourenço et al., 2013)] and a ‘Coalescent exponential growth’ as the tree prior. The remaining parameters were left as the default settings. We ran two independent runs of 50 million generations each with 10% burn-in, and convergence was checked in TRACER v1.7 (Rambaut et al., 2018). Next, log and tree files from each independent run were combined to produce a final dataset of 10,000 topologies by using LOGCOMBINER v2.5.0 (Bouckaert et al., 2019); this dataset was subsequently used to produce a maximum credibility tree in TREEANNOTATOR v2.5.0 (Bouckaert et al., 2019). Finally, the obtained ultrametric topology was visualized in FigTree v1.4.0 (Rambaut and Drummond, 2012).

Since the establishment/colonization age of a populations in a territory is directly related with their current levels of diversity and adaptive potential (Nei et al., 1975), consequently definning its priority for conservation, we implemented a spatiotemporal reconstruction of the ancestral area, dispersion routes (biogeographical portals) and colonization patterns of K. scorpioides between island and continent. A phylogeographic analysis using a Bayesian approach and the Relaxed Random Walk model, was run on BEAST v.1.8.2 (Drummond et al., 2012), in the CIPRES platform (Miller et al., 2010). We used the same dataset (intron + SNPs) for this approach, as described above. This method infers an evolutionary history in a continuum landscape based on the DNA sequences and the geographic coordinates (latitude and longitude of each sample), generating a genealogy and estimating the ancestral location of the internal nodes, considering uncertainties of the topology (Lemey et al., 2010). The substitution model selected was GTR+G, previously obtained from JModelTest (Posada, 2008). No ambiguities were used, since allele phasing has minimal impact on phylogenetic reconstruction from targeted nuclear genomic sequences (Kates et al., 2018). They were linked sites model and for the clock model, but not for the trees among loci. We used the uncorrelated log-normal relaxed clock model and coalescent exponential growth for the tree prior. For the temporal calibration of the phylogenetic tree, we used a substitution rate of 2.7x10-4 mutations per site per million year (Lourenço et al., 2013). Two independent runs where done, with a chain length of 50 million, sampled every 5,000 steps and a burn-in of 10%. Chain convergence, stationary behavior of the chains and Effective Sample Size (ESS) > 300, were examined in TRACER v.1.7.0 (Rambaut et al., 2018). Next, log and tree files from each independent run were combined to produce a final dataset of 10,000 topologies by using LOGCOMBINER v2.5.0 (Bouckaert et al., 2019); this dataset was subsequently used to produce a maximum credibility tree in TREEANNOTATOR v2.5.0 (Bouckaert et al., 2019). To generate spatiotemporal projections of the genealogies in keyhole markup language (kml), the software SpreaD3 v0.9.6 (Bielejec et al., 2016) was used, results that then were accessed and visualized on Google Maps (https://www.google.com/intl/es/earth/).

Considering that ESUs are defined as populations or groups of populations demonstrating significant divergence in allele frequencies (Dizon et al., 1992) or that are substantially reproductively isolated from other conspecific population units (Waples, 1991); several approaches were executed to test these principles using the SNPs dataset. First, we used STRUCTURE v.2.3.4 (Pritchard et al., 2000) to infer population genetic structure and admixture patterns among populations to test reproductive isolation, allowing us to examine the existence of biological groups within our sample. We conducted 20 independent runs for K = 1–6 by using 1,000,000 MCMC steps, a burn-in of 100,000 iterations and considering an admixture model and correlated frequencies. The α values and profile of posterior probabilities were examined to assess the convergence between independent runs. STRUCTURE Harvester (Earl and vonHoldt, 2012) was used to extract the Q values from each run for each K value, and then summarized in the program CLUMPP 1.1.2 (Jakobsson and Rosenberg, 2007). Finally, DISTRUCT v.1.1 (Rosenberg, 2004) was used to plot and visualize the obtained Q-matrix containing the ancestry probability of each individual in each predefined population. The ‘Evanno’ method (Evanno et al., 2005) and the Puechmille method (Puechmaille, 2016) were used to infer the most likely number of biological groups, using the web-based software StructureSelector (Li and Liu, 2018).

Patterns of genetic structure were further explored, performing a principal component analysis (PCA) and a discriminant analysis of principal component (DAPC) using the R- package ADEGENET (Jombart, 2008). These are dimensionality reduction methods used to transform and project data points onto fewer orthogonal axes that can explain the greatest amount of variance in the genetic data maximizing (DAPC) or not (PCA) the differences, helping with ESUs delimitation.

In addition, the genetic structure among geographical samples was calculated by the standardized pairwise F_ST (Meirmans, 2006 ), with 10,000 permutations and bootstraps, included in the Genepop v.1.1.7 R- package: “Population Genetic Data Analysis Using Genepop” (https://www.r-project.org, http://kimura.univ-montp2.fr/~rousset/Genepop.htm). When signals of genetic admixture between populations were observed through STRUCTURE or PCA/DAPC analyses, we quantified the gene flow in terms of the effective number of migrants per generation (Nm) between those populations (Muniz et al., 2018), using the program Migrate v3.6.11 (Beerli, 2006). For this purpose, we implemented the default model of free variation (*) for Theta and M parameters between populations, using exponential priors. The search strategy used static heating scheme with four chains at different temperatures (1.0, 1.5, 3.0, 1000000.0) and four replicates of one long chain with 50 million generations, discarding a burn-in of 5 million and sampling every 500 generations. To verify the convergence of parameter estimates, we examined the autocorrelation (<0.3), the effective sample size (ESS>500) and the correct distribution of the probability density plots for estimates of the parameters. To convert the scaled parameter “M” to biological data in terms of number of effective migrants (Nm), we assumed a conservative generation time of 28.2 years, estimated for Kinosternon members (Iverson et al., 2013); and a substitution rate for the SNPs of 2.7 × 10-4 mutations per site per million year (Lourenço et al., 2013).

In addition, to corroborate current gene flow among populations breaking the consistently congruent genealogical phylogenies (Avise, 1994), we analyzed the 10,000 genealogies from the above phylogenetic reconstruction, in DensiTree (Bouckaert et al., 2019), showing all variation over the Bayesian probability density of topologies among individuals from all populations.

Our results support the presence of two currently accepted subspecies of K. scorpioides in Colombia, one (K. s. albogulare) on San Andrés island and a second one (K. s. scorpioides) on the continent (Morales-Betancourt et al., 2015). However, in the calibrated phylogenetic reconstruction, STRUCTURE, PCA, DAPC and F_ST analyses, three additional differentiated groups were detected within K. s. scorpioides. In this analysis, clear differentiation and substantial reproductive isolation was found among the Caribbean coast (northwest, trans-Andean), Amazon and Orinoco regions (east, cis-Andean).

Configuration and diversification age for K. scorpioides populations in this study, occurred from ~ 8 to 1.9 Ma (late Miocene-Pliocene-Pleistocene), coinciding with the greatest orogenic events described by Hoorn and Wesselingh (2011), Hoorn et al. (2010) and Albert et al. (2018) between ~10 - 2 Ma, shaping the diversity in the Neotropics: the uplift of the eastern Andes (10 – 2 Ma), the Panamá Isthmus/arch submergence and emergence (9 - 3 Ma) and Vaupés Arch final formation (8 - 5 Ma).

As expected, the temporal origin, distribution and establishment of those differentiated populations through Colombia, inferred in the Relaxed Random Walk reconstruction analysis, suggested the ancestral area of Kinosternon scorpioides in a trans-Andean region in the northern South American continent (>8 Ma), with later colonization of the insular region (San Andres) about 1.72 Ma. There are at least two hypotheses that try to explain the arrival of K. S. albogulare to the island of San Andrés. One hypothesis suggests a natural origin of the population, via colonization of the island from the Central American continent in times of low sea level during the Miocene (Iturralde-Vinent, 1999; Vargas-Cuervo, 2004; Kirby et al., 2008). This hypothesis is supported by our data. However, an alternative hypothesis exists, suggesting the introduction of these turtles from the Central American populations by indigenous communities and traders as recently as the 1600s (Dunn and Saxe, 1950 in: Forero-Medina and Castaño-Mora, 2011).

According to O’Dea et al. (2016), between 9 and 6 Ma, rates of paleobathymetric change in sedimentary sequences reveal significant deepening across the Panama Arch, with posterior formation of the Isthmus of Panama sensu stricto around 2.8 Ma, allowing for faunal dispersion between South America and Central America (including its insular system). This was accompanied by a progressive and drastic decrease of eustatic sea-levels from 20 to 60 m (~2 – 1 Ma), coinciding with the establishing of K. scorpiodes in San Andres, and later in the Caribbean coast (1.04 Ma). Such results agree with the findings by Iverson et al. (2013), supporting a Mesoamerican origin of the family Kinosternidae, with the divergence of K. s. albogulare around ~ 8 - 4 Ma (in the late Miocene) and a later divergence of K. s. scorpioides around 6 – 2.5 Ma in the Pliocene.

However, an unexpected colonization pattern of cis-Andean populations was evidenced, mediated by more than one ancestor colonizing in independent events the north and south of cis-Andean populations (non-monophyletic pattern), through two different biogeographical portals. The first event occurred through Las Cruces Pass (Northern Central Andes) with later establishment in Amazon region; while the second one occurred through Mérida Cordillera at Táchira Depression (Northern Eastern Andes), with later expansion to the eastern lowlands of the Orinoco region. Both Las Cruces Pass and Táchira Depression, have been identified as biogeographical portals that, when closed, have generated allopatric speciation between trans- and cis-Andean ecoregions (Hazzi et al., 2018); as we are observing here for Kinosternon.

In fact, the differentiation of K. s. scorpioides in continental Colombia, giving origin to Orinoco cis-Andean population and Caribbean trans-Andean populations, could be related to vicariance events due to the uplift of the northern Andes and the Mérida Cordillera in the late Pliocene (5 - 3.4 Ma) (Diaz de Gamero, 1996; Audemard and Audemard, 2002). Most studies addressing this topic have identified the uplift of the northern Andes as a major driver for Neotropical diversification in a scenario consistent with allopatric speciation, wherein the complex topography of the Andes isolated populations on both sides of this barrier thus restricting gene flow (Brumfield and Capparella, 1996; Antonelli et al., 2009; Hoorn et al., 2010). Similar patterns of cis-Andean and trans-Andean population divergence resulting in speciation has been found in many vertebrates, including birds (Sedano and Burns, 2010; Arbeláez-Cortés, 2020), mammals (Patterson, 2020) and amphibians (Pirani et al., 2020).

However, we found divergence times shorter than 5-3.4 Ma for Caribbean Vs. Orinoco populations (1.9 Ma) and a very close genetic proximity (F_ST and variance in PCA/DAPC), suggesting that the northern eastern Andes uplift was not a definitive impermeable barrier at 3.4 Ma for K. scorpioides. Salgado-Roa et al. (2018) found similar divergences time between cis- and trans-Andean spider Gasteracantha cancriformis in Colombia [2.13 Ma], explained as a pattern of divergence in the presence of gene flow (because of its own dispersal ability), facilitated by remnants of low-elevation corridors along the Andes. In fact, they showed apparent evidence of admixture (STRUCTURE), suggesting no current reproductive isolation. However, Migrate-n analysis demonstrated that current gene flow does not actually exist between Caribbean and Orinoco populations. Therefore, one of the most plausible explanations for this admixture pattern is an isolation with retention of ancestral polymorphisms that can still be seen today, perhaps explained by the very slow generation time of the species preventing the rapid loss or fixation of new mutations by drift.

The Vaupés Arch fragmented the transcontinental paleo-Amazonas-Orinoco River about ~ 5 - 10 Ma (Hoorn et al., 1995), originating the final onset of the Amazon River during the Pleistocene (Reviewed in Albert et al., 2018), generating 61.2% exclusive Amazonian aquatic fauna and 16.6% exclusive Orinocan aquatic fauna. However, 22.2% are still shared between basins until the present, probably by the persistent connections between the Orinoco and Amazon basins via the Casiquiare Canal and other tributaries (Winemiller et al., 2008). However, these connections between Orinoco and Amazon are established through black water systems with extreme physical-chemical conditions, becoming a zoogeographic filter or barrier for many species that are not adapted to these conditions (Winemiller et al., 2008). This could be the case for K. scopioides, which inhabits whitewater bodies of the Amazon and Orinoco, so a Casiquiare-mediated connection would be unlikely for it at present, as reflected by the phylogeographical and genetic structure analyses.

All the characteristics that define an ESU are evidenced for each of the four sampled Colombian K. scorpioides populations. At present, we consider that our data provides support for the designation of what appear to be four Evolutionarily Significant Units (ESUs). Here, San Andres, Caribbean, Orinoco and Amazon populations show significant divergence of allele frequencies at nuclear loci (substantial reproductive isolation) and reciprocal monophyly (Moritz, 1994; Funk et al., 2012), that warrant separate management or priority for conservation because of high genetic distinctiveness and independent evolutionary trajectories (Funk et al., 2012).

The designation of four ESUs in Colombia, one for the population in San Andrés Island, and three in continental Colombia, one trans-Andean, in the Caribbean region, and two cis-Andean, in the Amazon and Orinoco regions, implies considering each of these as separate management units. Therefore, there is the need for separate management strategies for each of these units. For that reason, our recommendation to the environmental authorities in the country would be to work on independent population trend evaluations, as well as independent extinction risk analyses for each unit, aimed towards new risk categorization for these turtles in continental and insular Colombia.

Regarding the designation of the Amazon population as an ESU based on a small number of samples with a comparatively smaller number of SNPs and a high amount of missing data, we went deeper into the analyses of the 383 SNPs (and no missing data) dataset considered for this group and we found that 80 of these represented diagnostic sites, separating the Amazonian samples from all other samples from all other regions considered in our analyses. The presence of these diagnostic sites provides decisive support for this designation. Similar support for ESU or species designations for other reptile genus, for example geckos belonging in the genus Hemydactilus was based on a smaller dataset of 129 SNPs with no missing data (Leaché et al., 2014).

Population genetic theory predicts that recent colonization can result in reduced genetic diversity, known as founder effects; this in turn may lower the evolutionary potential and ultimately the ability of colonizing populations to establish and proliferate (Shirk et al., 2014). Therefore, special attention and additional studies to establish the population genetic diversity are needed for the cis-Andean ESUs, especially the one in the Amazon region, in terms of conservation priority, considering that it is the youngest population and recently established (accordingly Random Walk results).

According to the results, we believe that we present enough genomic information to support the designation of at least four new biological entities (species) within the genus, by upgrading K. s. albogulare to species level and designing the cis-Andean and trans-Andean K. s. scorpioides evidenced in our analyses as separate species. A species is a group of interbreeding natural populations that is reproductively isolated from other such groups [Biological Species Concept; (Mayr and Ashlock, 1991)], or, are defined as segments (subgroups) of population-level lineages [General Lineage Concept of Species (de Queiroz, 1998)]. Evidence supporting reproductive isolation and population-level lineages are presented here within the genus Kinosternon. Therefore, a full taxonomic revision of the genus is now recommended, as suggested by Spinks et al. (2014), including not only molecular data but also morphological information, as well as geographical range-wide taxon sampling, similar to the work done recently for the genus Chelus (Vargas-Ramírez et al., 2020).