Are Geckos Special in Sex Determination? Independently Evolved Differentiated ZZ/ZW Sex Chromosomes in Carphodactylid Geckos.

Amniotes possess astonishing variability in sex determination ranging from environmental sex determination (ESD) to genotypic sex determination (GSD) with highly differentiated sex chromosomes. Geckos are one of the few amniote groups with substantial variability in sex determination. What makes them special in this respect? We hypothesized that the extraordinary variability of sex determination in geckos can be explained by two alternatives: 1) unusual lability of sex determination, predicting that the current GSD systems were recently formed and are prone to turnovers; and 2) independent transitions from the ancestral ESD to later stable GSD, which assumes that geckos possessed ancestrally ESD, but once sex chromosomes emerged, they remain stable in the long term. Here, based on genomic data, we document that the differentiated ZZ/ZW sex chromosomes evolved within carphodactylid geckos independently from other gekkotan lineages and remained stable in the genera Nephrurus, Underwoodisaurus, and Saltuarius for at least 15 Myr and potentially up to 45 Myr. These results together with evidence for the stability of sex chromosomes in other gekkotan lineages support more our second hypothesis suggesting that geckos do not dramatically differ from the evolutionary transitions in sex determination observed in the majority of the amniote lineages.

Significance

Geckos, the species-rich group of lizards, show an extensive variety of sex determination systems, leading authors of earlier studies to the assumption that their sex determination systems are unstable and prone to turnovers. For the first time, we uncovered the gene content and investigated the sex chromosome homology of carphodactylid geckos, and we reviewed in a wider phylogenetic content the evolution of sex determination in geckos. We conclude that the observed variability on the sex determination systems in geckos can be explained by multiple and independent transitions from the ancestral environmental sex determination to later stable genotypic sex determination, with sex chromosomes comparably stable in the long term to other reptilian lineages.

Introduction

Sex determination is the process that decides whether an organism will develop as a male or a female. As such, it is of key importance at both individual and population levels since it can also influence population sex ratio. Amniote vertebrates possess two major sex determination systems: environmental sex determination (ESD) and genotypic sex determination (GSD). In species with ESD, the sex is dependent on environmental conditions and there are no systematic differences in genomes between males and females. Temperature is the most commonly studied environmental factor influencing sex under ESD, which is then frequently called temperature-dependent sex determination (TSD) (Bull 1983; Korpelainen 1990; Rhen and Schroeder 2010). In species with GSD, the sex of the developing embryo is determined at conception by sex-specific genotypes, that is, by the combination of sex chromosomes. Within GSD, there are two main sex determination systems: male heterogamety with XX/XY sex chromosomes and female heterogamety with ZZ/ZW sex chromosomes.

Some phylogenetic reconstructions of the evolution of sex determination suggest that ESD might be ancestral for squamate reptiles and even amniotes (Pokorná and Kratochvíl 2009; Gamble et al. 2015; Johnson Pokorná and Kratochvíl 2016; Straková et al. 2020; critically discussed in Kratochvíl et al. 2021). According to these models, GSD evolved independently multiple times either from ESD (i.e., the emergence of sex chromosomes) or GSD (i.e., sex chromosome turnover) ancestors. The emergence of highly differentiated sex chromosomes seems to act as an “evolutionary trap” stabilizing GSD systems in a long term (Pokorná and Kratochvíl 2009; Johnson Pokorná and Kratochvíl 2016). According to this hypothesis, transitions from GSD to ESD should be rare, as they would be connected with the loss of advantages stemming from the presence of sex chromosomes. In amniotes, the long-term stability of sex chromosomes was recently documented by molecular evidence in mammals, birds, softshell turtles, iguanas, caenophidian snakes, anguimorphan lizards, skinks, and lacertid lizards (Mank and Ellegren 2007; Veyrunes et al. 2008; Rovatsos, Altmanová, Johnson Pokorná, et al. 2014; Rovatsos, Altmanová, Pokorná, et al. 2014; Rovatsos, Pokorná, et al. 2014; Rovatsos, Vukić, et al. 2015; Rovatsos et al. 2017; Rovatsos, Rehák, et al. 2019; Altmanová et al. 2018; Kostmann et al. 2021). Transitions from GSD to ESD seem to be extremely rare within GSD lineages in amniotes (Harlow 2004; Gamble et al. 2015; Johnson Pokorná and Kratochvíl 2016; Rovatsos, Vukić, et al. 2019) and were never documented in several other GSD lineages, such as in amphibians or insects (Beukeboom and Perrin 2014).

