Characterization of Complete Mitochondrial Genome and Phylogeny of Sepia lycidas (Sepioidea, Sepiidae)

Sepia lycidas (Sepioidea, Sepiidae) is a demersal and neritic species (Dong, 1988), which mainly occurs throughout the East China Sea (including areas near southern Japan), South China Sea, and nearby regions of the Indian and West Pacific Oceans (Carpenter and Niem 1999; Okutani, 2005). Known as the kisslip cuttlefish, it is a common and large cuttlefish (38 cm adult mantle size), and has high nutritional and medicinal value (Dong, 1988). Extensive studies of S. lycidas have been focused on breeding, morphological, and biological characteristics (Natsukari and Tashiro, 2009; Nagai et al., 2001; Lucky et al., 2012), but preliminary studies about the phylogenetic analysis were based on single gene fragment (Montserrat et al., 2010; Wen et al., 2017). We sequenced the complete mt genome of S. lycidas (Genbank accession number: KJ162574) to assess its genomic structure and provide a more robust estimate of its phylogenetic relationships.

Most metazoan species possess a compact, circular mt genome, which varies in size from 15 to 20 kb that typically contain 37 genes, including 13 protein coding genes, two ribosomal RNA (rRNA) genes, and 22 transfer RNA (tRNA) genes necessary for translation of the proteins encoded by the mtDNA (Boore, 1999; Cheng et al., 2012). MtDNA has been extensively used for studying phylogenetic and evolutionary relationships among animal species, due to its maternal inheritance, rapid evolutionary rate, and lack of genetic recombination (Zheng et al., 2004; Cheng et al., 2013; He et al., 2016; Liu et al., 2016). Partial sequences of mtDNA genes, such as cytochrome oxidase I (COI), cytochrome oxidase III (COIII), and 16S rRNA, have proved to be important tools in intraspecific and interspecific phylogenetic studies of Cephalopoda and other mollusks (Akasaki et al., 2006; He et al., 2016; Mao et al., 2016; Wang et al., 2016). Compared to partial mt genes, complete mtDNA sequences can uncover gene rearrangements and other variation at the genome level for all phyla, and are especially useful because they have sufficient interspecies sequence variability for resolving species-level relationships (Cheng et al., 2013). The levels of variation of different genes vary significantly, making different genes useful for phylogenetic analysis at different taxonomic levels (Zheng et al., 2015).

The ability of different cephalopods to efficiently identify members of their own species is very important for breeding. Zheng (2001) showed that there were significant differences in radular structure among nine species of cuttlefish and suggested that radular morphology is informative for their taxonomy. This suggestion has been reinforced by other studies (Ma et al., 2016). Radular morphology could be used as a basis for classification of the cuttlefishes in the coastal areas of China. Morphological traits can be difficult to apply for classification of cephalopods due to phenotypic plasticity, overlapping geographical distributions and notable differences during growth and development. In addition, cephalopods are found throughout variable marine habitats and have diverged in response to a variety of ecological pressures (Boyle and Boletzky, 1996; Hanlon and Messenger, 1996; Allen et al., 2014). An understanding is needed of the phylogenetic relationships among genera of Sepiidae, especially when studying the differences among species within different groups.

Wakabayashi et al. (2012) studied phylogenetic relationships in Family Ommastrephidae, based on mtDNA 16S rRNA and COI genes. The relationships among subspecies in the resulting maximum parsimony trees were consistent with neighbor joining trees produced using allozyme data (Yokawa, 1994). However, these relationships differed from those inferred from analyses of morphological characters (Roeleveld, 1988). The first determination of the whole genome sequence of the human mt genome was made in 1981 (Anderson et al., 1981). Since then, genome-level sequencing has provided many new insights into species relationships. Both molecular and morphological analyses are needed to provide the most informative classification within cuttlefish and among higher taxonomic ranks.

Specimens of S. lycidas were collected from the Zhanjiang fishing ground (GuangDong Province, 2011.8). The sample was preserved at -80°C until use. Total DNA was extracted from wrist muscle tissue of S. lycidas with the conventional phenol-extraction method (Sam, 1989).

