Abstract

In this study, we report, for the first time, mitochondrial genome of Narcine timlei (Bloch & Schneider, 1801) and its phylogenetic relationships within the order Torpediniformes. Narcine timlei is a medium-sized ray that occurs in nearshore waters of the Indo-Pacific and is classified as “‘vulnerable” category on the IUCN Red List. The mitogenome is assembled from short Illumina reads (150 bp paired end reads). It is 17,964 bp long and includes 13 protein-coding genes (PCGs), 22 tRNA genes, and 2 rRNA genes. The gene order, size, and nucleotide composition are largely consistent with mitogenomic characteristics of previously reported other Narcine spp. The slightly larger mitogenome length of N. timlei than other Narcine spp. may be due to the presence of a putative control region of 1,916 bp with three tandem repeats. Phylogenetic reconstruction using concatenated PCGs (n = 13) of 9 Torpediniformes based on maximum likelihood and Bayesian inference analysis revealed identical topologies. The tree showed two main clades: one clade containing members of the family Narcinidae and the second sister clade consisting of the families Narkidae and Torpedinidae. Our result supports the monophyletic nature of Narcinidae based on mtDNA. The information obtained in this study will contribute to a better understanding of the population genetics, phylogenetic analysis, conservation, and evolutionary biology research of N. timlei.

1. Introduction

With nearly 650 species in 4 orders and 23 families, the superorder Batoidea (rays) forms one of the most speciose groups of the subclass Elasmobranchii [1]. Phylogenetically, they are sister to the superorder Selachimorpha (sharks) [2, 3]. One of the 4 orders of Batoidea is Torpediniformes, which is grouped into 5 families with approximately 68 valid species (Eschmeyer’s Catalog of Fishes assessed on 15th December 2022). They are commonly referred as electric rays due to their ability to generate electrical discharges to stun prey and to defend themselves [4]. Electric rays play an important role in the benthic ecosystem as they are the predators that feed on the diverse invertebrates and small fishes; however, their contribution in benthic dynamic is least known [5, 6].

Previous studies on the phylogenetic positioning of Torpediniformes suggest that they are a sister group to other orders of batoids based on synapomorphies [7], ribosomal genes, and karyological structures [8, 9]. However, molecular phylogenetic studies using different nuclear and mitochondrial genes have shown that the Rajiformes are related to other orders of Batoids and Torpediniformes are related to the order Myliobatiformes [8]. Previous phylogenetic studies based on a few molecular markers, mostly cytochrome oxidase I (COXI) and/or NADH dehydrogenase 2 (ND2), recognized 4 subfamilies within Torpediniformes, but the relationships among them remain confusing [10, 11]. Morphological characters indicate that the families are monophyletic, i.e., separate clades for Platyrhinidae, Narkidae, Narcinidae, Hypnidae, and Torpedinidae [7]. However, molecular phylogenetics using ND2 markers suggests polyphyly of the genus Narcine [2].

In the last decade, mitochondrial DNA (mtDNA) has been widely sequenced to elucidate phylogenetic relationships among taxa because it provides finer taxonomic resolution, especially in cartilaginous fishes [11]. The number of complete mitochondrial genome sequences for elasmobranchs is increasing, mainly due to the reduced cost of sequencing and ease of bioinformatic data analysis, which improves our phylogenetic understanding of fishes [12, 13]. Based on the NCBI database (checked on 20th Jan, 2022), more than 175 complete or partial mtDNA are available for elasmobranchs. However, mtDNA is still scarcely available for fishes of the order Torpediniformes. Complete mtDNA has been reported only for 5 of the 68 valid species: Narcine entemedor [14], N. bancroftii and N. brasiliensis [11], Torpedo marmorata [2], and Narke japonica (GenBank accession: MZ417389.1), and 3 species have partial mtDNA, i.e., N. tasmaniensis, Typhlonarke aysoni, and Tetronarce macneilli [3]. Gaitán-Espitia et al. [11] established the complete phylogeny of Torpediniformes based on 6 mitogenomes, suggesting that the individual orders of Batoidea formed a separate clade, i.e., monophyletic in nature, and that Torpediniformes belong to a group that includes the orders Myliobatiformes, Pristiformes, and Rajiformes. They also showed the genus Narcine is monophyletic, contradicting the earlier report of polyphyly. The recent mtDNA phylogeny encompassing all elasmobranchs suggests that Torpediniformes and Rajiformes form a sister clade, albeit with low support node values [12, 13, 15]. It is worth noting that these recent mtDNA phylogeny studies had one or a few representatives of the order Torpediniformes. Therefore, it is important to generate mtDNA for more species of Torpediniformes to clarify their phylogenetic position.