An alternative hypothesis postulates that the ancestral sex determination system in amniotes was GSD (Ezaz et al. 2017; Singchat et al. 2018) and that GSD to ESD transitions occurred regularly during the amniote evolution (Ezaz et al. 2017; Pennell et al. 2018). However, many of the reported GSD to ESD transitions suggested within lacertids, skinks, varanids, and chameleons appeared to be based on an erroneous assignment of GSD species as ESD (Rovatsos, Johnson Pokorná, et al. 2015; Hill et al. 2018; Nielsen et al. 2018; Iannucci et al. 2019; Rovatsos, Vukić, et al. 2019; Cornejo-Páramo et al. 2020; Kostmann et al. 2021) and it seems that many major clades of amniotes fixed their GSD quite a long time ago (reviewed in Kratochvíl et al. 2021). The variability in sex determination found in some of these lineages, such as snakes (Rovatsos, Vukić, et al. 2015; Gamble et al. 2017; Augstenová et al. 2018), chameleons (Rovatsos, Johnson Pokorná, et al. 2015; Nielsen et al. 2018; Rovatsos, Vukić, et al. 2019), or iguanas (Acosta et al. 2019; Nielsen, Guzmán-Méndez, et al. 2019) concerns turnovers of sex chromosomes, that is, transitions within GSD. Lineages with the well-supported cooccurrence of both GSD and ESD are of particular interest to explore the mechanisms driving the evolution of sex determination systems; however, it seems that they are only a few of them among amniotes: the turtles (Valenzuela and Adams 2011; Bista and Valenzuela 2020), the dragon lizards (Agamidae) (Ezaz et al. 2009; Pokorná and Kratochvíl 2009), and the geckos (Pokorná and Kratochvíl 2009; Gamble 2010; Gamble et al. 2015). In turtles, GSD seems to have evolved at least five times independently from the ancestral ESD, which would support the evolutionary trap models (Johnson Pokorná and Kratochvíl 2016; Mazzoleni et al. 2020). Some authors reconstructed GSD to ESD transitions deeply in turtle history (Valenzuela and Adams 2011; Montiel et al. 2016; Bista and Valenzuela 2020). Nevertheless, there seems to be an agreement on the evolutionary reconstructions in the cryptodiran turtles, where the ancestral ESD turned into GSD four times and GSD stayed stable then. The situation is even more complicated and far from resolved in dragon lizards and geckos, which makes them excellent systems for testing evolutionary scenarios on sex determination.

Geckos (Gekkota) are an ancient and species-rich group of squamates, with more than 2,000 currently described species (Uetz et al. 2020). Their basal split was dated approximately 57–180 Ma (Kumar et al. 2017). Geckos are currently classified into seven families: Carphodactylidae, Diplodactylidae, Eublepharidae, Gekkonidae, Phyllodactylidae, Pygopodidae, and Sphaerodactylidae (Gamble et al. 2008; Pyron 2013). Although sex determination has been studied only in a minority of species (∼3–4%; Pokorná and Kratochvíl 2009; Gamble et al. 2015; Kratochvíl et al. 2021), it is evident that geckos as a whole show extreme variability in sex-determining systems, ranging from ESD to GSD with both male and female heterogamety (Pokorná and Kratochvíl 2009; Pokorná et al. 2014; Gamble et al. 2015). Based on the phylogenetic analyses, it was suggested that sex chromosomes evolved within amniotes around 40 times (Johnson Pokorná and Kratochvíl 2016), out of it around 20 times in geckos (Gekkota) (Gamble et al. 2015).

What makes geckos among amniotes so special in sex determination? We suggest that the variability observed in geckos could be explained by two hypotheses:

Extreme lability of sex determination: according to this hypothesis, geckos are prone to turnovers in sex determination, ESD evolved several times within this group from GSD, and turnovers of sex chromosomes were frequent as well. This hypothesis suggests that geckos are different from the most other amniote groups in the stability of sex determination systems and are closer to fish and amphibian lineages with rapid turnovers in sex determination (Ross et al. 2009; Myosho et al. 2015; Gammerdinger and Kocher 2018; Jeffries et al. 2018; Böhne et al. 2019).

