Intraspecific variations in Cyt b and D-loop sequences of Testudine species, Lissemys punctata from south Karnataka.

The freshwater Testudine species have gained importance in recent years, as most of their population is threatened due to exploitation for delicacy and pet trade. In this regard, Lissemys punctata, a freshwater terrapin, predominantly distributed in Asian countries has gained its significance for the study. A pilot study report on mitochondrial markers (Cyt b and D-loop) conducted on L. punctata species from southern Karnataka, India was presented in this investigation. A complete region spanning 1.14 kb and ∼1 kb was amplified by HotStart PCR and sequenced by Sanger sequencing. The Cyt b sequence revealed 85 substitution sites, no indels and 17 parsimony informative sites, whereas D-loop showed 189 variable sites, 51 parsimony informative sites with 5' functional domains TAS, CSB-F, CSBs (1, 2, 3) preceding tandem repeat at 3' end. Current data highlights the intraspecific variations in these target regions and variations validated using suitable evolutionary models points out that the overall point mutations observed in the region are transitions leading to no structural and functional alterations. The mitochondrial data generated uncover the genetic diversity within species and conservationist can utilize the data to estimate the effective population size or for forensic identification of animal or its seizures during unlawful trade activities.

Introduction

Lissemys punctata is a cryptodiran omnivorous freshwater turtle that belongs to order Testudines, family Trionychidae. The Reptile Database [1] provides comprehensive information regarding species taxonomic position, distribution, number of species, and their conservation status. Trionychidae turtles (soft shelled turtles) are the preferred group of freshwater turtles for delicacy in Asia, due to its low bone to body ratio [2]. According to Bhupathy et al. [3] L. punctata is being exploited in mass for their meat and eggs in recent years at various parts of India. TRAFFIC, a wildlife trade monitoring network (Wildlife Institute) of India and IUCN (2000) conservation action report reveals that L. punctata is threatened due to high trade activity (CITES Appendix II). However, IUCN (Red list category) has categorized it as Low risk/least concern, which needs up gradation. It is possible only through usage of quick and consistent molecular tool for unbiased identification of the species during illegal trade activity.

DNA based forensic identification is gaining significance in current era, where morphological identification of species is no longer perceptible, such as in animal derived products (carapace powder, fragments of shells, meat etc). In this context, mitochondrial DNA and its loci have made remarkable contributions in turtle studies. Pereira [4] has explained the vertebrate mitogenome organization and the salient features of mtDNA that made it a popular marker in various genetic applications from past two decades. Out of several mtDNA loci studied among vertebrate, the Cyt b gene and regulatory non coding D-loop are the loci that show high resolving power for species identification, due to its high genetic variability [5].

Cyt b is an electron transport chain protein and ROS generator by its function, transmembrane in location show specific sequence variability that can be used to determine relationship within families and genera [6]. Praschag et al. [7] identified three distinct clade of phylogenetic relationship using Cyt b along with other mitochondrial loci in L. punctata subspecies from India and Srilanka. Cyt b sequence is used as valuable locus on mtDNA in forensic sciences for wildlife species identification [8], [9]. D-loop is the predominant regulatory most variable region of mtDNA, used to identify genetically discrete populations, foraging ground, and nesting behavior of turtles, intraspecific variability, and phylogenetic relationship [10], [11], [12], [13]. D-loop mutation has promoted longevity in centenarians [14] and occurrence of tandem repeats (VNTRs) is regarded as valuable molecular marker in turtle species studies [15].

Current mtDNA study in L. punctata was started by specifically extracting Cyt b and D-loop sequence from complete mitogenome sequence of L. punctata available at NCBI with Accession no. NC_012414.1. The mitogenome was composed of 13 protein coding genes, 22 tRNAs, 2 rRNAs with single D-loop/ control region as that of other vertebrates [4]. Partial Cyt b sequences for L. punctata from various states of India (Andra Pradesh, Goa, Gujarat, Karnataka (Mangalore district), Kerala, Maharashtra, Odisha and Tamil Nadu) with respective Accession numbers at NCBI (FR850622, FR850625, FR850626, FR850631, FR850632, FR850635, FR850637, and FR850642) have been submitted by different authors [7]. In this regard, current study is the first to report mitochondrial sequence data of L. punctata from Mysore districts of Karnataka, India. This work focused on intraspecific variations in Cyt b and D-loop region, which would be cumulative to the existing knowledge and growing database of mitogenome. This provides an imperative groundwork for future conservational studies on L. punctata species from south Karnataka, India.

Material and methods