We used 10 pairs of primers that amplify contiguous, overlapping segments of the complete mt genome of S. lycidas (Table I). The primers were designed from the complete mt genome sequence of S. officinalis (NC007895, Akasaki et al., 2006), S. japonica (NC028731, Zheng et al., 2015) and S. pharaonis (NC021146, Wang and Wu, 2014). The Long-PCR reaction volume was 50 μl, containing 31.5 μl sterile deionized water, 5.0 μl 10× LA PCR Buffer (Mg2+ plus), 8.0 μl dNTPs (2.5 mM each), 1 μl of each primer (25 pmol/ml), 0.5 μl LA Taq DNA polymerase (Takara), and 3 μl DNA template, with the following reaction parameters: 2 min high-denaturation (93°C), followed by 20 cycles of 10 s denaturation (92°C), 30 s annealing (53°C), and a 10 min extension (68°C), with a 7 min final extension (68°C).

The PCR products were electrophoresed on a 1% agarose gel and stained with ethidium bromide for band characterization via ultraviolet transillumination. The products were recovered and concentrated using a Gel Extraction Kit (CW2302, CWBIO) according to the manufacturer’s protocol.

The PCR products were ligated into pMD 18-T vectors and cloned using a TOPO TA Cloning Kit (Invitrogen, USA) according to the manufacturer’s instructions. Each cloned DNA fragment was sequenced according to the manufacturer’s protocol with TOPO vector inner primers T7 and M13R (Invitrogen, USA). Labeled fragments were analyzed on an ABI Prism 3730 DNA Sequencer (Applied Biosystems, USA).