Multiple independent transitions from the ancestral ESD to GSD: this hypothesis assumes that geckos might possess ancestrally ESD (Pokorná and Kratochvíl 2009; Gamble et al. 2015), which was inherited by several extant species, resulting in homologous ESD even among species from distinct families. GSD and hence sex chromosomes would then emerged several times independently from this ancestral ESD in parallel and could have stayed then in each clade with independently evolved GSD stable in a long term. This hypothesis assumes unequal transitions rates between ESD and GSD, and suggests that the situation in geckos is not different from many other amniote lineages, maintaining stable sex chromosomes for 40–165 Myr (Mank and Ellegren 2007; Veyrunes et al. 2008; Rovatsos, Altmanová, Johnson Pokorná, et al. 2014; Rovatsos, Altmanová, Pokorná, et al. 2014; Rovatsos, Pokorná, et al. 2014; Rovatsos, Vukić, et al. 2015; Rovatsos et al. 2017; Altmanová et al. 2018; Rovatsos, Rehák, et al. 2019; Cornejo-Páramo et al. 2020; Kostmann et al. 2021). According to this alternative, the geckos would be a precise analogy for instance to the cryptodiran turtles, where the common ancestor of the extant radiations probably had ESD, and GSD evolved independently several times and remained evolutionary stable then (Johnson Pokorná and Kratochvíl 2016; Straková et al. 2020), just the diversity would be larger in the highly diversified gekkotan lineage.

These hypotheses give different predictions concerning stability and homology of sex determination and sex chromosomes within geckos. To test them, we need a systematic survey of the homology of sex chromosomes based on knowledge of gene content of sex chromosomes across gecko phylogeny. Unfortunately, the current knowledge in this respect is restricted to few gecko lineages (Kawai et al. 2009; Nielsen, Daza, et al. 2019; Rovatsos, Farkačová, et al. 2019; Pensabene et al. 2020; Kratochvíl et al. 2021).

The family Carphodactylidae with the informative phylogenetic position was largely neglected in this respect. Up to now, sex determination is known only in a single species, Underwoodisaurus milii, where female heterogamety with differentiated ZZ/ZW sex chromosomes was revealed by molecular cytogenetic methods (Pokorná et al. 2014). ESD was reported for members of the genus Nephrurus (Brown 2012), but without any rigorous testing. The hypothesis on the multiple independent transitions from the ancestral ESD to GSD in geckos predicts that ESD in this genus is either not supported, or that GSD species of carphodactylids phylogenetically separated by ESD should have nonhomologous sex chromosomes. In order to expand our knowledge on sex determination in geckos and to test this specific prediction, here, we sequenced the whole genome from a male and a female individual of U. milii and a male and a female individual of Saltuarius cornutus, with the aim to uncover gene content of their sex chromosomes by comparative genome coverage analysis (e.g., Vicoso and Bachtrog 2011; Vicoso et al. 2013; Picard et al. 2018), and to validate sex linkage of a subset of genes revealed from the bioinformatic analysis by quantitative real-time PCR (qPCR) (e.g., Rovatsos, Altmanová, Johnson Pokorná, et al. 2014; Rovatsos, Altmanová, Pokorná, et al. 2014; Rovatsos, Vukić, et al. 2015; Rovatsos, Vukić, et al. 2016). The sex linkage of genes was further tested by qPCR in other species from the family Carphodactylidae to explore the homology of sex chromosomes across the phylogenetic spectrum of this gecko lineage and to test homology and stability of GSD across gekkotan lizards.

Results

Identification of Z-Specific Loci by Genome Coverage Analysis

Illumina reads were successfully mapped to 20,156 genes in both sexes (supplementary tables S1, S2, Supplementary Material online). We identified genes as candidate Z-specific, if they fulfilled two criteria: 1) female to male ratio for read coverage between 0.35 and 0.65, and 2) absence of SNPs in at least 80% of the exons per gene. In S. cornutus, we identified 86 candidate Z-specific genes with known homology to the chicken genome: 22 genes have orthologs linked to GGA17, 24 genes to GGA22, and 19 genes to GGA24, whereas 21 genes have orthologs scattered to 12 other chromosomes (fig. 1). In U. milii, we identified 147 candidate Z-specific genes with known homology to the chicken genome: 56 genes have orthologs linked to GGA10, whereas the remaining 91 genes have orthologs scattered to 24 other chromosomes (including six genes with orthologs to GGA17) (fig. 1 and supplementary table S2, Supplementary Material online).

Fig. 1