Blood samples from L. punctata (named Lp1 to Lp8) were collected during the field work period from nearby freshwater bodies at Mysore district of Karnataka, India. Permission for collecting blood samples was obtained from Principal Chief Conservator of Forest (Wildlife), Bangalore, Karnataka, India vide letter No. D/WL/CR/149/2010 and PS/WL/CR/21/2013. Without sacrificing or anesthetizing the animal, about 0.3 mL venous blood was drawn from the hind limb femoral vein using lithium/sodium heparin coated BD vacutainers [16]. The animals were left to their habitat after medicating the spot using Betadine solution. Whole blood was immediately stored in 5 mL of Longmire lysis buffer (100 mM TrisHCl pH 8, 100 mM EDTA, 10 mM NaCl and 0.5% SDS) [17] that is recommended for field and room temperature storage of blood samples for future genomic DNA extraction. Genomic DNA extraction followed conventional phenol–chloroform method with proteinase K digestion [18]. Finally, the pellets were resuspended in 1x TE buffer and stored in −20 °C for future genetic analysis.

Set of primers were designed for amplifying Cyt b and D-loop region using specific regions from reference sequence (NC_012414.1) with the help of Primer 3 online tool [19] and presented in Tables 1 and 2. Hotstart PCR was performed at room temperature using JumpStart Taq Ready mix (Sigma-Aldrich Co. LLC, Bangalore, India) with total reaction volume of 50 µL consisting components according to manufacturer’s protocol. Amplifications were performed in Applied BiosystemsVeriti™ 96 well thermal cycler with initial denaturation at 95 °C −1 min, denaturation at 94 °C −30 sec, respective annealing temperatures for Cyt b and D-loop primers (Tables 1 and 2) for 1 min, extension at 72 °C −2 min for 45 cycles and final extensions at 72 °C −10 min and hold at 4 °C. Amplified products were inspected in 2% low melting agarose gel and amplicon size was assessed by running prestained 50 bp DNA ladder (MolBio™- HiMedia, Mysore, India).

Table 1

List of Cyt b gene primers designed using Primer3 online tool.

Oligo (5′-3′)	Length	Tm (°C)	Ta (°C)	Amplicon (bps)
F-GCAACAAATCTACGAAAACATCAC	24	57	55	226
R-CGTATTGTACGTCTCGGGTG	20	58
F-GCCAACGGAGCATCACTATT	20	58.3	55	298
R-AGTGGAAGGTGAAGAATCGGT	21	59
F-CACGAAACTGGATCAAATAACC	22	58.4	55	178
R-TGGCTGGTGAGAAGTTGTCT	20	58.4
F-CCAATAACCCAAACACTATTCTGAT	25	57	55	162
R-AGGCTGGAGAGTGGTATGAG	20	58
F-TAACATTCCGCCCAATAACC	20	59.6	55	180
R-TTTGTTCTCGATTAGGCTGGA	21	59.8
F-TACAATGAATTTGAGGTGGCTTC	23	60	55	200
R-TTTGTACGAGAAGTATGGGTGGA	23	60
F-ACCCAAACATACTTGGAGACC	21	57	55	400
R-TAATGGAGTATTTTGTTCTCGATTAG	26	57
F-CACTACTCACCAAATACTATAACAGCA	27	58	55	155
R-CCGTAGTAAATTCCTCGTCCA	21	59

Table 2

List of D-loop primers designed using Primer3 online tool.

Oligo (5′-3′)	Length	Tm (°C)	Ta (°C)	Amplicon (bps)
F-TCCGCTAGCATATCACCTAT	20	58.5	55	247
R-CCTGAAACTGGTAATGGTGT	20	58.9
F-AGGCCCATTGATAGCTGGAG	20	59	55	201
R-CGGGCCTGAAGACAGAAAGA	20	59
F-CCCATTGATAGCTGGAGGAC	20	59	55	217
R-TCGGCAGACATCAGTTATGC	20	59
F-CATTCGTTCAAGTTGCTTGC	20	59	56	197
R-TTGGGGTTTGACGAGGATTA	20	60
F-CATTCGTTCAAGTTGCTTGC	20	59	57	462
R-GTTGTGATGTCCAAGACATAAAGG	24	59
F-CCCAAAGCCGGAATTTTTA	19	59	55	345
R-AGCTATCAATGGGCCTGAAA	20	59

PCR products were submitted for Sanger Sequencing at Chromous Biotech Pvt. Ltd, Bangalore, India. The products were column purified using Chromous PCR clean-up kit and standard protocol was followed to sequence both strands in Applied Biosystem 3500xl genetic analyzer. Further obtained fasta sequences were aligned using MUSCLE online tool. They were further edited and analyzed using BioEdit v7.2.5 [20] and MEGA v6.06 [21] user friendly software.

