Complete sequence and characterization of the Mobula tarapacana (Sicklefin Devilray) mitochondrial genome and its phylogenetic implications

2.1 Sample collection and mitochondrial DNA extraction

Mobula tarapacana tissue sample was collected from Nagapattinam (10°45′04′' N and 79° 50′46′' E) fish landing centre, Tamil Nadu, India. The collected tissue sample was preserved in liquid nitrogen and carried to the laboratory. The mtDNA was extracted using QIAamp DNA Micro Kit (Cat.No: 56304, Germany) followed by manufacturer's instructions. The purity of DNA was analysed using NanoDrop 2000 spectrophotometer (Thermo fisher scientific, USA). Also, the quantity of the DNA was examined by Qubit fluorometer (Thermo Fisher Scientific, USA).

2.2 Construction of 2 × 150 NextSeq500 shotgun library

The pure DNA was extracted using TruSeq Nano DNA Library Prep Kit (Illumina, San Diego, California, USA) of the paired-end sequencing library. DNA fragmentation was carried out using Covaris M220 with the mean fragment dissemination of 350 bp. Covaris dissolution results in dsDNA fragments with 3′ or 5′ overhangs. The fragmented overhangs are subjected to blunt ends with the help of an end-repair mix. The 3′ overhangs were removed during the 3′ to 5′ exonuclease activity and also polymerase activity of the 5′ to 3′ overhangs occupies the 5′ overhangs followed by adapter association with the fragments. This scenario provides a low rate of fragment formation. Size selection of the fragments was done using AMPure XP beads. PCR amplification with Index primer was carried out to the size-selected fragments. The ends of the DNA fragments were tied up with the indexing adapters as they were prepared for assimilation onto a flow cell.

2.3 Cluster generation and sequencing

NextSeq500 was loaded with PE Illumina libraries for cluster formation and further sequencing. The Qubit Concentration and the Mean peak value of Agilent TapeStation 4200 were also loaded with other library inputs. Forward and backward direction sequencing was possible through paired-end sequencing on the NextSeq500 platform. In the paired-end flow cell, samples were tied together with the oligonucleotide adapters using kit reagents. These oligos were responsible for selective cleavages on the forward strands and re-synthesis of the complimentary strand throughout the process. The backward sequencing was done in a counter direction and the opposite end of the fragment with the reverse strand of the previous sequence.

2.4 Sequence assembly and annotation

The high-quality pair-end reads acquired from Illumina NextSeq 500 platform were reference-based assembled by SeqMan NGen (DNASTAR Inc., Madison, WI, USA). Sequence alignment was described by NCBI nucleotide database with blastn tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi). PCR amplification was done with the help of standard method (Sun et al., 2016). The characteristic of the assembled mitogenome of M. tarapacana was assessed by MitoAnnotator (http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html) (Iwasaki et al., 2013).

2.5 Mitogenome analysis