Log2-transformed female to male ratios of DNA-seq read coverage per gene in Saltuarius cornutus and Underwoodisaurus milii. The Z-specific genes have half female to male read coverage ratio (log2-transformed ratios of ∼−1.00) than autosomal and pseudoautosomal genes (log2-transformed ratios of ∼0.00). The position of gene orthologs in chicken chromosomes are illustrated. Genomic regions with Z-specific genes in S. cornutus (GGA17, GGA22, GGA24) and U. milii (GGA10, GGA17) are indicated by arrows.

Validation of Z-Specific Loci by qPCR and Homology of Sex Chromosomes by qPCR

We tested by qPCR a data set of 4 genes linked to GGA4 (maml3, mbnl3, noct, zdhhc9), 1 to GGA5 (bmf), 12 genes to GGA10 (aen, ap4e1, cep152, chrna5, cilp, dennd4a, fan1, herc1, kif7, mrpl46, rhcg, vps13c), 1 to GGA15 (derl3), 3 to GGA17 (adamts13, st6galnac4, tor2a), 2 to GGA22 (adam9, aebp1), and 1 to GGA24 (snx19) in S. cornutus and U. milii (fig. 2 and supplementary table S3, Supplementary Material online). The qPCR results confirmed the results of the comparative genome coverage analysis. The Z-specific region in S. cornutus consists of genes with orthologs linked to GGA17, GGA22, and GGA24, whereas the Z-specific region in U. milii has gene content with homologs in GGA10 (aen, ap4e1, cep152, chrna5, cilp, dennd4a, fan1, herc1, kif7, mrpl46, rhcg, vps13c) and GGA17 (adamts13, st6galnac4, tor2a). The genes linked to sex chromosomes in the geckos of the genus Paroedura (GGA4: maml3, mbnl3, noct, zdhhc9; GGA15: derl3) (Rovatsos, Farkačová, et al. 2019) and in pygopodid geckos (GGA5: bmf) (Rovatsos et al. 2021) are autosomal or pseudoautosomal in both S. cornutus and U. milii (fig. 2 and supplementary table S3, Supplementary Material online).

Fig. 2

Gene dose ratios between sexes of autosomal (∼1.00) and Z-specific genes (∼0.5), across ten carphodactylids. Missing bars indicate that the specific gene was not successfully amplified by qPCR in the given species. Phylogenetic branching patterns are according to Pyron et al. (2013). All data are presented in supplementary table S3, Supplementary Material online.

The same set of genes was tested by qPCR across five species from the genus Nephrurus and three additional species from the genus Saltuarius. The qPCR test revealed that the species from the genera Underwoodisaurus and Nephrurus share similar gene content in their Z chromosomes, consisting mainly of genes with orthologs to GGA10 and GGA17 (fig. 2 and supplementary table S3, Supplementary Material online). Notably, there is variability on the Z-specificity of the genes from GGA17 region (adamts13, st6galnac4, tor2a), which occasionally appear to be autosomal or pseudoautosomal in species of the genus Nephrurus (fig. 2 and supplementary table S3, Supplementary Material online). On the other hand, all four species of the genus Saltuarius show similar Z-specific gene content (GGA17: adamts13, st6galnac4, tor2a; GGA22: adam9, aebp1; GGA24: snx19). Genes linked to the sex chromosomes of pygopodids (GGA5: bmf) (Rovatsos et al. 2021) and geckos of the genus Paroedura (GGA4: maml3, mbnl3, noct, zdhhc9; GGA15: derl3) (Rovatsos, Farkačová, et al. 2019) are autosomal or pseudoautosomal in all tested carphodactylid species (fig. 2 and supplementary table S3, Supplementary Material online).

Discussion

Gene Content of Sex Chromosomes in Carphodactylids

We uncovered the Z-specific gene content in two species of carphodactylid geckos: the thick-tailed gecko U. milii and the northern leaf-tailed gecko S. cornutus. Sex chromosomes of U. milii seem to share genes with GGA10 and GGA17, whereas sex chromosomes of S. cornutus contain genes from GGA17, GGA22, and GGA24. The tested species from the genera Underwoodisaurus and Nephrurus share homologous ZZ/ZW sex chromosomes, but the tested loci Z-specific in members of these genera are autosomal or linked to poorly differentiated regions of the sex chromosomes in the carphodactylid species from the genus Saltuarius (fig. 2). Notably, the genes Z-specific in several gecko species from the genus Paroedura (Rovatsos, Farkačová, et al. 2019), and X-specific in the pygopodid geckos (Rovatsos et al. 2021), seem not to be Z-specific in the tested carphodactylids (fig. 2 and supplementary tables S2 and S3, Supplementary Material online), indicating that sex chromosomes in the carphodactylid geckos evolved independently from the sex chromosomes of other gekkotan lineages.

