Introduction

Illegal trade emerges as a prominent driver of global biodiversity loss in vertebrates, affecting approximately 5 579 species, including around 1 147 reptiles primarily traded as pets and products (^{Scheffers et al., 2019}). Recent estimates revealed that a staggering 421 million wildlife individuals listed under CITES fell victim to global trafficking between 1998 and 2018, with a disproportionately severe impact on developing nations facing economic challenges (^{Liew et al., 2021}). This illicit trade generates significant revenue, around $ 3.2 billion annually in the U.S., fueling criminal organizations and threatens both biodiversity and global biosecurity (^{Rosen & Smith, 2010}; ^{Tow et al., 2021}).

Turtles (order Testudines) are significantly impacted by illegal trade, resulting in substantial reductions in their populations, with 51 % of species classified as threatened (^{Böhm et al., 2013}; ^{Lin et al., 2021}). The primary contributors to this decline include habitat degradation, as well as overexploitation for food, pets, and traditional medicine (^{Stanford et al., 2020}). Concerningly, for marine turtles, an estimated 44 000 individuals are traded anually (^{Senko et al., 2022}). However, studies in species like the hawksbill sea turtles (Eretmochelys imbricata) suggest that this figure may be underestimated (^{Miller et al., 2019}). The illegal trade of continental turtles is also pervasive, exemplified by the Indian Star Tortoises (Geochelone elegans), with records from a single “trade hub” in India of approximately 50 000 individuals over 40 years (^{D’Cruze et al., 2015}). This widespread trade not only poses a significant threat to wild populations but also increases the potential risk of zoonotic disease outbreaks, as seen in cases of salmonellosis linked to pet turtles (^{Montague et al., 2022}).

Colombia is home of 33 turtle species; however, a concerning 37 % of them fall under a conservation threat category (^{Morales-Betancourt et al., 2015}). This problem is significantly driven by illegal trade, as turtles are extensively exploited for their meat and eggs, kept as pets, used in traditional medicine, and their parts are utilized in handicrafts (^{Sollund, 2017}). Among the 28 species of tortoises and freshwater turtles, 23 are exploited by local communities, and 25 are illegally traded (^{Páez et al., 2012}), Particularly impacting the genera Chelonoidis, Trachemys, Kinosternon, Podocnemis, and Rhinoclemmys (^{Arroyave-Bermudez et al., 2014}; ^{Bernal-Restrepo, 2021}).

The Colombian slider turtle (Trachemys venusta callirostris) is the most traded turtle in the country with over a million individuals extracted annually from just one locality (^{Bock et al., 2015}; ^{Mendivelso-Gamboa & Montenegro, 2007}). This turtle is classified as vulnerable (VU), primarily attributed to overexploitation of the population and habitat degradation (^{Bock et al., 2015}). The red-footed tortoise (Chelonoidis carbonarius) ranks as the second-third most traded turtle in Colombia and the foremost in some central regions (^{Bernal-Restrepo, 2021}; ^{Suárez-Giorgi, 2016}). Nationally classified as vulnerable (VU), its main threats include the illegal pet trade, driven by the belief that its possession brings good luck, and to a lesser extent, its meat consumption (^{Medina-Rangel, 2015}). Rhinoclemmys turtles constitute the fifth most traded genus in the country, with species like the Colombian wood turtle (Rhinoclemmys melanosterna) and the Maracaibo wood turtle (Rhinoclemmys diademata) classified as Near Threatened (NT) and Endangered (EN), respectively (^{Morales-Betancourt et al., 2015}; ^{Arroyave-Bermudez et al., 2014}). These turtles are hunted in their natural habitats for food consumption in various localities, as pets, and to produce handicrafts made with their shells (^{Echeverri-García et al., 2012}).

The extensive extraction and subsequent seizure of turtles demand an adequate management of captive individuals. Information regarding the Colombian slider turtle reveals that approximately 70 % of seized turtles exhibit at least one clinical problem, with a significant 56 % succumbing during the rehabilitation process (^{Castro-Cortés et al., 2022}). For successfully rehabilitated individuals, the optimal goal is their release back into the wild. However, the complex trade networks, which link international markets and extraction of turtles from remote regions to supply areas with high demand, complicates identifying their geographic origins (^{Arroyave-Bermudez et al., 2014}). As a result, current protocols release seized individuals without distinguishing their original populations, placing them in general areas within the species’ distribution range.

To contribute to the conservation of turtles in Colombia, ex-situ conservation projects have been proposed, aiming at establishing hatchery programs (^{Morales-Betancourt et al., 2015}). This vital conservation strategy for turtles aims to strengthen wild populations, while preserving their genetic diversity, and has been proven effective in recovering populations at risk of extinction (^{Barbanti et al., 2019}; ^{Martins et al., 2021}). However, in some of the breeding stocks used in reptiles ex-situ conservation in Colombia, precise information regarding their origin and genetic background is lacking, limiting the implementation of these conservation strategies (^{Saldarriaga-Gómez et al., 2023}).

Several turtle species in the country exhibit intraspecific phylogeographic patterns with individuals from different geographic locations representing distinct genetic or evolutionary lineages (^{Gallego-García et al., 2023}; ^{Vargas-Ramírez et al., 2013}). Therefore, translocating individuals without knowledge of their geographic origin can lead to unwanted outcomes, as mixing populations from distant genetic clusters risks homogenizing natural genetic diversity, outbreeding depression, and reduction of the offspring fitness (^{Frankham et al., 2011}; ^{Gippoliti et al., 2018}; ^{Oklander et al., 2020}). Reports of potential outbreeding depression have surfaced in cases of hybridization between sea turtle species, resulting in reduced reproductive success (^{Arantes et al., 2020}). Furthermore, a preliminary and tentative study identified inbreeding-outbreeding tension in hawksbill turtle populations (^{Phillips et al., 2017}). Recognizing that this risk is crucial in the conservation programs for turtle populations, the identification of geographic origin for translocation and release scenarios is mandatory (^{Alacs et al., 2007}).

