The first complete mitogenome of the South China deep‐sea giant isopod Bathynomus sp. (Crustacea: Isopoda: Cirolanidae) allows insights into the early mitogenomic evolution of isopods

Abstract In this study, the complete mitochondrial (mt) genome sequence of the South China deep‐sea giant isopod Bathynomus sp. was determined, and this study is the first to explore in detail the mt genome of a deep‐sea member of the order Isopoda. This species belongs to the genus Bathynomus, the members of which are saprophagous residents of the deep‐sea benthic environment; based on their large size, Bathynomus is included in the “supergiant group” of isopods. The mt genome of Bathynomus sp. is 14,965 bp in length and consists of 13 protein‐coding genes, two ribosomal RNA genes, only 18 transfer RNA genes, and a noncoding control region 362 bp in length, which is the smallest control region discovered in Isopoda to date. Although the overall genome organization is typical for metazoans, the mt genome of Bathynomus sp. shows a number of derived characters, such as an inversion of 10 genes when compared to the pancrustacean ground pattern. Rearrangements in some genes (e.g., cob, trnT, nad5, and trnF) are shared by nearly all isopod mt genomes analyzed thus far, and when compared to the putative isopod ground pattern, five rearrangements were found in Bathynomus sp. Two tRNAs exhibit modified secondary structures: The TΨC arm is absent from trnQ, and trnC lacks the DHU. Within the class Malacostraca, trnC arm loss is only found in other isopods. Phylogenetic analysis revealed that Bathynomus sp. (Cymothoida) and Sphaeroma serratum (Sphaeromatidea) form a single clade, although it is unclear whether Cymothoida is monophyletic or paraphyletic. Moreover, the evolutionary rate of Bathynomus sp. (dN/dS [nonsynonymous mutational rate/synonymous mutational rate] = 0.0705) is the slowest measured to date among Cymothoida, which may be associated with its relatively constant deep‐sea environment. Overall, our results may provide useful information for understanding the evolution of deep‐sea Isopoda species.


