Tracing Genetic Divergence and Phylogeographic Patterns of Gekko gecko Linnaeus, 1758 (Squamata: Gekkonidae) Across Southeast Asia Using RAG1 Sequence

TL;DR

The study reveals genetic diversity in tokay geckos across Southeast Asia, suggesting multiple hidden species within what was thought to be one.

Contribution

The study identifies three distinct genetic lineages in Gekko gecko using RAG1 sequences, supporting the need for integrative taxonomy.

Findings

Analysis of RAG1 sequences revealed 20 haplotypes and strong regional genetic structuring in Gekko gecko.

Three distinct evolutionary lineages (A–C) were identified, supported by species delimitation tests.

Color morphs (red- and black-spotted) do not correspond to genetic differentiation at the RAG1 locus.

Abstract

The tokay gecko (Gekko gecko) is a colorful lizard found in Southeast Asia and is widely used in traditional medicine. In this study, we analyzed DNA from tokay geckos in Thailand, Laos, and Cambodia to understand their genetic differences. We found clear genetic variation between populations and discovered three main genetic lineages. These results suggest that what is thought to be one species may include hidden, distinct groups. This information is important for protecting the species and managing its trade.

The tokay gecko (Gekko gecko) is a widely distributed lizard species in Southeast Asia, with significant importance in traditional medicine and the pet trade. Previous studies using mitochondrial DNA sequences revealed extensive genetic variation across its range, indicating the presence of distinct evolutionary lineages. In this study, we assessed the nuclear genetic variation and phylogenetic pattern of G. gecko using the recombination activating gene 1 (RAG1). We analyzed 105 RAG1 sequences from 16 localities across Thailand, Laos, and Cambodia, along with additional sequences from GenBank. Sequence analysis revealed 20 variable sites and 20 haplotypes (TgR1–TgR20). Haplotype network and phylogenetic analyses revealed strong regional structuring and at least three distinct evolutionary lineages (A–C), supported by the species delimitation test (PTP). Both red- and black-spotted morphs were present in different clades, indicating that external coloration does not correspond to genetic differentiation at this locus. Our results support the presence of distinct evolutionary lineages in G. gecko and emphasize the importance of integrative taxonomy for accurate species delimitation. These findings have implications for conservation, sustainable management, and regulation of international trade in this commercially exploited species.

Figures

Full text

1. Introduction

The tokay gecko (Gekko gecko), a large nocturnal lizard species in the family Gekkonidae, is widely distributed across tropical and subtropical Asia, including much of Southeast Asia. Additionally, it has been introduced to several regions outside its native range, such as the United States, Mexico, Puerto Rico, Martinique, and Hong Kong [1]. Beyond its ecological role, G. gecko is highly valued in traditional medicine across Asia, where it is believed to possess therapeutic properties against asthma, tuberculosis, cancer, diabetes, skin disorders, and even HIV/AIDS [2]. Demand is intense: over 1.2 million dried individuals are exported annually from Java [3], and an estimated 2–5 million individuals are exported from Thailand each year [4]. Additional sourcing occurs in Cambodia and Laos [5,6], with most trade directed toward China, Hong Kong, and Malaysia for use in traditional Chinese medicine [7].

Two main color morphs of G. gecko are recognized, namely the black-spotted form, occurring mainly in Guangxi and Yunnan, southern China, and the red-spotted form, which is widespread across Asian countries [8]. Morphological, behavioral, and genetic differences between these two morphs have been reported, including differences in coloration [9], vocalizations [10], karyotypes [8], RAPD profiles [11], allozymes [5], and both nuclear and mitochondrial DNA markers [6,9,12,13].