The characterization of the mt genome of S. lycidas was performed in comparison with other cuttlefish mt genomes (Akasaki et al., 2006; Yokobori et al., 2007; Wang et al., 2014). MEGA 6.06 software (Tamura et al., 2007) was used for sequence splicing and analyses. Concatenated nucleotide sequences were loaded into DNAStar 7.1 statistical software packages to analyze the full mtDNA sequence, nucleotide content, codon usage, and amino acids. Tandem Repeats Finder (http://tandem.bu.edu/trf/trf.html; Benson, 1999; Huang et al., 2017) was used to analyze sequences that were rich with AT-repeat regions. Finally, we used the online mapping software OGDraw v1.2 (http://ogdraw.mpimp-golm.mpg.de/index.shtml) to diagram the S. lycidas mtDNA circular structure. Protein coding genes were analyzed by ORF Finder using the invertebrate mitochondrial code. Base composition and codon usage were calculated with DNAStar software.

All 13 protein-coding gene sequences of mtDNA were used to estimate phylogenetic relationships among 37 species in Cephalopoda, using maximum likelihood (ML) and Bayesian inference (BI). The mtDNA protein-coding regions were concatenated and aligned using the ClustalW algorithm with default parameters implemented in MEGA 6.06 (Tamura et al., 2013). In order to discuss phylogenetic relationships within Sepiidae and with other families of Cephalopoda (Ommastrephidae, Loliginidae, Vampyroteuthidae, and Octopodidae) thirty-seven species in Genbank were included (Table II). Allonautilus scrobiculatus (NC. 026997.1) and Nautilus macromphalus (NC. 007980.0) were used as outgroup taxa. The ML tree was evaluated with 1000 bootstrap replicates in MEGA 6.06.

MrBayes ver 3.1.2 (Huelsenbeck et al., 2001) was used for Bayesian analyses under the Mtzoa-F + C4 model to partition different amino acids positions. Two runs with four chains of MCMC iterations were performed for 100,000 generations, sampling trees every 100 generations.

The complete mt genome of S. lycidas is 16,228 bp in length and consists of 13 protein-coding genes, 22 transfer RNA genes (tRNA), 2 ribosomal RNA genes (rRNA), and 2 long-noncoding regions (CR). The structural organization of the complete mt genome is given in Figure 1 and Table III. The base composition is 35.8% T, 14.8% C, 41.3% A, and 8.1% G. As seen in other invertebrate mt genomes, A+T content (77.1%; Table IV) is much higher than the G+C content.

The 13 protein-coding genes of S. lycidas can be classified into two categories: COIII, ND2, COI, COIII, ATP8, ATP6, and ND3 are encoded by the light strand, whereas ND1, ND5, ND4, ND4L, Cyt b, and ND6 are encoded by the heavy strand (Table III). The overall A+T content of the protein genes is 76.2%, with the highest in ATP8 (83.1%), and the lowest in COI (70.1%). AT and CG skew of the 13 protein-coding genes is 0.087 and -0.294, respectively (Table IV).

Most of the protein-coding genes start with an ATG initiation codon, except for ND4, which starts with ATA. Ten protein-coding genes (COI-III, ND1-5, ATP8, and ND4L) use TAA as the termination codon, while the other three (ATP6, Cyt b, and ND6) share the termination codon TAG. The pattern of codon usage of S. lycidas mtDNA is shown in Table V. The most frequently used amino acids are Leu (15.7%), followed by Ile (10.2%), Phe (8.5%), Met (7.9%), Ser (7.3%), and Gly (6.0%). The least frequently used amino acid is Gln (1.6%). The most frequently used codon in S. lycidas and other cuttlefish is UUA (appearing 456 times).

A total of 22 transfer RNA genes are interspersed in the S. lycidas mt genome and range from 63 bp to 94 bp, with a base composition of 39.2% T, 12.6% C, 39.2% A,

and 9.0% G (Table IV). In addition, tRNASer and tRNALeu both appear twice; the other 18 types of tRNA appear only once. Nineteen of the tRNA are encoded by the light strand and three are encoded by the heavy strand. The overall AT and CG skew of the 22 tRNA genes is 0 and -0.167, respectively. The CG skew of tRNAArg (CGR), tRNAAsp (GAD), and tRNASer (AGS) is zero.

S. lycidas has two ribosomal RNA genes, 12S rRNA (981 bp) and 16S rRNA (1274 bp). They are located between tRNALeu and tRNACys genes, and separated by the tRNAVal gene. The base composition of the two rRNAs gene sequences is 42.0% A, 37.6% T, 13.7% C, and 6.7% G. The overall AT and CG skew of the rRNA genes is 0.055 and -0.343. The A+T content of 12S rRNA and 16S rRNA is 80.0% and 79.3%, respectively (Table IV).

The complete mt genome of S. lycidas has two long, noncoding regions that are both control regions, regulating replication and transcription. Their overall base composition is rich in A and T (Jing et al., 2016). Nucleotide composition across the non-coding regions is 42.0% A, 37.6% T, 13.7% C, and 6.7% G. The A+T content is 79.6%. AT and CG skew is 0.055 and -0.343, respectively (Table IV).

Two non-coding regions are both 577 bp in length and contain termination-associated sequences and conserved sequence blocks: one is located between tRNAGly and tRNAAsn, another is located between tRNAGlu and the protein-coding gene COIII. The control region in S. lycidas exhibits the typical tripartite structure with an extended termination-associated sequence domain (ETAS), central conserved sequence block domain (CD), and conserved sequence block domain (CSB). The possible function of the LNCR as the control region is deduced from the potential stem and loop structures (Tomita et al., 2002). The termination-associated sequence (TAS) domains, which are thought to act as a signal for the termination of heavy strand elongation, were identified in ETAS. This is a hypervariable domain that may be useful for analyzing interspecies variation within S. lycidas (Southern et al., 1988).

Total codon number: 3736.

Many systematic and population genetic studies have been based on genetic markers in the mt genomes at either the nucleotide or amino acid level (Zou et al., 2011). Phylogenetic relationships inferred from ML and Bayesian analyses were basically consistent with each other (Figs. 2 and 3) and with the existing morphological classification.

The use of long DNA sequences can help to resolve major phylogenetic relationships and provide resolution of closely related species. In our study, the phylogenetic trees were inferred using protein sequences. In both phylogenetic trees, the ten sampled species of Sepiidae were divided into three clades. S. lycidas forms a clade with S. pharaonis, S. aculeata and S. esculenta.

Both the ML and BI analyses divided Decapodiformes into three parts: Oegopsida, Myopsida, and Sepiida; the Oegopsida consisted of Ommastrephidae, Architeuthidae, Bathyteuthidae, and Enoploteuthidae; the Myopsida only had one family and the Loliginidae contained eight species. Octobrachia can be divided into two suborders: Cirrina and Incirrina. However, there are no complete mt genomes available for Cirrina, and so our analyses only include Incirrina. Octopus conispadiceus is separated from other Octopodidae in both trees (Tanner et al., 2017).

The position of Ommastrephes bartramii differed between the ML and BI trees. In the ML tree, O. bartramii was sister to Sthenoteuthis oualaniensis, whereas in the BI tree O. bartramii was sister to the remainder of Ommastrephidae (Figs. 2 and 3).

The mt genomes of Ommastrephidae and Enoploteuthidae contain six duplicated genes, including COI, COII, COIII, ATP 6 and 8, and tRNAAsp (Fig. 3). These genes occur between 12S RNA and 16S RNA (Yokobori et al., 2007). Each gene in the mt genomes for these taxa possesses specific functions. The level of variation between duplicate genes is low, and every gene copy contains several insertions and deletions (Kawashima et al., 2013). These duplicated genes likely have lost their function, but the secondary structures and anticodon positions of functional genes would be disrupted without them (Eda et al., 2010). These duplications are characteristic of all members of Ommatrephida. It is likely that the octopus-type mt genome is most similar to the ancestral state, this state being maintained from at least the Cephalopoda ancestor to the common ancestor of Oegopsida, Myopsida and Sepiolida (Kawashima et al., 2013).

The previous classification for cephalopoda was based solely on morphological traits (Chen et al., 2009). The three sampled members of the family Ommastrephidae (Dosidicus gigas, O. bartramii, and S. oualaniensis) belong to the subfamily Ommastrephinae and were monophyletic in the ML tree. O. bartramii is sister to S. oualaniensis in our ML tree, but is placed more distantly in our BI tree and in an ML tree base on COI (Wakabayashi et al., 2012; Tanner et al., 2017).

The largest difference between our results and some of the previous studies is in the Octopodidae. Vampyroteuthis infernalis was separated from the remainder of Octopodidae in our phylogenies, as was also found by Allcock et al. (2011). However, other phylogenetic analyses of the mt genomes of V. infernalis and other Octopodidae suggested a very close relationship (Kawashima et al., 2013). In both the ML and BI trees, Amphioctopus aegina and O. ocellatus are sister species despite belonging to different genera in the traditional taxonomy (Chen et al., 2009). Additional data from the nuclear genome are needed to further test the placement of these species.

The relationships that we inferred within family Sepiidae (Figs. 2 and 3) are consistent with the traditional morphological classification (Chen et al., 2009). We found a close relationship between S. apama and S. latimanus, as was also seen by Akasaki et al. (2006). However, our results differ from some previous studies in the position of S. officinalis and the relationships among S. esculenta, S. aculeata, S. lycidas, and S. pharaonis. Yokobori et al. (2007) produced phylogenies based on the DNA sequences of 16S, 12S, and COI, and showed a further relationship between S. lycidas and S. esculenta. In comparison, our trees place S. lycidas sister to an S. aculeata + S. pharaonis clade, which was same as the ML phylogenetic relationships based on nucleotide and amino acid sequence data (Zhang et al., 2015; Groth et al., 2015; Strugnell et al., 2017). Our MP and BI analyses of amino acid sequences allowed inferences of phylogenetic relationships of different families with results that were consistent with Kawashima et al. (2013).

To fully resolve relationships and provide a robust classification, the valuable and extensive information available in mt genomes should be combined with other sources of data, such as nuclear genes and morphological and ecological characteristics.

This work was funded by the project of Hong Kong, the Macao and Taiwan science and technology cooperation project (2014DFT30120), Zhejiang Provincial Natural Science Foundation of China (LY14C190002), National Natural Science Foundation of China (41576131, 41406138), National Spark Program (2015GA700014), Zhejiang Provincial Natural Scicence Foundation of China (LY17C190006) and Zhejiang science and technology Program (2015F50055).

Authors have declared no conflict of interest.