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

Baoying Guo*1, Yu Chen1, Chuan Zhang2, Zhenming Lv1, Kaida Xu3, Hongling Ping4 and Huilai Shi4

1National Engineering Research Center of Maricultural Facilities of China, College of Marine Science and Technology, Zhejiang Ocean University, Zhoushan 316022, PR China

2Zhoushan Fisheries Research Institute, Zhoushan 316000, PR China

3Scientific Observing and Experimental Station of Fishery Resources for Key Fishing Grounds, MOA, Key Laboratory of Sustainable Utilization of Technology Research for Fisheries Resources of Zhejiang Province, Zhoushan 316021, China

4Marine Fisheries Research Institute of Zhejiang, Zhejiang Province Key Lab of Mariculture and Enhancement, Zhoushan 316021, PR China

ABSTRACT

We sequenced the complete mitochondrial (mt) genome of the kisslip cuttlefish, Sepia lycidas, made the comparison with the mt genomes of other cuttlefishes, and constructed phylogenetic trees estimating their relationships. The genome was 16,228 bp and contained 13 protein-coding genes, 22 transfer RNA genes, 2 ribosomal RNA genes, and 2 long-noncoding regions [both in the control regions (CR)]. The composition and order of genes in S. lycidas were similar to those of most other invertebrates. The overall base composition of S. lycidas is 35.8% T, 14.8% C, 41.3% A, and 8.1% G, with extremely high A+T content (77.1%). Both control regions contain termination-associated sequences and conserved sequence blocks. Maximum likelihood and Bayesian methods were used to build phylogenetic trees based on protein coding mtDNA genes of 37 cuttlefish species. S. lycidas have a close relationship with S. pharaonis, S. aculeata, and S. esculenta. This result confirmed the relationships of S. lycidas as being similar to the traditional taxonomy. This study plays an important role in the investigation of phylogenetic relationships, taxonomic resolution, and phylogeography for Sepiidae species.

Article Information

Received 07 July 2017

Revised 30 August 2017

Accepted 23 September 2017

Available online 25 June 2018

Authors’ Contribution

BG and ZL conceived the idea. BG, YC and CZ designed the experiments and wrote the article. KX, HP and HS provided the specimens. CZ conducted the experimental work. YC analyzed the sequence data and submitted to the Genbank.

Key words

Cephalopoda, Sepiida, Sepia lycidas, Mitochondrial genome, Phylogeny.

DOI: http://dx.doi.org/10.17582/journal.pjz/2018.50.4.1497.1508

* Corresponding authors: guobaobao2000@126.com

0030-9923/2018/0004-1497 $ 9.00/0

Introduction

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.

Table I.- Primers used for amplification of the Sepia lycidas mitochondrial genome.

Primer pair	Primer sequence (5′-3′)	Region	Size (bp)
1 F	GTRGGWATAGAYGTWGATACACGAGCYTATT	2934-8397	5463
1 R	TGTGCCAGCATYYGCGGTT
2 F	TAAAGAATAATAGGGTATCTAATCCTAGT	8187-9640	1454
2 R	TTTATAACATCAATATAAYCCGTTCA
3 F	GAGGCMTTTAACTGTTAATTAAAT	9548-11260	1713
3 R	AATTGCVGGDATAACTAAAAGAGC
4 F	GATATTATTATTACWCCYAATTGACT	10964-12598	1635
4 R	AATTGCAGGDTCWATRATTTTAGC
5 F	AAWGAYTTAAYATCMCTTTGWCC	12416-14063	1648
5 R	TGARAATTTTATTCCBGCTAAYCC
6 F	GGAATWGAACGTAAAATWGCAT	13968-15463	1496
6 R	GCTAARWHTTWAAGCTATTGGGTTC
7 F	GGRATWGCYGAWACTAAATTAG	14361-15772	1412
7 R	ATAWGCTMRAGGGATGTTTGAGAG
8 F	GGGTATGAACCCAATAGCTT	15433-16172, 1-808	1548
8 R	GCTTAAATTCGGCCACTTAAT
9 F	TCHACYTTYTTTGTAGCTACAGG	577-2087	1511
9 R	RTGRTTWGTKGAGAAWARTCATCG
10 F	CCWAAYWTAACYAAACAAATAACTTG	1400-3041	1642
10 R	WGAWCCATARATAGTRGCTAATCA