In mainland Southeast Asia, studies using multi-locus enzyme electrophoresis and mitochondrial DNA sequences have demonstrated considerable intra- and inter-population diversity, with evidence for multiple genetic clades [5,6]. For example, six haplogroups (G1–G6) from various localities in Southeast Asia have been identified based on mitochondrial tRNA-Gln/tRNA-Met/partial ND2 sequence, which were grouped into five clades (A–E) [6]. Genetic breaks corresponded to major mountain ranges (e.g., Phetchabun, Dong Paya Yen, Phi Pan Nam) and river barriers (e.g., Mekong), suggesting that geographic barriers strongly influence lineage divergence [5,6]. These findings demonstrate distinct evolutionary lineages and possible species-level distinctions, underscoring the need for integrative taxonomy.

While mitochondrial DNA provides high resolution for detecting recent divergence, it represents only the maternal lineage. Nuclear markers provide complementary biparental inheritance and are critical for resolving deeper evolutionary history. The recombination activating gene 1 (RAG1) is a single-copy nuclear coding gene, which encodes an essential enzyme for initiating V(D)J (variable, diversity, and joining gene segments) recombination during lymphocyte development. It has been proven useful in vertebrate phylogeny and in exploring population structure and divergence in fishes and reptiles due to its moderate mutation rate, low paralogy, biparental inheritance, and conserved coding regions [14]. Previous studies have successfully applied RAG1 in reptiles, including G. gecko [1,12].

In this study, we analyze RAG1 sequences from G. gecko across Thailand, Lao PDR, and Cambodia, together with additional GenBank records, to investigate genetic divergence, population structure, and phylogenetic pattern. Specifically, we test for phylogeographic structuring, evaluate congruence with mitochondrial studies, and assess potential cryptic species using multiple delimitation methods. These findings aim to clarify the evolutionary history of G. gecko and provide a framework for conservation and sustainable trade management.

2. Materials and Methods

2.1. Sample Collection

Tokay gecko samples from 13 localities in Thailand, Laos, and Cambodia were obtained from a previous study [6]. An additional 23 samples were newly collected from three localities in Thailand and Laos (Table 1 and Figure 1). Geckos were captured using the fishing pole method [15]. Prior to tissue collection, individuals were anesthetized with isoflurane [16] to minimize stress and pain. Each gecko was placed in an induction chamber containing a cotton pad soaked with 3–5% isoflurane. Anesthesia was confirmed by the absence of the righting reflex and reduced response to tactile stimuli. A tail-tip sample (~1–2 mm) was excised aseptically and preserved in 80% ethanol. After recovery, all geckos were released back into their natural habitats. Samples were transported to the Walai Rukhavej Botanical Research Institute, Mahasarakham University, and stored at room temperature until DNA extraction.

2.2. DNA Extraction and PCR Amplification

Genomic DNA was extracted from tail tissue using the E.Z.N.A.^®^ Tissue DNA kit (Omega bio-tek, Norcross, GA, USA) following the manufacturer’s protocol. An estimated 700 bp fragment of RAG1 gene was amplified using primers RAG1-F (5′-CCA GAG GAA GTT CAG CAG TGT C-3′) and RAG1-R (5′-GCT TCC AAC TCA TCA GCT TGT C-3′) [12]. PCR conditions were: 95 °C for 3 min, followed by 35 cycles of 94 °C for 35 s, 57 °C for 45 s, and 72 °C for 90 s, with a final extension step at 72 °C for 5 min. Amplicons were purified with an E.Z.N.A.^®^ Gel Purification Kit (Omega bio-tek, Norcross, GA, USA). The purified PCR products were sent for DNA sequencing at ATGC Co; Ltd., (Pathum Thani, Thailand).

2.3. DNA Sequence Analysis