Notably, few genes from GGA17 genomic region appear to be Z-specific in both Underwoodisaurus + Nephrurus and Saltuarius, according to both the comparative coverage analysis and the qPCR test for sex linkage. It is possible that the region with homologs linked to GGA17 represents the ancestral carphodactylid sex chromosomes and the regions homologous to GGA10, and GGA22 and GGA24 were later added independently in Underwoodisaurus + Nephrurus and Saltuarius, respectively, into differentiated parts of their ZZ/ZW sex chromosomes. Alternatively, these two lineages evolved sex chromosomes independently and the fragment homologous to a part of GGA17 was a part of the ancestral sex chromosomes in one of these lineages and it was more recently added into the differentiated part of sex chromosomes of the other lineage. We are not able to differentiate between these alternatives due to the small number of shared genes on sex chromosomes of these lineages. The identification of the sex-determining locus in both carphodactylid lineages would be decisive, as homologous sex determination systems should share a sex-determining gene.

The knowledge of sex-linked genes can be used in the qPCR-based method of molecular sexing in members of the genera Underwoodisaurus, Nephrurus, and Saltuarius, as was previously developed for anguimorphan reptiles, caenophidian snakes, iguanas, lacertids, skinks, and trionychid turtles (Rovatsos, Pokorná, Altmanová, et al. 2014; Rovatsos, Vukić, et al. 2015; Rovatsos, Vukić, et al. 2016; Rovatsos et al. 2017; Rovatsos, Rehák, et al. 2019; Rovatsos, Vukić, et al. 2019; Kostmann et al. 2021). Such molecular sexing method can be important for instance in breeding projects and developmental studies requiring knowledge of the sex of embryos.

Evolution of Sex Determination Systems in Geckos

Our survey on sex determination across gekkotan lizards revealed that sex determination is currently known in 63 species of geckos (fig. 3). GSD was reported in species from all seven gecko families by detecting balanced hatchling sex ratios across a range of incubation temperatures, the presence of sex chromosomes (cytogenetic methods), or sex-specific genomic regions (next-generation sequencing methodologies). Nevertheless, species with well-supported ESD are rare. In fact, ESD was well documented only in the genera Eublepharis and Hemitheconyx from the small family Eublepharidae (Pokorná et al. 2010), Phelsuma from the family Gekkonidae and Tarentola from the family Phyllodactylidae (Nettmann and Rykena 1985; Viets et al. 1994; Gamble et al. 2011). ESD was reported for more species after detecting unbalanced hatchlings sex ratios across a range of incubation temperatures, which is the traditional method to detect ESD, but many of these studies are questionable, as they either had a small sample size or did not cover a wide range of temperatures (Gamble et al. 2015). In addition, the hatchling sex ratio may vary even in GSD species, due to factors such as sex-biased embryonic mortality (Eiby et al. 2008), maternal hormones in ovum (Merkling et al. 2018; Firman 2020), and temperature-induced sex reversals (Holleley et al. 2015; Ehl et al. 2017), which makes the identification of ESD species even more complicated. For example, ESD was previously reported in the diplodactylid gecko Correlophus ciliatus due to the higher proportion of females at low incubation temperatures (Seipp and Henkel 2000; de Vosjoli et al. 2003), but GSD with ZZ/ZW sex chromosomes was recently revealed by restriction site-associated DNA sequencing (RAD-seq) (Gamble et al. 2015). Similarly, ESD was reported based on incubation experiments covering three constant temperatures in Gekko japonicus (Tokunaga 1985; Ding et al. 2012), but heteromorphic XX/XY sex chromosomes were reported in the same species by Yoshida and Itoh (1974). Temperature-dependent sex ratio in G. japonicus cannot be fully explained by sex-specific mortality (Tokunaga 1985; Ding et al. 2012). However, the intermediate temperature gave rise to nearly equal sex ratio and none of the tested conditions led to all-female offspring, which is an unusual pattern for an ESD species (Viets et al. 1993; Whiteley et al. 2021). The presence of sex chromosomes together with sex ratio dependent on incubation temperatures could be explained by a thermally induced sex reversal well-documented in the bearded dragon Pogona vitticeps (Holleley et al. 2015; Ehl et al. 2017) and other agamids (Whiteley et al. 2021). In any case, the situation in G. japonicus as well as in the whole genus Gekko deserves further research. Last but not least, ESD was previously reported in the geckos of the genus Nephrurus (Brown 2012), but the identification of sex chromosomes in the current study undoubtedly support the presence of GSD, thus, further expanding the growing list of erroneous reports of ESD in squamate reptiles (Rovatsos, Johnson Pokorná, et al. 2015; Rovatsos, Vukić, et al. 2019; Kostmann et al. 2021).

