Central American Trachemys revisited: New sampling questions current understanding of taxonomy and distribution (Testudines: Emydidae)

Using 3226-bp-long mtDNA sequences and five nuclear loci (Cmos, ODC, R35, Rag1, Rag2, together 3409 bp), we examine genetic differentiation and relationships of Central American slider turtles ( Trachemys grayi , T. venusta ). Our investigation also included samples from taxa endemic to North America ( T. gaigeae , T. scripta ), the Antilles ( T. decorata , T. decussata , T. stejnegeri , T. ter-rapen )


Introduction
Slider turtles (Trachemys Agassiz, 1857) represent one of the most speciose and most widely distributed genera of turtles. Currently, 11-13 species plus 9-11 subspecies are recognized that are distributed from the Great Lakes region of North America across Mexico and Central America to northern South America, with additional isolated occurrences in northeastern Brazil (Maranhão) and the Rio de la Plata region of Brazil, Argentina, and Uruguay. Four further species with three additional subspecies occur on the Antillean islands (TTWG 2021). The taxonomy and distribution of Mexican and Central American slider turtles are insufficiently known, with many taxonomic changes in recent decades (Moll and Legler 1971;Legler 1990;McCord et al. 2010;Fritz et al. 2012;Legler and Vogt 2013;Parham et al. 2015;Seidel and Ernst 2017;Vargas-Ramírez et al. 2017). Until today, the species boundaries and the validity of several taxa are contentious (see the reviews in TTWG 2017TTWG , 2021. In the present study, we follow the taxonomy of the latest edition of the checklist of the Turtle Taxonomy Working Group (TTWG 2021), but treat T. venusta callirostris and T. v. chichiriviche as subspecies of T. venusta (and not as representing a distinct species T. callirostris with two subspecies) and T. dorbigni adiutrix as a subspecies of T. dorbigni (and not as a full species), as proposed by Fritz et al. (2012).
Based on morphology, several subspecies of slider turtles were described in the past decade from southern Mexico and Central America (McCord et al. 2010). However, the validity of one of these subspecies (Trachemys venusta uhrigi; Fig. 1) was soon challenged based on mitochondrial and nuclear DNA sequence data (Fritz et al. 2012). Other authors (e.g., Legler and Vogt 2013) completely ignored the descriptions by McCord et al. (2010). Only in a few cases are the distribution ranges of the currently recognized Mexican and Central American taxa supported by genetically verified data. In particular, the putative geographic distributions of the subspecies described by McCord et al. (2010) were largely specula-tive and have later been tentatively revised (TTWG 2017(TTWG , 2021 using the patchy sampling of a molecular study (Fritz et al. 2012). Following the currently accepted understanding (TTWG 2017(TTWG , 2021, there are two distinct species occurring in the Pacific and Caribbean drainages of Central America, T. grayi and T. venusta, respectively. However, near Acapulco de Juárez (Guerrero, Mexico), i.e., along the Pacific versant, slider turtles from the Caribbean (T. venusta) are thought to be introduced (Parham et al. 2015).
To shed new light on the distribution and taxonomy of Central American slider turtles, we supplemented the genetic data of Fritz et al. (2012) and Vargas-Ramírez et al. (2017) with DNA sequences from material collected by James R. McCranie in Honduras and by Rodney Scott and Raúl Fournier in Costa Rica until 2012. Furthermore, to complete the taxon sampling from Mexico and the Antilles, we sequenced some specimens from the collection of the Museum of Zoology, Senckenberg Dresden, using aDNA approaches. We also sequenced two museum specimens from Honduras, among them the holotype of T. v. uhrigi McCord, Joseph-Ouni, Hagen & Blanck, 2010 from the Florida Museum of Natural History, Gainesville, Florida. Our sampling also includes additional material from the region of Acapulco de Juárez to re-examine the status of those populations.

Materials and Methods
We supplemented previously published mitochondrial and nuclear DNA sequences of Trachemys and some related taxa (Fritz et al. 2012;Vargas-Ramírez et al. 2017;Vamberger et al. 2020) with newly generated homologous data for 40 additional samples (Table S1).
According to the state of preservation of the samples, we used different workflows. For 10 blood and tissue samples stored at -80°C, we obtained sequences of the same mitochondrial and nuclear genes as in Fritz et al. (2012) and Vargas-Ramírez et al. (2017) using Sanger sequencing (mtDNA: 12S, ND4L, ND4, cyt b plus adjacent tRNA-Thr, together 3226 bp; nuclear loci: Cmos, ODC, R35, Rag1, Rag2, together approx. 3400 bp). Seventeen additional samples were extracted DNA which had been frozen for approximately 10 years at -20°C. For these samples, only a few microliters of high-molecular-weight DNA were available. We processed these DNA samples and an additional sample from a roadkill from Grand Cayman using a Next Generation Sequencing (NGS) approach, including in-solution hybridization capture. Twelve further samples were tissues from museum specimens preserved 30-50 years ago. We subjected these to the same general NGS approach, but worked in a cleanroom facility that is physically separated from the main laboratory to avoid contamination with foreign DNA. Using the described NGS workflows, we obtained for each specimen a complete or near-complete mitochondrial genome (mt-genome) and DNA sequences for the same five nuclear loci as for the Sanger-sequenced samples. Details on the laboratory procedures as well as subsequent data handling are explained in the Supplementary Information (Supplementary Material 1).
For phylogenetic analyses, we used the above-mentioned mitochondrial genes, which were already utilized in previous studies (Fritz et al. 2012;Vargas-Ramírez et al. 2017;Vamberger et al. 2020), and the Maximum Likelihood (ML) and Bayesian Inference approaches implemented in RAxML 8.0.0 (Stamatakis 2014) and MrBayes 3.2.6 (Ronquist et al. 2012). Software settings and further details are specified in the Supplementary Information (Supplementary Material 1).
In addition, we used the concatenated nuclear DNA sequences to create a phylogenetic splits network in the program SplitsTree4 v4.18.3 (Huson and Bryant 2006) using uncorrected p distances and the NeighborNet algorithm. We also performed Principal Components Analyses (PCAs) using the R package Adegenet (Jombart 2008) on the same nuclear dataset. For PCAs, we utilized only polymorphic sites and coded sequence data as alternate alleles (1 or 0). Missing data and ambiguous sites were imputed with the mean allele frequency for that locus, and all sites were then centered and scaled. To create graphical summaries of the genetic variation within taxa in the PCAs, we delineated inertia ellipses around the cloud of points using the default ellipse size of 1.5 in Adegenet.

DNA sequences from fresh and collection material
For museum specimens, the roadkill, and previously fro zen DNA samples, complete (D 39069, D 39078, D 41608) or almost complete mt-genomes were assembled (remaining 27 samples; Table S1). The mt-genomes had a length of approximately 16,600 bp. The near-complete mt-genomes only missed part of the 5ˈ-end of the control region (see Supplementary Material 1 for the individual mapping results). The arrangement of the 13 protein-coding genes, 2 rRNA genes, the control region, and 22 tRNA genes corresponds to that of a typical vertebrate mt-genome, e.g., Homo sapiens Linnaeus, 1758. For the remaining samples, the same mitochondrial genes (12S, ND4L, ND4, cyt b plus adjacent tRNA-Thr; together 3226 bp) were Sanger-sequenced as in Fritz et al. (2012) and Vargas-Ramírez et al. (2017). The sequences of the nuclear loci obtained either with Sanger sequencing or NGS approaches had the following lengths: Cmos -563 bp, ODC -621 bp, R35 -974 bp, Rag1 -614 bp, and Rag2 -637 bp. DNA sequences generated for the present study are available under the European Nucleotide Archive (ENA) project accession number PRJEB60683; individual ENA accession numbers for genes and mt-genomes are listed in Table S1.

Mitochondrial phylogeny
Trachemys constitutes a maximally supported clade that is sister to Graptemys. Trachemys contains in a polytomy three major, maximally supported clades, in Figure 2 from top to bottom, clade A including T. venusta, T. grayi, T. medemi, and T. dorbigni from Mexico, Central America and South America, clade B corresponding to the Antillean species T. decussata, T. terrapen, T. decorata, and T. stejnegeri, and clade C containing the North American species T. scripta and T. gaigeae.
Within clade A, which is in the focus of the present investigation, T. venusta, T. grayi, T. medemi, and T. dorbigni constitute distinct, well-supported clades. Trachemys venusta and T. grayi are sister taxa, though only with weak support, and constitute together the sister clade of another clade containing T. medemi and T. dorbigni. However, not all currently recognized subspecies of these species are well resolved. Only the two subspecies of T. dor-  bigni are clearly distinct and reciprocally monophyletic. In contrast, both within T. venusta and T. grayi, mtDNA divergences are shallow and the individual subspecies, if being reciprocally monophyletic at all, represent clades with very short basal branches. Within T. venusta, the relationships of the individual subspecies are particularly weakly resolved. A sequence of T. v. cataspila is sister to all remaining taxa. The remaining taxa are placed in a weakly resolved polytomy. Two sequences assigned to T. v. venusta sensu lato (s. l.) from the putatively introduced populations around Acapulco de Juárez (Guerrero; Parham et al. 2015) are paraphyletic with respect to the remaining subspecies of T. venusta. Trachemys venusta callirostris and T. v. chichiriviche and two further clades corresponding to T. v. venusta s. l. (Fig. 2, clades 1 and 2) are monophyletic, but with short basal branches. Clade 1 contains mainly sequences from the putative distribution range of T. v. uhrigi (Honduras), including the holotype, but also from Acapulco and Guatemala. Clade 2 contains two sequences from Veracruz (T. v. venusta sensu stricto) and further sequences from Acapulco. The phylogenetic relationships of these clades are weakly resolved, except for the well-supported sister group relationship of the two South American taxa T. v. callirostris and T. v. chichiriviche. In contrast, sequences of the three subspecies of T. grayi are reciprocally monophyletic.
Our new sequences from Costa Rica cluster together with previously published sequences (Fritz et al. 2012) of T. g. panamensis and T. g. emolli. The only locality on the Pacific side (Terraba Sierpe) represents T. g. panamensis, whereas the remaining three Costa Rican sites (Barra del Colorado, Caño Negro, Gandoca) correspond to T. g. emolli. Remarkably, Barra del Colorado and Gandoca are situated on the Caribbean side of Costa Rica (Fig. 3), where T. v. uhrigi is expected to occur (TTWG 2017, 2021).

Nuclear DNA
Our Splitstree analysis revealed entangled relationships among the individual Trachemys taxa, showing that the studied nuclear loci only incompletely resolve the relationships within Trachemys. Yet, three clusters are visible that represent two clearly distinct groupings (termed "North" and "South" in Fig. 4). One cluster contains the Antillean species T. decorata, T. decussata, T. stejnegeri, and T. terrapen and the three subspecies of the North American T. scripta, i.e., T. s. scripta, T. s. elegans, and T. s. troostii. The only available sample of T. gaigeae is distinct from this cluster, but closer to it than to the remaining samples, and represents a distinct cluster. Trachemys gaigeae and the Antillean and North American species together constitute the grouping "North." The third cluster is equal to the grouping "South" and contains all studied species from Mexico, Central and South America (T. dorbigni, T. grayi, T. medemi, T. venusta). Except for T. v. venusta s. l. and one sample of T. g. panamensis, South American taxa (T. d. dorbigni, T. d. adiutrix, T. medemi, T. v. callirostris, T. v. chichiriviche) are within this cluster distinct from Central American and Mexican taxa (T. g. grayi, T. g. emolli, T. g. panamensis, T. v. cataspila).
When the nuclear sequences are analyzed in PCAs, the described genetic differentiation is also reflected. We present here PCAs for the grouping "South" and for T. grayi and T. venusta. When all species of the grouping "South" are processed together (Fig. 5)   cies of T. dorbigni and T. medemi are removed, T. grayi and T. venusta are more clearly separated (Fig. 6).

Discussion
Our data provide new insights into the distribution of slider turtles in southern Central America. First, we present evidence that Trachemys grayi is not restricted to the Pacific versant of southern Central America as previously thought (TTWG 2017(TTWG , 2021. Mitochondrial and nuclear DNA sequences (Figs 2 and 4) indicate that T. g. emolli also occurs in Costa Rica along the Caribbean coast (Fig. 3). This is supported by photographs of the sampled turtles ( Fig. 7A-D), which morphologically match T. g. emolli from the Lake Nicaragua and Lake Managua drainage in Nicaragua, Honduras, and Costa Rica (cf. Legler 1990;McCranie et al. 2013). The Caribbean populations of T. g. emolli are connected to the previously known distribution range via the Río San Juan, from where the subspecies was known only from the inland section close to Lake Nicaragua (Legler 1990). Our samples from Barra del Colorado are from the mouth of the Río Colorado, a tributary of the Río San Juan. Second, we provide the first genetically substantiated record of T. g. panamensis (Fig. 7E-F) for the Pacific side of Costa Rica. Based on geographic considerations, the occurrence of this taxon along the Pacific coast of Costa Rica was already suggested (TTWG 2017(TTWG , 2021. Third, all slider turtles along the Caribbean versant of southern Central America are currently identified as T. venusta (TTWG 2017(TTWG , 2021, but our records of T. g. emolli from Caribbean Costa Rica imply that the identity of the slider populations in adjacent Nicaragua and Panama should be reexamined. For Caribbean Costa Rica, photographic records from Tortuguero (Moll 1994(Moll , 2010 suggest that T. g. emolli ranges farther southwards. The male figured in Moll (1994) shows a completely broken wide postorbital stripe and the female figured in Moll (2010) has a wide bilobate postorbital stripe, both matching the variation known in T. g. emolli. In contrast, a photograph in Legler (1990) shows a slider turtle from Puerto Cabezas (Caribbean coast of Nicaragua) with a thin continuous postorbital stripe resembling our Figure 1 of T. v. uhrigi, and two turtles from Juan Mina near Colón (Caribbean side of Panama) figured in Moll and Legler (1971: p. 5) have similar narrow head stripes 1 , suggesting that T. g. emolli occurs in neither locality. Without genetic evidence, we can only speculate about the identity of these turtles, even though we wish to highlight that the specimens from Juan Mina could represent T. g. panamensis. Our mitochondrial phylogeny of slider turtles (Fig. 2) is in agreement with two previous publications (Fritz et al. 2012;Vargas-Ramírez et al. 2017) and reveals a well-supported clade for the Central and South American taxa, another one for the four Antillean species, and a third one for the North American species T. gaigeae and T. scripta. These three well-supported clades occur in an unresolved polytomy. Within the Central and South American taxa, T. medemi and the two subspecies of T. dorbigni (T. d. dorbigni and T. d. adiutrix) are the sister clade of T. venusta and T. grayi.
Our nuclear DNA sequences reveal less differentiation than mtDNA, but indicate a clear differentiation between northern and southern taxa. It is noteworthy that the North American T. scripta subspecies (T. s. scripta, T. s. elegans, T. s. troostii) cluster in our Splitstree analysis together with the Antillean species T. decorata, T. decussata, T. stejnegeri, and T. terrapen (Fig. 4). This is in line with shared morphological and ethological traits: In all of these taxa, males have greatly elongated foreclaws playing an important role in the innate stereotypic courtship behavior, the so-called "titillation behavior," while other Trachemys taxa lack long claws in males and this complicated courtship behavior (Fritz 1990(Fritz , 1991Seidel and Ernst 2017). In addition, in both T. scripta and Antillean taxa, the posterior marginal scutes are medially notched, giving the posterior shell margin a serrated appearance, whereas in Mexican, Central American, and South American Trachemys species the marginal scutes lack a pronounced notch, resulting in a more or less smooth posterior rim of the shell. Trachemys gaigeae is the only geographically northern Trachemys species without long claws in males and without pronounced marginal notches. Its distinct placement in our Splitstree (Fig. 4) and the conflicting close sister group relationship of T. gaigeae and T. scripta in our mitochondrial phylogeny (Fig. 2) suggest either a distinct evolutionary trajectory for T. gaigeae or old hybridization with T. scripta and mitochondrial capture-with both possibilities not being mutually exclusive and calling for further research.
Our nuclear markers incompletely resolve the relationships among Antillean species and T. scripta, although the Antillean taxa are in the Splitstree distinct from T. scripta (Fig. 4). For the southern grouping, relationships are even more entangled. There is an incomplete separation between taxa mainly from South America (T. v. callirostris, T. v. chichiriviche, T. medemi, T. d. dorbigni, T. d. adiutrix) versus taxa mainly from Central America and Mexico (T. g. grayi, T. g. emolli, T. g. panamensis, T. v. cataspila), with sequences of T. venusta and T. g. panamensis occurring in both clusters.
For inferring fine-scale relationships among individual Trachemys taxa, more nuclear genomic markers would be needed. To get a complete picture for the genus, also the Mexican taxa that were not included in the present study need to be considered (T. hartwegi, T. nebulosa, T. ornata, T. taylori, T. venusta iversoni, T. yaquia). Based on 15 nuclear loci and a more comprehensive taxon sampling, Thomson et al. (2021) recently provided a global phylogeny for the majority of extant turtle species. Their nuclear phylogeny suggests a sister group relationship between the North American T. scripta and the Antillean Trachemys species, supporting our Splitstree results. However, another sample of T. scripta clusters in their phylogeny remotely with a sample labeled as "Trachemys adiutrix" (= T. d. adiutrix from Maranhão, Brazil) and these two samples together represent the sister clade of Mexican, Central and South American taxa plus T. gaigeae. In addition, a sample labeled as T. venusta unexpectedly appears as sister taxon of T. grayi and not of T. v. callirostris, and T. taylori, a species from northern Mexico not included in the present study, is the sister taxon of T. dorbigni from the Rio de la Plata region of South America. This is unexpected and difficult to explain. These issues suggest that samples might have been misidentified or that an outdated taxonomy was used (T. g. panamensis, now assigned to T. grayi [Fritz et al. 2012;TTWG 2017TTWG , 2021, was formerly classified as a subspecies of T. venusta [McCord et al. 2010;TTWG 2014]). Further inconsistencies are visible in the tree presented in the supplementary information of Thomson et al. (2021), which comprises more samples. It is beyond the scope of the present study to go into details here, but these confusing results highlight the need for further research. Necessary prerequisites for enhancing our knowledge will be not only more informative nuclear genomic markers, but also samples of unambiguous taxonomic identity and known provenance.
Our analyses of nuclear and mitochondrial DNA sequences (Figs 2, 4-6) indicate that T. v. callirostris and T. v. chichiriviche are weakly differentiated from Central American and Mexican subspecies of T. venusta. This supports their conspecificity with T. venusta. With regard to mtDNA sequences, the differentiation within T. venusta is generally weak, and sequences identified as T. v. venusta s. l. from Acapulco are paraphyletic with respect to the South American T. v. callirostris and T. v. chichiriviche plus two clades with sequences of T. venusta from Mexico, Guatemala, and Honduras.
Some of our T. venusta originate from the region of Acapulco, a major tourist destination. The populations in this region are thought to be introduced from elsewhere (Parham et al. 2015). This is supported by our results, because sequences from Acapulco appear in both clades of T. venusta and in the basal polytomy containing all sequences of T. venusta except for one representative of T. v. cataspila. This suggests multiple introductions from different source regions.
One of the two mtDNA clades within T. v. venusta s. l. (clade 1 in Fig. 2) contains a sequence of the holotype of T. v. uhrigi. The validity of this subspecies is debated. Using the same genetic markers as in the present study, Fritz et al. (2012) proposed that T. v. uhrigi should be synonymized with T. v. venusta because they found a sample from Honduras, identified as T. v. uhrigi, only negligibly differentiated from another sample identified as T. v. venusta from Guatemala. The synonymy of T. v. uhrigi and T. v. venusta was supported by McCranie et al. (2013), who revealed little variation among many Trachemys taxa using approximately 750-bp-long mtDNA sequences. McCranie et al. (2013) also argued that the description of T. v. uhrigi was only based on coloration and pattern differences, and that these are more variable than suggested in the original description (McCord et al. 2010). This argument was later followed by Parham et al. (2015), but not by the Turtle Taxonomy Working Group (TTWG 2017(TTWG , 2021, who continued to treat T. v. uhrigi as a valid subspecies. Using more samples, our present mitochondrial phylogeny (Fig. 2) shows that there are two clades plus some unassigned individuals within T. v. venusta s. l.. The type specimen of T. v. uhrigi is placed in the same clade as other samples from Honduras, but also from Guatemala (and from the introduced populations around Acapulco; clade 1 in Fig. 2). The second clade within T. v. venusta s. l. (clade 2 in Fig. 2) comprises sequences from Veracruz and the region of Acapulco. This suggests that all T. venusta from Guatemala and Honduras represent the genetic lineage described as T. v. uhrigi and that Fritz et al. (2012) and later authors erred. This would be in agreement with the tentative assignment of slider turtles from Guatemala and Honduras to T. v. uhrigi in the last edition of the turtle checklist of the Turtle Taxonomy Working Group (TTWG 2021; see also our Fig. 3). However, the situation is not straightforward because the mitochondrial differentiation within T. venusta is generally shallow and our nuclear markers cannot resolve it. Specifically, in our mitochondrial phylogeny (Fig. 2), two other currently recognized subspecies of T. venusta, T. v. callirostris and T. v. chichiriviche, are reciprocally monophyletic sister clades and placed together with the two clades identified as T. v. venusta s. l. in a polytomy that also contains the two unassigned samples of the Acapulco population of T. venusta (D 39071 and D 39077).
We cannot disentangle this intricate situation and refrain from further taxonomic conclusions. However, it should be noted that the subspecific diversity within T. venusta (and of T. grayi with similarly shallow divergences) could be overestimated. These species have conspicuous color patterns, in particular in juveniles. This allows differentiating turtles from different regions and these populations are currently identified as distinct subspecies (TTWG 2021). For instance, the allopatric T. v. callirostris and T. v. chichiriviche from Colombia and Venezuela are morphologically very similar, but differ in the number and shape of their unique snout spots (Pritchard and Trebbau 1984).
For morphologically complex groups of mollusks, beetles, and butterflies, the number of recognized taxa tends to be inflated compared to less conspicuous taxa (Páll-Gergely et al. 2019). It seems likely that a similar taxonomic artifact contributes to the many recognized species and subspecies in the colorful and conspicuously patterned slider turtles and their kin. For instance, Praschag et al. (2017) found in the closely related genus Graptemys less genetically differentiated taxa than currently recognized, and the genetic results presented by Vamberger et al. (2020) for T. scripta suggest that the traditionally recognized three subspecies (T. s. scripta, T. s. elegans, T. s. troostii) merely differ in coloration and pattern and might not even warrant taxonomic distinction in the face of virtually non-existent genetic differences. A similar result was obtained using genomic datasets for T. s. scripta, T. s. elegans, and T. s. troostii (Parham et al. 2020). The hypothesis of taxonomic inflation should be extended to Central and South American slider turtles and should be tested using broader geographic sampling, a morphological assessment, and additional genetic markers.