In cases where precise geographic origin information is unavailable, the use of tools such as genetic markers has been proposed to address this limitation (^{Oklander et al., 2020}). This approach employing mitochondrial markers has proven effective in determining the geographic origin of Seized individuals and products in turtles as the Hawksbill turtle (Eretmochelys imbricata) in Asia and Australia (^{LaCasella et al., 2021}), the Star tortoises (Geochelone elegans) in Singapore (^{Gaur et al., 2006}), and the Sea turtle (Chelonia mydas) in China, (^{Gaillard et al., 2021}). Additionally, this strategy has proven efficient in determining the origin of captive individuals in Colombia, such as the Orinoco matamata turtle (Chelus orinocensis) (^{Lasso et al., 2018}), and the savanna side-necked turtle (Podocnemis vogli) (^{Cárdenas-Barrantes et al., 2024}). However, similar studies are lacking for several highly trafficked turtle species in the country.

In this study, we applied mitochondrial sequencing to estimate the origins of captive individuals from three of Colombia’s most trafficked turtle species: Chelonoidis carbonarius, Trachemys venusta callirostris, and Rhinoclemmys melanosterna. This approach aims to provide new tools to enhance release protocols, mitigate potential negative impacts on wild populations, and offer crucial insights into the characterization of trade networks.

Materials and methods

Reference genetic databases: In this study, we established a reference genetic database integrating secondary information from previous studies and incorporating new sequences obtained from tissue samples from individuals with known origin. These new samples were sourced from the Banco de ADN y Tejidos de la Biodiversidad (BTBC), stored at the Instituto de Genética, Universidad Nacional de Colombia (SMT 1).

For the Red-footed tortoise (C. carbonarius), we used 27 published sequences of the Cytochrome b (CYTB) gene from a previous study and sequenced 162 new samples from previously unsampled localities. For the Rhinoclemmys turtles, we employed 75 CYTB sequences from an earlier study and added 28 tissue samples from individuals collected from unsampled areas. For the Colombian slider turtle (T. v. callirostris), we obtained 58 Control Region (CR) sequences.

Sample collection from individuals of unknown origin: Tissue samples lacking information on geographic origin were obtained from both seized individuals and breeding stock turtles used in ex-situ conservation efforts. The seized individuals were sourced from the Wild Animal Rescue and Rehabilitation Unit (URRAS for its acronym in Spanish) of the Universidad Nacional de Colombia and the Wild Fauna Reception Center of the Regional Autonomous Corporation of Cundinamarca (CAR for its acronym in Spanish). Simultaneously, samples from breeding stock turtles were sourced from the Estación de Biología Tropical Roberto Franco (EBTRF) ex-situ research and conservation program at the Universidad Nacional de Colombia in Villavicencio, Meta (SMT 1).

In the analysis of the Red-footed tortoise, a total of 77 samples of unknown origin were examined, comprising 73 from the URRAS-CAR and four from the EBTRF. For the Colombian wood turtle and Maracaibo wood turtle, samples from 38 individuals were included, with 11 from the URRAS-CAR and 27 from the EBTRF. Additionally, for the freshwater Colombian slider turtle, 42 samples were collected, consisting of 32 from the URRAS-CAR and 10 from the EBTRF.

DNA Extraction, PCR and sequencing: Whole DNA was extracted from tissue samples using the NucleoSpin Tissue Kit from MACHEREY-NAGEL, Germany, following the manufacturer’s instructions. DNA concentration and purity levels (measured by 280-260 and 230-260 ratios) were quantified with an EzDrop 1 000 Micro-Volume Spectrophotometer (Bleu-Ray). For the amplification of a fragment ranging from 719 to 723 bp of the CYTB gene for C. carbonarius and Rhinoclemmys turtles, the primers mt-c-For2 (TGAGG VCARATATCATTYTGAG) and mt-f-na (GCRAATARRAAGTATCATTCTGG) (^{Fritz et al., 2006}) were used. A fragment of 638 bp for the CR gene of T. v. callirostris was obtained using the primers Tca_CRL (CCAGCTCCCAAAGCTGAGAT) and Tca_CRH (GTGGCTATTCGGACTGGTGT) (^{Balcero-Deaquiz, 2022}). The PCR solution consisted of 25 μl, with 30-100 ng of DNA (0.5-5 μl), 12.5 μl of OneTaq® 2X Master Mix with Standard Buffer (Biolabs), 1.25 μl of the forward primer, 1.25 μl of the reverse primers, and adjusted with ultrapure water. PCR reactions were carried out using a Mastercycler PRO S 6325 Thermal Cycler (Eppendorf). The protocol included a denaturation step at 94 °C for 5 minutes, followed by 35-40 cycles of denaturation at 94°C for 45 seconds, annealing at 50-60 °C for 52 seconds, extension at 72 °C for 80 seconds, and a final elongation cycle at 72 °C for 10 minutes. Amplification products were verified using 1% agarose gel electrophoresis stained with Z-vision. PCR products were purified using the ammonium acetate protocol of Bensch et al. (2 000). The purified products were sequenced using Sanger sequencing technology with an ABI 3 130 sequencer (Applied Biosystems).