In the present study, we report, for the first time, the mitochondrial genome sequence of Narcine timlei (Bloch & Schneider, 1801). This species belongs to the family Narcinidae, commonly known as spotted numbfishes. It is a medium-sized ray with large oval/shovel-shaped discs, stout tails, and a naked body (without dermal denticles) [16, 17]. They are known to occur in nearshore waters of the Indo-Pacific ranging from Pakistan to southern China [18]. Their IUCN conservation status was recently changed from data deficient to “vulnerable” [19], yet they are common bycatch batoids in mechanized and artisanal fisheries on the southeast coast of India [20] (authors’ per. obs.). We characterized the mitogenome organization of N. timlei and compared it to other available Torpediniformes to examine the evolutionary relationship within the order.

2. Materials and Methods

2.1. Specimen Collection

The specimen of spotted numbfish N. timlei was collected in November 2021 during our routine survey at the Covelong Fish Landing Center (12°47′31″N; 80°15′04″E) to determine the diversity of catches and bycatch. Covelong fisher folks engage in artisanal fishing, mainly using gillnets, and bottom gillnets at depths of 0–20 m within 5–7 km of shore [21]. Collected specimens were cleaned and photographed in the field before being taken to the laboratory for the detailed study of morphological and meristic characters. Specimens were identified using standard keys and descriptions [16, 17].

2.2. DNA Extraction, Library Preparation, and Sequencing

Total genomic DNA was extracted using Omega Bio-tek E.Z.N.A. Blood and Tissue DNA Kit, as described in [22], and treated with RNase (Promega Corp, USA). The intactness of the DNA was checked by 1% agarose gel electrophoresis. Quantification was performed using the Qubit™ dsDNA BR assay kit (Catalog: Q32853, Thermo Fisher Scientific), and measurements were performed in the Qubit 3.0 Fluorometer (Thermo Fisher Scientific).

After ensuring the quality of genomic DNA, whole genome sequencing libraries were prepared using the NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina (Catalog: E7805S, New England Biolabs). Briefly, 500 ng of DNA was enzymatically fragmented using a fragmentation reagent by targeting 275 bp to 475 bp in size. DNA fragments were subjected to end repair to convert them into blunt ends. The 3′-5′ exonuclease activity of the end repair mixture removes the 3′ overhangs and the polymerase activity fills the 5′ overhangs. The fragments with blunt ends were adenylated by adding a single “A” nucleotide to the 3′ ends. Loop adapters were ligated to the adenylated fragments and cleaved with the uracil-specific excision reagent (USER) enzyme. Size selection was performed according to the manufacturer’s protocol with the addition of AMPure XP beads (catalog: A63881, Beckman Coulter) to achieve a final library size of 400–600 bp. In addition, DNA was amplified by 6 PCR cycles with the addition of NEBNext Ultra II Q5 Master Mix and “NEBNext® Multiplex Oligos for Illumina” to facilitate multiplexing during sequencing. The amplified products were then purified with 0.9X AMPure XP beads (Beckman Coulter), and the final DNA library was eluted in 15 μl of 0.1X TE buffer. Library concentration was determined using Qubit 3 fluorometer, and quality was assessed using the Agilent D1000 Screen Tape System. Paired end sequencing (2 ∗ 150 bp) was performed using Illumina NovoSeq 6000 (Illumina Inc., USA).

2.3. Mitochondrial Genome Assembly and Annotation