| INTRODUCTION
Mitochondria, which are regarded as relicts of bacterial or alphabacterial endosymbionts incorporated into an early eukaryotic cell, have their own genetic material and genetic system. Since the discovery of intraorganellar DNA in chick embryo mitochondria (Nass & Nass, 1963), determination of the structure of mitochondrial genomes (mitogenomes) has become an important aspect of genome research because of the important data it can provide for phylogenetic analyses.
In general, the metazoan mitogenome is a circular, double-stranded DNA molecule approximately 12-20 kb in length that consists of the same 37 genes: 13 protein-coding subunits, two ribosomal RNAs, and 22 transfer RNAs (Wolstenholme, 1992). Due to the different rates of evolutionary change among different segments of a given mitogenome (Helfenbein, Fourcade, Vanjani, & Boore, 2004;Liebers, de Knijff, & Helbig, 2004), the gene order (Roehrdanz, Degrugillier, & Black, 2002), and the RNA secondary structure (Macey, Schulte, & Larson, 2000), these sequences can provide large datasets for phylogenetic analyses at different levels and can also serve as ideal models of gene rearrangement and genome evolution.
Since Clary and Wolstenholme (1985) sequenced the mitogenome of Drosophila yakuba in 1985, 1,160 Arthropoda mitogenomes have been determined. However, only 203 mitogenomes within the class Crustacea have been determined, and the number of isopod mitogenomes currently available is even lower, with only three complete mitogenomes (one each for Limnoriidae, Ligiidae, and Phreatoicidae) and another eight sequences that are almost complete (one each for Sphaeromatidae, Chaetiliidae, Cirolanidae, Asellidae, Armadillidiidae, Idoteidae, Cylisticidaeand, and Trachelipodidae; see Table 1). Among crustacean mitogenomes, most share the ancestral pancrustacean (crustacean + hexapod) gene order that shows only a trnL-UUR translocation relative to the ancestral arthropod arrangement found in the horseshoe crab Limulus polyphemus (Lavrov, Boore, & Brown, 2000) or present only minor tRNA translocations (Yang & Yang, 2008).
Approximately 77% of the ocean floor and 60% of our planet's surface are covered with deep-sea habitats that remain largely unexplored. In these deep-sea environments, a lack of sunlight, extremely high pressures, and low oxygen levels prevent the formation of typical biological assemblages due to the lack of photosynthesis, which is responsible for biosphere primary production. Nonetheless, the deep sea is a species-rich habitat, as has been well documented since marine biologists began to extensively sample the bathyal and abyssal depths (Grassle, 1989;Hessler & Sanders, 1967;Wolff, 1977). For example, Isopoda comprises a highly diverse and species-rich group of crustaceans, with many members living in the abyssal benthos in all oceans (Hessler, Wilson, & Thistle, 1979;Wolff, 1962), and the existence of such biological communities has produced a profound change in our perception of deep-sea life (Van Dover, German, Speer, Parson, & Vrijenhoek, 2002). In addition, mitochondria are the energy metabolism centers of the cell because more than 95% of cellular energy is generated by mitochondria through oxidative phosphorylation (OXPHOS). Mitochondrial-encoded OXPHOS genes may therefore evolve under selection due to metabolic requirements and display evidence of adaptive evolution in mammals, birds, and fishes (Shen, Shi, Sun, & Zhang, 2009;Sun, Shen, Irwin, & Zhang, 2011).
Because of their large size, Bathynomus spp. (Crustacea: Isopoda: Cirolanidae) are classified into the "supergiant group" of isopods (Lowry & Dempsey, 2006). These animals are important scavengers in the deep-sea benthic environment, from the gloomy sublittoral zone, at a depth of 170 m, to the dark of the bathypelagic zone at 2,140 m, and they are often found at depths between 365 and 730 m (Holthuis & Mikulka, 1972). Sankar et al. (2011) also mentioned these species as the first recorded deep-sea isopods in the waters off India. However, no complete mitogenome data are available to date for any deep-sea isopod, even though studying the early mitogenomic evolution of deep-sea isopods using mt DNA fragments is facilitated by effective Crustacea-specific versatile primers (Crandall & Fitzpatrick, 1996;Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994;Merritt et al., 1998).
In this study, we report the first complete mitogenome sequence of the Bathynomus sp. mitogenome, consisting of the same 13 protein-coding genes and two rRNAs found in most metazoans but only containing 18 tRNAs, with tRNA-H, tRNA-I, tRNA-L1, and tRNA-S1 missing from the typical structure. Moreover, the genes show surprising rearrangements compared to the ancestral pancrustacean order and the isopod ground pattern. Phylogenetic analyses indicate that Bathynomus sp. (Cymothoida) and Sphaeroma serratum (Sphaeromatidea) form one clade and that the evolutionary rate of Bathynomus sp. (dN/dS = 0.0705) is the slowest measured to date within Cymothoida. Our results may provide useful information for understanding the evolution of deep-sea isopod species.

| MATERIALS AND METHODS
Experiments were performed in accordance with the recommendations of the Ethics Committee of the Institute of Hydrobiology, Chinese Academy of Sciences. These policies were enacted according to the Chinese Association for Laboratory Animal Sciences and the Institutional Animal Care and Use Committee protocols.  S1). All specimens were preserved in 99% ethanol until DNA extraction. Species-level morphological identification was performed based on original morphological descriptions, locality data, and additional information about the giant deep-sea scavenger genus Bathynomus (Lowry & Dempsey, 2006). Total genomic DNA was extracted from one or several pleopods using the EZNA ® Tissue DNA Kit (OMEGA, Wuhan, China) following the manufacturer's instructions.
To facilitate the subsequent long PCR, we chose candidates with high melting temperatures (Table S1). Long PCRs were performed in a volume of 25 μl containing 12.75 μl of sterilized ultrapure water, 2.5 μl of 10× Takara LA PCR buffer (including MgCl 2 ), 4 μl of Takara dNTP mixture (2.5 mmol/L each), 0.25 μl of Takara LA Taq polymerase (5 units/ μl), 1 μl of primer mixture (10 mmol/L each), and 4 μl of DNA template (100 ng/μl). The corresponding thermal cycler protocol consisted of an initial denaturation step (94°C, 1 min), followed by 30 cycles of denaturation (98°C, 10 s), annealing and extension (64-68°C, 12 min), and ending with another extension (72°C, 10 min). The long PCR products were sequenced using a primer-walking strategy. Considering the lack of data regarding genome arrangement for this species and to obtain precise genomic sequences, five gene-specific primer pairs were designed to cover any gaps and inaccurate sequencing.
All of the PCR products were visualized on 1.2% low-melting agarose gels stained with ethidium bromide. The products were then purified and sequenced using an ABI3730XL sequencing system, and the primers are described in Table S1.