Identification of phylogeographic groups: To investigate whether the studied species exhibit a phylogeographic structure that allows for separation into geographically delimited genetic groups, we employed 1) a phylogenetic approach using Bayesian and Maximum Likelihood algorithms, and 2) haplotype network reconstruction. Initially, we aligned all sequences within each genetic reference database (comprising samples of known geographic origin) using the MUSCLE function available in MEGA 11. The Bayesian species tree was constructed using BEAST 2.6.7 and BEAUTI 2 software. The most appropriate substitution model was determined using JmodelTest 2.1.1 software based on the Bayesian information criterion (BIC). Two Markov Chain Monte Carlo (MCMC) runs were conducted separately. Each MCMC included 50 000 000 generations with samples collected every 1 000 generations. The results of both MCMC runs were combined using LogCombiner 2, ensuring adequate sampling processes and convergence. The maximum credibility tree (20 % burn-in) and posterior probabilities (PP) were obtained using TreeAnnotator 2.6.6. The maximum likelihood analysis was performed using IQ-tree version 1.6.12, and the optimal substitution model was assessed using the ModelFinder function. A consensus tree with bootstrap support (100 000 bootstrap replicates) was obtained. The maximum credibility tree and consensus tree were visualized and edited using FigTree 1.4.4 software.

For the haplotypes analisys, sequences were collapsed into haplotypes using DnaSP version 6.12.03. The Haplotype networks were constructed and visualized in PopART version 1.7 utilizing the TCS algorithm (^{Cruz et al., 2021}). The relationship between haplotypes and geographic origin was depicted on a haplotype map, also created using PopART. To achieve a more precise delimitation of genetic groups, we employed the Bayesian Analysis of Population Structure (BAPS) algorithm in the R environment (RhierBAPS) using default parameters (^{Tonkin-Hill et al., 2018}).

Genetic distances and structure: For a quantitative measure of the genetic differences between phylogeographic/genetic groups, we computed p-uncorrected genetic distances using MEGA 11 software with 1 000 bootstrap repetitions (^{Santos-Rodrigues et al., 2016}). To assess the statistical significance of differences between the revealed evolutionary lineages, we performed an analysis of molecular variance (AMOVA) and differentiation coefficients (F _ST) using the Arlequin 3.5.2.2 software.

Approximation of the probable geographic origin in captive individuals: Sequences from both seized and breeding stock parental individuals were jointly analyzed with the genetic reference databases to determine their most likely origin. These analyses were based on the association of each seized individual’s haplotype with the previously identified phylogroups (^{Gaur et al., 2006}; ^{LaCasella et al., 2021}). To provide further insight into this association, we performed a repeat haplotype analysis in PopArt and RhierBAPS analysis, incorporating samples from the captive individuals.

Results

Red-footed tortoise (Chelonoidis carbonarius): The phylogenetic analyses unveiled a distinctive geographic pattern within this turtle species. Initial separation is evident between samples from the Cis-Andean region and those from the Trans-Andean region in the Northwest of the continent (Fig. 1). Within the Cis-Andean group, further subdivision is observed between samples in the Southwest of Paraguay (CH-PA) and those in the Northeast part of South America. The latter group forms two clades, comprising samples from the Eastern region in French Guiana (CH-FG) and those from the Orinoco Llanos of Colombia (CH-OL). In the Trans-Andean group (CH-TA), four subclades are discerned, although they do not constitute monophyletic groups based on geographic origins.

Fig. 1 Bayesian inference tree of the CYTB gene for the Red-footed tortoise. The corresponding support values are shown as posterior probability (PP) (left) and bootstrap support (right). 

The haplotype network and RhierBAPS analysis (Fig. 2A, Fig. 2C) further differentiate the samples from Paraguay (CH-PA), French Guiana (CH-FG), and the Orinoco Llanos (CH-OL) as distinct genetic groups. Additionally, the analysis subdivides the samples from the Trans-Andean region into four separate groups. While these clusters exhibit some geographic separation, the presence of shared haplotypes across biogeographic regions underscores the complexity of genetic relationships. Despite the absence of strict geographic correspondence, notable phylogeographic patterns emerge within these groups. The genetic group CH-TA1 is predominantly found in the East Caribe savannas and the middle Magdalena inter-Andean Valleys. CH-TA2 is highly dominant in the Guajira Peninsula in the Northwestern part of Colombia. Haplotypes of CH-TA3 are more frequent in the Sinú valleys. Lastly, CH-TA-4 is associated with a wide range of ecosystems in the North to Northwest of the country, encompassing the Magdalena River Basin, the Sinú valleys, and the Northern Pacific plains. Notably, samples from the Northwestern Pacific plains exhibit multiple genetic group associations, with 40 % of them displaying a unique haplotype (TA3-6).

Fig. 2 A. Haplotype network of the Red-footed tortoise constructed using sequences of individuals with known geographic origins. B. Haplotype network of Red-footed tortoise including sequences from seized individuals in Colombia. Hatch marks indicate the number of mutations separating each haplotype. Colored boxes represent the haplotypes grouped under a population according to RhierBAPS Analysis. C. Haplotype map plot for the Red-footed tortoise. The size of the circles indicates the number of samples with a respective haplotype. In set photography of an adult female Chelonoidis carbonarius from the EBTRF. Photo by Juan Manuel Vargas Ramírez. 

To assess the significance of differences between the studied regions (French Guyana, Paraguay, Orinoco Llanos, Magdalena River Basin, Middle Magdalena Inter-Andean Valleys, East Caribe savannas, Sinú Valleys, Pacific Plains, and Guajira Peninsula), we conducted an Analysis of Molecular Variance (AMOVA) and an F _ST test. The analysis yielded a highly significant value (p < 0.01), with 83 % of the observed variance occurring between populations. Pairwise comparisons further confirmed significant differences (p < 0.05) between all regions, except for the Magdalena River Basin vs. East Caribe savannas (p = 0.25) and the Middle Magdalena valleys vs. East Caribe savannas (p = 0.25). Consequently, the Magdalena valleys and East Caribe savannas were considered a single genetic group.

