A Comprehensive Description and Evolutionary Analysis of Testudines Mitochondrial Genomes

With a long and successful evolutionary history, turtles are one of the most specialized and ancient groups of reptiles. It appears in the carboniferous, at the same time as the mighty dinosaurs, which disappeared at the end of the cretaceous period. However, tortoises and turtles have survived and multiplied so far. Therefore, tortoises and turtles have played an important role in the systematic evolution of vertebrates (Zhang et al., 2007). There are not many species of turtles and some species have become rare or even endangered due to the changes in the ecological environment, the destruction of human pet market trade, the use of food and medicinal materials and other factors. Therefore, it is of great theoretical and practical significance to study turtle species in terms of the protection of biological diversity and the sustainable utilization of turtle resources. According to statistics, there are about 13 families, 89 genera and more than 270 species of tortoises in the world (Zhou et al., 2004), among which Batagurinae are the largest family. There are 23 genera in Batagurinae, which include Batagur, Callagur, Chinemys, Hieremys, Kachuga, Ocadia, Malayemys, Hardella, Morenia, Geoclemys, Orlitia, Siebenrockiella, Annamemys, Geoemyda, Mauremys, Melanochelys, Notochelys, Sacalia, Rhinoclemmys, Cuora, Pyxidea, Cyclemys and Hoesemys; more than 60 species are included; they are mainly distributed in Europe, North Africa and Asia, with another genus distributed in central and South America.
       
Previous work has been done on morphology (Gauthier et al., 1988; DeBraga and Rieppel, 2010) and chromosomes (Hedges and Poling, 1999) have been studied, but it is difficult to reach a consistent conclusion.
       
The complete mitochondrial genome sequence, critical information (Guangxin, 2015) for the study of genome evolution and species phylogeny, remains unavailable for most of the Geoemydidae. We determined one complete mitochondrial genome of Geoemyda spengleri and examined the relationships among 23 Testudines species based on phylogenetic analyses of mitochondrial genome. We generated complete nucleotide sequence for the mitochondrial genomes and determined the mitochondrial genomic structure. Finally, we conducted phylogenetic analyses based on the mitochondrial sequence data with the main aim of investigating the phylogenetic relationships within Geoemydidae.

The experiment was conducted in the engineering center of Zhejiang Ocean University from 2018-6 to 2019-5.
 
Sample collection and DNA extraction
 
Tissue sample of G. spengleri was obtained from a mountain forest in Guangxi, China (21.7°N,108.0°E). We took the muscle tissue from the severed tail and disinfected the turtle’s wound with alcohol. After the turtle has recovered, we released it into its habitat.The part of the tail muscle was used to determine the complete mitochondrial sequence. Both morphological and PCR-based (COI gene) were used for species identification; the tail muscle was sampled and preserved in 95% ethanol. Total DNA was extracted from the muscle tissue using the proteinase K method (Sambrook and Russell, 2001) and kept at -20°C for PCR amplification. The rest of the tail sample was stored in the Marine Biological Museum of Zhejiang Ocean University, under the registration codes: 1508504.
 
PCR amplifications and sequencing
 
Sixteen pairs of degenerate primers were designed (Supplementary Table 1), based upon the reported complete mitochondrial genome sequences for Mauremys mutica and Sacalia quadriocellata, using the program: 95°C pre-denaturing for 2 min, followed by 35 cycles of 94°C for 40 s, 51°C -58°C for 45 s and 72°C for 1 min and a final extension at 72^oC for 10 min. Each reaction contained 100 ng template DNA, 2.5 μL 10× Buffer (TaKaRa, Dalian, China), 2 μL MgCl 2 (2.5 mol/L), 1.5 μL of each dNTP, 0.25 μL of each primer (25 mol/L) and 1 U Taq DNA polymerase (5 U/μL, TaKaRa). PCR products were electrophoresed on a 1.0% agarose gel and DNA fragments of intended sizes were recovered using a Gel Extract Purification Kit (TaKaRa, Dalian, China). Then the cleaned PCR products were sequenced in both directions with an ABI3730 automated sequencer (Invitrogen Biotechnology Co., Ltd, USA). To ensure the reliability of the data, we amplified and confirmed multiplicity until almost all the sequencing chromatograms displayed single peak and the resultant sequences were as the expected length.
 

 
Sequence analysis
 