A total of 16,107,264 reads were generated, and the quality of the data was checked using FastQC [23] and MultiQC [24]. Low-quality reads (Phred score < 30) and adapter sequences were removed using fastp [25]. After quality filtering, the reads were assembled into contigs using Megahit v.1.1.3 [26] with kmer sizes 21, 49, 77, 105, 133, and141. Contigs of less than 200 bp were removed from the assembly. The final assembled mitogenome of 17,964 bp was obtained and subjected to BLAST homology against the NCBI nucleotide database. In addition, annotations were performed with MitoAnnotator [27] using the genetic code of vertebrate mitochondria. Mitogenome visualization was performed with the CGView server [28] using the composite FASTA sequence and map file from the output of MitoAnnotator. Codon usages and relative synonymous codon usages (RSCU) for each protein-coding gene (PCGs) were predicted in the Codon Usage web server (https://www.bioinformatics.org/sms2/codon_usage.html) and MEGA X [29] using the vertebrate mitochondrial code. tRNA genes were identified using ARWEN software [30] implemented in the MITOS web server [31], and secondary structure was predicted using tRNAscan-SE v.2.0 [32]. The putative control region (POR) was analyzed for the presence of repeats using the Tandem Repeat Finder v.4.09 web server (https://tandem.bu.edu/trf/trf.html).

2.4. Phylogenetic Analysis

The phylogenetic position of N. timlei among other species of Torpediniformes was investigated. The assembled mitogenome of N. timlei, 8 other members of Torpediniformes, and Gymnura poecilura (Table S1) were used for mitophylogenetic analysis, performed using the MitoPhAST pipeline [33]. G. poecilura, which belongs to the order Myliobatiformes, was selected as an outgroup. The MitoPhAST pipeline extracts the nucleotide sequence for 13 PCGs from each of the 10 GenBank files of mitogenomes, aligns each gene with MAFT [34] and TranslatorX [35], trims it with Gblocks [36] to remove ambiguously aligned regions, and concatenates it into supermatrices with FASconCAT-G [37]. The best-fitting substitution models were selected for each partition using ProtTest [38]. The best model for the current dataset was mtMAM + I + G4 for ATP6, ND5, ND3, ND4L, ND4, ATP8, and ND2; mtMAM + I + G4 for COX1, COX2, COX3, ND1, and CYTB; and mtZOA + I for ND6. The rate gamma and rate invariable for ATP6, ND5, ND3, ND4L, ND4, ATP8, and ND2 were 0.823 and 0.246, respectively, and for COX1, COX2, COX3, ND1, and CYTB were 0.640 and 0.412, respectively. The rate of invariance for ND6 was 0.431. Supermatrices along with partition information were used to perform maximum likelihood (ML) phylogenetic analysis by IQ-TREE [39]. The robustness of the ML tree was analyzed by reiterating the observed data with an ultrafast bootstrap approximation for 1000 generations [40]. In addition, gene order information was also obtained for comparative analysis. We also performed phylogenetic analysis using Bayesian inference (BI) in MrBayes [41]. The analysis was performed for 1,00,000 generations (as the standard deviation of split frequencies of <0.005 was achieved), every 100th tree was sampled from the MCMC analysis, and a consensus tree was obtained after discarding the first 25% of the sampled trees. Support for the nodes in the BI tree was obtained by the posterior probability values.

3. Results and Discussion

3.1. Mitogenome Organization

The mitogenome of Narcine timlei was successfully sequenced and assembled, and it was deposited in the NCBI GenBank under the accession number OM404361. The size of the assembled mitogenome was 17,964 bp, which is the expected size range for batoids [13]. However, the size is slightly longer than the previously published mitogenome of other Narcine spp. (Table S2), e.g., 17081 bp in N. entemedor, 16971 bp in N. bancroftii, and 16997 bp in N. brasiliensis [11, 14]. The mitogenome of N. timlei encodes typical mitochondrial DNA genes of metazoans, including 13 protein-coding genes (PCGs) (COX1, COX2, COX3, CYTB, ND1, ND2, ND3, ND4, ND4L, ND5, ND6, ATP6, and ATP8), small and large ribosomal RNAs, and a complete set of 22 tRNAs (Table 1, Figure 1). With the exception of the ND6 gene, all PCGs were transcribed from the heavy strands (H). These PCGs began with the common start codon ATG, with the exception of COX1, which began with the GTG codon. Most PCGs terminated with a complete codon (TAA/TAG/AGA), whereas incomplete termination was observed at ND4 (T). The incomplete termination at T could be extended to TAA through polyadenylation of the 3′ end of the mRNA at the posttranscriptional level, a common phenomenon in the metazoan mitogenome [24].

The base composition of the mtDNA was in the following order: A (36.2%) > T (29.2%), C (22.7%), and G (11.9%), with a tendency towards A + T content. The A + T bias was also observed in all PCGs. The mtDNA showed a significant AC bias (skew_AT = 0.11 and skew_GC = −0.31), indicating a greater abundance of A than T and C than G (Table 2). Similar skewness was also found in the complete genomes of other Narcine spp. (Table S2) indicating a common pattern in this genus.

The A + T bias increases the AT-rich codons in codon usage, which appears to be a common pattern in most vertebrates [42]. The most frequently used codons were ATTIle (5.77%), CTALeu (5.11%), MetATA (4.43%), TTALeu (3.99%), and ThrACA (3.99%), followed by others (Table 3).

The two ribosomal RNAs (large, 16S rRNA, and small, 12S rRNA) were transcribed from the H-strand. 12S rRNA consisted of 944 bp and was located between tRNA-Phe and tRNA-Val. 16S rRNA consisted of 1663 bp and was located between tRNA-Val and tRNA-Leu. Both rRNA genes had a positive AT skew (∼0.20) and a negative GC skew (∼0.1). Of the 22 transfer RNA genes identified, 8 were transcribed from the L-strand and the remaining from the H-strand. Their size ranges from 67 to 75 bp and exhibits a typical cloverleaf secondary structure, except for one tRNA-Ser that contained a simple loop without a D-arm (Figure S1), similar to many metazoan mitogenomes [43].

Gene order, size, and nucleotide composition were consistent with mitogenomic features of previously reported Narcine spp. [11, 14] (Figure S2 and Table S2).

Between the tRNA-Pro and tRNA-Phe genes, we found a putative control region (PCR) of 1916 bp, comparatively longer than in Torpediniformes, which ranges between 1060 and 1328 bp. The difference could be due to the insertion and/or tandem repeats in the control regions [13]. The base composition of the PCR was 31.1% for G, 15.2% for C, 35.2% for T, and 18.4% for T with the negative GC skew (−0.39) and positive AT skew (0.25). We found three repeats in this region, first from 270 to 360 bp with period size 47, second from 1403 to 1427 bp with period size 10, and third from 1714 to 1758 bp with period size 22 (Table S3 and Figure S3). The larger size of the control region might be the reason for the larger mitogenome size of N. timlei compared to other Narcine spp. [13].

3.2. Phylogenetic Reconstruction

MtDNA sequences are considered to have enough phylogenetic information to reveal relationships in fishes because they show small and stable changes over a long period of time and are better than the phylogeny of a single gene or two concatenated genes [43]. We used mtDNA of 9 species representing 3 families of order Torpediniformes, which is by far the most for any mtDNA phylogenetic study on Torpediniformes. Kousteni et al. [13] took 3 species of 2 families and Amaral et al. [12] took 1 species of Torpediniformes in elasmobranch mtDNA phylogeny. The most complete mtDNA phylogeny of Torpediniformes to date was established by Gaitán–Espitia et al. [11], with 6 species from 3 families. In the present study, phylogenetic reconstruction using ML and BI analyses revealed identical topologies with similar branch lengths. We obtained two main clades: one consisted of Narcinidae, while the second consisted of Narkidae and Torpedinidae (Figure 2). Within the family Narcinidae, N. tasmaniensis diverged early from other species in the geological time scale. In addition, N. timlei branched off and formed a separate subclade containing N. enetmedor, N. brasiliensis, and N. bancroftii. The nodes and internodes of the Narcinidae clade were supported by high bootstrap and posterior probability values. Our result supports the monophyletic hypothesis of the family Narcinidae based on the mitogenome [19], in contrast to previous studies that used the ND2 gene phylogeny and suggested polyphyly of the Narcinidae [35]. The earlier studies suggest that the Narcinidae are monophyletic only with the inclusion of the Narkidae [7]. It has also been suggested that some narkids are derived members of the Narcinidae based on comparative anatomy [15] and some genera such as Narcine are sister to Torpedinidae and Hypnidae, while genus Discopyge is sister to Benthobatis and Typhlonarke [45]. The inclusion of Narke japonica in the phylogenetic tree suggests that N. japonica branched early from Typhlonarke ayosni, a sister genus of the Torpedinidae, although nodal support for these branches is lower (<50%). Excluding N. japonica, Torpedo has been reported to split early from other families [10, 11]. Our analysis suggests that tree topologies and interrelationships among the members of the order Torpediniformes have changed with the inclusion of additional species. Therefore, it is necessary to obtain the complete mtDNA of more species to achieve a more accurate phylogenetic resolution within the order.

Data Availability

Mitogenome generated in the present study is submitted to NCBI GenBank under the accession number OM404361.

Disclosure

A preprint has previously been published on Research Square [46], under reference no. 1804785.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors are thankful to the management of Sathyabama Institute of Science and Technology, Chennai, for providing the necessary facilities and funding to carry out the research work.

Supplementary Materials

Figure S1: secondary structure of tRNAs in the mitochondrial genome of N. timlei predicted by tRNAscan-SE v2.0. Figure S2: gene order and patterns in N. timlei and other Torpediniformes. Figure S3: sequence of mitochondrial putative control region highlighting tandem repeats. Table S1: species names and mitochondrial genome sequences used for phylogenetic analysis. Table S2: structural characteristics of mitochondrial genomes of species used in this study. Table S3: description of tandem repeats found in the control region of the assembled mitogenomes. (Supplementary Materials)