Genetic distances between French Guyana, Paraguay, Orinoco Llanos, and Northern groups ranged between 2-3 %, whereas the distance within the Northern groups was much lower, ranging from 0.06-0.2 % (SMT 2).

Given these findings, we evaluated the utility of mitochondrial sequencing to estimate the most probable origin in 73 turtles rescued from illegal trade and four individuals from the ex-situ conservation program of the EBTRF (Fig. 2B). The analyses revealed that two of the seized specimens were misidentified and corresponded to the yellow-footed tortoise (Chelonoidis denticulata). Among the seized individuals, 75 % (50) of the samples were assigned to the Trans Andean group, while only 25 % (17) were assigned to the Cis-Andean Orinoco Llanos (Fig. 2B). Within the Trans Andean region, there appeared to be a higher impact on the Guajira Peninsula, as 27 % (18) of the samples were assigned to CH-TA2. Additionally, there was a significant impact on the East Caribe Savannas and the Middle Magdalena Inter-Andean Valleys, with 23 % (16) assigned to CH-TA1. Conversely, the Western ecosystems seemed to be less exploited, with 14 % (9) assigned to CH-TA3 (associated with the North to Northwestern ecoregions) and 11 % (7) to CH-TA4. Within the CH-TA4 population (associated with the Sinú valleys), one sample was associated with the haplotype found exclusively in the Pacific plains (TA4-17). Regarding the four samples from the ex-situ conservation program, three of them were assigned to the Orinoco Llanos clade, and one to the CH-TA-2 group, most likely originating from the Guajira Peninsula.

Colombian wood turtle (Rhinoclemmys melanosterna): The phylogenetic analysis revealed several geographically delimited groups within this species and a non-monophyletic relationship between species within the R. melanosterna complex (Fig. 3). An initial separation is observed between the Western and Eastern clades. The Western clade comprises R. melanosterna turtles from the Pacific plains to Sinú plains, as well as R. funerea samples from Central America. On the other hand, the Eastern clade includes R. melanosterna turtles in Northeast Colombia, R. diademata from Venezuela, and R. punctularia from French Guyana.

Fig. 3 Bayesian inference tree of the CYTB gene for the Colombian wood turtle. The corresponding support values are shown as posterior probability (PP) (left) and bootstrap support (right). 

The Western clade was further subdivided into samples of R. funerea from Central America and R. melanosterna in West Colombia. The latter was also subdivided into samples from the Northwest and those in the Southwest. The Northwest clade showed further subdivision into two genetic groups. The first one was predominantly delimited to the Northwest area (RH-NW), encompassing samples from the North of the Pacific Plains (Antioquia and North-Chocó) and those in the North of the Sinú Valleys. The second one grouped samples from the Southern part of the Sinú Valleys (RH-SSV) close to the Mompos Depression. On the other hand, in the Southwest clade, two genetic groups were observed. The first one included samples from the delta of the San Juan River in Chocó (RH-RSJ), and the other extended from the central to the Southern part of the Pacific Plains (RH-PPCS), ranging from South Chocó and Valle del Cauca in Colombia to the Northwest of Ecuador.

In the Eastern clade, we observed a separation into two groups. The first one corresponded to samples from the Magdalena Valleys (RH-MV), and the second one included turtles in the East. The latter clade was subsequently divided into a group comprising samples from the R. punctularia in French Guyana (RH-RP) and an aggregation of R. melanosterna from the East Caribe savannas of Colombia (RH-EC) and R. diademata from Venezuela (RH-RD).

The haplotype and RhierBAPS analysis (Fig. 4A, Fig. 4C) predicted eight genetic groups within this turtle. These groups align with the previously described clades in the phylogenetic analysis. However, the RH-EC and RH-RD clades were grouped into a single genetic group (RH-E). Nonetheless, they possess distinctive haplotypes. In the phylogenetic and haplotype analyses, one sample from the North of the Sinú Valleys has a genotype associated with the Southern Sinú genetic group (RH-SSV), and one sample from the central part of the Pacific plain was grouped under the Northwest group (RH-NW). Apart from these exceptions, the geographic regions exhibited a strong genetic identity, with no shared haplotypes among them.

Fig. 4 A. Haplotype network of the Colombian wood turtle constructed using sequences from individuals with known geographic origins. B. Haplotype network of Colombian wood turtle including sequences from seized individuals in Colombia. Hatch marks indicate the number of mutations separating each haplotype. Colored boxes represent the haplotypes grouped under a population according to RhierBAPS Analysis. C. Haplotype map plot for the Colombian wood turtle. The size of the circles indicates the number of samples with a respective haplotype. In set photography of an adult female Rhinoclemmys melanosterna from the EBTRF. Photo by Juan Manuel Vargas Ramírez. 

We explored the genetic differences between the North-West, Pacific Plains-Río San Juan, Pacific Plains Central-South, Magdalena Valleys, East Caribe Savannas, and East Colombia-Venezuela (R. diademata) geographic regions using the AMOVA test. We excluded the R. punctularia and R. funera groups from the analysis, as they were represented by only three and two samples, respectively, which could lead to biases. Additionally, these groups already showed a noticeable separation in the haplotype analysis. Regarding R. diademata, given that the reference database only included three samples, we supplemented it with two additional sequences from seized individuals. These individuals were identified as belonging to this species based on reliable morphological differences particularly the coloration patterns. These analyses confirmed the genetic differences between geographic regions (p < 0.001) with 94 % of the variance found between populations. The pairwise analysis revealed a significant difference between all the regions (p < 0.05). The genetic distances between the regions ranged from 0.3-3.5 %, being lower when compared to geographically close groups, while higher when comparing groups in distant localities (SMT 2).

