Evolutionary Relations of Hexanchiformes Deep-Sea Sharks Elucidated by Whole Mitochondrial Genome Sequences

Hexanchiformes is regarded as a monophyletic taxon, but the morphological and genetic relationships between the five extant species within the order are still uncertain. In this study, we determined the whole mitochondrial DNA (mtDNA) sequences of seven sharks including representatives of the five Hexanchiformes, one squaliform, and one carcharhiniform and inferred the phylogenetic relationships among those species and 12 other Chondrichthyes (cartilaginous fishes) species for which the complete mitogenome is available. The monophyly of Hexanchiformes and its close relation with all other Squaliformes sharks were strongly supported by likelihood and Bayesian phylogenetic analysis of 13,749 aligned nucleotides of 13 protein coding genes and two rRNA genes that were derived from the whole mDNA sequences of the 19 species. The phylogeny suggested that Hexanchiformes is in the superorder Squalomorphi, Chlamydoselachus anguineus (frilled shark) is the sister species to all other Hexanchiformes, and the relations within Hexanchiformes are well resolved as Chlamydoselachus, (Notorynchus, (Heptranchias, (Hexanchus griseus, H. nakamurai))). Based on our phylogeny, we discussed evolutionary scenarios of the jaw suspension mechanism and gill slit numbers that are significant features in the sharks.


Introduction
The subdivision Selachii or modern sharks, along with skates and rays, comprises the subclass Neoselachii within the class Chondrichthyes or cartilaginous fishes. Chondrichthyans, including the Neoselachii, chimaeroids, and several fossil forms are defined as jawed fish with skeletons made of prismatic cartilage rather than bone and pelvic claspers in males.
Among the sharks, Hexanchiformes is regarded as an ancient order of sharks with just five extant species that are characterized by having only one dorsal fin, either six or seven gill clefts, and no nictitating membrane in the eyes [2]. The Hexanchiformes is usually divided into two families, the Chlamydoselachidae (Chlamydoselachus anguineus (C. anguineus)) and the Hexanchidae (Hexanchus griseus (H. griseus), Hexanchus nakamurai (H. nakamurai), Notorynchus cepedianus (N. cepedianus), and Heptranchias perlo (H. perlo)) with the latter family also known as "cow sharks. " The frilled shark, C. anguineus, is very different from the cow sharks, and its own order of Chlamydoselachiformes was proposed [3]. However, derived features (e.g., the extra gill arch and moreheart valve rows) shared with other Hexanchiformes support its retention within the Hexanchiformes. A third family, Notorynchidae, was also proposed for the Notorynchus species because of morphological and behavioral differences from the other members of the family Hexanchidae [3]. Interestingly, the tooth structure and composition of one of the Hexanchiformes, C. anguineus, is similar to that of the stem-group fossil shark Cladoselache sp., although such features are not observed in the other Hexanchiformes species [4]. Therefore, the accurate placement of Hexanchiformes is essential to understand the evolution of morphology in sharks. However, the lack of available DNA sequence data for most shark species and orders remains a major limitation to obtaining reliable results in molecular phylogenetic studies.
The mitochondrial DNA (mtDNA) has been one of the most widely used molecular markers for diversity and phylogenetic studies in animals because of its size, maternal mode of inheritance, high rate of mutation, and simple genomic structure [5]. Although mtDNA sequences have proved valuable in determining phylogenetic relationships, the choice of a gene as a molecular marker and clock in phylogeny is also important [6,7]. Recent phylogenetic studies in different taxa suggest that full-length mitochondrial genomic sequences provide an improved resolution for reconstructing a robust phylogeny and for molecular dating of divergence events within a phylogeny [8]. So far, of all the shark species, the complete mtDNA sequences were determined in the subclass Neoselachii of some species, and the sequences were used for elucidating interrelationships between sharks and bony fishes and between sharks and rays [9][10][11][12][13][14]. Therefore, there were no analyses of the relationships between species within cartilaginous fish orders using whole mtDNA sequences, although intrarelationships were estimated by using partial mtDNA and nuclear DNA sequences [6,7,10,[15][16][17][18][19][20].
In order to obtain more sequence data for different shark species and to allow accurate placement within the order Hexanchiformes, we chose to determine the complete mtDNA sequences of seven shark species including five species of the order Hexanchiformes, along with Somniosus pacificus (S. pacificus), which is a member of the order Squaliformes and Pseudotriakis microdon (P. microdon), which is a member of the order Carcharhiniformes. We then analysed the phylogenetic relationships among the five Hexanchiformes species and between the Hexanchiformes and ten other shark species using the complete mitochondrial genomic sequences of 19 cartilaginous fish species. On the basis of the results of our phylogenetic analysis, we propose new scenarios for the evolution of the jaw suspension mechanism and the number of gill clefts that are significant features in sharks.   Figure 1). In this paper, we use the species name Hexanchus nakamurai [21] instead of its synonym Hexanchus vitulus [22]. The mtDNAs were isolated from the muscle or spleen tissue using the mtDNA Extractor CT Kit (Wako Pure Chemical Industries, Ltd., Osaka, Japan) or standard phenol-chloroform method [23]. The quantity and quality of the isolated DNA samples were measured and estimated, respectively, by the spectral absorbance of the DNA at 260 nm and 280 nm. For PCR amplification of the cytochrome b (CYTB) region, the 10 L amplification reaction contained 10 ng of mtDNA, 1.0 unit of Ex Taq polymerase (TaKaRa Shuzo, Otsu, Japan), 1x PCR buffer, 2.5 mM MgCl 2 , 2 mM of each dNTP, and 0.5 M of each primer. The cycling parameters were as follows: an initial denaturation of 96 ∘ C for 3 min followed by 30 cycles of 96 ∘ C for 30 sec, 50 ∘ C for 30 sec, and 72 ∘ C for 1 min and followed by a final cycle of 72 ∘ C for 4 min. For long-ranged PCR amplifications, the 20 L amplification reaction contained 10 ng of mtDNA, 0.4 unit of KOD-FX polymerase (TOYOBO, Osaka, Japan), 2x PCR buffer, 2.5 mM MgCl 2 , 400 M of each dNTP, and 0.5 M of each primer. The cycling parameters were as follows: an initial denaturation of 94 ∘ C for 2 min followed by 35 cycles of 98 ∘ C for 10 sec, 60 or 68 ∘ C for 30 sec, and 68 ∘ C for 2 to 20 min. PCR reactions were performed by using the thermal cycler GeneAmp PCR system 9700 (Applied Biosystems, CA, USA). The longranged PCR size was 6.5 kb on average and ranged from 1,711 bp to 16,800 bp (Supplementary table 1).