Results and discussion

As per Parson et al. [22] the complete Cyt b (1140 bp) and D-loop region (∼1000 bp) were identified solely by comparing with database mitogenome sequence of L. punctata (NC_012414.1). All the sequences subjected to preferential BLAST with L. punctata, scored highest BLAST hits (i.e., E-value less than equal to 0) provides first level confirmation to rule out the chances of 'numts' (nuclear copies of mitochondrial origin) [23]. The common but, not universal phenomenon of transposition of mtDNA fragments into nuclear genome is referred to 'numts' [23]. These are referred as common contaminant encountered, while using genomic DNA from blood sample as a source for mtDNA studies in lower vertebrates [23], [24]. Sorenson and Quinn [23] suggested hints to recognize and avoid 'numts' from actual mtDNA sequence. It includes avoid using universal primers designed for the taxa, instead suggest to use newly designed primers using reference sequence available [23]. Accordingly, all the primers in the current study were newly designed choosing particular region of interest from the reference sequence. Also, the presence of unusual substitutions or stop codons [23] in the sequenced fragments was verified to confirm the absence of 'numts'. According to Spinks and Shaffer [25] the presence of multiple peaks in the sequenced chromatogram indicates 'numts'. The chromatogram obtained in the current study showed clear single peaks and no mixed peaks were observed.

Intraspecific Cyt b sequence variations

Multiple sequence alignment of Cyt b sequence with its reference (Fig. 1) revealed that the open reading frame begins with 'ATG' (Methionine) and terminates with 'TAA' was similar to other reported Trionychidae species in database, except Pelodiscus sinensis from Republic of Korea (AY962573.1) reported 'ATT' start codon for Cyt b [26]. The AT rich (60%) nucleotide composition of Cyt b gene in current study reflects the characteristic vertebrate mitogenome composition [27]. 'TGA' (or UGA) believed to be a universal termination codon for nuclear genes appeared in the current Cyt b sequence stretch but, here it codes for Tryptophan as in all vertebrate mtDNA [28].

Fig. 1

Multiple sequence alignment of L. punctata Cyt b sequences with that of reference (BioEdit v7.2.5). Dots represent consensus nucleotides and Nucleotide (A, T, G, C) denotes variations.

Sequence identity for Cyt b calculated in BioEdit v7.2.5 for pairwise alignment of all L. punctata samples showed 98.6 to 99.9% identity between the individuals. The percentage identity is analogous to Shen et al. [29] report for intraspecific variations that is usually less than 1–2% and do not exceed 5%. The 1140 bp Cyt b sequence analyzed in MEGA v6.06 [21] (Fig. 1) revealed 85 substitution sites, which includes high transition (55) when compared to transversion (28), only 2 sites showed both type of mutations, no indels and 17 parsimony informative sites. With respect to amino acid changes observed after translating the Cyt b nucleotide sequence, there were 45 – 3rd codon transitions and 6 – 1st codon transitions together 51 synonymous substitutions and only 22 non-synonymous substitutions. Avise et al. [30] reported that protein coding genes will have synonymous substitutions higher than non-synonymous substitutions; this appears to be true with our data.

Further, the Cyt b sequences were individually subjected to CDD at NCBI [31] for annotation of protein specific domains. The search parameters included selection of CDD v3.16 database, with no low complexity filter and E-value threshold of 0.01. To increase the consistency of domain search rescue borderline hits and suppress weak overlapping hit parameters were also selected to perform live search at CDD. Out of 500 maximum hits, the best five hits are depicted in Fig. 2, which showed E-value < 0.01. All 8 query sequences showed specific hits with CYTB (MTH00119) a member of cl27766/QcrB protein superfamily representing proteins involved in energy production and conversion of taxon Sauropsida (snakes, beaked reptiles and turtles) (Fig. 2). The conserved functional domains of Cyt b gene like hemebL binding site, hemebH binding site and Qi binding sites spans between nucleotide residues 52nd to 621st in current data (Fig. 2). In spite of nucleotide variations observed in Cyt b sequence the presence of functional conserved domains highlights its role in electron transport system, thereby confirm the protein itself. The Cyt b sequences of L. punctata from the present study can be accessed in future at NCBI using Accession numbers (KY946735, KY946736, KY946737, KY946738, KY946739, KY946740, KY946741 and KY946742) respectively.

Fig. 2