We used the CYTB reference to identify the most probable origin for 11 seized individuals and 27 turtles from the ex-situ conservation program of the EBTRF. Of the seized individuals, nine corresponded to R. melanosterna, and seven of them were assigned to the East clade, with six associated with the Magdalena Valley (RH-MV) group, and one to the East group (RH-E) having the genotype of East Colombia (RH-EC, E-1). According to the haplotypes, those of the Magdalena most likely come from the middle Magdalena-inter Andean Valley. Only two of these samples were associated with the West clade, both assigned to the Central to South Pacific Plains (RH-PPCS) group. The other two seized turtles corresponded to R. diademata, and as expected, were associated with the East group (RH-E). However, both presented new haplotypes, with one of them (E-6) showing a noticeable degree of separation from previously reported haplotypes.

In contrast, all the Colombian wood turtles being used in the EBTRF conservation program were assigned to groups on the Western clade (Fig. 4B). 44 % (12) of the samples were assigned to the Central to South Pacific Plains (RH-PPCS) group, 37 % (10) to the Southern Sinú plains (RH-SSV), 11 % (3) to the Pacific Plains-Río San Juan (RH-RSJ), and 7 % (2) to the North-West group (RH-NW).

Freshwater Colombian slider turtle (Trachemys venusta callirostris): The phylogenetic analysis (Fig. 5) shows two separate clades within the Colombian territory. The first one grouped turtles from the Northern region of the country across the Caribe lowlands (TC-N), while the second one corresponded to turtles from the inter-Andean valley at the middle Magdalena plains area (TC-IA). The haplotype and RhBAPS analysis (Fig. 6A, Fig. 6C) also indicated these two genetic groups. The AMOVA test showed a significant genetic difference between these two regions (p < 0.001) with 81 % of the variance found between populations. The genetic distance between these two groups was relatively low, at 0.6 %.

Fig. 5 Bayesian inference tree of the CR gene tree for the Colombian slider turtle. The corresponding support values are shown as posterior probability (PP) (left) and bootstrap support (right). 

Fig. 6 A. Haplotype network of the Freshwater Colombian slider turtle constructed using sequences of individuals with known geographic origins. B. Haplotype network of the Freshwater Colombian slider turtle including sequences from seized individuals in Colombia. Hatch marks indicate the number of mutations separating each haplotype. Colored boxes represent the haplotypes grouped under a population according to RhierBAPS Analysis. C. Haplotype map plot for the Freshwater Colombian slider turtle. The size of the circles indicates the number of samples with a respective haplotype. In set photography of an adult female Trachemys venusta callirostris from the EBTRF. Photo by Juan Manuel Vargas Ramírez. 

Using this information, we estimated the provenance of 32 seized turtles and 10 individuals from the EBTRF conservation program. All the individuals of seized turtles were assigned to the Northern group (Fig. 6C). Also, most of the turtles from the EBTRF (nine individuals) were assigned to the Northern group, while only one sample had the inter-Andean genotype. Table 1 presents a summary of the results from the haplotype and genetic differences analysis, along with the phylogeographic groups and proposed independent conservation units.

Discussion

Our results confirmed the existence of intraspecific genetic diversity and structure in two continental freshwater turtles and one tortoise from Colombia. Furthermore, we observed that captive individuals from confiscations and an ex-situ conservation program originate from multiple genetic/geographic groups. This emphasizes the significance of a thorough release process for captive individuals, and the importance of such programs in preserving the species’ genetic diversity. Colombian environmental law recognizes the risk of genetic contamination to wildlife populations during release processes and advocates for the use of genetic tests to assess genetic variability before the release of seized individuals (^{Choperena-Palencia & Mancera-Rodríguez, 2016}; ^{MADS Resolución 2064, 2010}). Nevertheless, the insufficient understanding of the phylogeography and population genetics of traded turtle species, coupled with the absence of standardized and validated laboratory and analytical protocols, has restricted the application of genetic analysis to inform the release of individuals. Consequently, a significant portion of reintroduction and release processes in the country occurs without reliable information about the animals’ geographic origin, attracting significant criticism (^{Jiménez & Cadena, 2004}).

Documented efforts to enhance translocation protocols for seized or captive wildlife in the country are relatively scarce. In the case of turtles, a study on the Orinoco matamata turtle (Chelus orinocensis) utilized mitochondrial COI and CR regions to determine the most probable geographic origin of seized individuals. The findings revealed that, despite being rescued in Leticia-Amazonas, they originated from a distant locality in the Orinoco Llanos, allowing for the proper release of the individuals to their original geographic region (^{Lasso et al., 2018}). Additionally, a study on the savanna side-necked turtle (Podocnemis vogli) employed microsatellite markers to estimate the origin of 26 captive individuals (^{Cárdenas-Barrantes et al., 2024}).

Internationally, there are more instances of the use of forensic DNA analysis to study the trade of turtles and improve translocation processes. A study utilizing CYTB allowed for the approximation of geographic origin in the endangered European pond turtles (Emys orbicularis) retrieved from recovery centers in Spain, Portugal, and Morocco (^{Velo-Antón et al., 2021}). In Asia, CYTB and CR genes were employed to address the geographic origin of seized Star tortoises (Geochelone elegans), enabling their repatriation to their original populations (^{Gaur et al., 2006}). Similarly, this approach has been used in Hawksbill turtles (Eretmochelys imbricata) and green turtles (Chelonia mydas), allowing for the identification of hunting hotspots (^{Joseph et al., 2019}; ^{LaCasella et al., 2021}). These previous works and our results highlight the importance of genetic analysis in the management of seized turtles.