Materials and Methods

Sample and DNA extraction

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).

Polymerase chain reaction (PCR)

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).

Table II.- The complete mitochondrial genome sequence in this study.

Species	Length(bp)	Taxonomic status		Accession No.
Bathyteuthis abyssicola	20,075	Bathyteuthidae	Oegopsida	NC_016423
Doryteuthis opalescens	17,370	Loliginidae	Myopsida	KP336703
Loligo bleekeri	17,211	Loliginidae	Myopsida	NC_002507
Loliolus beka	17,483	Loliginidae	Myopsida	NC_028034
Loliolus uyii	17,134	Loliginidae	Myopsida	NC_026724
Sepioteuthis lessoniana	16,694	Loliginidae	Myopsida	KM878671
Uroteuthis chinensis	17,353	Loliginidae	Myopsida	NC_028189
Uroteuthis duvauceli	17,413	Loliginidae	Myopsida	NC_027729
Uroteuthis edulis	17,360	Loliginidae	Myopsida	NC_017746
Allonautilus scrobiculatus	16,132	Nautilidae	Nautiloidea	NC_026997
Nautilus macromphalus	16,258	Nautilidae	Nautiloidea	NC_007980
Amphioctopus aegina	15,545	Octopodidae	Octopoda	NC_029702
Cistopus chinensis	15,706	Octopodidae	Octopoda	KF017606
Cistopus taiwanicus	15,793	Octopodidae	Octopoda	NC_023257
Octopus bimaculatus	16,084	Octopodidae	Octopoda	KT581981
Octopus bimaculoides	15,733	Octopodidae	Octopoda	NC_029723
Octopus conispadiceus	16,027	Octopodidae	Octopoda	NC_029747
Octopus minor	15,974	Octopodidae	Octopoda	NC_015896
Octopus ocellatus	15,979	Octopodidae	Octopoda	NC_007896
Octopus vulgaris	15,744	Octopodidae	Octopoda	NC_006353
Architeuthis dux	20,333	Architeuthidae	Oegopsida	KC701764
Watasenia scintillans	20,089	Enoploteuthidae	Oegopsida	KJ845633
Dosidicus gigas	20,324	Ommastrephidae	Oegopsida	NC_009734
Illex argentinus	20,278	Ommastrephidae	Oegopsida	NC_026908
Ommastrephes bartramii	20,308	Ommastrephidae	Oegopsida	NC_020348
Sthenoteuthis oualaniensis	20,306	Ommastrephidae	Oegopsida	NC_010636
Todarodes pacificus	20,254	Ommastrephidae	Oegopsida	NC_006354
Sepia aculeata	16,219	Sepiidae	Sepiida	KF690633
Sepia apama	16,184	Sepiidae	Sepiida	NC_022466
Sepia esculenta	16,199	Sepiidae	Sepiida	NC_009690
Sepia latimanus	16,225	Sepiidae	Sepiida	NC_022467
Sepia lycidas	16,228	Sepiidae	Sepiida	KJ162574
Sepia officinalis	16,163	Sepiidae	Sepiida	NC_007895
Sepia pharaonis	16,208	Sepiidae	Sepiida	NC_021146
Sepiella inermis	16,191	Sepiidae	Sepiida	KF040369
Sepiella japonica	16,170	Sepiidae	Sepiida	NC_028731
Vampyroteuthis infernalis	15,617	Vampyroteuthida	Vampyromorphida	NC_009689

Table III.- Mitochondrial genome characteristics of S. lycidas.