The Sequencher TM (GeneCode, AnnArbor, MI, USA) program was used for editing and assembling the contiguous, overlapping sequences. The most DNA sequence data were characterized using BLAST searches at NCBI. Sequence data of G. spengleri was analyzed with EditSeq (DNASTAR). The locations of protein-coding, rRNA and tRNA genes were identified by tRNA Scan-SE1.21 (http://lowelab.ucsc.edu/tRNAScan-SE) and SQUEIN (ver.5.35) program as well as comparisons with the corresponding sequences of other known turtles. Besides, analyses of the control region were performed with DNAsis program. The complete mtDNA sequence of G. spengleri was deposited in GenBank under accession number KU641028.
       
For phylogenetic analyses of major turtle lineages, a total of 23 turtle mtDNA sequences available in the database were chosen (Supplementary Table 2). In addition, Crocodylia and Sauria two clades of diapsid reptile, were used as the outgroup, based on the alliance between turtles, crocodilians and saurians. One typical and uncontroversial crocodilians, Mecistops cataphractus (# MT554033) and one typical and uncontroversial saurians Basiliscus vittatus (#NC 012829), were used here. We aligned the nucleotide data for the 13 protein-coding genes resulting in 11,244 nucleotide sites. Phylogenetic trees were constructed by maximum-likelihood (ML) using MEGA7.0.26.
 

Genome features
 
The mt genome of G. spengleri was 17,448 bp long and conformed to other consensus vertebrate mitochondrial form (Fig 2). It consisted of 13 protein-coding, 2 rRNA, 22 tRNA genes and one control region (D-loop), all of which were similar in length to their counterparts in other turtles (Supplementary Table 1). There were few or small noncoding intergenic spacers, whereas a 26 bp intervening sequence existed between tRNAAsn and tRNACys genes, a 16 bp intervening sequence existed between ND4 and tRNAHis genes,and three partially overlapping, i.e., ATP8-ATP6, ND4L-ND4 and ND5-ND6 sharing 22, 7 and 5 nucleotides, respectively, were found (Table 1).
 

       
Totally, it encodes 3,775 amino acids. The most frequently used amino acids was Leu (16.2%), Thr (9.2%), Ile (8.1%) and Ser (7.5%), while the least common amino acids were Cys (0.9%), Asp (1.8%), Arg (1.8%) and Lys (2.3%) (Fig 1A). Relative synonymous codon usage (RSCU) values for the third positions of the 13 PCGs was shown in Fig 1B. The usage of both two-fold and four-fold degenerate codons was biased toward the use of codons abundant in T or A, in accord with other turtles.
 

       
The nucleotide compositions of the 18 mitochondrial genomes were shown in Table 2. The base composition of the major coding strand of G. spengleri mtDNA was A, 33.7%; G, 13.1%; C, 25.6% and T, 27.6%, showing high average A%+T% and low G% contents as in most other vertebrates. Strand asymmetry of nucleotide composition is usually described by AT and GC skews. The mitochondrial genome AT skews of these chelonian species ranged from 0.098 to 0.136, while the GC skews were all negative, arranging from -0.375 to -0.320. Most of the G. spengleri mt genome protein-coding genes had a start codon of ATG, ND2 start with an ATA initiation codon, COI start with a GTG initiation codon, while ND6 uses CCT as the initiation codon. Protein-coding genes of six terminated with TAA, three ended with TAG, two ended with CAT, one ended with AGG and the remainder of the reading frames had an incomplete termination of a single nucleotide T (Table 1), where the post-transcriptional polyadenylation could produce a TAA stop codon. Besides, we discovered an extra base “G” at a specific position in the ND3 gene of G. spengleri.
 

 

 
Ribosomal and transfer RNA genes  
 
16S rRNA is located between the tRNAIle and tRNAVal genes and 12S rRNA is located between the tRNAVal and tRNATrp  genes. The lengths of 12S rRNA and 16S rRNA genes of G. spengleri are 964bp and 1681bp, respectively. All the mitochondrial genomes of the 18 chelonian species contain 22 tRNAs, 14 of which are encoded by the H-strand and 8 encoded by the L-strand. Among the 22 tRNAs, two forms of tRNALeu (UUR and CUN) and two forms of tRNASer (UUR and CUN) were observed in all chelonian species, ranging in size from 62 to 76 bp and in few of their end sequences overlapping with neighboring tRNAs or rRNAs. All tRNAs except tRNASer (AGY) could be folded into the typical clover-leaf secondary structure.
 
Control region (CR)
 
The major noncoding CR of G. spengleri was 1,937 bp long, flanked by tRNAPro and tRNAPhe genes (Table 1). In this 67.99% A+T rich region, several conserved motifs including three conserved sequence blocks (CSB1, CSB2 and CSB3) were identified. The most striking feature of this control region was the presence of longer tandem repeats of (ATATTATTATATTATTATATATC) n  (n=3) at the 3’ end, downstream to the CSB3. Furthermore, it was followed by a AT-enriched microsatellite sequence flanking the tRNAPhe gene (Supplementary Fig 1).
 

 
Phylogenetic analysis
 
The ML trees are illustrated in Fig 3. From the resultant ML tree, the ingroup species were divided into two major clades: the Pleurodira (Chelodina longicollis, Pelusios castaneus, Podocnemis unifilis) and an assemblage of 20 cryptodires. Within the Cryptodira, 9 geoemydids formed a clade with a bootstrap value of 100%. The phylogenetic tree shows that Geoemydidae was monophyletic in origin and the genera and species in the family were gradually clustered. There were four clades in the family: 1. The C. dentata-S. quadriocellata assemblage and N. platynota formed a clade. 2. C. mouhotii and M. annamensis formed a clade. 3. Heosemys annandalii was a separate clade. 4. G. spengleri and the B. trivittata-P. sylhetensis assemblage forms a clade. In general, the Geoemydidae was the sister group to the monophyly of turtle (M. emys) of Testudinidae and was then clustered with the D. coriacea-N. depressus assemblage fllowed by a clade comprising of P. megacephalum and three representative turtles of Emydidae. The assemblage of the above species was in order clustered with M. temminckii followed by K. leucostomum and P. cantorii. Finally, the Carettochelyidae comprising of C. insculpta was sister group to all remaining Cryptodira used in this study.
 

       
Nine controversial geoemydid genera (nine species) were studied, based on morphological characters, Bramble (Bramble, 1974) hypothesized that Cyclemys, Pyxidea and Cuora formed a closely related assemblage (Cyclemys group) and originated from a Heosemys-like ancestor. In our resultant trees, Cyclemys, Sacalia and Notochelys formed a monophyletic group, Cuora and Mauremys formed another monophyletic group. These two monophyletic groups and their common ancestor formed a new monophyletic group and the intimate relationships between the common ancestor of all the descendants of this new monophyletic group and Heosemys were well supported (100% BP under ML), which seems to be consistent with their morphological characters.
 
As far the relationships of the Geoemydidae, Testudinidae and Emydidae, some authors suggested a closer affinity of the geoemydid turtles with the emydids (Ernst and Barbour, 1989; Iverson, 1992) and the others suggested the Geoemydidae and Testudinidae originated from a common ancestor and sister-relationship of the two families inferred from morphological and molecular data (e.g., Shaffer et al., 1997; Hirayama, 1984; Honda et al., 2002; Lee, 2019).Our ML analyses supported the latter view that the geoemydid turtles have closer affinities to the testudids than to emydids with high statistical supports (Fig 3).
       
Based on morphological and partial mtDNA data, the Chelydridae and Platysternidae were considered as sister taxa (Gaffney and Meylan, 1988; Gaffney et al., 1991; Shaffer et al., 1997). In contrast, most other workers elucidated that Platysternon was more similar to emydids and testudinoids based on morphological evidence (Whetstone,1978), karyotype characteristics (Haiduk and Bickham,1982) and molecular data (Parham et al., 2006). The present ML analysis  implied that Platysternon is a sister taxon to the Emydidae (C. picta, M. terrapin and T. scripta) with a strong support of 98%. Obviously, Platysternon had more intimate relationship with emydidae than Chelydridae according to the phylogram. But some workers agreed with placing Platysternon within Testudinoidea (Wu, 2004). Therefore, further verifications are definitely needed to determine their exact phylogenetic relationships.
       
In this study, the distribution maps of 9 species of the Geoemydidae were drawn to illustrate the species distribution of 9 genera of the Geoemydidae. Looking at the habitat distribution of G. spengleri, the B. trivittata and the Pangshura sylhetensis, we found something interesting (Fig 4). All three species are found in countries or regions with abundant water resources, but there is no overlap in their habitats. In addition, these three species are only distributed in the southern hemisphere and are mostly located in the area between 8 degrees and 28 degrees North latitude.
 

           
Several other species of the family have a wider range of habitats in terms of latitude than these three species. The latitudinal zonality is mainly manifested by the regular change of climate, soil, organisms and their environment from the equator to the poles (Das, 2018). Because the latitude distribution of these three species is more similar than that of the other turtles and the longitude distribution of the habitat of B. trivittata and G. spengleri is more similar than that of P. sylhetensis, we guess that G. spengleri and B. trivittata are more closely related. This is also consistent with the results shown in the phylogenetic tree.

In summary, our mtDNA genome sequence data provided a new stock of useful information to review the relationships among major turtle lineages and verify preliminary hypotheses to turtle evolution further. Obviously, larger number of sample groups and more gene information are needed towards getting a definitive understanding of the grand evolution problem.

We thank many anonymous reviewers for providing helpful suggestions and the technical editors for critically reviewing of the manuscript. This study was jointly supported by the Fundamental Research Funds for Zhejiang Provincial Universities and Research Institutes [No. 2019J00051, No. 2019J00052], the National Natural Science Fund [No. 31700322] and the Natural Science Foundation of Zhejiang Province [No. LQ16C030003, No. LY18C190010], the project of science and technology of Zhoushan [No. 2019C21016].