| Gene annotation and sequence analysis
Overlapping fragments obtained by sequencing were edited and aligned using BioEdit v7.0.9.0 (Hall, 1999). The MITOS webserver (Bernt et al., 2013) was used to annotate the Bathynomus sp. mitogenome, and protein-coding and ribosomal RNA genes were rechecked by aligning them with publicly available mitogenomes of 11 isopod species (Table 1). The reliable identification of tRNA genes is generally not trivial, and the most common approach is to search for base pairings that conform to a typical tRNA clover-leaf structure. Therefore, tRNAs were initially re-detected using two computer programs: tRNAscan-SE version 1.21 (Lowe & Eddy, 1997) and ARWEN 1.2.3.c (Laslett & Canback, 2008). Confirmation was performed by blast searches of all tRNAs in MitoZoa 2.0 (de Meo et al., 2012).
CGView (Stothard & Wishart, 2005) was used for a circular display of the Bathynomus sp. mitogenome, which was modified manually. The complete genome sequence has been submitted to NCBI (GenBank: KU057374). The nucleotide composition was analyzed with MEGA 5 (Tamura et al., 2011). Strand skew values were calculated according to the formulae given by Perna and Kocher (1995): AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C), where A, T, C, G are the four bases.

| Phylogenetic analysis and evolutionary rate estimation
We utilized the 11 other determined isopod mitogenomes for phylogenetic analyses (Table 1) Maximum likelihood analyses of the nucleotide and amino acid alignments were assembled in PhyML 3.0 (Guindon & Gascuel, 2003), with 1,000 replicates performed under the models GTR + G + I and MtArt + G + F, respectively. Bayesian analyses of the nucleotide and amino acid alignments were carried out using MrBayes 3.1.2 (Ronquist & Huelsenbeck, 2003) under the models GTR + G + I and MtRev + I + G + F (see above), respectively, with 2,000,000 generations in two runs of eight chains each.
To estimate the evolutionary rate, standard branch models calculated with the CODEML program of PAML 4.6 (Yang, 2007) were used, and an adjusted Chi-square test (Storey & Tibshirani, 2003) was applied for testing p values.

| Mitogenome organization
The mitogenome of the South China deep-sea giant isopod Bathynomus sp.
is a 14,965-bp circular molecule ( Figure 1)  Incomplete determination with a partial cob sequence, and a few tRNAs and the control region were not sequenced.
c Incomplete determination, with a few tRNAs and the control region not sequenced.
of A% = 26.6, G% = 24.3, T% = 32.1, and C% = 17.0. The genome has an overall AT content of 58.7%, which appears to be low for isopods (typical range: 54.4%-71.2%). The composition is skewed away from G in favor of C (the GC skew is +0.177) but is almost balanced for A and T (the AT skew is −0.094). This feature is well conserved among isopods (Table 1).
This species has the smallest complete mitogenome found in isopods (Table 1) thus far. The genome contains the same 13 protein-coding genes and two ribosomal RNAs found in most metazoans. However, it exhibits incomplete tRNA-encoding capacity (18 tRNAs instead of the usual 22, see Table 2). Twenty-one mt genes are transcribed from one strand (the plus strand) and the remaining 12 from the other (the minus strand). A total of 817 noncoding intergenic nucleotides were found, with the largest continuous region (362 bp, AT% = 60.2) located between trnA and cob. Due to its location and AT richness, we predict that this part of the genome is the mt control region and that it likely contains the origins of replication and regulatory elements for transcription.
However, this predicted control region is considerably smaller than that in other isopods (Table 1) and shows no sequence similarity with any previously sequenced Isopoda control region (data not shown).

| Protein-coding genes and ribosomal RNAs
With regard to protein-coding genes, ten (cox1-3, atp8, atp6, nad1-3, and nad5-6) are encoded by the plus strand and the remaining three (cob, nad4, and nad4L) the minus strand ( Figure 1, Table 2). This orientation is shared by all isopod mitogenomes sequenced to date, except for nad1. The 13 protein-coding genes appear to start with the codon ATN, which is typical for metazoan mitogenomes (Wolstenholme, 1992), and TAN is the termination codon for 11.
The rrnL and rrnS genes of Bathynomus sp. are 1,181 bp (AT% = 60.5) and 742 bp (AT% = 59.7) in length, respectively. These lengths are typical for crustaceans, whereas the AT contents are slightly lower than those of other isopods (Table 1).