All 553 bp RAG1 sequences generated in this study were aligned using the ClustalW program version 2.0 [17], and the variable sites between haplotypes were compared in the BioEdit program version 7.2.5 [18]. Molecular diversity indices and haplotype data analysis were calculated using the DnaSp program version 5 [19]. The genetic difference between populations was calculated based on p-distance and Kimura 2-parameter (K2P) distance [20] using the program MEGA version 11 [21]. A minimum-spanning haplotype network was constructed in the PopART version 1.7 (https://popart.maths.otago.ac.nz/, accessed on 30 September 2025) based on a median-joining network [22] using all sequences generated in this study. Analysis of Molecular Variance (AMOVA) was conducted using the Arlequin program version 3.5.2.2 [23]. Furthermore, to investigate the genetic structure across geographical barriers within the study area in this study (Thailand, Laos, and Cambodia), a Spatial Analysis of Molecular Variance (SAMOVA), to cluster the RAG1 sequences into genetically and geographically homogeneous groups, was conducted in SAMOVA version 1.0 [24]. The analyses were performed for K = 2–8 with 1000 simulated annealing steps from 100 random starting conditions. Due to incomplete coordinate information, the RAG1 sequences of G. gecko from NCBI GenBank were excluded.

2.4. Phylogenetic Tree Analysis

Phylogenetic relationships among G. gecko were constructed using Bayesian inference (BI), Maximum likelihood (ML), and Neighbor joining (NJ) methods based on RAG1 sequences. The best-fit evolutionary model (HKY + G) was identified under the corrected Akaike information criterion (AIC) using MrModeltest 2.4 [25]. For Bayesian Inference (BI), the nucleotide substitution model was set in Mesquite v.4.01 [26] with nst = 2 and rates = gamma. The analysis was then performed in MrBayes v.3.2.7 [27] using two parallel runs of four Markov chains for 30 million generations, sampling every 1000 generations. Convergence was assessed in Tracer 1.7.2 [28], with the first 25% of trees were discarded as burn-in. Posterior probabilities were estimated from a 50% majority-rule consensus tree. ML and NJ analyses were conducted in MEGA XI [21] with 1000 bootstrap replicates, using the HKY + G model and Kimura 2-parameter (K2P), respectively. To provide phylogenetic context, representative sequences from related gecko taxa (Cyrtodactylus, Dravidogecko, Agaura, Tropiocolotes, and Hemidactylus) were included, with Cnemaspis tucdupensis designated as the outgroup.

2.5. Species Delimitation Analysis

Species boundaries were evaluated using three single-locus species delimitation methods in this study on G. gecko RAG1 sequences. Automatic Barcode Gap Discovery (ABGD) [29], Assemble Species by Automatic Partitioning (ASAP) [30], and the tree-based Poisson tree (PTP) [31]. The PTP method was applied to RAG1-based phylogeny using the web server (https://mptp.h-its.org/#/tree/; accessed on 28 July 2025). ABGD and ASAP analyses were implemented in the bioinformatic toolkit iTaxoTools v0.1 [32]. For ABGD, the Kimura (K80) substitution model was employed, with the default maximum and minimum intraspecific distances (Pmax = 0.1, Pmin = 0.001). The barcode gap width was set to 1.5. A non-recursive partition was selected, with a prior maximum intraspecific divergence (p) = 2.8 × 10^−3^. For ASAP, the partition with the lowest ASAP score and an appropriate threshold distance (dT) was chosen under the Kimura (K80) model, using a default transition/transversion ratio of 2.0 and a probability of 0.01.

3. Results

3.1. Phylogenetic Tree and Species Delimitation

Phylogenetic analyses based on RAG1 sequences consistently recovered three well-supported clades (A–C) with species boundaries further evaluated using the Poisson Tree Processes (PTP) method (Figure 2). Clade A contained the majority of individuals, including the samples from Thailand, Myanmar, China, and introduced populations in the USA. Clade B consists of a sample of black-spotted from China and represents a highly supported, genetically distinct group. Clade C comprises individual samples from Thailand, Vietnam, China, Malaysia, and the USA. Clades A and C include both red- and black-spotted tokay geckos, suggesting that coloration alone does not reflect deep genetic divergence.

The genetic distance-based species delimitation analyses using ABGD and ASAP revealed a non-recursive partition of two genetic clades and a recursive partition of three genetic clades, respectively. The barcode gap analysis from both ABGD and ASAP methods indicated multiple narrow intervals within the range of approximately 0.01–0.00%. In ABGD, partitions 1 and 2, with prior maximum intraspecific divergence of p = 1.0 × 10^−3^ and p = 1.7 × 10^−3^, respectively, produced over-splitting results with 20 genetic groups. Therefore, partition 3 (p = 2.80 × 10^−3^), which yielded two genetic groups, was considered the most appropriate for further interpretation. Unlike ABGD, which relies on manually defined prior P, ASAP ranks partitions using a scoring system. The best partition, with the lowest ASAP score (ASAP score = 2.00; dT = 0.0036), identified three genetic groups. However, this result was not fully consistent with the phylogenetic tree analyses (ML, NJ, and BI) (Figure 2). In contrast, the tree-based species delimitation method (PTP) inferred three genetic clades based on branching rates, and its partition pattern was the most congruent with phylogenetic analysis. Therefore, the PTP-inferred clades were the putative species boundaries for subsequent analyses.

3.2. Haplotype Network

A total of 20 haplotypes (TgR1–TgR20) were identified from the RAG1 sequences based on 20 variable nucleotide sites. The median-joining haplotype network revealed three major haplogroups (Figure 3), which corresponded to the species delimitation (PTP) analysis (Figure 2). Haplotypes TgR1 and TgR2 represented central, widely distributed haplotypes. Among them, TgR1 was the most prevalent and geographically widespread, occurring in samples from Laos, Thailand (north, northeastern, and south), and Cambodia. Several other haplotypes (e.g., TgR4, TgR5, TgR6, TgR10, TgR13, TgR15, TgR18–TgR20) radiated from TgR1, indicating localized mutations or population differentiation. Conversely, TgR2 was predominantly composed of sequences retrieved from GenBank, particularly from the USA, Indonesia, and Thailand (Phuket), along with contributions from other Southeast Asian countries (e.g., Vietnam, Philippines, and Malaysia). Several unique haplotypes (e.g., TgR7 and TgR14) were also derived from TgR2, suggesting a divergent clade from that found in mainland Southeast Asia. Notably, TgR11 was shared among Thailand (northeastern), Vietnam, and the Philippines, indicating potential genetic connectivity among these regions. In contrast, TgR17 formed a separate and distinct clade (Clade B).

3.3. Spatial Analysis of Molecular Variance (SAMOVA)

The SAMOVA analysis based on RAG1 data revealed clear genetic structuring of G. gecko across mainland and island populations. The proportion of variation among groups (FCT) was relatively high, ranging from 0.69525 at K = 2 to 0.75897 at K = 8 (Table S2), and all F-statistics for each independent run were significant (p < 0.01). Although the maximum FCT occurred at K = 8, this value indicated possible over-fragmentation. Instead, genetic structure reached an optimal plateau at K = 5 (FCT = 0.73889, p < 0.001), where the variation among populations within groups (FSC) sharply declined (FSC = 0.00115, p < 0.001) before turning negative at K = 6 (FSC = −0.06172, p < 0.001) (Table S2).

At K = 5, the analysis identified five geographically coherent genetic clusters (Figure 4). Group A was confined to SK in northeastern Thailand, corresponding to the Central Highlands. Group B encompassed populations from Laos (VV, SV, and CS) and northeastern Thailand (KS), corresponding mainly to the Central Highlands and Annamite Range. Group C was restricted to NR in northeastern Thailand, while Group D comprised populations from Cambodia (CS) and northeastern Thailand (SR). Finally, Group E represented the isolated population from Laos (VT), northern Thailand (CR, LP, LN, and UT), and southern Thailand (RG, PG, and PN). The distribution of these groups broadly matched major physiographic regions, with mountain ranges and river systems (e.g., the Mekong River) serving as likely barriers to gene flow.

3.4. RAG1 Sequence Variation

A total of 105 RAG1 sequences, each with a length of 553 bp, were analyzed across 16 geographical localities examined in this study and deposited in GenBank under accession numbers PX432629–PX432691, along with 6 additional localities retrieved from GenBank. Sequence comparison revealed 20 variable nucleotide sites (3.6%), with 13 variable amino acid positions after translation (Table S1). The number of segregating sites (S) varied from 5 to 13. The number of haplotypes (h) ranged from 6 to 13 among populations, with several populations exhibiting unique haplotypes. Haplotype diversity (Hd) ranged from 0.396 ± 0.075 to 0.865 ± 0.030, with an average of 0.772 ± 0.039, and nucleotide diversity (π) values ranged from 0.0009 ± 0.0002 to 0.0030 ± 0.0003, with an average of 0.0049 ± 0.0003 (Table 2).

3.5. Genetic Differences and Genetic Structure

Pairwise genetic distances among G. gecko populations from the three genetic clades (A–C) were calculated using both the Kimura 2-Parameter (K2P) model (lower triangle) and uncorrected p-distance (upper triangle) based on RAG1 sequences (Table 3). The lowest genetic distances were observed between clade C and the other clades (p-distance = 0.0082; K2P = 0.0082–0.0083), whereas the highest genetic distances were found between clade A and clade B (p-distance = 0.0110; K2P = 0.0111).

Analysis of molecular variance (AMOVA) based on RAG1 sequences revealed pronounced genetic structuring among G. gecko populations across different geographic scales (Table 4). When populations were grouped into different countries, namely Thailand, Laos, Cambodia, and USA, 5.26% of the variation was found among groups (FCT = 0.05257, p > 0.05), 56.39% among populations within groups (FSC = 0.53569, p < 0.01), and 48.87% within populations, resulting in an overall FST of 0.51128 (p < 0.01). When populations were grouped according to the SAMOVA-defined clusters (Groups A–E), the proportion of genetic variation among groups increased substantially to 73.89%, with a high and significant fixation index (FCT = 0.73889, p < 0.01), suggesting strong genetic structuring among these groups. In contrast, only 0.03% of the variation occurred among populations within groups, as indicated by a low fixation index (FSC = 0.00115, p < 0.01), while 26.08% of the variation was found within populations (FST = 0.73919, p < 0.01). These results demonstrate strong genetic structuring and restricted gene flow among the SAMOVA-defined groups of G. gecko, highlighting substantial genetic divergence among geographically separated populations.

4. Discussion

This study reveals significant nuclear genetic divergence and population structure in G. gecko across Southeast Asia using the RAG1 gene. The previous studies suggest the presence of multiple evolutionary lineages and possible cryptic species of G. gecko in Southeast Asia [5,6]. These findings complement and extend earlier mitochondrial studies that uncovered deep population structure and highlighted the influence of geographic barriers in shaping genetic diversity [6]. The mitochondrial ND2 study by Saijuntha et al. [6] analyzed 166 red-spotted tokay geckos from 23 localities and identified six major haplogroups (G1–G6), which were associated with five clades (A–E). Strong phylogeographic structuring was observed, with genetic breaks corresponding to mountain ranges such as Phetchabun, Dong Paya Yen, and Phi Pan Nam, and rivers like the Mekong. These natural features appear to have restricted gene flow, leading to regional lineage divergence. Our nuclear RAG1 results parallel these findings, revealing high genetic differences, especially among Thailand, Laos, and Cambodia, and identifying distinct clades with limited haplotype sharing.

While mitochondrial markers such as ND2 and Cyt-b often provide higher resolution for detecting recent divergence, they are maternally inherited and may not reflect the full evolutionary history of a species [6,9,12,13]. Nuclear markers like RAG1 provide independent biparental information and are critical for resolving deeper phylogenetic relationships and confirming species boundaries [1,12]. In our study, the combination of RAG1 haplotype network, phylogenetic tree, and species delimitation (PTP) analyses confirmed regional divergence and supported the presence of three major RAG1 clades, partially corresponding to the mitochondrial clades found in previous studies. Additional support for this pattern comes from Fieldsend et al. [1], who reported two distinct mito-nuclear lineages in G. gecko populations introduced to Florida, matching lineages from the native range. These lineages showed reproductive isolation, suggesting they may represent separate species. Likewise, Wang et al. [12] and others have noted consistent nuclear-mitochondrial differentiation between red- and black-spotted morphotypes, adding weight to the hypothesis of cryptic speciation.

However, in this study, the RAG1 gene did not clearly distinguish between the black- and red-spotted morphotypes of G. gecko, suggesting low genetic differentiation at this nuclear locus between these color forms. The black-spotted form is generally smaller and darker, whereas the red-spotted form is larger with red markings. Although morphological, karyotypic, and genetic differences have led some studies to propose their recognition as separate species, their taxonomic status remains controversial [33]. This suggests that the observed color polymorphism may not be strongly associated with genetic divergence but may instead reflect phenotypic plasticity, recent divergence, ongoing gene flow, or selection related to ecological conditions. Previous research has also indicated that evolutionary transitions to color polymorphism are often associated with increased diversification rates in lizards [34]. Furthermore, a recent study showed that the black- and red-spotted forms of G. gecko occupy distinct ecological niches, with niche differentiation primarily influenced by climatic factors such as isothermality, temperature seasonality, and precipitation during the warmest quarter. These results suggest that ecological divergence, together with morphological variation, may contribute to maintaining their separate distributions and potential taxonomic distinction [33].

Nuclear markers such as RAG1, which evolve more slowly than mitochondrial genes, often lack resolution for detecting recent or shallow divergences, particularly within morphologically variable but genetically cohesive lineages. Interestingly, the presence of a single black-spotted individual from China [12] assigned to clade B, distinct from other sampled individuals, may indicate localized differentiation, introgression from a divergent lineage, or even a cryptic species. This unique placement suggests that while most RAG1 sequences are highly conserved across morphotypes, some degree of regional or population-specific divergence may exist. However, limited sample sizes for both black- and red-spotted G. gecko in certain populations reduce the resolution of our molecular analyses. Consequently, these results should be interpreted with caution, and expanded sampling across broader geographic ranges is needed to better assess genetic differentiation and the potential influence of color polymorphism on population structure.

Furthermore, the tokay gecko’s heavy exploitation for traditional medicine and the pet trade underscores the importance of understanding its genetic structure. Millions of individuals are harvested annually from Southeast Asia [3], raising concerns about the loss of unique genetic lineages. If G. gecko comprises a complex of cryptic species, unregulated harvest may disproportionately impact certain lineages and drive local extinctions. Moreover, international trade regulations, including CITES listings, may need to consider lineage-level distinctions for effective protection. In light of these findings, we recommend that future studies use multi-locus nuclear markers, broader geographic sampling (particularly in peninsular Thailand and Malaysia), and morphological/ecological comparisons to clarify species limits. Conservation strategies should also consider the presence of genetically distinct populations when developing region-specific management plans.

5. Conclusions

This study highlights significant genetic divergence and population structuring of G. gecko across Southeast Asia based on nuclear RAG1 sequence analysis. The identification of three major genetic lineages, supported by robust statistical approaches, such as AMOVA, SAMOVA, phylogenetic tree, and species delimitation analyses, points to distinct evolutionary lineages within this widely distributed species. Our findings underscore the need for an integrative taxonomic approach combining nuclear and mitochondrial markers, morphology, and ecological data to accurately assess species boundaries. These results contribute valuable information for conservation planning and regulatory frameworks, particularly in light of ongoing habitat loss and intense commercial harvesting of G. gecko.