Gene	Position		Size (bp)		Codon		Intergenic nucleotidesa	Strand
Gene	From	To	Nucleotide	Amino acid	Initiation	Stop	Intergenic nucleotidesa	Strand
CO3	1	780	780	259	ATG	TAA	0	L
tRNALys	791	858	68				10	L
tRNAAla	857	923	67				-2	L
tRNAArg	923	989	67				-1	L
tRNASer	993	1058	66				3	L
ND2	1059	2090	1032	343	ATG	TAA	0	L
CO1	2068	3600	1533	510	ATG	TAA	-23	L
CO2	3603	4289	687	228	ATG	TAA	2	L
ATP8	4354	4548	195	64	ATG	TAA	64	L
ATP6	4511	5203	693	230	ATG	TAG	-38	L
tRNAPhe	5230	5293	64				26	L
ND1	5295	6233	939	312	ATG	TAA	1	H
tRNALeu	6234	6306	73				0	L
tRNALeu	6307	6381	75				0	L
16S rRNA	6382	7655	1274				0	L
tRNAVal	7656	7726	71				0	L
12S rRNA	7727	8707	981				0	L
tRNACys	8708	8770	63				0	L
tRNATyr	8769	8832	64				-2	L
tRNAGln	8834	8898	65				1	L
tRNAGly	8905	8977	73				6	L
Control region	8982	9558	577				4	L
tRNAAsn	9561	9626	66				2	H
tRNAIle	9628	9692	65				1	L
ND3	9693	10,046	354	117	ATG	TAA	0	L
tRNAAsp	10,052	10,117	66				5	H
ND5	10,129	11,850	1722	573	ATG	TAA	11	H
tRNAHis	11,851	11,915	65				0	L
ND4	11,914	13,275	1362	453	ATA	TAA	-2	H
ND4L	13,272	13,568	297	98	ATG	TAA	-4	H
tRNAThr	13,575	13,639	65				6	L
tRNASer	13,640	13,705	66				0	H
Cyt b	13,705	14,844	1140	379	ATG	TAG	-1	H
ND6	14,837	15,349	513	170	ATG	TAG	-8	H
tRNAPro	15,351	15,418	68				1	L
tRNAMet	15,419	15,486	68				0	L
tRNATrp	15,489	15,555	67				2	L
tRNAGlu	15,557	15,650	94				1	L
Control region	15,652	16,228	577				1	L

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.

Cloning and sequencing of the PCR products

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).

Sequence analysis

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.

Phylogenetic analysis

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.

Results

Genome organization and composition

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.

Protein-coding genes and codon usage

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).

Transfer RNA genes and RNA genes

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,

Table IV.- Nucleotide composition in mitochondrial genome of S. lycidas.

Species		Size (bp)	A%	T%	C%	G%	A+T%	C+G%
S. lycidas	Whole genome	16,228	41.3	35.8	14.8	8.1	77.1	22.9
	13 PCGs	11,174	41.4	34.8	15.4	8.4	76.2	23.8
	tRNA	1,501	39.2	39.2	12.6	9.0	78.4	21.6
	rRNA	2,255	42.0	37.6	13.7	6.7	79.6	20.4
	2 LNCRs	1,154	41.6	36.0	14.4	8.0	77.6	22.4

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).

Non-coding regions

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).

Table V.- Frequency and count for genetic codons and codon usage in Sepia lycidas mitochondrial genome.