| Transfer RNAs
We analyzed the entire genome sequence of Bathynomus sp. and successfully identified 18 tRNA genes based on their potential secondary structures ( Figure 2) using the MITOS webserver. Moreover, tRNAscan-SE was used to recheck five genes and ARWEN for 17 genes.
The 18 tRNA genes are spread throughout the entire genome and are located on both strands (Figure 2, Table 2). Sixteen of these genes display a common clover-leaf secondary structure, and the remaining two exhibit a t-shaped secondary structure.

| Phylogenetic analysis and evolutionary rate estimation
We examined both the nucleotide and amino acid sequences of protein-coding genes in intraorder phylogenetic analyses. For each dataset, ML and Bayesian analyses resulted in the same tree topology ( Figure 5). The evolutionary rate was estimated using the nucleotide sequences of 13 protein-coding genes and shown in F I G U R E 2 Putative secondary structures for the 18 transfer RNAs of the Bathynomus sp. mitogenome. All 18 structures were generated using the MITOS webserver and verified using tRNAscan-SE and ARWEN. Most tRNAs feature a standard clover-leaf structure. Exceptions: The TΨC arm is absent from trnQ, and the DHU arm is absent from trnC 4 | DISCUSSION

| Transfer RNAs
The DHU arm of trnC and the TΨC arm of trnQ are completely missing ( Figure 2). Loss of the DHU arm from tRNA-C has also been observed in other isopods, such as Ligia oceanica (Kilpert & Podsiadlowski, 2006) and Eophreatoicus sp.-14 (Kilpert & Podsiadlowski, 2010). This feature might be a putative autapomorphy of isopoda, as it is absent in other malacostracan crustaceans. In general, tRNAs with a U at the wobble position (the first position) of the anticodon recognize either fourfold degenerate or NNR codons, whereas those with a G at this position only recognize NNY codons (Yang & Yang, 2008). All tRNAs of Bathynomus sp. mitogenome obey this rule.
Despite extensive efforts to identify secondary structure in noncoding regions, the genes trnH, trnI, trnL1, and trnS1 were not found in the mitogenomic sequence. Considering the absence of trnQ in the L. oceanica mitogenome (Kilpert & Podsiadlowski, 2006), this is not an exception among isopods. Indeed, a lack of these genes is also common in arthropods, such as the absence of trnQ in Aleurodicus dugesii, trnS in Schizaphis graminum (Thao, Baumann, & Baumann, 2004), and trnD in Centruroides limpidus (Davila, Pinero, Bustos, Cevallos, & F I G U R E 3 Comparison of mitochondrial gene arrangements in 12 species in the order Isopoda. In addition, the pancrustacean ground pattern (putative ancestral state) is provided. The mt complete genomes of Bathynomus sp., Limnoria quadripunctata, Ligia oceanica, and Eophreatoicus sp.-14 are available; the other mt genomes are incomplete and are missing some tRNAs and the control region (CR). All tRNAs are designated by single letters (except L1, L2, S1, and S2 for trnL-CUN, trnL-UUR, trnS-AGN, and trnS-UCN, respectively). Colored genes denote translocated genes in comparison with the pancrustacean ground pattern. Red indicates genes that share a derived position in isopods. Uniquely derived gene positions of individual species are depicted in yellow Davila, 2005). tRNA gene deficiencies have often been observed in protozoans, fungi, algae, plants, and lower metazoans (Schneider & Marechal-Drouard, 2000); in these cases, imported nuclear DNAencoded tRNAs compensate for the lack of mt tRNA. Furthermore, marsupial mt tRNAs exhibit interesting patterns. Janke and Paabo identified in Didelphis virginiana a pseudogene similar to trnD with an anticodon of GCC instead of the usual GUC (Janke & Paabo, 1993); the authors found that the cytosine is changed to uridine during posttranscriptional RNA editing. Based on its nonfunctional secondary structure, Dorner, Altmann, Paabo, and Morl (2001) discovered a trnK pseudogene in the same marsupial species, with compensation due to import of cytoplasmic tRNA rather than RNA editing. As there appears to be no candidate template for RNA editing in the case of Bathynomus sp., import of tRNA is a plausible explanation. In addition, some other explanations also may be plausible explanation: (1) specific bias in codon usage. Indeed, a recent work has succeeded in replacing seven codons in Escherichia coli (at least part of it), so that the genome can be decoded on 57 triplets rather than 64 (Thiele et al., 2012). It could be that some codons are missing, so that some tRNAs would not be necessary; (2) a hypothesis could be that some of the present tRNA would be edited, so that they would recognize both different amino acids and have different anticodons. This is not entirely impossible, as it is known, in some cases, that the specification of the amino acid loaded on the tRNA uses the anticodon as a signature. In a way, this would require much less genes from the host (just editing) than create a whole RNA translocation across the double membrane of mitochondria.