Fig. 3

Overview of the current knowledge on the sex determination systems across the seven gecko families: Carphodactylidae (A), Diplodactylidae (B), Pygopodidae (C), Eublepharidae (D), Sphaerodactylidae (E), Phyllodactylidae (F), and Gekkonidae (G). Phylogeny according to Rocha et al. (2010), Pyron et al. (2013), and Keating et al. (2021). The methods applied to uncover the sex determination system and the sex chromosome homology to the chicken genome (wherever known) are presented. Note that in Gecko japonicus (highlighted with *), GSD and ESD were contradictory reported by independent studies. Data originate from current and previous studies (Yoshida and Itoh 1974; Branch 1980; Wagner 1980; Solleder and Schmid 1984; Nettmann and Rykena 1985; Tokunaga 1985; McBee et al. 1987; Osadnik 1987; Olmo and Signorino 2005; Kawai et al. 2009; Pokorná et al. 2010; Trifonov et al. 2011; Ding et al. 2012; Matsubara et al. 2013; Koubová et al. 2014; Pokorná et al. 2014; Kasai et al. 2019; Pensabene et al. 2020; King 1978, 1987; Moritz 1984a, 1984b, 1990; Ota et al. 1992, 2001; Viets et al. 1993, 1994; Schmid et al. 2014a, 2014b; Gamble et al. 2015, 2018; Rovatsos et al. 2016a, 2019, 2021; Nielsen et al. 2019b; Keating et al. 2020, 2021).

Our major research question was whether geckos have for amniotes extremely labile sex determination, or whether the pattern is consistent with multiple emergences of GSD from the ancestral ESD. This ambitious question cannot be resolved by the current data, as we still lack the knowledge on the gene content of sex chromosomes in many gekkotan GSD lineages and systematic multiple-species tests on the homology of GSD are still scarce. Moreover, the partial gene content of sex chromosomes is known only in eight gekkotan lineages (fig. 3): in ten species of the carphodactylids (this study), in the limbless pygopodids (Rovatsos et al. 2021), in two species of the eublepharid genus Coleonyx (Pensabene et al. 2020), in the sphaerodactylid genus Aristelliger (Gamble et al. 2015; Keating et al. 2020), in the phyllodactylid species Phyllodactylus wirshingi (Nielsen, Daza, et al. 2019), in the gekkonid Gekko hokouensis (Kawai et al. 2009) and Cyrtodactylus pharbaungensis (Keating et al. 2021), and in six species of the genus Paroedura (Gekkonidae) (Rovatsos, Farkačová, et al. 2019). However, it is notable that the sex chromosomes of the studied gekkotan lineages are homologous to different genomic regions with the exception of 1) Nephrurus-Underwoodisaurus and Cyrtodactylus pharbaungensis, with sex chromosomes partially homologous to GGA10; and 2) Phyllodactylus wirshingi and Gekko hokouensis, with sex chromosomes homologous to GGAZ (reviewed in fig. 3). However, the members of these groups are mutually phylogenetically distant and separated by numerous lineages with other sex determination systems (fig. 3). Therefore, we assume that these groups of geckos probably evolved their sex chromosomes from autosomes by convergence, that is, by independent cooption of either the GGA10 or the GGAZ genomic region for the role of sex chromosomes, rather than that they inherited these sex chromosomes directly from their common ancestors, who shared the same sex chromosomes (reviewed in Kratochvíl et al. 2021). This hypothesis should be tested during further research on gekkotan sex determination.

