Geckos are among the most species-rich reptile groups and the sister clade to all other lizards and snakes. Geckos possess a suite of distinctive characteristics, including adhesive digits, nocturnal activity, hard, calcareous eggshells, and a lack of eyelids. However, one gecko clade, the Eublepharidae, appears to be the exception to most of these ‘rules’ and lacks adhesive toe pads, has eyelids, and lays eggs with soft, leathery eggshells. These differences make eublepharids an important component of any investigation into the underlying genomic innovations contributing to the distinctive phenotypes in ‘typical’ geckos.
We report high-depth genome sequencing, assembly, and annotation for a male leopard gecko, Eublepharis macularius (Eublepharidae). Illumina sequence data were generated from seven insert libraries (ranging from 170 to 20 kb), representing a raw sequencing depth of 136X from 303 Gb of data, reduced to 84X and 187 Gb after filtering. The assembled genome of 2.02 Gb was close to the 2.23 Gb estimated by k-mer analysis. Scaffold and contig N50 sizes of 664 and 20 kb, respectively, were comparable to the previously published Gekko japonicus genome. Repetitive elements accounted for 42 % of the genome. Gene annotation yielded 24,755 protein-coding genes, of which 93 % were functionally annotated. CEGMA and BUSCO assessment showed that our assembly captured 91 % (225 of 248) of the core eukaryotic genes, and 76 % of vertebrate universal single-copy orthologs.
Assembly of the leopard gecko genome provides a valuable resource for future comparative genomic studies of geckos and other squamate reptiles.
Genomic DNA was extracted from the tail tissue of a male leopard gecko (Eublepharis macularius: NCBI taxonomy ID 481883; specimen ID TG1477) (Fig. 1). All tissues were collected in accordance with University of Minnesota animal use protocols 0810A50001 and 1108A03545. This animal was captive born from 30+ generations of inbreeding of a strain originating from animals of Indian origin at the Woodland Park Zoo (Seattle) and imports from Pakistan at the National Zoo (Washington, DC) . A total of seven paired-end libraries with a gradient insert size ranging from 170 to 20 kb were constructed and sequenced on an Illumina HiSeq 2000 platform according to the manufacturer’s instructions (Illumina, San Diego, California, USA). For long insert size libraries (2, 5, 10 and 20 kb), the sequenced read length was 49 bp, while for short insert size libraries (170, 500 and 800 bp), the sequenced read lengths were 100 and 150 bp (Table 1). A total of 303 Gb (136X) raw sequences were eventually obtained (Table 1). Before assembly, strict quality control was performed for raw reads using SOAPfilter, a software application in the SOAPdenovo package , which included removing low-quality reads and duplicate reads arising from PCR amplification during library construction. Sequencing errors were corrected using the k-mer frequency method in SOAPec (version 2.02) . After filtering and correction, 187 Gb (84X) high-quality sequences were obtained for genome assembly (Table 1).
We first performed a 17-mer analysis  to estimate the leopard gecko genome size using 54 Gb clean sequences from 170 and 500 bp insert size libraries. Briefly, reads were divided into sliding short sequences of 17 bp, overlapping by 16 bp, with the exception of the first base pair. The count distribution of 17-mers followed a Poisson distribution (Additional file 1). The genome size was estimated as 2.23 Gb for E. macularius by dividing the total number of 17-mers by the peak of distribution (Table 2).
We then assembled a high-quality leopard gecko genome using SOAPdenovo (version 2.0)  in three steps: contig construction, scaffolding, and gap filling. In the contig construction step, SOAPdenovo was used to a de Bruijn graph by dividing high-quality reads from short insert libraries into kmers in which paired-end information was ignored, and kmers were then merged, tips clipped, bubbles merged, and low coverage links removed. Next, contigs displaying unambiguous connections in de Bruijn graphs were collected. A series of kmer lengths were tested and a 33-mer was selected to generate a contig assembly with the longest N50 value. In the scaffolding step, reads from both small and large insert libraries were mapped to contig sequences to construct scaffolds using distance information from read pairs, with the requirement that at least three read pairs were used to form a reliable connection between two contigs. To close intra-scaffold gaps (the gap filling step), overlapping paired-end reads from the 170 bp insert library were first connected using COPE , then Kgf  was employed to close gaps using these connected reads together with reads from other short insert size libraries. An additional local assembly for reads with one end of a read pair uniquely aligned to a contig and the other end located within the gap was performed using GapCloser . The end result was a leopard gecko genome assembly with a total length of 2.0 Gb and scaffold and contig N50s of 664 and 20 kb, respectively, which is comparable to the previously reported Gekko japonicus genome assembly (Table 3) . Comparison of assembly N50s for the leopard gecko genome with eleven previously published reptile genomes (Anolis carolinensis , Python molurus bivittatus , Ophiophagus hannah , Alligator sinensis [8, 9], Alligator mississippiensis, Gavialis gangeticus, Crocodylus porosus , Chelonia mydas, Pelodiscus sinensis , Pogona vitticeps , and Chrysemys picta bellii ) further confirmed that our results were of comparable or better quality (Table 4).
Estimation of genome completeness
We evaluated the completeness of the assembly using CEGMA  and BUSCO , which quantitatively assess genome completeness using evolutionarily informed expectations of gene content. CEGMA assessment showed that our assembly captured 225 (91 %) of the 248 ultra-conserved core eukaryotic genes, of which 210 (85 %) were complete. BUSCO analysis showed that 58 and 18 % of the 3023 expected vertebrata genes were identified as complete and fragmented, respectively, while 24 % were considered missing in the assembly. Both assessment methods showed that our assembly was more complete than the previously reported Gekko japonicus genome assembly (Tables 5 and 6).
We combined a homology-based and de novo method to identify transposable elements (TEs) and other repetitive elements in the leopard gecko genome. Using the homology-based method, we identified known TEs using RepeatMasker  to search against the Repbase TE library (RepBase21.01)  and RepeatProteinMask within the RepeatMasker package to search against the TE protein database. In the de novo method, we first constructed a de novo leopard gecko repeat library using RepeatModeler (http://www.repeatmasker.org/RepeatModeler.html, version 1.0.5) and Piler , and the de novo TE library was subsequently used by RepeatMasker to annotate repeats in the leopard gecko genome. Finally, we used TRF  to predict tandem repeats, with the following parameters: Match = 2, Mismatch = 7, Delta = 7, PM = 80, PI = 10, Minscore = 50. Overall, we identified a total of 851 Mb of non-redundant, repetitive sequences, accounting for 42 % of the leopard gecko genome. The most predominant elements were long interspersed nuclear elements (LINEs), which accounted for 30 % of all TE sequences and 13 % of the genome (Table 7).
We combined homology-based, de novo, and transcriptome-based methods to predict protein-coding genes in the leopard gecko genome.
In the homology-based methods, we downloaded the gene sets of Taeniopygia guttata, Homo sapiens, Anolis carolinensis, Pelodiscus sinensis and Xenopus tropicalis from the Ensembl database (release-73). We first aligned these homologous protein sequences to the leopard gecko genome assembly using TBLASTN with an E-value cutoff of 1e-5, and linked the BLAST hits into candidate gene loci with GenBlastA . We then extracted genomic sequences of candidate loci, together with 3 kb flanking sequences, using GeneWise  to determine gene models. Finally, we filtered pseudogenes that had only one exon with frame errors, as these loci were probably derived from retrotransposition.
In the de novo method, we randomly selected 1000 leopard gecko genes with intact open reading frames (ORFs) and the highest GeneWise score from the homology-based gene set to train the Augustus  gene prediction tool with default parameters. Augustus was then used to perform a de novo gene prediction on repeat-masked genome sequences. Gene models with incomplete ORFs and small genes with a protein-coding length <150 bp were filtered out. Finally, a BLASTP search of predicted genes was performed against the SwissProt database . Genes with matches to SwissProt proteins containing any one of the following keywords were filtered: transpose, transposon, retro-transposon, retrovirus, retrotransposon, reverse transcriptase, transposase, and retroviral.
Transcriptome-based gene prediction was then performed using leopard gecko RNA-seq data from liver, salivary gland, scent gland, and skin tissues obtained from the NCBI database (accession number SRR629643, ERR216315, ERR216316, ERR216322, ERR216325, ERR216304 and ERR216306) . Tophat (v1.3.3) was used to align the RNA-seq reads against the leopard gecko genome assembly to identify splice junctions, and cufflinks (v2.2.1) was used to assemble transcripts using the aligned RNA-seq reads .
Finally, the results of homology-, de novo-, and transcriptome-based analyses were merged to yield a non-redundant reference gene set based on a priority order of transcriptome-based evidence > homology-based evidence > de novo-based evidence. We employed an in-house annotation pipeline to merge the gene data as follows:
A Markov model was estimated with 1000 high-quality genes, which were previously used to train Augustus, using the trainGlimmerHMM tool included in the GlimmerHMM software package . The coding potential of each transcript assembled from the transcriptome data was then identified using the Markov model. Transcripts with complete ORFs were extracted and multiple isoforms from the same locus were collapsed by retaining the longest ORF.
These non-redundant ORFs were then integrated with homology-based gene models to form the core gene set using a custom script. If a gene model with a higher priority overlapped with a model with a lower priority (overlapping length >100 bp), the latter was removed. If two gene models with the same priority overlapped, the one with a longer ORF was preferred.
Homology-based gene models not supported by transcriptome-based evidence but supported by homologous evidence from at least two species were added to the core gene set.
De novo-based gene models not supported by homology-based and transcriptome-based evidence were added to the core gene set where significant hits (BLASTP E-value <1e-5) for non-transposon proteins in the SwissProt database were obtained.
As a result of these steps, a total of 24,755 non-redundant protein-coding genes were annotated in the leopard gecko genome assembly.
Functional annotation of protein-coding genes
We assigned names to 93.59 % of all leopard gecko protein-coding genes by searching against the function databases TrEMBL and SwissProt  using BLASTP (Table 8). We then searched the leopard gecko protein sequences against the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database  using BLASTP to identify molecular pathways that the genes might be involved in. Protein domains and motifs were annotated using InterProScan (version 5.16)  using seven different models (Profilescan, blastprodom, HmmSmart, HmmPanther, HmmPfam, FPrintScan and PatternScan). This revealed that 20,958 of the predicted leopard gecko proteins had conserved functional motifs. We also obtained 1028 Gene Ontology (GO)  terms that were assigned to 15,873 leopard gecko proteins from the corresponding InterPro entry.
Liu B, Yuan J, Yiu S-M, Li Z, Xie Y, Chen Y, Shi Y, Zhang H, Li Y, Lam T-W. COPE: an accurate k-mer-based pair-end reads connection tool to facilitate genome assembly. Bioinformatics. 2012;28(22):2870–4.
Alfoldi J, Di Palma F, Grabherr M, Williams C, Kong L, Mauceli E, Russell P, Lowe CB, Glor RE, Jaffe JD, et al. The genome of the green anole lizard and a comparative analysis with birds and mammals. Nature. 2011;477(7366):587–91.
Castoe TA, de Koning AP, Hall KT, Card DC, Schield DR, Fujita MK, Ruggiero RP, Degner JF, Daza JM, Gu W, et al. The Burmese python genome reveals the molecular basis for extreme adaptation in snakes. Proc Natl Acad Sci U S A. 2013;110(51):20645–50.
Vonk FJ, Casewell NR, Henkel CV, Heimberg AM, Jansen HJ, McCleary RJ, Kerkkamp HM, Vos RA, Guerreiro I, Calvete JJ, et al. The king cobra genome reveals dynamic gene evolution and adaptation in the snake venom system. Proc Natl Acad Sci U S A. 2013;110(51):20651–6.
Wan QH, Pan SK, Hu L, Zhu Y, Xu PW, Xia JQ, Chen H, He GY, He J, Ni XW, et al. Genome analysis and signature discovery for diving and sensory properties of the endangered Chinese alligator. Cell Res. 2013;23(9):1091–105.
Chen H, He G, Hu L, Pan S, Wan Q, Xia J, Xu P, Zhu Y, He J, Ni X et al. Genomic data of the Chinese alligator (Alligator sinensis). Gigascience Database. 2014. http://dx.doi.org/10.5524/100077.
Green RE, Braun EL, Armstrong J, Earl D, Nguyen N, Hickey G, Vandewege MW, St John JA, Capella-Gutierrez S, Castoe TA, et al. Three crocodilian genomes reveal ancestral patterns of evolution among archosaurs. Science. 2014;346(6215):1254449.
Wang Z, Pascual-Anaya J, Zadissa A, Li W, Niimura Y, Huang Z, Li C, White S, Xiong Z, Fang D, et al. The draft genomes of soft-shell turtle and green sea turtle yield insights into the development and evolution of the turtle-specific body plan. Nat Genet. 2013;45(6):701–6.
Georges A, Li Q, Lian J, O’Meally D, Deakin J, Wang Z, Zhang P, Fujita M, Patel HR, Holleley CE, et al. High-coverage sequencing and annotated assembly of the genome of the Australian dragon lizard Pogona vitticeps. Gigascience. 2015;4:45.
Shaffer HB, Minx P, Warren DE, Shedlock AM, Thomson RC, Valenzuela N, Abramyan J, Amemiya CT, Badenhorst D, Biggar KK, et al. The western painted turtle genome, a model for the evolution of extreme physiological adaptations in a slowly evolving lineage. Genome Biol. 2013;14(3):1.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.
Xiong Z, Li F, Li Q, Zhou L, Gamble T, Zheng J, Kui L, Li C, Li S, Yang H et al. Supporting data for “Draft genome of the leopard gecko, Eublepharis macularius”. GigaScience Database 2016. http://dx.doi.org/10.5524/100246.
This work was funded by the China National GeneBank. This research was supported by the Genome 10 k (G10k) project. We thank the faculty and staff in the BGI-Shenzhen, who contributed to the sequencing of the leopard gecko genome, and R. Tremper for providing experimental animals.
Availability of supporting data
Supporting datasets are available at GigaDB . Raw sequencing reads have been deposited in the SRA (Sequence Read Archive) database under SRA ID SRA451060 and Bioproject ID PRJNA339626.
GZ and QL conceived and supervised the project. TG provided the leopard gecko samples. ZX, FL, LZ and JZ performed genome assembly, repeat annotation, gene annotation and gene function annotation. LK, CL and SL provided materials and advice. ZX and FL drafted the manuscript. QL revised the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Authors and Affiliations
College of Forensic Medicine, Xi’an Jiaotong University, Xi’an, Shaanxi, 710061, China
Zijun Xiong & Shengbin Li
China National GeneBank, BGI–Shenzhen, Shenzhen, 518083, China
Zijun Xiong, Fang Li, Qiye Li, Long Zhou, Jiao Zheng, Cai Li & Guojie Zhang
State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences (CAS), Kunming, Yunnan, 650223, China
Zijun Xiong, Qiye Li, Ling Kui & Guojie Zhang
Centre for GeoGenetics, Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, 1350, Copenhagen, Denmark
Department of Biological Sciences, Marquette University, Milwaukee, WI, 53201, USA
Centre for Social Evolution, Department of Biology, University of Copenhagen, Universitetsparken 15, DK-2100, Copenhagen, Denmark
BGI-Shenzhen, Shenzhen, 518083, China
James D. Watson Institute of Genome Sciences, Hangzhou, 310058, China
Frequency distribution of 17-mer analysis. 17-mers are counted from a subset of paired-end reads from 170 bp and 500 bp libraries. The peak depth is 21X. The total number of 17-mers present in this subset is 46,813,180,882. The genome size, estimated by dividing the total number of 17-mer by the peak depth, is 2.229 Gb. (PDF 179 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.