| Genome rearrangement
Compared to the gene orders of other isopod species, the Bathynomus sp. mitogenome has undergone significant changes in gene arrangement. First, all 11 isopod mt genomes differ in gene order, but, except for Limnoria quadripunctata, most of the differences are limited to the position of one or a few tRNA genes. Second, the gene order of the isopod mt genomes can be clearly distinguished from the pancrustacean ground pattern. Inferred changes in the Bathynomus sp. genome involve tRNA genes, protein-coding genes, rRNA genes, and the mt control region (CR). Kilpert et al. showed that isopods share a derived order of genes, with individual species experiencing modifications. Using the similarities common to most isopods, a most parsimonious hypothesis for the mt gene order of the isopod ground pattern was inferred, which includes a unique arrangement of nad1, trnL1, rrnS, CR, trnS1, cob, trnT, nad5, and trnF (Kilpert et al., 2012). For Bathynomus sp., cob, trnT, nad5, and trnF fit the isopod ground pattern, whereas nad1, trnL1, and rrnS strongly diverge from this pattern.
We also found that the mitogenome of Bathynomus sp. harbors significant alterations in gene arrangement compared to the gene order of the putative isopod ground pattern described by Kilpert et al. (2012). In total, we identified rearrangements in this species (Figure 4): One rearrangement involves protein-coding genes, and the remaining are tRNA translocations. Relative to its ancestral position, trnR is at a position upstream of nad3. Glyptonotus cf. Antarcticus (Kilpert et al., 2012), Cylisticus convexus, and Trachelipus rathkei (Chandler et al., 2015) share this rearrangement. The positions of trnE and trnW are changed to positions downstream of rrnS, which has not been found in other isopods to date. Finally, the nad1 gene is moved to a position downstream of nad3 and from the minus strand to the plus strand and this translocation is novel feature found only in isopods (Kilpert & Podsiadlowski, 2010;Kilpert et al., 2012;Marcade et al., 2007).

| Phylogenetic analysis and evolutionary rate estimation
Indeed, the nucleotide and amino acid trees are nearly the same, except for the position of the Asellota clade (Asellus aquaticus) and the Phreatoicidea clade . In the case of the nucleotide dataset, the Asellota clade is placed at the base of Isopoda ( Figure 5a) and is well supported (ML/BPP = 100/1.00). However, in the amino acid tree, the Phreatoicidea clade is placed at the base of Isopoda (Figure 5b), which is in agreement with the findings of previous research (Kilpert et al., 2012) and is also well supported (ML/BPP = 100/1.00). This discrepancy may be due to the fact that the two clades are represented by only one species each. analyses of the two dataset ( Figure 5). However, as research on the phylogenetic relationships within this infraorder is rare (Kilpert et al., 2012), intensive sampling and analysis of a greater number of species are necessary to determine the phylogenetic relationships among members of Isopoda.
We found that the evolutionary rate of Bathynomus sp. is the slowest measured to date within Cymothoida (Table 3), which means that this species may have experienced nonsynonymous mutations harmful to its survival and that its mt genes have been under strong purifying selection (Yang, 2007). This slow evolutionary rate may be associated with the relatively constant deep-sea environment of this species. For F I G U R E 5 Phylogenetic trees showing relationships among isopods based on nucleotide (a) and amino acid datasets (b). Six decapod species served as outgroups.
Branch lengths and topologies are derived from Bayesian analyses. Numbers next to nodes specify bootstrap percentages from the maximum likelihood (ML) analysis plus Bayesian posterior probabilities (BPP). Specifically, <50 indicates bootstrap percentages from ML values below 50% example, Bathynomus giganteus was already present as early as 160 million years ago, and a previous study showed that B. giganteus specimens from Australia, Mexico, and India are almost exactly the same, possibly because of their nearly identical environments (Parker, 2003).
This phenomenon has been found in many high-altitude habitats species (Shen et al., 2009;Sun et al., 2011).
This study is the first determination of the mitogenome of a deepsea member of the order Isopoda, an effort that not only increases sampling of Isopoda but also more importantly provides useful information for understanding the evolution of deep-sea isopods.

ACKNOWLEDGMENTS
The specimens used in the study were collected by the "ROV Lander