Codon /1000(C)	Codon /1000(C)	Codon /1000(C)	Codon /1000(C)
UUU 280(74.9)	UCU 84(22.5)	UAU 153(41.0)	UGU 61(16.3)
Phe	Ser	Tyr	Cys
UUC 38(10.2)	UCC 15(4.0)	UAC 24(6.4)	UGC 7(1.9)
Phe	Ser	Tyr	Cys
UUA 456(122.1)	UCA 94(25.2)	UAA 10(2.7)	UGA 81(21.7)
Leu	Ser	Stop	Trp
UUG 47(12.6)	UCG 4(1.1)	UAG 3(0.8)	UGG 17(4.6)
Leu	Ser	Stop	Trp
CUU 39(10.4)	CCU 64(1.7)	CAU 70(18.7)	CGU 23(6.1)
Leu	Pro	His	Arg
CUC 4(1.1)	CCC 4(1.1)	CAC 18(4.8)	CGC 5(1.3)
Leu	Pro	His	Arg
CUA 39(10.4)	CCA 49(13.1)	CAA 55(14.7)	CGA 19(5.1)
Leu	Pro	Gln	Arg
CUG 0	CCG 1(0.3)	CAG 6(1.6)	CGG 6(1.6)
Leu	Pro	Gln	Arg
AUU 355(95.0)	ACU 75(20.1)	AAU 162(43.4)	AGU 64(17.1)
Ile	Thr	Asn	Ser
AUC 27(7.2)	ACC 12(3.2)	AAC 21(5.6)	AGC 18(4.8)
Ile	Thr	Asn	Ser
AUA 260(69.6)	ACA 64(17.1)	AAA 91(24.4)	AGA 47(12.6)
Met	Thr	Lys	Ser
AUG 35(.4)	ACG 0	AAG 9(2.4)	AGG 16(4.3)
Met	Thr	Lys	Ser
GUU 90(24.1)	GCU 75(20.1)	GAU 60(16.0)	GGU 128(34.3)
Val	Ala	Asp	Gly
GUC 2(0.5)	GCC 15(4.0)	GAC 10(2.7)	GGC 5(1.3)
Val	Ala	Asp	Gly
GUA 97(26.0)	GCA 40(10.7)	GAA 73(19.5)	GGA 60(16.1)
Val	Ala	Glu	Gly
GUG 11(2.9)	GCG 3(0.8)	GAG 15(4.0)	GGG 31(8.3)
Val	Ala	Glu	Gly

Total codon number: 3736.

Phylogenetic status of S. lycidas

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).

Discussion

Gene arrangement

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).

Phylogenetic analysis

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.

Acknowledgments

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).

Statement of conflict of interest

Authors have declared no conflict of interest.

References

Akasaki, T., Nikaido, M., Tsuchiya, K., Segawa, S., Hasegawa, M. and Okada, N., 2006. Extensive mitochondrial gene arrangements in coleoid Cephalopoda and their phylogenetic implications. Mol. Phylogen. Evolut., 38: 648-658. https://doi.org/10.1016/j.ympev.2005.10.018

Allcock, A.L., Cooke, I.R. and Strugnell, J.M., 2011. What can the mitochondrial genome reveal about higher-level phylogeny of the molluscan class Cephalopoda? Zool. J. Linnean Soc., 161: 573–586. https://doi.org/10.1111/j.1096-3642.2010.00656.x

Allen, J.J., Bell, G.R., Kuzirian, A.M., Velankar, S.S. and Hanlon, R.T., 2014. Comparative morphology of changeable skin papillae in octopus and cuttlefish. J. Morphol., 275: 371-390. https://doi.org/10.1002/jmor.20221

Anderson, S., Bankier, A.T., Barrell, B.G., De, B.M., Coulson, A.R., Drouin, J., Eperon, I.C., Nierlich, D.P., Roe, B.A. and Sanger, F., 1981. Sequence and organization of the human mitochondrial genome. Nature, 290: 457-465. https://doi.org/10.1038/290457a0

Benson, G., 1999. Tandem repeats finder: a program to analyze DNA sequences. Nucl. Acids Res., 27: 573-580. https://doi.org/10.1093/nar/27.2.573

Boore, J.L., 1999. Animal mitochondrial genomes. Nucl. Acids Res., 27: 1767-1780. https://doi.org/10.1093/nar/27.8.1767

Boyle, P.R. and Boletzky, S.V., 1996. Cephalopod populations, definition and dynamics. Phil. Trans. R. Soc. B, 351: 985-1002. https://doi.org/10.1098/rstb.1996.0089

Carpenter, K.E. and Niem, V.H., 1999. The living marine resources of the western central pacific. batoid fishes, chimaeras and bony fishes (Elopidae to Linophrynidae). Food and Agriculture Organization of the United Nations.

Chen, X., Liu, B. and Wang, R.G., 2009. Cephalopods of the World. China Ocean Press, Beijing.