The structural modification of the identified transfer RNAs (tRNAs) was predicted using the tRNAmod platform (https://webs.iiitd.edu.in/raghava/trnamod/index.html). OGDRAW 1.2 platform was used for generating the circular mitochondrial genome of M. tarapacana (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) (Lohse et al., 2013). MEGA 6.0 software was used in the analysis of nucleotide composition and codon usage bias (Tamura et al., 2013). The bioinformatic tool DnaSp 5.10.01 were used to calculate the ratio of synonymous and non-synonymous codon usage between the protein coding genes of different species (Rozas, 2009). The following formula AT-skew = (A − T)/ (A + T) and GC-skew = (G − C)/ (G + C) were used for calculating the A + T composition of whole Mitogenome of protein-coding genes (PCGs), tRNA, rRNA and control regions. (Junqueira et al., 2004).

2.6 Phylogenetic analysis

The phylogenetic analysis was performed on the Mobula species were studied using Maximum Likelihood methods, with the dataset of 21 species containing the nucleic acid and concatenated amino acid sequences. Bootstrap values were calculated using 1,000 replicates to assess node support. Evolutionary analyses were done using MEGA X (Kumar et al., 2018).

3.1 Deciphering the mitochondrial genome

An isolated mitochondrial DNA was sequenced with Illumina NextSeq500 platform. Sequencing has resulted in 7,800,506 reads, with a minimum read length of 300 bp. The high-quality data of 2.3 GB was generated and processed to decode the mitochondrial genome. The sequencing analysis has identified the circular mitochondrial genome of Mobula tarapacana with a size of 15,686 bp, and submitted it to GenBank (accession number MH669414.1). An online tool MitoAnnotator were used for annotate the M. tarapacana mitogenome (Iwasaki et al., 2013). In total of 37 genes were annotated which includes 13 protein-coding genes, 22 tRNA genes, and 2 rRNA genes, respectively (Fig. 1). However, the similar mitogenome pattern has been reported in other Mobula species like Mobula thurstoni, Mobula mobular, Mobula japonica. These results are in strongly supported with fish mitogenome (Miya and Nishida, 1999).

Close

Fig. 1

Map of the Mobula tarapacana mitochondrial genome. The innermost circle of the images represents GC% per every 5 bp of the mitogenome; the darker lines are, the higher their GC%.

Export to PPT

The nucleotide composition of M. tarapacana was observed and A + T content was found to be 60.74% of the whole mitogenome. The A + T content of PCGs, tRNA, and rRNA constitutes 60.30%, 62.76%, and 61.33% of the whole mitogenome respectively. The positive AT skewness (0.02) was observed in the whole mitogenome which is lower than other Mobula species representing the presence of higher adenine (As) compared with thymine (Ts). The negative GC skewness (-0.32) represents the presence of more cytosine (Cs) compared with guanine (Gs), which is similar to other Mobula species (Table 1). Previous studies reported that the base composition bias was playing the significant part in the transcription and replication process (Chang and Clayton, 1986).

The bases of 2 gene regions overlapped with their neighbouring genes, from 2 −5 bp in size. The predominant overlapping region observed among ND4L and ND4 with a size of 5 bp. The intergenic spacer nucleotides were observed in 18 places ranging from 1 − 144 bp and the longest spacing sequence of 144 bp was detected between tRNA- Lys and ATP6 (Table 2).

3.2 Protein-coding genes

There are 13 protein-coding genes obtained from the mitogenome of M. tarapacana with the length of 11, 236 bp, which was shorter than other Mobula species except for Mobula japanica, though it covers 60.30% of A + T content. The start codon ATG was used as a canonical putative initiative for all PCGs, except COX1 and ATP8 where they use GTG and CCT as a start codon respectively. These features of initial start and stop codons are generally ascertained in all elasmobranchs (Liu et al., 2013). This was common among all vertebrate mitogenomes (Yue et al., 2016). Among 13 PCGs, 12 of them were determined by heavy strand (ND1, ND2, COXI, COXII, ATP6, ATP8, COXIII, ND3, ND4L, ND4, ND6, Cytb) and ND6 was stated by light strand. The protein-coding genes ND1, COXI, Cytb, COXIII use TAA as termination codon, whereas ND3, ND4L, ND5 use AAT as a termination codon. Likewise, COXII, ATP6, ATP8, ND4, and ND6 use TTT, TCT, CTA, CCT, and TAG respectively. ND2 was encoded by an incomplete termination codon TA (Table 2). The same decoding pattern was observed in the vertebrate mitogenome and extended to a complete TAA termination codon via post-transcription polyadenylation (Pindaro et al., 2016).

The predominant amino acids predicted in the protein-coding genes (PCGs) of M. tarapacana were leucine, proline, and serine in the range of 12.51%, 8.24%, and 10.34% respectively (Table 3). The number of codon usage within the PCGs of M. tarapacana was compared with other Mobula species, and it revealed that leucine, serine, and proline were higher in all species (Figs. 2 & 3).

Close

Fig. 2

Comparison of codon usage within the mitochondrial genome of members of the Mobula species (Mobula tarapacana, Mobula eregoodootenkee, Mobula thurstoni, Mobula mobular, Mobula japonica, Mobula munkiana and Mobula kuhlii).

Export to PPT

Close

Fig. 3

Codon distribution in members of seven species in the Mobula. CDspT = codons per thousand codons.

Export to PPT

3.3 Transfer RNAs (tRNA) and ribosomal RNAs (rRNA)

The tRNA act as a vital role in adaptor molecule during the protein synthesis. The tRNAs of M. tarapacana was 1552 bp in length with 62.76% of A + T content, which includes 32.41% A, 30.35% T, 19.72% G and 17.53% C. AT skewness of tRNA of M. tarapacana was 0.03 whereas the GC skewness was 0.06 which was similar with some Mobula species. The rRNA of the M. tarapacana mitochondrial genome was 2661 bp in length and covers 61.33% of A + T content. Here the As (34.87%) are higher than Ts (26.46%) resulting in a positive AT skewness of 0.14 which was lower than other Mobula species. Similarly, Gs (17.06%) was lower than Cs (21.61%) resulting in a negative GC skewness of 0.12 which is similar to M. thurstoni and M. kuhlii (Table 1). There are 22 tRNA molecules predicted in M. tarapacana mitochondrial genome which is about 70–80 nucleotides long and exhibiting a characteristic cloverleaf structure. The Post-transcriptional modifications carried out in rRNA, mRNA, and tRNA molecules leads to a certain stability, translational efficiency, and fidelity of tRNA. The tRNAmod is a tool that is used to identify the modifications carried out in the tRNA structure. Using this tool on M. tarapacana, it was observed that 20 tRNAs were showing the possibility of modifications in their structure, rather 2 of them were not undergoing any modifications (Fig. 4). Codon-anticodon interactions were highly influenced based on the position of modification, which is usually near the wobble position. This property was well preserved in eukaryotes that resulted in regulation of translational efficiency, frameshifting, and maintenance (Ranjan and Rodnina, 2016; Tuorto and Lyko, 2016). The stability of tRNA depends on the changes occurring within the central tRNA structure, which may result in tRNA degradation and pool discrepancy (Lorenz et al., 2017).

Close

Fig. 4

Schematic representation of tRNA modification in mitogenome of M. tarapacana.

Export to PPT

3.4 Codon usage bias and the control region

The control region of M. tarapacana, M. eregoodootenkee, M. munkiana, and M. kuhlii were not defined, while a few Mobula species have the control region information of their mitochondrial genome (Table 2). The previous finding revealed, the mitochondrial DNA act as a tool for studying various classes of fish including teleost and elasmobranch. They were focused on the control region (CR) as a marker for studying intraspecific variations. The control region exhibits variation in many vertebrates comprising teleosts (Aquadro & Greenberg, 1983), Humans (Lee et al., 1995), and birds (Wenink et al., 1994). Martin et al. (1992) reported the evolutionary rates of elasmobranchs in mitochondrial genes are minimum compared with other vertebrates. In contrast, it is difficult to find out the difference between diverse populations with commonly used CR sequencing approaches. The present study focused on the mitogenome markers except for the control region for the phylogenetic evaluation of M. tarapacana. The Relative synonymous codon usage concerning amino acid utilization of the mitochondrial genome of M. tarapacana and their closely related species were exhibited in Table 3 and Fig. 5. Relative synonymous codon usage (RSCU) analysis of PCGs in M. tarapacana revealed that the amino acid codons encoding Leu, Pro, Ser, Stp, Lys, Thr, Asn, and Ile were the most recurrent ones and those encoding Asp and Lys were rare. The usage of hydrophobic amino acid codons is comparatively higher in the vertebrate mitogenome than general nuclear genome (Satoh et al., 2010). This indicates that the genomic region closer to CR is highly used, hence exhibiting high translational efficiency. This can result in a favorable translation in the vertebrate mitogenome.

Close

Fig. 5

Relative Synonymous Codon Usage (RSCU) of the mitochondrial genome of seven species in the Mobula. Codon families are plotted on the x-axis. Codons indicated above the bar are not present in the mitogenome.

Export to PPT

3.5 The nucleotide substitution of PCGs

Two homologous PCGs of closely related species were compared to find out the evolutionary behavior of similar species. The Ka/Ks ratio was the key to evaluate this evolutionary behavior. The Value of Ka = Ma/Na, the number of nonsynonymous substitutions observed divided by the total number of nonsynonymous changes in the entire gene sequences of two species. Similarly, Ks = Ms/Ns, the number of synonymous substitutions observed divided by the total number of synonymous changes in the entire gene sequences of two species. Hence, the Ka/Ks ratio reflects the rate of evolutionary changes that occurred throughout the sequence and also reflects the selective pressure on the evolution of the organism. Ka/Ks > 1, Ka = Ks and Ka/Ks < 1 are the possible cases and correspond to positive, neutral, and negative selection respectively (Li et al., 2015). The sequence divergence of Mobula mtDNA (Mobula mobular, Mobula tarapacana, and Mobula thurstoni) was determined by Ka and Ks substitution rates. The Ka/Ks value of 13 PCGs is varying from 0.008 (COI) to 1.9 (ATPase8) (Fig. 6).

Close

Fig. 6

Ka/Ks ratios for the 13 mitochondrial protein-coding genes among the reference Mobula tarapacana, Mobula thurstoni and Mobula mobular.

Export to PPT

Here the nonsynonymous substitution rate is faster than the synonymous substitutions taken place. The % ratio of variable sites of MTA/MTH were maximum in ATPase6 and ATPase8 and minimum in the COI gene, signifying that ATPase6 and ATPase8 were broadly distributed and COI was selectively distributed in the mitochondrial proteins. In M. thurstoni and M. mobular, the Ka/Ks ratio was minimum in all 13 protein-coding genes compared to M. tarapacana, indicating that the MMO and MTH had a closer phylogenetic relationship than M. tarapacana.

3.6 Phylogenetic analysis

Phylogenetic investigation of mtDNA is a broadly utilized for outlining the relationship between species (Moritz, 1994) considering its quick pace of succession advancement and fast ancestry arranging comparative with the nuclear genetic material (Avise, 1989). Nevertheless, the data obtained from mtDNA could provide negligible information when small regions are studied. Taxonomy is recognized to be the foundation of the understanding of biodiversity and evolutionary behavior. Phylogenetic analysis on the other hand is used to compare and study the similarities and dissimilarities within a family. Such studies have been challenging due to poor fossil record, misidentification of organism alongside lack of mitochondrial sequence data (Hinojosa-Alvarez et al., 2015). There exist three methods to predict phylogeny, namely maximum parsimony, maximum likelihood, and distance-based phylogenetic methods. Among the three, the maximum likelihood (ML) method yields the closest match and widely accepted among vertebrate genome analysis. The lengthier the branch, causing multiple problems in the Phylogenetic tree due to the nucleotide changes leads to extensive evolutionary patterns (Broughton et al., 2001).

The phylogenetic data from the whole mitochondrial genome data analysis of 21 species and 7 genera in the order: Myliobatiformes. The Maximum likelihood phylogenies of the Mobula species were compatible with the same species. All the clades are in agreement with the bootstrap values and Bayesian posterior probability. There were three clades separated among the Mobula species. The sequence of our strain Mobula tarapacana CK02, MH669414 is closely related to Mobula kuhlii MKU MG10KM364987 and Mobula eregoodootenkee NC 025954. Furthermore, Mobula mobular KT203434 and Mobula japonica JX392983 formed a cluster, and Mobula munkiana MMU MG01KM364990 and Mobula thurstoni NC 037,219 clustered in a separate clade. Given that the sister relationship between Mobula tarapacana CK02MH669414 is closely related to Mobula kuhlii MKU MG10KM364987 based on nucleotide sequence (Fig. 7). Family Mobulidae comprises two genera namely Mobula and Manta. The Mobula genus is holds nine recognized species, whereas the Manta genus has two species. Based on both morphological data and molecular data these species are categorized into 3 clades. The first clade has larger Mobulid species like, M. birostris, M. alfredi, M. tarapacana, M. mobular, and M. japanica. The other two clades include smaller species, three each. M. eregoodootenkee, M. kuhlii, and M. thrustoni belong to the second clade, and M. mukiana, M. hypostoma, and M. rochebrunei belong to the third clade (Muller and Henle, 1841; Aschliman, 2011; Poortvliet and Hoarau, 2013; William et al., 2018).

Close

Fig. 7

Phylogenetic trees of Mobula tarapacana relationships from the complete genome. The concatenated amino acid datasets and evolutionary history was inferred using the maximum likelihood method. The evolutionary distances were computed using the Poisson correction method and are in the units of the number of amino acid substitutions per site. This analysis involved 21 amino acid sequences. Evolutionary analyses were conducted in MEGA X.

Export to PPT

The genetic divergence study reveals genetically indistinguishable pairs in Mobulidae species. Cephalic lode length, tooth morphology, and bronchial filter plate were used as a key to distinguish species of Mobulidae. Morphology along with molecular data analysis shows the indistinguishable nature within a pair of species. M. mobular and M. japanica are morphologically similar except for the difference in tooth morphology, yet further detailed studies revealed the similarities of the species pair at the molecular level (Adnet et al., 2012). M. kuhlii and M. eregoodootenkee share the same evolutionary independent unit based on Phylogenetic analysis (White et al., 2017). Similarly, the molecular data reveals the similarity between M. hypostoma and M. rochebrunei with 100% bootstrap (White et al., 2017). The nuclear exon data reveals a strong link between M. mukiana and M. hypostoma, where the width of the lower tooth band is the only distinguishing morphological features present (William et al., 2018). Morphological features like pectoral fins, tail, and structure beneath and above the first dorsal fin are the key characters to distinguish between interspecies. These morphological studies revealed that these key features remain small in M. tarapacana, compared with M. mobular, M. kuhlii and M. thurstoni (Paig-Tran et al., 2013).