Graphical summary of conserved domains found in Cyt b sequences of current study along with list of domain hits from CDD-NCBI.

Intraspecific D-loop/control region sequence variations

The length of D-loop region flanked by tRNA- Pro and Phen, amplified in the current study range from 999 to 1008 bp and the alignment is as shown in Fig. 3. The overall nucleotide composition consist 33–34% (A), 30–32% (T), 11–12% (G) and 21–22% (C), highlighting AT rich vertebrate mtDNA composition [27]. D-loop sequences (∼1000 bp) analyzed in MEGA v6.06 [21] for intraspecific variations revealed 189 variable sites, 51 parsimony informative sites and only 89 conserved sites. Though D-loop is the most variable, fast evolving part of mitochondrial genome [32], they too possess functional conserved domains preceding VNTR region like TAS, CD (CSB-F) and CSBs (1, 2 and 3) as found in mammals [33], freshwater turtles of order- Geoemydidae [32] and Trionychidae [34]. These functional sites have been identified in the current data and their intraspecific variations are depicted in Table 3.

Fig. 3

Multiple sequence alignment of L. punctata D-loop with that of reference (BioEdit v7.2.5). Dots represent consensus nucleotides, (-) represent gaps and Nucleotides (A, T, G, C) denotes variations.

Table 3

Intraspecies variation of D-loop conserved regions after alignment from BioEdit v7.2.5.

	TAS (36–40)	RC of TAS (26–30)	CSBF (254–266)	CSB1 (595–613)	CSB2 (680–698)	CSB3 (729–747)
NC_012414.1	TACAT	ATGTA	AGAAATAAGCATC	TTAATGCTAGATAGACATA	TTAAACCCCCCTACCCCCC	TCGTCAAACCCCAAAATCC
Lp1	*****	*****	*************	******T***GAG******	C*******T*********T	*******************
Lp2	*****	*****	*************	**********GAG******	C*******T*********T	*******************
Lp3	*****	*****	*************	**********GAG******	C*******T*********T	*******************
Lp4	*****	*****	*************	**********GAG******	C*******T*********T	*******************
Lp5	*****	*****	*************	**********GAG******	C*******T*********T	*******************
Lp6	*****	*****	*************	******T***GAG******	C*******T*********T	*******************
Lp7	*****	*****	*C***********	******T***GAG******	C*******T*********T	*******************
Lp8	*****	*****	*************	******T***GAG******	C*******T*********T	*******************

Xiong et al. [34] reported that TAS domain is characteristic of both pleurodiran and cryptodiran turtles with the sequence 5′-TACAT-3′ and its reverse complementary (RC) sequence 5′-ATGTA-3′ near 5′ region involved in termination of H strand synthesis via stable hairpin loop formation. Among L. punctata samples studied 'TACAT' and 'ATGTA' was observed close to 5′ end of D-loop region (Table 3). According to Xiong et al. [34] the conserved domain (CD) with conserved sequence block (CSB-F) is the only found CD unit among soft shelled turtles. It is characterized by 'AGAAATAAGCATC' sequence. In current data also the region is conserved with the exception of only 1 intraspecific variation - transversion (G to C-2nd residue) (Table 3). CSB-F has been reported in all different classes of vertebrates along with other units (B, C, D, E) signifying its role in H strand replication across the vertebrate lineage [34], [35].

Both Xiong et al. [34] and Zhang et al. [35] observed conserved sequence block CSB1, CSB2, CSB3 in CSB domain in turtles of family Trionychidae and Geoemydidae. These blocks were found along with some regulative regions, H-strand replication origin sites, and transcription promoters for H/L-strand genes [35]. Wang et al. [36] reported that among the three, CSB1 6 bp motif 'GACATA' is conserved from fishes to mammals, birds and also in turtles [34]. CSB1 in soft shelled turtles is identified by 'TTAATGCTAGATAGACATA' sequence [34]. Present data revealed mutations upstream to conserved 6 bp motif i.e., Lp1, Lp6, Lp7 and Lp8 showed transition at position 7 (C to T) and all 8 samples showed 2 transitions (A to G) at position 11, 13, 1 transversion (T to A) at position 12 in CSB1 sequence (Table 3). Xiong et al. [34] reported that a typical CSB2 is characterized by two (C)6 series separated by 'TA', with sequence being 'TT(A)3(C)6TA(C)6′. Our data deviates by 3 transitions in this sequence resulting in 'CT(A)3(C)3T(C)2TA(C)5T' (Table 3). CSB3 is represented by 'TCGTC(A)3(C)4(A)4TCC'sequence showed no intraspecific variations in our data which coincides with Xiong et al. [34] report on interspecies data for the region (Table 3). Wang et al. [36] also reported absence of CSB2 and CSB3 in pleurodiran turtles, birds and some mammals.