Red-footed tortoise (Chelonoidis carbonarius): In agreement with previous research, our study revealed a substantial phylogeographic structure in the red-footed tortoise across its distribution range. We identified four main clades geographically delimited to the Southwest in Paraguay, an Eastern clade in French Guiana, the Orinoco Llanos of Colombia, and the Trans-Andean region. The original work that explored the CYTB genetic structure in the species also demonstrated these four major clades and reported an additional one in Brazil (^{Vargas-Ramírez et al., 2010}). Subsequently, a study utilizing RADseq analysis (30 327 neutral SNPs) to compare the genetic structure of the Orinoco Llanos and Trans-Andean groups revealed significant differences between them, providing evidence of allopatric reproductive isolation, local adaptation, and ecological divergence (^{Gallego-García et al., 2023}). The multilocus SNPs genomic analysis revealed moderately differentiated clusters within the Trans-Andean region, proposing four management units: the Caribbean Coast (Magdalena River Delta to the Sinú Valleys), the Darien Mountain Forest (Pacific plains), the Rainforest of the Magdalena River Valley (East Caribe savannas), and the Tropical Dry Forest in the inter-Andean valleys (Middle Magdalena inter-Andean Valleys) (^{Gallego-García et al., 2023}). Our study also detected a significant, albeit less pronounced structure among ecoregions in the Trans-Andean region of Colombia. However, we did not observe a significant difference between the East Caribe savannas and the Middle Magdalena inter-Andean Valleys. Additionally, our analysis indicated that samples from the Magdalena River basin were more closely related to those from the East Caribe savannas than to the Sinú valley. Furthermore, with the inclusion of samples from the Guajira Peninsula, we identified a significant structure in this region. Given the higher resolution and more comprehensive information provided by RADseq data, we recommend maintaining the proposed geographic delimitations for management units as suggested by Gallego-García et al. 2023, while incorporating the Guajira Peninsula as a new separate management unit.

The observed patterns of defined phylogeographic groups with the admixture of individuals in the Trans Andean region have been linked to the recent natural and/or human-mediated dispersion of the species (^{Gallego-García et al., 2023}). This is one of the most traded reptiles in the country and unlike other turtles; it is primarily sought after as a pet rather than for use in food products (^{Echeverry-Alcendra, 2019}). Consequently, this species is susceptible to being artificially translocated and colonizing areas far from their original habitat, either through escape or intentional release by their owners.

The analysis of seized individuals revealed that red-footed tortoises recovered in the central part of Colombia originate from several regions and genetically distinct groups along the country. Both the Trans Andean and Orinoco Llanos were identified as sources for traded tortoises, with a higher number of specimens assigned to the former (75 % vs. 25 %, respectively). Within the Trans Andean clade, the Northeastern territory of the Guajira Peninsula and East Caribe Savannas to Middle Magdalena inter-Andean Valleys were the most affected. Nevertheless, the ecoregions of the Northwest also appear to be used as hunting sites, but in lower proportions.

The extraction of individuals from various regions of the country explains the high number of reported seizures of this turtle, making it one of the two most traded turtles in Colombia (^{Arroyave-Bermudez et al., 2014}). It has been documented that hunting and exploitation of this species in the Orinoco Llanos are common, and trade networks in the area include both local consumption and transportation outside the country to Venezuela (^{Medina-Rangel, 2015}). Thus, not all poached individuals are transported to the interior of Colombia, and although we found a lower number of individuals coming from the Orinoco Llanos compared to those from the Trans-Andean region, the severity of hunting in this region should not be underestimated.

A previous study on the trade networks of the species, relying on information about the seizure place and reported origin collected by environmental entities in Colombia, suggested that Chelonoidis turtles seized in Cundinamarca were primarily (or exclusive) harvested from the Orinoco Llanos, and that this region served as an “intermediate or passage area” for their mobilization to Northern areas like Antioquia and Santander (^{Arroyave-Bermudez et al., 2014}). However, our results contradict this hypothesis, as only a fraction of the turtles seized in Cundinamarca originated from the Orinoco Llanos. The majority came from multiple areas of the Trans Andean Coastal regions. It is also noteworthy that the Red-footed tortoise is the most seized reptile in the central region of Colombia (Cundinamarca and Boyacá), representing up to 65 % of the seized reptiles (^{Suárez-Giorgi, 2016}). The substantial number of traded individuals obtained from various regions across the country, implies that central regions of Colombia, such as Cundinamarca, may function as the final destination or a distribution hub for multiple trade networks in the country. Moreover, the discrepancies between studies may suggest that information collected during the rescue process is unreliable, possibly due to traders concealing or falsifying details about the trade networks. However, these differences may also be attributed to the varying timeframes of the studies, reflecting shifting trends in consumption and hunting preferences across different sampling periods.

Our results emphasize the importance of conducting release processes of this turtle with proper genetic analysis, given its pronounced phylogeographic structure, the diverse origin of seized individuals, and the lack of reported morphological differences between genetic groups. This highlights the necessity for a revision of the current release protocols, as the widespread release of seized individuals without information about their geographic origins is a common practice (^{Echeverry-Alcendra, 2019}; ^{Medina-Rangel, 2015}).

Our data also emphasizes the significance of genetic evaluation in ex-situ conservation programs. In the case of the four individuals from the EBTRF, three of them originate from the Orinoco Llanos, and one is most likely from the Guajira Peninsula. Consequently, we recommend segregating of Cis Andean individuals to prevent unintended hybridization. Moreover, we suggest expanding the representation of individuals from the Trans Andean region, encompassing a broader range of ecoregions and genetic groups to ensure the preservation of the genetic diversity within the species. Additionally, increasing the number of animals in this conservation strategy is crucial, as a minimum of 15 founders is recommended to minimize the loss of genetic diversity (^{Witzenberger & Hochkirch, 2011}).