Genomic Sequencing Strategy and Sequence Analysis.
The PCR products were subjected to complete and bidirectional shotgun sequencing with an average 7.2x coverage which was sufficient for assembly and analysis of the entire sequence using previously established procedures [26] and direct sequencing using PCR primers as sequencing primers. DNA sequencing was performed by the cycle sequencing method using AmpliTaq-DNA polymerase FS and the fluorescently labeled BigDye terminators in a GeneAmp PCR system 9700 (Applied Biosystems, Foster City, CA, USA). A 3130xl Genetic Analyzer was used for automated fluorescent sequencing (Applied Biosystems, Foster City, CA, USA). Individual sequences were minimally edited to remove vector sequences and assembled into contigs using the Sequencher 4.2 software (Gene Codes CO., MI, USA). Remaining gaps or ambiguous nucleotides were determined by the direct sequencing of PCR products obtained with appropriate PCR primers or by nucleotide sequence determination of shotgun clones.
Nucleotide similarities between sequences were calculated by the "Search Homology" tool of GENETYX-MAC ver. 12.0 (Software Development Co. Ltd., Tokyo, Japan), and those with nucleotide sequences in GenBank/EMBL/ DDBJ were searched by BLAST program (http://www.ncbi .nlm.nih.gov/BLAST/). The newly determined mtDNA sequences were annotated by comparison with known mtDNA sequence information of other shark species.
Phylogenetic trees were constructed using the modelbased maximum likelihood (ML) and Bayesian inference (BI: MrBayes Ver. 3.2.1 [28]) methods and the distance-based neighbour-joining (NJ) method (MEGA5). For the ML analyses we used "Find best DNA/protein models (ML)" program of MEGA5 to estimate the most likely model of sequence evolution including third codon sites. Based on maximum likelihood values and the Akaike information criterion (AIC) [29], the TN93+G+I model was selected as the most likely model (ln = −118970.445, AIC = 238464.763 for Align Set 1, and ln = −110086.677, AIC = 220253.367 for Align Set 2) for a nucleotide based ML tree using 10,000 ML-bootstrap replicates. The neighbor-joining tree was constructed using distances corrected according to the Kimura 2-parameter model with 1.0 gamma parameters [30] and assessed using 10,000 bootstrap replicates. For the BI analyses, we used MrAIC Ver.

Biological and Genetic
Information. The biological and genetic features of the seven shark species sampled for this study are shown in Table 1. Six or seven gill clefts and one dorsal fin were observed in the Hexanchiformes species, whereas five gill clefts and two dorsal fins were observed in S. pacificus and P. microdon and in all non-Hexanchiformes and Neoselachii, except some skates and myliobatoid rays in which dorsal fins were presumably lost secondarily [32]. Our preliminary sequence analysis of the mtDNA Cytb or 12S genes for six of the seven species, excluding H. nakamurai, using a BLAST search (GenBank) primarily showed the highest nucleotide identities (97.0% to 100%) with previously published nucleotide sequences of the target species and secondarily showed lower nucleotide identities (80.2% to 93.1%) with the nucleotide sequences of different species (Table 1). Therefore, this analysis, on the basis of biological and/or genetic features, helped us to cross-check and confirm that all our collected samples were identified correctly as the targeted species.

Phylogeny of the Hexanchiformes Using Complete mtDNA
Sequences. Figure 2 shows an ML phylogenetic tree constructed by using the TN93+G+I model of the ML method. The tree was constructed using the 13,749 bp nucleotide   alignment (Align Set 1) of the 13 protein coding genes and two rRNA genes from 19 species with C. monstrosa selected as the outgroup. In the aligned sequence 48.7% nucleotides (6,693 sites) were constant sites. The similar topology, supported clades, and bootstrap support values or posterior probability were shown by the GTR+G+I model of the BI method and the Kimura 2-parameter model of the NJ method using the same nucleotide alignment as for the ML method (Supplementary figure 1). Moreover, even if we set three Batoidea species as outgroups, the same topology, supported clades, and bootstrap support values or posterior probability were shown by the TN93+G+I model within the ML method, the GTR+G+I model within the BI method and the Kimura 2-parameter model within the NJ method using 13,784 nucleotide alignment (Align Set 2) of the 13 protein coding genes and two rRNA genes from 18 species (Supplementary figure 2). Our phylogenetic analyses of nucleotide sequences using the ML, BI, and NJ methods suggest that the Selachii is divided into the two superorders, Squalomorphi and Galeomorphi. The Squalomorphi superorder includes the orders Hexanchiformes and Squaliformes, and the Galeomorphi superorder includes the orders Orectolobiformes, Lamniformes, and Carcharhiniformes ( Figure 2). In Figure 2 tree, the Hexanchiformes clade is monophyletic and separates from the Squaliformes clade in the Squalomorphi order that is separated from the Galeomorphi lineage. The Hexanchiformes and Squaliformes lineages showed clades with extremely high bootstrap support values (100) and posterior probabilities (100) (Figure 2, Supplementary figures 1 and 2). The mitogenomic data also show good resolution within the Galeomorphi with Orectolobiformes separated from "Lamniformes plus Carcharhiniformes" with high bootstrap support values (61∼93) and posterior probabilities (100) (Figure 2, Supplementary figures 1 and 2).
Within the Hexanchiformes, C. anguineus is sister of all the others with the longest branch, so, assuming clock like evolution, it could be the oldest extant shark lineage. The N. cepedianus lineage emerged next. By comparison, the terminal branch lengths of H. griseus, H. nakamurai, and H. perlo are relatively short.

Phylogeny of Cartilaginous Fishes.
In the case of Batoidea, Carcharhiniformes, Hexanchiformes, and Squaliformes, our phylogeny supports, in a number of respects, the previous findings derived from partial mitochondrial genomes, genes COX1, NADH2, NADH4, CYTB, 12S, 16S, and/or tRNAs and nuclear genome genes 5.8S, 18S, 28S, and RAG1 [6,7,[15][16][17][18][19]35]. Recently, Naylor et al. [7] published the most comprehensive phylogeny of sharks using a total of 595 shark species representing eight orders and 159 genera and 56 families, but mainly using a single mitochondrial gene 8 BioMed Research International (NADH2) as the molecular marker. Although we did not have any examples of the Echinorhiniformes, Pristiophoriformes, and Squatiniformes in our study, the remaining orders within Squalomorphi were generally similar to the relationships reported by Vélez-Zuazo [6] and by Naylor et al. [7] of Hexanchiformes and Squaliformes within Squalomorphi and Lamniformes, Orectolobiformes, and Carcharhiniformes within Galeomorphi. Recent myological studies on Hexanchiformes also support the inclusion of the Chlamydoselachidae and Hexanchidae in the Squalomorphi [36]. Our analysis, though limited to just one lamniform representative, found strong support for Lamniformes as the sister order of Carcharhiniformes [2,19] instead of the Orectolobiformes [6,35]. In addition, our phylogeny is inconsistent with some previously reported morphology-based phylogenies such as the hypnosqualean hypothesis that places the batoids within sharks [2,37], the tooth structure-based tree that places Hexanchiformes outside of batoids [38,39], the jaw protrusion, and feeding-based trees [40] and a Bayesian analysis based on the CYTB gene that supported the position of Hexanchiformes as the sister of all other shark orders [20]. The full-length mitochondrial genomic sequences provide strong statistical support and an improved resolution for reconstructing a robust phylogeny in the cartilaginous fish [8]. Therefore, our phylogeny appears to be a reasonable reconstruction of the evolutionary process dividing the Selachii into the two superorders, Squalomorphi and Galeomorphi. However, the divergence time for a node of Squaliformes and Hexanchiformes was estimated to be 115.4 Mya (109.9 Mya∼123.7 Mya). This is not largely consistent with the known fossil record of Hexanchiformes (190 Mya) [33]. In this regard, in case of setting the divergence time of Hexanchiformes for 190 Mya, the divergence times for each node of Selachii/Batoidea and Chlamydoselachidae were estimated as 351.3 Mya and 134.9 Mya, respectively. The divergence time of Chlamydoselachidae is consistent with the known fossil record of Chlamydoselachidae (85 Mya) [33] but that of Selachii/Batoidea is not consistent with the previously estimated divergence time of Selachii and Batoidea (213.4 Mya (203.3∼228.8 Mya in the 95% confidence intervals)) [24]. In this study we used only 15 Selachii and three Batoidea species for the phylogenetic analysis, but some Selachii species that show ambiguous classification such as Echinorhiniformes, Pristiophoriformes, and Squatiniformes were not included in this study. Therefore, detailed phylogenetic analysis based on full-length mitochondrial genomic sequences using additional Selachii species that are thought to be diverged on evolutionary important positions, will be necessary for estimation of the precise divergence times in future.

Phylogenetic Relationships in Hexanchiformes.
Vélez-Zuazo and Agnarsson [6] reported that the Hexanchidae within Hexanchiformes was paraphyletic because it also contained the only species of Chlamydoselachidae. Their Bayesian analysis of Hexanchiformes, in using only a small portion (15%) of a single nucleotide sequence composed of five genes, COX1, NADH2, CYTB, 16S, and Rag1, may be compromised by too little data. They showed low support with a posterior probability of 45% at the node splitting the C. anguineus and N. maculates taxa. In addition, although they found that N. maculates was a sister group to the N. cepedianus, N. maculates is in fact a synonym of N. cepedianus and therefore the same species [6]. In our analysis, N. cepedianus is clearly the sister of the Heptranchias-Hexanchus clade with which it forms a clade separate from Chlamydoselachus with all nodes strongly supported (Figure 2). The sequencing of the complete mitochondrial genome of N. maculates, if it is a separate species from N. cepedianus, should help to resolve this issue.
The frilled shark, C. anguineus, in our phylogenetic analysis, was found to be the sister group to all the other Hexanchiformes, similar to the results on their anatomical studies [2,3,36]. Some systematists have proposed a separate order for the frilled shark [3]. Recently, a new Chlamydoselachidae species (Chlamydoselachus africana; C. africana) was discovered around southern Africa [41]. Although it is not known if the mtDNA of C. africana will group phylogenetically with C. anguineus, our phylogeny supports the retention of C. anguineus within the Hexanchiformes rather than within its own order. The frilled shark, C. anguineus, is placed in its own family, Chlamydoselachidae, because of its many unusual features such as elongated and eel-like body; its low dorsal fin; blunt snout, long jaws that are narrower at the tip than at the corners; terminal mouth, similar upper and lower teeth with three prong-like cusps; and gill clefts with frilly margins and the first gill slit continuous across the throat [42].
The position of N. cepedianus in our phylogenetic tree favours its own family name of Notorynchidae rather than being placed into the Hexanchidae family with the one Heptranchias species and the two Hexanchus species. Incidentally, H. nakamurai [21] in our study corresponds to H. viulus [22] in the Vélez-Zuazo and Agnarsson study [6], as the latter species name is junior to the former name [42,43]. However, Naylor et al. [44] suggest that both H. nakamurai and H. viulus may actually represent valid species.

New Insights on the Morphological Evolution of Sharks.
Our phylogenetic analysis based on mitochondrial genomic sequences is useful for comparing morphological features to phylogenetic relationships among the sharks. For example, the Hexanchiformes species have six or seven gill clefts and one dorsal fin, whereas most of other Selachii species have five or six gill clefts and mostly two dorsal fins (Table 1). Most fossils of Agnathans and some fossils of acanthodian (stem Chondrichthyes or stem osteichthyes) support the presence of multiple gill clefts [45], and the sea lamprey, which is in the class Petromyzontida, has seven gill clefts. However, the rabbitfish, which is in the subclass Holocephali, is from the sister clade of sharks and Batoidea, and it has one gill cleft. In this regard, the phylogenetic tree suggests that the multiple gill clefts have been maintained by species-specific increases and decreases in both the Petromyzontida and Neoselachii lineages (Figure 3(a)), but holocephalians could easily have reduced gill slits to one.
The phylogenetic placement of Hexanchiformes in our study suggests a clearer scenario for the evolution of the upper jaw suspension in sharks. Upper jaws (palatoquadrate) of sharks are "suspended" from the cranium by the hyomandibular cartilage and several articulations between the cranium and upper jaw (Figure 3(b)) and shark jaws evolved with a general trend towards jaw shortening, increased kinesis of upper jaw suspension, and protrusion [46]. Many researchers have focused on the evolution of the jaw suspension mechanism in order to clarify the phylogenetic relationship within sharks [47][48][49]. Because the jaw suspension mechanism is strongly coupled with the jaw protrusion capability, the evolution of jaw protrusion behaviour has been reconstructed based on the evolution of jaw suspension mechanisms [46]. However, the position of Hexanchiformes as the sister to Squaliformes (Figure 3) suggests that the Hexanchiformes represents an evolved state of the jaw suspension mechanism independent of the other shark orders. It can be expected that investigation of fossils from non-Neoselachii Chondrichthyes such as the hybodonts and cladodonts will help to elucidate the polarization of the evolution of the jaw suspension mechanism.
The Hexanchiformes species have two characteristic articulations between their cranium and palatoquadrate. One is the orbital articulation, which is between the orbital process of the palatoquadrate and the orbital walls of the cranium, and the other is the postorbital articulation, which is between the otic process of the palatoquadrate and the postorbital process of the cranium [47,48]. Due to the postorbital articulation, jaw protrusion capabilities of the Hexanchiformes species are strongly restricted. On the other hand, Squaliformes species lack the postorbital articulations, and Galeomorphi species lack both the orbital articulations and the postorbital articulation [47,48]. Because the postorbital and/or orbital articulations are absent from these two clades, in contrast to Hexanchiformes species, they have a greater capability of jaw protrusion, flexibility, and maneuverability [46]. Extensive upper jaw protrusion in modern sharks was found to allow faster closure of jaws to gouge or cut smaller pieces of prey to fit into the mouth [50]. Therefore, based on the differences in jaw suspension mechanisms among Hexanchiformes, Squaliformes, and Galeomorphi, the following evolutionary scenario can be proposed. First, the orbital process and postorbital articulation were missing from the common ancestor of all shark orders. Second, the orbital articulation was gained at the base of Squaliformes + Hexanchiformes clade, and this modification restricted the evolved sharks to protrude their jaws. Third, the postorbital articulation was gained at the base of Hexanchiformes clade, and this modification might have further restricted the Hexanchiformes from the capability of upper jaw protrusion.
Previous research suggested that the presence of the orbital process in the palatoquadrate is one of the strong morphological characteristics indicating the monophyly of Hexanchiformes and Squaliformes. Thus, Maisey named the Hexanchiformes + Squaliformes clade as "orbitostylic sharks" [25,[47][48][49]. Our result suggests that the presence of the orbital process is an evolved character of the Hexanchiformes + Squaliformes clades and not a plesiomorphic character of the whole shark clade (Figure 3(a)).
In summary, we sequenced and analysed the complete mtDNA sequences of five Hexanchiformes species. Our phylogeny and the known morphological features of sharks resolved interrelationships of major Hexanchiformes species. Further insights into phylogeny of the mtDNA sequences were provided by comparative analyses using other shark and nonshark species. A similar approach using the whole mtDNA genome of sharks in the other orders should help to resolve the intraspecific and interspecific relationships within Chondrichthyes (cartilaginous fishes). The Hexanchiformes are observed mainly in deep-sea areas (<1000 m) all over the world. However, they rarely occur in the continental shelf shallower than 200 m [51]. This may be a relic of a past behavioral habit, when their ancestors had inhabited the coastal shelf [52,53]. Although some Hexanchiformes specific morphological features such as tooth, gill cleft numbers, jaw suspension, and no nictitating membrane in the eyes. have been reported so far, the other features such as distribution of the living areas, physiology, reproduction and genetic diversity are unknown. Therefore, it will be necessary to compare between those features and phylogenetic relationships derived from nucleotide sequences to comprehensively understand evolution of the sharks including Hexanchiformes.