Variable number tandem repeats (VNTRs) as the name suggest is the most variable part both in length and type of repeats, a characteristic feature of most vertebrate D-loop region [35]. The distribution of VNTRs is random but, still occurs at certain hypervariable sites or domains at 5′ or 3′ region of D-loop [35], [37]. The dynamic VNTRs decide the length heteroplasmy observed in mitogenome across species lineage, originated by strand slippage or replication mispairing [35]. The D-loop sequences of current study were submitted to Tandem Repeat Finder (TRF) tool [38], revealed only one type of tandem repeat 'AT' with copy number range 23.5 to 24.66 at 3′ end. The similar observation is reported by Xiong et al. [34] for L. punctata species. The D-loop sequences of L. punctata obtained in the current study can be accessed from NCBI database with following accession numbers (KY946743, KY946744, KY946745, KY946746, KY946747, KY946748, KY946749 and KY946750).

Nucleotide substitution analysis for Cyt b and D-loop sequences

The Cyt b and D-loop sequences were subjected to best fit model test in MEGA v6.06 [21]. This resulted in selection of TN93+G [39] (Table 4) and HKY+G+I [40] (Table 5) models respectively, based on lowest AIC corrected criterion values [41]. The nucleotide frequencies used by these models are A = 33.27%, T = 26.91%, C = 28.75%, G = 11.07% for Cyt b and A = 34.31%, T = 31.63%, G = 12.09%, C = 21.97% for D-loop sequence. The maximum Log likelihood values for substitution matrix under these respective model were calculated to be −2042.852 (Cyt b) and −2280.411 (D-loop).

Table 4

Maximum Likelihood estimation of substitution matrix of Cyt b gene of L. punctata based on TN93+G model performed in MEGA v6.06.

	A	T/U	C	G
A	–	4.38	4.68	4.64
T/U	5.41	–	25.25	1.80
C	5.41	23.63	–	1.80
G	13.93	4.38	4.68	–

Tamura-Nei (1993) model discrete Gamma distribution (+G) (5categories; parameter = 0.1571). Lowest AICc (Akaike Information Criterion corrected) value was 4130.998. Gaps and missing data were eliminated. Rates of transitions in bold and transversions in italics.

Table 5

Maximum Likelihood estimation of substitution matrix for D-loop sequence of L. punctata performed in MEGA v6.06.

HKY (G + I) model
	A	T	C	G
A	–	5.77	4.00	7.68
T	6.26	–	13.96	2.20
C	6.26	20.10	–	2.20
G	21.80	5.77	4.00	–

HKY- Hasegawa-Kishino-Yano with [+G] discrete Gamma distribution (5 categories, parameter = 0.6307) and 44.35% [+I] evolutionarily invariable sites were allowed by variation rate model. Lowest AICc (Akaike Information Criterion corrected) value was 4602.931. Rate of transitions in bold and tranvsersions in italics.

The transition/transversion rate ratio observed between Cyt b and D-loop region of L. punctata in the current study are 2.08 and 1.56; this may be due to difference in mutation accumulation, mutation rate and frequency between the regions [42]. The ratios are not comparable with that of other species as the frequency varies between species, classes, order, and genera [43]. Also Keller et al. [44] pointed that transition/transversion bias is not ubiquitous to all vertebrates and invertebrates, rather it is species specific. Cyt b sequence in the study showed transitions and transversions leading to synonymous substitutions of amino acid in the conserved domain and other regions of protein, which in turn impart no change in protein folding, structure and activity. In contrast to this, complete D-loop sequence showed more transversions but, maintain high transitions in their conserved domains so that their regulatory function is not disturbed. Probably, these transversions are not harmful to the organism as the region is noncoding.

Conclusions

The data uncover the intraspecific variations in mitochondrial Cyt b and D-loop regions of the L. punctata. This helps to understand the genetic structure of L. punctata dwelling in the geographical location mentioned in the article. Hence from the analysis, it is clear that mitochondrial protein coding gene Cyt b showed high synonymous substitutions which indicate natural selection favors transitions in protein coding genes so as to maintain its quaternary structure necessary for functional constancy in turtles as well. In case of complete noncoding D-loop region transversions outnumbered transition but, with respect to functional domains transition were more indicating that they are necessary for its regulatory function.

Conflict of interest

The authors declare that there is no conflict of interest.