In the country, there are also zoo-breeding farms for the species that are mandated by law to release part of the breeding animals as part of the repopulation quotas, intending to improve the conservation status of the species (^{Echeverry-Alcendra, 2019}). However, our results indicate that these actions may pose a risk of genetic contamination to the receiving populations. Therefore, genetic evaluation of animals in zoo-breeding farms is also necessary to meet the conservation expectations of these sustainable use initiatives.

Colombian wood turtle (Rhinoclemmys melanosterna): We identified a significant phylogeographic structure in the R. melanosterna complex, comprising eight distinctive genetic groups. Our results align with the phylogeographic analysis conducted by ^{Vargas-Ramírez et al. (2013)}, but the inclusion of new samples has led to some changes in the delimitations proposed by the authors. Our findings support the separation of R. funerea from Centro America (RH-RF) and R. punctularia (RH-RP) from French Guyana as distinctive genetic groups. Additionally, we concur with the grouping of R. diademata and R. melanosterna in the East Caribe savannas as a single genetic group (RH-E), albeit with distinctive haplotypes. Furthermore, our results confirm the presence of a genetic group in the South of the Sinú plains (RH-SSP), near the Mompos Depression.

In the original study, two phylogroups were proposed for the Northwest, one in the North of the Pacific plains extending from Panama to the North of the Sinú plains in Colombia, and a second in the Urabá-Antioquia region. However, the second group was represented by a single sample (^{Vargas-Ramírez et al., 2013}). In our analysis, we included more samples from Urabá in Antioquia, as well as samples from North Chocó, and found a single genetic group for the Northwest of South America, spanning from the North of the Pacific plains (including the Urabá-Antioquia region, North Chocó, and Panama) to the North of the Sinú plains (RH-NW).

In a previous study, a genetic group was identified in the Southwest Pacific plains and hypothesized to be distributed throughout the Chocó region (^{Vargas-Ramírez et al., 2013}). While we also detected this group, our findings show that it is restricted to the San Juan River delta (RH-SJR) rather than the entire Chocó region. The original analysis also described a genetic group extending from Valle del Cauca to Northwest Ecuador (^{Vargas-Ramírez et al., 2013}), but our results indicate that this group is also present in the Southern part of the Chocó region (CSPP).

In the Eastern region of Colombia, ^{Vargas-Ramírez et al. (2013)} described a genetic group delimited to the Middle Magdalena Inter-Andean Valley. Our analysis also supports this group and expands its range to the Magdalena River basin, suggesting that this clade is distributed across the low and middle Magdalena valley.

Given the notable phylogenetic structure in this turtle, we successfully identified the most probable origin of 11 seized individuals. Our findings indicate that turtles rescued from illegal trade have diverse origins, with a more pronounced impact on the Eastern part of the country, particularly in the Middle Magdalena inter Andean valley, which is closer to the interior of the country (Cundinamarca) where the individuals were seized. However, two individuals were assigned to a markedly distant area in the central-South Pacific plains, indicating that commercialization from distant trade networks also occurs. The two samples from the Maracaibo wood turtle (Rhinoclemmys diademata) were, as expected, associated with the East group but presented new haplotypes. The reference database only had a few samples of this turtle from Venezuela, potentially explaining the reported new haplotypes, as the seized individuals most likely come from the Colombian populations. Considering the species’ distribution range, they most likely come from the Eastern part of Norte de Santander, which is its only area of distribution in Colombia (^{Armesto et al., 2014}).

The strong phylogeographic structure of natural populations and heterogeneity in the origin of seized individuals underscore the necessity for the release process to involve proper genetic analysis or reliable information about the geographic origin of the animals. In this turtle species, differences in head stripe patterns between geographic regions have been reported and could potentially serve as a simpler and more cost-effective method for determine their origin (^{Vargas-Ramírez et al., 2013}). However, the reliability of this method has not been tested, and these types of morphological traits are not always clear across different individuals, especially in juveniles.

The trade networks of this species in Colombia remain largely unknown, and the extent of distribution routes, particularly for the pet trade, has not been characterized (^{Echeverri-García et al., 2012}). This lack of information limits the development of proper management plans to address the illegal trade of this turtle, a notable problem, given that seized reports of Rhinoclemmys turtles have been drastically increasing (^{Arroyave-Bermudez et al., 2014}). Our analysis unveiled at least four trade networks for this turtle, involving the supply of individuals to the central part of the country from their extraction areas in the Magdalena valleys, East Caribe plains (in Cesar), Norte de Santander, and the central-South Pacific plains. This data aligns with information regarding the harvesting of this turtle in the communities from the Magdalena valley and Pacific Plains (^{Echeverri-García et al., 2012}). However, commercialization of this turtle is relatively low in the central region of Colombia compared to Chelonoidis and Trachemys turtles (^{Suárez-Giorgi, 2016}), while in other markets, like Buenaventura and Cali, this is one of the most traded turtles (^{Corredor et al., 2007}). Thus, replicating this analysis in other regions of the country would be valuable to gain a better understanding of the trade networks.

The analysis of turtles from the EBTRF revealed that this conservation program has representation from the four Western genetic groups, while there was no representation from the Eastern clades. Thus, the separation of the turtles according to their assigned phylogeographic group is necessary to avoid the generation of hybrids with potentially lower environmental adaptability (^{Barbanti et al., 2019}; ^{Weeks et al., 2011}). We also recommend including founders from the Eastern clades to formulate a more comprehensive conservation plan that preserves the genetic diversity of the country.

