Characterization of five complete Cyrtodactylus mitogenome structures reveals low structural diversity and conservation of repeated sequences in the lineage

Specimen collection and DNA extraction

Individual specimens of each of the five species (C. peguensis, C. tigroides, C. thirakhupti, C. auribalteatus, and C. chanhomeae), taken from Nakhon Ratchasima Zoo during 2012–2014, were stored in 100% ethanol and detailed information is presented in Table 1. Individual species were classified based on their morphology (Boulenger, 1893; Bauer, Sumontha & Pauwels, 2003; Pauwels et al., 2004; Sumontha, Panitvong & Deein, 2010). A piece of muscle clipped from each sample was collected to provide a source of DNA. Whole genomic DNA was extracted following the standard salting-out protocol as described previously (Supikamolseni et al., 2015). Animal care and all experimental procedures were approved by the Animal Experiment Committee, Kasetsart University, Thailand (approval no. ACKU-SCI-022), and conducted in accordance with the Regulations on Animal Experiments at Kasetsart University.

Complete mitogenome sequencing

Complete mitogenome sequences were obtained using a PCR (polymerase chain reaction)-based strategy to amplify overlapping mitochondrial fragments. Forty-six PCR primer pairs were designed based on five squamate reptile mitogenome sequences: Gekko chinensis (NC_027191.1), G. gecko (NC_007627.1), G. japonicas (NC_028035.1), G. swinhonis (NC_018050.1), and Hemidactylus bowringii (NC_025938.1). Nineteen primer pairs were also taken from Kumazawa & Endo (2004) (Table S1). PCR amplification was performed using 20 µl of 1×ThermoPol buffer containing 1.5 mM MgCl₂, 0.2 mM dNTPs, 5.0 µM of primers, 0.5 U of Taq polymerase (Apsalagen Co. Ltd., Bangkok, Thailand), and 25 ng of genomic DNA. PCR conditions were as follows: initial denaturation at 94 °C for 3 min, followed by 35 cycles of 94 °C for 30 s, 45–60 °C for 30 s, 72 °C for 1 min 30 s, and then a final extension at 72 °C for 5 min. PCR products were detected by electrophoresis on 1% agarose gels. No multiple bands were found in any of the PCR products. PCR products were purified using FavorPrep GEL/PCR Purification Mini Kit (Favorgen Biotech Corp., Ping-Tung, Taiwan), and nucleotide sequences of the DNA fragments were determined using a 3500 Genetic Analyzer (Life Technology, California, USA). BLASTn and BLASTx programs (http://blast.ncbi.nlm.nih.gov/Blast.cgi) were used to search for nucleotide sequences in the National Center for Biotechnology Information database to confirm the identity of the DNA fragments amplified in the present study.

Nucleotide sequence annotation and analysis

Sequence assembly was performed to combine all overlapping PCR fragments into one contig strand using SeqScape software v 2.5 (Life Technology, California, USA) and re-checked carefully by visual inspection to ensure accuracy of the variable sites and avoid the problems of repeated sequences identified by the program in the CR. Gene structures of all mitogenomes were annotated using web-based MITOS (Bernt et al., 2013) with manual inspections to compare DNA or amino acid sequences with known sequences from several gecko lizards. ExPASy-translate tool (Gasteiger et al., 2003) was used to characterize the nucleotide sequences of encoded genes. To identify mitochondrial tRNA genes, the nucleotide sequences were searched for regions which could form characteristic secondary structures using tRNA Scan-SE 2.0 by parameter sequence source: vertebrate mitochondrial genome search mode (http://lowelab.ucsc.edu/tRNAscan-SE), the RNAfold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) and UNFOLD (http://unafold.rna.albany.edu/) under fold algorithms as basic options is minimum free energy (MFE) and partition function (Pereira et al., 2008). Overlapping regions and intergenic spacers between genes were counted manually. The CR domains comprising the extended termination associated sequence (ETAS) domain, central conserved region (CCR) domain, and conserved sequence block (CSB) domain were characterized following Douzery & Randi (1997). Tandem repeated sequences in CR were identified using the program ‘Tandem Repeats Finder’ with default parameters (Benson, 1999). The most thermodynamically stable putative secondary structures of VNTR were determined using the UNFOLD and RNAfold web servers. All nucleotide sequences were then deposited in the DNA Data Bank of Japan (DDBJ) (Table 1). Base composition and codon usage were analyzed using the default parameters of Molecular Evolutionary Genetics Analysis 7 (MEGA7) software (Center for Evolutionary Functional Genomics, The Biodesign Institute, Tempe, AZ, USA; Kumar, Stecher & Tamura, 2016). The A +T% and G +C% values, and the AT-skew and GC-skew (Perna & Kocher, 1995) for the H-strand were also calculated for five Cyrtodactylus mitogenomes, another 25 gecko lizard mitogenomes, and Iguana iguana (NC_002793.1), Plestiodon egregious (NC_000888.1), and Varanus salvator (NC_010974.1) mitogenomes as the outgroup. Values obtained were subsequently mapped as a scatter plot.

Phylogenetic placement

Multiple sequence alignment of concatenated heavy-strand encoded protein-coding genes was performed on 33 sequences for 30 gecko lizards with the outgroup (I. iguana, P.  egregious, and V. salvator) using the default parameters of MEGA 7 software. Additionally, the ND6 gene was used separately to perform alignment from phylogenetic analyses of the datasets as it was encoded by the L-strand and possessed base composition bias (Dong & Kumazawa, 2005; Wang et al., 2009). The dataset of sequence lengths of concatenated heavy-strand encoded protein-coding genes ranged from 10,766 to 10,899 bp, while ND6 sequence lengths ranged from 510 to 546 bp. All unalignable and gap-containing sites were carefully removed and trimmed from the datasets. Datasets of aligned concatenated heavy-strand encoded protein-coding genes and ND6 fragments showed 10,813 and 548 nucleotides, respectively. Level of sequence divergence between species was estimated using uncorrected pairwise distances (p-distances) as implemented in MEGA7. We also performed multiple sequence alignment comprising gap-containing sites (as insertion and deletion) using the GUIDANCE2 Server (Sela et al., 2015). Datasets of aligned concatenated heavy-strand encoded protein-coding genes and ND6 fragments showed 11,322 and 607 nucleotides, respectively. Phylogenetic analysis was performed using Bayesian inference (BI) and maximum likelihood (ML). The best-fit model of DNA substitution was determined for each gene using Kakusan4 (Tanabe, 2011). The GTR+I+G model gave the best-fit of each locus for BI and ML analyses. BI analysis was performed with MrBayes v 3.2.6 (Huelsenbeck & Ronquist, 2001). The Markov chain Monte Carlo process was used to run four chains simultaneously for one million generations. Log likelihood and parameter values were accessed through Tracer version 1.6 (Rambaut et al., 2015). After the log-likelihood value stabilized, a sampling procedure was performed every 100 generations to obtain 10,000 trees and a majority-rule consensus tree with average branch lengths was generated. All sample points prior to attaining convergence were discarded as burn-in, and the Bayesian posterior probability in the sampled tree population was calculated as a percentage. ML analyses for all datasets were performed using RAxML (Stamatakis, 2014), followed by executing 10 runs of random additional sequences generating 1,000 rapid bootstrap replicates. By combining the definite phylogenetic relationships of different species with the distribution pattern of mitogenome types in gecko lizards, evolutionary processes of the mitogenome structure were determined. Mitochondrial gene rearrangements with similar components and order were classified as identical. All mitogenome types were visualized by linearized organization and drawn by OrganellarGenomeDRAW v 1.2 (OGDraw) (Lohse et al., 2013).

Mitochondrial gene rearrangement of Cyrtodactylus species in gecko lizards

Generally, organization of the 37 genes and the CR of most gecko lizards tend to be conserved in squamate reptiles (Zhou et al., 2006; Kumazawa et al., 2014; Kim et al., 2015; Hao, Ping & Zhang, 2016). However, various rearrangement patterns of mitogenomes have been found in snakes (Yan, Li & Zhou, 2008), varanid lizards (Amer & Kumazawa, 2008), and acrodont lizards (Okajima & Kumazawa, 2010) as a consequence of the shuffling of tRNA gene clusters, translocations and/or duplications of genes, and gene loss (Anderson et al., 1981; Fujita, Boore & Moritz, 2007; Kumazawa et al., 2014; Li et al., 2015). Comparison of gecko lizard mitogenomes revealed six types of mitochondrial gene rearrangements as I, II, III, IV, V, and VI based on available gecko lizard mitogenome sequences. Type I is distributed in most gecko lizards (80%) including Cyrtodactylus species. Types II and III are found in Tropiocolotes tripolitanus and Stenodactylus petrii, respectively. Type II was derived from the Type I by duplication of tRNA-Gln and ND4, whereas Type III showed multiple copies of O_L within the WANCY cluster. Type IV was observed in Uroplatus fimbriatus, resulting from the translocation of tRNA-Glu from 5′ proximal ND 6 to 5′ proximal tRNA-Pro (Fig. 2). Deletion of tRNA-Glu was found in Uroplatus ebenaui as Type V, while duplication of 5′ proximal ND4- to 5′ proximal tRNA-Ile was observed in Heteronotia binoei as Type VI. Our phylogenetic placement indicated that gecko lizard mitogenome Types II–VI were found at the terminus of the tree in Gekkonidae, suggesting that structural variation of mitogenomes has occurred independently in gekkonids. Rare changes of mitochondrial gene rearrangement in gecko lizard lineages may indicate the possibility of homoplasy-free datasets. However, many family-level taxa are represented by a few species; thus, more gecko lizard mitogenome sequences are required to better understand the evolutionary history in this lineage.

Dynamics of CR in Cyrtodactylus

Although the CRs of the five Cyrtodactylus species were mainly composed of high A+T%, distribution of AT was not homogenous among the three domains. The AT content of the CSB domain was lower than the ETAS and CCR domains (Table 4). The size of CRs varied from 1,445 to 1,728 bp among Cyrtodactylus species and tallied with those of gecko lizards. This feature was mainly caused by the duplication of tandem repeat at VNTR and frequently occurs in ETAS and/or CSB domains (Larizza et al., 2002). In the ETAS domain, different copy numbers of the 75-bp box tandem repeat were found in most Cyrtodactylus with at least 76% similarities, except for C. auribalteatus containing 225-bp box. These results collectively suggest that the 75-bp box tandem repeat was present in ancestral mitogenomes before divergence of Cyrtodactylus species. Here, the 75-bp box is designated as 75-bp box family. Variable copy numbers have resulted from independent duplication in each species, a factor commonly observed in other vertebrate groups (Zhou et al., 2006; Kumazawa et al., 2014; Starostová & Musilová, 2016).

Comparison of 75-bp and 225-bp boxes among Cyrtodactylus species revealed 69% similarities. The 225-bp box is likely to contain two 75-bp boxes tandemly and one 75-bp box-like (Fig. 3). This result suggests that the presence of 225-bp box resulted from duplication of 75-bp boxes with nucleotide substitution, followed by homogenization and fixation in C. auribalteatus after it diverged from Cyrtodactylus. Interestingly, both the sequence of 75-bp and 225-bp boxes could form stable secondary structures that might be responsible for increasing possible misalignment (H-strand and L-strand), involving rearrangement in CR during DNA replication of mitochondrial DNA (Broughton & Dowling, 1997; Pereira et al., 2008). Secondary structures of VNTR are believed to be termination sequences in the replication process in fish (Terencio et al., 2013). Although the function of these secondary structures remains completely unknown, they were probably retained in Cyrtodactylus under selective pressure. Structural and functional studies are required to explain this molecular mechanism. To investigate the distribution of the 75-bp box family in gecko lizards, a comparative search (BLAST) was conducted for similar sequences housed in the GenBank. A tandem repeat with 55–75% similarities was found in most gecko lizards, even for Eublepharidae and Pygopodidae as the basal group. This result suggests that the tandem repeat of 75-bp box family was acquired in the genome of the common ancestor of gecko lizards and evolved gradually through independent nucleotide mutation in each lineage.

The CCR domain is highly conserved among gecko lizards compared to ETAS and CSB domains. The CSB-F block characterized the CCR domain in Cyrtodactylus and different levels of similarity (59.26–96.15%) were observed when compared with other gecko lizards. However, two blocks (CSB-E and/or CSB-D) observed in other vertebrates were not identified in Cyrtodactylus and other gecko lizards (Anderson et al., 1981; Broughton & Dowling, 1994; Lee et al., 1995). This result suggests plausible plasticity across vertebrates concerning the CSB-F, CSB-E, and CSB-D blocks. The CSB domain in Cyrtodactylus and gecko lizards is characterized by the three conserved blocks CSB-1, CSB-2, and CSB-3 which involve regulation of replication and transcription (Li et al., 2013; Kumazawa et al., 2014). These three blocks have also been identified in vertebrate mitogenomes (Liu, 2002; Brehm et al., 2003; Zhang et al., 2009; Xiong et al., 2010). The sequence of CSB-3 blocks is the most highly conserved in Cyrtodactylus and other gecko lizards, whereas the CSB-1 block showed the most variability. Repeated sequences as minisatellites and/or microsatellites (TA)_n or (CA)_n were found between CSB-1 and CSB-2 or between CSB-3 and tRNA-Phe in Cyrtodactylus and other gecko lizards (Table S8). Minisatellites were found in four gecko lizards (Hemitheconyx caudicinctus, Gekko japonicas, U. fimbriatus, and U. ebenaui) which did not contain tandem repeats in the ETAS domain. This result suggests that repeated sequences are detected in at least ETAS or CSB domains of CR. Thus, VNTRs can be used as molecular markers to provide phylogenetic information in population genetics, species identification, genetic diversity, and conservation (Zardoya & Meyer, 1998; Zhang et al., 2009). However, VNTRs differ in length, copy number, and base composition even in very closely related species, and may be shared by distantly related species but not by close ones, such as in C. auribalteatus. This result suggests that the presence/absence of repeated sequences might not be a reliable phylogenetic marker.