All gekkotan lineages with knowledge on partial gene content separated by ESD seem to have clearly nonhomologous sex chromosomes (fig. 3). The phylogenetic reconstruction of the evolution of sex determination in geckos does not provide strong support for a GSD to ESD transition. In the same context, GSD seems to be evolutionary stable in long term in geckos, comparable to other amniote lineages (for a review on the age of sex chromosomes, see Kratochvíl et al. 2021). Such long-term stability was recently documented in pygopodids (30–72 Myr) (Rovatsos et al. 2021), the geckos of the genus Paroedura (62–90 Myr) (Rovatsos, Farkačová, et al. 2019), the eublepharids Coleonyx elegans and Coleonyx mitratus (34 Myr) (Pensabene et al. 2020), and the sphaerodactylids of the genus Aristelliger (Keating et al. 2020). Even in carphodactylids, the genera Nephrurus and Underwoodisaurus share homologous sex chromosomes for at least 15–36 Myr and the genus Saltuarius for 16 Myr, and if these two lineages share a homologous sex determination system, it will be dated 29–46 Myr old. In fact, any potential instability of sex determination is restricted in three genera: Gekko, Hemidactylus, and Cyrtodactylus, but even in these lineages, it seems that particular sex chromosome systems can be several dozens of million years old (Keating et al. 2021) and still far more stable in comparison to nonamniotic lineages, such as the true frogs of the genus Rana, where at least six sex chromosome turnovers were documented, resulting to novel sex determination systems with age less than 10 Myr (Jeffries et al. 2018), and the frogs of the genus Xenopus (Roco et al. 2015).

In conclusion, we believe that geckos perhaps do not have for amniotes extreme lability of sex determination systems, but instead, the variability in their sex determination can be explained by their high species diversity, high age of the whole lineage, and potentially ancestral ESD, prone to switch to GSD. Geckos are often viewed as a single homogenous group of lizards by nonspecialists, but the extant gekkotan families represent very old radiations, with the basal splits estimated to 57–180 Ma. Each of the gekkotan families is thus comparable in the age to other amniote lineages with the variability in sex determination systems, such as snakes (75–165 Ma) and cryptodiran turtles (97–250 Ma). In addition, the past imprecise records on sex determination systems in various lineages, driven by misidentification of ESD, artificially increased the observed variability of sex determination in geckos. It seems that instead of large lability, many gekkotan lineages show long-term stability of sex determination systems. In the current contribution, we explicitly formulated these working hypotheses and show that the emerging, still incomplete picture in this important diversified group seems to support the latter possibility, but we call for further experimental testing.

Materials and Methods

Studied Material and DNA Isolation

We collected blood and/or tip of tail samples from both sexes in ten species of carphodactylid geckos, obtained from the pet trade (supplementary table S1, Supplementary Material online). DNA was isolated by the DNeasy Blood and Tissue Kit (Qiagen), using the manufacturer’s protocol. DNA concentration and quality were measured by NanoDrop ND‐2000 spectrophotometer (Thermo Fisher Scientific Inc.). All methods were carried out in accordance with relevant Czech and EU guidelines and regulations (ethical committee permit MSMT-8604/2019-7), and by researchers accredited for animal experimental design (accreditations L.K.: CZ02535, M.R.: CZ03540).

DNA-Seq and Coverage Analysis in U. milii and S. cornutus

Genomic DNA from one male and one female of the thick-tailed gecko U. milii and the northern leaf-tailed gecko S. cornutus were sequenced at high coverage by Novogene (Cambridge, UK) on Illumina HiSeq2500 platform, with 350 base pairs (bp) pair-end option (DNA-seq). The raw Illumina reads are deposited in NCBI database under the BioProject PRJNA701686. Adapters and low-quality bases from raw reads were trimmed by Trimmomatic (Bolger et al. 2014) and reads shorter than 50 bp were removed. Trimmed reads were checked in FASTQC (Andrews 2010) and MULTIQC (Ewels et al. 2016).