Freshwater Colombian slider turtle (Trachemys venusta callirostris): Unlike previous examples, the freshwater Colombian slider turtle exhibited a less pronounced phylogeographic structure with only two differentiable groups. The first group (TC-N) was in the North of Colombia across the Caribbean plains (including the Magdalena and Sinú valleys) and the Guajira Peninsula. The second group corresponded to turtles from the middle Magdalena inter-Andean valleys (TC-IA). These were the same groups reported in the original study on the species used as a reference database. The author of this population genetic study additionally used microsatellite markers, retrieving also these two populations (^{Balcero-Deaquiz, 2022}). This is consistent with a study using allozyme analysis in turtles from the Caribbean plains at the Mompos Depression, which found low genetic variability and a lack of structuration (^{Martínez et al., 2007}). Thus, there is no evidence of additional structuration in the country.

All 32 seized turtles evaluated were assigned to the North group, suggesting stronger exploitation of the Northern coastal plains populations. Previous reports have also suggested a trade network for the species with poaching in the coastal region towards the interior of the country in the Antioquia, Quindío, and Santander departments (^{Arroyave-Bermudez et al., 2014}). Our results support the hypothesis of an extraction center in the Caribbean coastal region and suggest that part of these harvested turtles also supplies the demand in the most central regions of Cundinamarca department.

Our analysis suggests that the liberation process of individuals from the central region of Colombia should focus on the coastal region of the country to avoid genetic contamination in the inter-Andean populations. However, we recommend further evaluations with larger sample sizes to better assess the use (if any) of the inter-Andean group in the trade network. Considering that this species is the most exploited turtle in Colombia, with millions of individuals harvested yearly (^{Restrepo et al., 2014}), a precise evaluation of trade patterns may require a higher number of samples. This characterization is important as there are not reported morphological characteristic that allow the differentiation between the two genetic groups.

It is important to note that the pattern found in the central part of the country (Cundinamarca) may differ in other areas with an even higher demand for the species, such as Antioquia and Santander (^{Arroyave-Bermudez et al., 2014}). Considering that the Inter-Andean clade is distributed within the Santander and Antioquia territories, it may be exploited to supply the internal demand in these areas. Hence, genetic analysis of seized samples from other territories is required to better understand the trade networks and define the most adequate release protocols for the species.

Regarding the EBTRF samples, nine out of 10 also came from the Northern coastal clade, while only one of them was associated with the middle Magdalena Inter-Andean region. Thus, we suggest separating this individual to prevent unwanted crossbreeding. Additionally, we recommend increasing the number of founders from the Inter-Andean region to establish a program preserving both lineages within the country.

General remarks: Our study significantly contributes to the understanding of the phylogeographic structure of turtles of conservation concern within the Colombian territory. We emphasize the importance of developing intraspecific conservation plans for these turtles, utilizing the proposed phylogeographic groups as independent conservation units. Additionally, our findings underscore the critical need to improve the management practices for captive individuals, particularly in translocation and release processes, to prevent genetic contamination and potential outbreeding depression in wild populations. Consequently, genotyping seized individuals and founders from ex-situ conservation programs should be a mandatory step before their release or use in breeding programs.

While this study focused on three turtle species, it is crucial to acknowledge the documented evidence of genetic structure within Colombia for other turtles commonly involved in illegal trade, such as the scorpion mud turtle (Kinosternon scorpioides), matamata turtles (Chelus spp.), the side-necked turtle (Podocnemis vogli), and the Magdalena River turtle (Podocnemis lewyana) (^{Cárdenas-Barrantes et al., 2024}; ^{Hurtado-Gómez et al., 2024}; ^{Vargas-Ramírez et al., 2012}, ^{Vargas-Ramírez et al., 2020}). Therefore, obtaining information on the geographic origin should be a preliminary step in the release process for all turtle species, unless evidence demonstrating a lack of genetic structure in their populations is available.

This study validates the efficacy of mitochondrial markers for approximating the geographic origin of three of the most traded turtles in Colombia. We recommend the use of these tools to inform the management of seized individuals, contributing to more effective and secure release programs as a conservation strategy. Additionally, we propose employing these markers in the initial decision-making processes for managing and forming reproductive groups in ex-situ breeding programs. For this propose the level of classification presented different degrees of confiability among the studied species. For T. v. callirostris and R. melanosterna, the phylogeographic groups were more distinctly differentiated, resulting in higher confidence in their assignment. In contrast, for C. carbonarius, only the division between trans-Andean and cis-Andean groups was clearly defined, while differentiation within the trans-Andean regions was less evident (with shared haplotypes across regions), making the assignments less reliable.

Our study underscores the potential of molecular tools to reconstruct trade networks in Colombia and to verify current trafficking route hypotheses. However, a more accurate reconstruction would require extensive sampling across additional regions and the integration of other data sources to better identify poaching locations and transport routes. Our study also provides an initial understanding of the most vulnerable regions targeted by hunters. However, the phylogeographic analyses have relatively limited resolution, allowing us to identify only broad ecoregions in the country. To achieve a more precise determination of poaching areas, the use of additional tools would be required.

It is important to note that mitochondrial markers only provide information about the maternal genotype and do not identify individuals resulting from hybridization between lineages. Therefore, in cases where there are concerns about potential hybridization, the use of additional markers such as microsatellite or SNPs data is crucial. This is particularly relevant in ex-situ conservation programs and zoo breeding centers that lack proper documentation on the origin of founder animals over several generations. Additionally, the incorporation of additional markers offers insights into genetic diversity, inbreeding, heterozygosity, parental veracity, and relatedness; important information for determining the fitness of captive populations and their suitability for translocation and release processes (^{Barbanti et al., 2019}; ^{Saldarriaga-Gómez et al., 2023}).

Ethical statement: The authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgments section. A signed document has been filed in the journal archives.

See supplementary material