The differences in the gene copy numbers between sexes are expected to be proportional to the differences in coverage of the reads from DNA-seq in Illumina HiSeq platforms (Vicoso et al. 2013; Picard et al. 2018). Therefore, the Z-specific loci are expected to have half-read coverage in ZW females in comparison to ZZ males in ZZ/ZW sex determination systems with degenerated W chromosomes, whereas autosomal, pseudoautosomal, and poorly differentiated loci should have equal read coverage in both sexes. We mapped independently the trimmed Illumina reads from the male and the female specimen to a reference data set consisting of 170,981 exons, extracted from the G. japonicus genome project (Liu et al. 2015), using Geneious Prime 2020.0.4 (https://www.geneious.com). The average read coverage per gene was calculated in each specimen after filtering all exons with unexpectedly high or low coverage (3-fold difference from the average read coverage of each specimen). In S. cornutus, the average read coverage was 38× in male and 41× in female, whereas in U. milii, it was 31× in the male and 28× in the female (supplementary table S2, Supplementary Material online). Subsequently, we calculated the ratio of female to male read coverage for each gene, normalized to the total number of assembled reads per specimen (see Vicoso et al. 2013) (supplementary table S2, Supplementary Material online). We recently applied the same methodology to uncover the sex chromosome gene content in the Yucatán-banded gecko C. elegans (Pensabene et al. 2020), the common sandfish Scincus (Kostmann et al. 2021), and the Florida softshell turtle Apalone ferox (Rovatsos and Kratochvíl 2021).

In addition, Z-specific, single-copy genes are hemizygous in the heterogametic sex, therefore, such loci should not have single-nucleotide polymorphisms (SNPs) in our map-to-reference assembly from the female specimen. Therefore, we calculated the presence/absence of SNPs per exon and gene in the assembly of the females, in order to use this information for further validation of the Z-specific genes that were revealed from the comparative coverage analysis.

The chromosome level assemblies of the green anole Anolis carolinensis (Alföldi et al. 2011), the common wall lizard Podarcis muralis (Andrade et al. 2019), and the chicken Gallus gallus (Warren et al. 2017) were used to explore the gene homology for genome-wide cross-species comparisons (supplementary table S2, Supplementary Material online).

Primer Design and Validation of Z-Specific Genes by qPCR

We designed primers to amplify products of 120–200 bp in size in Primer‐Blast software (Ye et al. 2012) using Primer3 (Untergasser et al. 2012) from exons which were revealed to be Z-specific from the comparative coverage analysis either in U. milii or S. cornutus. Prior to qPCR, all primers were tested by standard PCR in DNA templates from U. milii and S. cornutus to select primer pairs that do not amplify secondary products.

We designed primers for 18 candidate Z-specific genes to be tested by qPCR (adam9, adamts13, aebp1, aen, ap4e1, cep152, chrna5, cilp, dennd4a, fan1, herc1, kif7, mrpl46, rhcg, snx19, st6galnac4, tor2a, vps13c) (supplementary table S3, Supplementary Material online). These genes have homologs to the genomic regions GGA10, GGA17, GGA22, and GGA24, which are involved on the sex chromosomes of either U. milii or S. cornutus (fig. 1 and supplementary tables S2 and S3, Supplementary Material online). In addition, we also tested by qPCR five genes homologous to GGA4 (bmf, maml3, mbnl3, zdhhc9) and to GGA15 (derl3) which were Z-linked in several species of geckos from the genus Paroedura (Rovatsos, Farkačová, et al. 2019) and a single gene homologous to GGA5 (noct), which is X-specific in pygopodid geckos (Rovatsos et al. 2021) (supplementary table S3, Supplementary Material online), in order to explore the potential homology or cooption of sex chromosomes among lineages of geckos. The genes mecom and rag1 were used for the normalization of the qPCR values (Rovatsos, Farkačová, et al. 2019) (supplementary table S3, Supplementary Material online).

We used a qPCR method to calculate the relative gene copy number variation between the male and female genome and to test the Z-specificity of the candidate Z-specific genes (Rovatsos, Altmanová, Pokorná, et al. 2014; Rovatsos, Pokorná, et al. 2014; Rovatsos, Vukić, et al. 2015; Rovatsos, Vukić, et al. 2016; Rovatsos et al. 2017; Nielsen, Guzmán-Méndez, et al. 2019; Rovatsos, Farkačová, et al. 2019; Rovatsos, Rehák, et al. 2019; Rovatsos, Vukić, et al. 2019). With the same reasoning as for the comparative genome coverage, males (ZZ) have double copies of Z-specific genes compared with females (ZW) in species with degenerated nonrecombining W chromosomes. Therefore, the female to male ratio (r) in gene copy number is expected to be 0.5 for the Z-specific genes, 1.0 for autosomal or pseudoautosomal autosomal genes, and 2.0 for the X-specific genes. For the qPCR test, we tested each primer pair and per specimen in triplicates. The test was performed using LightCycler II 480 (Roche Diagnostics). The cycler conditions and the protocol for the qPCR mix were previously described in Rovatsos, Altmanová, Pokorná, et al. (2014).

qPCR Test of Sex Chromosome Homology

The primers for the candidate Z-specific genes in S. cornutus and U. milii were subsequently tested by qPCR method to explore the homology of sex chromosomes in eight additional carphodactylid species from the genera Nephrurus and Saltuarius (supplementary table S3, Supplementary Material online).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.