Brown hagfish from the northwest and east coasts of Honshu, Japan are genetically different

ABSTRACT

The brown hagfish (Eptatretus atami) is one of several known hagfish species occurring in Japanese coastal waters. To date, there has been no research studying genetic polymorphisms in the species. In the present study, we analyzed differences in nucleotide sequences between two populations: one from Suruga Bay on the Pacific coast of Honshu, Japan, and the other from the Sea of Japan, off Akita on the northwest coast of Honshu. We sequenced part of the cytochrome oxidase subunit 1 gene (COX1) from the mitochondrial genome, and three G protein-coupled receptor genes from the nuclear genome. Phylogenetic networks of all four genes showed divergence between the two populations. Further, comparison of the COX1 data using a phylogenetic tree for a range of hagfish species indicated clear differences between the populations, suggesting that they differ at the species level. The numbers of their teeth, in particular of fused cusps (anterior/posterior multicusps), also supported these findings. Individuals of the Suruga Bay population had 3/3 fused cusps, as described for E. atami, whereas individuals of the Akita population had 3/2 fused cusps. These results suggest that the brown hagfish from the Sea of Japan, off the northwest coast of Honshu, is a distinct species from E. atami.

INTRODUCTION

Hagfish are jawless vertebrates which, along with lampreys, belong to the Cyclostomata. Over the past few years, several studies have investigated their phylogeny and taxonomy using nucleotide sequence data (e.g., Kuo et al., 2003, 2010; Chen et al., 2005; Fernholm et al., 2013; Zintzen et al., 2015). However, in this regard, little attention has been given to the species occurring in Japanese coastal waters. The inshore hagfish (Eptatretus burgeri), which is found along the Pacific coast of southern Japan, South Korea, and northern and northwestern Taiwan, is commonly used as a representative of hagfish in molecular studies.

The brown hagfish (E. atami, Dean, 1904) is another species that occurs in Japanese coastal waters, and it can be distinguished from the inshore hagfish based on the spaces between the gill pouches. Both species have six pairs of gill pouches, but whereas those of the inshore hagfish are well spaced and arranged in a linear pattern, those of the brown hagfish are closely spaced and sometimes not linearly arranged (Nakabo and Kai, 2013).

The brown hagfish is one of the most important hagfishes from the perspective of the fisheries industry. The grilled meat of this species is popular in parts of Japan and in Korea, and the skin is used for leather.

Brown hagfish are thought to inhabit the waters on both sides of the Japanese archipelago (the Pacific and the Sea of Japan), and are called “Kuro-nutaunagi” (kuro = black, nutaunagi = hagfish) in Japanese (Nakabo and Kai, 2013). McMillan and Wisner (2004) compared the morphologies of hagfish species of the Northwestern Pacific Ocean and reported that Paramyxine moki and P. walkeri are closely related to the brown hagfish in terms of morphology. They reported that P. moki is only found in Misaki, Sagami Bay, on the Pacific coast of Honshu. In contrast, P. walkeri is found both off Choshi, on the east (Pacific) coast of Honshu, and along the northwest coast, in the Sea of Japan, from Izumozaki to Niigata. Although McMillan and Wisner (2004) used Paramyxine as the genus name, molecular phylogenetic analyses have since indicated that this genus is indistinguishable from Eptatretus (Jansson et al., 1995; Kuo et al., 2003, 2010; Fernholm et al., 2013). We therefore use Eptatretus as the genus name for these species.

To date, no studies have attempted to identify DNA sequence polymorphisms specific to the brown hagfish. We therefore evaluated nucleotide sequence differences between two populations of this species occurring off Honshu: one from Suruga Bay, on the Pacific coast, and another from Akita, in the Sea of Japan, near Akita Prefecture on the northwest coast. As a mtDNA marker, we employed a part of the cytochrome oxidase subunit I gene (COX1) that is widely used for phylogenetic studies of hagfishes (e.g., Rasmussen et al., 1998; Steinke et al., 2009; April et al., 2011; Mabragaña et al., 2011; Fernholm et al., 2013; Zintzen et al., 2015). As nuclear gene markers, we used three G protein-coupled receptor (GPR) genes (GPR4, GPR27 and GPR63). Typical GPRs consist of approximately 400 amino acid residues and are encoded in a single exon. Because allelic length polymorphisms are generally rarer in exons than in introns and intergenic regions, we were able to obtain usable long-range nucleotide sequences for comparisons using a PCR-direct sequencing method. Thus, we selected the GPR genes as nuclear gene markers.

MATERIALS AND METHODS

Samples and DNA extraction

Thirteen brown hagfish individuals were captured at a depth of approximately 300 m in Suruga Bay, on the Pacific coast of Honshu, Japan in June 2011 (Supplementary Fig. S1). A further 16 individuals were captured at a depth of approximately 120 m in the Sea of Japan near Akita Prefecture (off Akita) on the northwest coast of Honshu in September 2013. We also used transcriptome sequence data from liver and gill samples of one individual from Suruga Bay (Suzuki et al., 2017). Genomic DNA was extracted from muscle tissues using the conventional sodium dodecyl sulfate lysis and phenol-chloroform method.

PCRs and sequencing

PCRs were performed in a 10-μl reaction volume containing 20–50 μg of template DNA, 1 × Ex Taq buffer, 0.2 mM dNTP mixture, 1.0 μM each of forward and reverse primers, and 1 U of TaKaRa Ex Taq (TaKaRa Bio). The list of primers used for PCRs and sequencing is given in Supplementary Table S1.

We used the transcriptome sequence data for the brown hagfish (Suzuki et al., 2017) as references for primer design. We checked that the three GPR genes had human orthologs and were located on different chromosomes (GPR4:19q13.3, GPR27:3p21-p14 and GPR63:6q16.1-q16.3) in humans. This implies that these genes were not located close to each other in the brown hagfish genome, as would be the case for tandem duplicated genes, but were instead located in sufficiently different genomic regions.

The PCR conditions consisted of 35 cycles of 10 s denaturation at 98 ℃, followed by 30 s primer annealing at 50–60 ℃ (50 ℃ for COX1, 55 ℃ for GPR4, and 60 ℃ for GPR27 and GPR63) and 1.5–2.5 min (1.5 min for GPR4 and GPR63, and 2.5 min for GPR27 and COX1) extension at 72 ℃. The PCR products were confirmed by 1.5% agarose gel electrophoresis and purified using a MonoFas DNA Purification Kit I (GL Science) or a High Pure PCR Product Purification Kit (Roche Diagnostics).

DNA sequencing was performed on the PCR products using the BigDye Terminator v3.1 Cycle Sequencing Kit and ABI PRISM 310/3130 Genetic Analyzers (Applied Biosystems). The Phred/Phrap software program (Ewing et al., 1998) was used for base-calling and assembly, and to obtain quality scores for the assembled data. Editing was performed using the Consed package (Gordon et al., 1998) to identify all low-quality bases and to check whether the assembly was correct, based on linking information. Both strands of the PCR products were sequenced, and two overlapping peaks for one site in the nuclear genes were used to indicate heterozygosity. The site was designated by following the Nomenclature Committee of the International Union of Biochemistry (NC-IUB) protocol.

Polymorphic and phylogenetic analyses

Haplotypes of GPR4, GPR27 and GPR63 were predicted using the PHASE 2.1.1 program (Stephens et al., 2001; Stephens and Scheet, 2005). Character-based phylogenetic networks for the haplotypes of each gene were constructed manually following the procedures of Bandelt (1994) and Saitou and Yamamoto (1997). DnaSP v5 software (Librado and Rozas, 2009) was used to estimate polymorphic parameters. The phylogenetic tree of COX1 was constructed using the neighbor-joining method (Saitou and Nei, 1987), and the Tamura and Nei (1993) model with MEGA7 (Kumar et al., 2016) software. The reliability of the inferred tree was tested using the bootstrap method with 1,000 replicates (Felsenstein, 1985).

Counts of teeth

The teeth were extracted from each individual and washed with ethanol. Using 10× and 20× loupe magnifiers, we counted the numbers of anterior (outer) fused cusps (multicusps), posterior (inner) multicusps, anterior unicusps (AUC) and posterior unicusps (PUC) for 12 of the Suruga Bay individuals and all 16 of the Akita individuals.

RESULTS

Nucleotide sequences of brown hagfish genes

We obtained nucleotide sequences of 735 bp for GPR4, 793 bp for GPR27, 627 bp for GPR63 and 1551 bp for COX1 in brown hagfish. Nucleotide sequence data were deposited in the DDBJ/EMBL/GenBank International Nucleotide Sequence Database, and given the accession numbers LC178903–LC179015. Predicted haplotypes of the three nuclear genes and haplotypes of COX1 are given in Supplementary Table S2. Levels of polymorphism for each gene in the two populations are summarized in Table 1. For COX1, the average nucleotide diversity (π) of the Suruga Bay population (0.0009) was approximately five times higher than that of the Akita population (0.0002). In contrast, for the three nuclear genes (GPR4, GPR27 and GPR63), there were no such large differences in π values between the two populations.

Table 1. Polymorphisms in brown hagfishes in three nuclear genes and one mitochondrial gene

Gene	Population	N	S	h	π	θW	DT
GPR4	Suruga Bay	24	10	8	0.0032	0.0036	−0.4473
735 bp	Akita	32	7	7	0.0017	0.0024	−0.8206
	Merged	56	18	15	0.0074	0.0053	1.2113
GPR27	Suruga Bay	24	2	3	0.0008	0.0007	0.3448
793 bp	Akita	32	3	5	0.0008	0.0009	−0.4159
	Merged	56	9	8	0.0037	0.0025	1.3311
GPR63	Suruga Bay	22	6	7	0.0023	0.0026	−0.3440
627 bp	Akita	32	7	7	0.0033	0.0028	0.5059
	Merged	54	16	14	0.0063	0.0056	0.3758
COX1	Suruga Bay	13	7	6	0.0009	0.0015	−1.3784
1551 bp	Akita	16	1	2	0.0002	0.0002	1.3090
	Merged	29	68	8	0.0208	0.0112	3.1682***

N: number of sequences; S: number of polymorphic sites; h: number of haplotypes; π: average number of nucleotide differences per site between two sequences; θ_W: nucleotide diversity based on the proportion of segregating sites in a sample; D_T: Tajima’s D (Tajima, 1989). *** indicates the significance level P < 0.001.

Overall, the phylogenetic networks for the three nuclear genes indicated splits between the two populations—four sites in GPR4 (207, 217, 337 and 717), four sites in GPR27 (133, 240, 393 and 744) and three sites in GPR63 (19, 67 and 622)—although there were also some incompatible sites in the phylogenetic network of GPR4 (Fig. 1). The COX1 sequence data, however, unequivocally indicated a divergence between the two populations (Fig. 2). The average numbers of nucleotide differences per site between populations (Dxy) were 0.0127, 0.0066 and 0.0097 for GPR4, GPR27 and GPR63, respectively. In contrast, Dxy for COX1 was 0.0400, approximately four times higher than that of the three nuclear genes combined (0.0096).

Fig. 1.

Phylogenetic networks of nucleotide sequence data for GPR4 (A), GPR27 (B) and GPR63 (C) in brown hagfish. The circle at each node indicates whether that haplotype occurred in the Suruga Bay (black) or Akita (gray) sample, and the size of the circle indicates the frequency of each haplotype, as shown on the lower right. The numbers on each branch are the positions of the nucleotide substitutions responsible for those branches, and those in bold denote nonsynonymous substitutions. Bold lines indicate the split between the two populations.

Fig. 2.

Phylogenetic network of nucleotide sequence data for COX1 in brown hagfish. Details are as for Fig. 1.

Numbers of teeth

All 12 of the Suruga Bay individuals whose teeth we counted had 3/3 fused cusps (three anterior and three posterior multicusps) and all 16 of the Akita individuals had 3/2 fused cusps (Supplementary Fig. S2). The Suruga Bay individuals had 8–9, 5–9 and 42–47 AUC, PUC and total cusps, respectively (one individual had only five PUC on the left), and the Akita individuals had 6–8, 7–9 and 37–42 (Fig. 3).

Fig. 3.

Numbers of anterior unicusps (AUC) (A) and posterior unicusps (PUC) (B) in brown hagfish from the Suruga Bay (black) and Akita (gray) populations.

The phylogenetic tree for COX1

Because COX1 nucleotide sequence data for other hagfish species were available in the DNA database, we constructed a phylogenetic tree for COX1 for the Eptatretus and Myxine species, to investigate the phylogenetic relationships of E. atami with other hagfish species (Fig. 4). The list of species and sequences used for the tree construction is given in Supplementary Table S3. Haplotypes from our two study populations did not form a cluster; instead, the haplotype COX1-s1 from the Suruga Bay population formed a cluster with E. cf. fernholmi, from Taiwan, and COX1-a1 from the Akita population was identical to that of E. sp. ‘Korea’.

Fig. 4.

Phylogenetic tree of nucleotide sequence data for COX1 in hagfish species, constructed using the neighbor-joining method (Saitou and Nei, 1987) and the Tamura and Nei (1993) model. Nucleotide sequence data are listed in Supplementary Table S3. A sequence from Rubicundus lopheliae was used as the outgroup. Scale bar: 0.02 nucleotide substitutions per site. Bootstrap values (percentages of 1000 replicates) of more than 90% are shown at the nodes. Brown hagfish sequences identified in this study are enclosed in rectangles. Haplotypes COX1-s1 and COX1-a1 in Fig. 2 were used as E. atami (Suruga Bay) and E. atami (Akita), respectively.

The average nucleotide difference per site between Eptatretus and Myxine species was estimated as 0.233, based on the phylogenetic tree for COX1 (Fig. 4). Kuraku and Kuratani (2006) estimated the divergence between these two genera to have occurred 60–90 million years ago (MYA). Thus, we estimated the evolutionary rate to be 1.29 × 10⁻⁹–1.94 × 10⁻⁹, using the formula d = 2λt, where d is nucleotide difference, λ is evolutionary rate and t is divergence time. Since the nucleotide difference per site between E. burgeri and our Akita sample of E. atami was estimated (from the phylogenetic tree) to be 0.084, based on our estimated evolutionary rate we estimated the divergence time between these two species to be 22–32 MYA. This corresponds to the period from the Aquitanian of the Miocene to the Rupelian of the Oligocene.

DISCUSSION

Interpretation of haplotype networks

In general terms, reticulations within a single phylogenetic network indicate incompatible phylogenetic information (or character states). If there is no incompatible phylogenetic information in the data, the phylogenetic network automatically becomes an unrooted phylogenetic tree. A phylogenetic network is thus a generalization of the concept of a phylogenetic tree (Kitano, 2012).

The phylogenetic networks for all three nuclear genes contained reticulations (Fig. 1), indicating that these sequence data contained incompatible phylogenetic information. For example, in Fig. 1B, the 306^th site divides GPR27-a2–GPR27-a4–GPR27-a5 from the rest of the haplotypes, and the 456^th site divides GPR27-a2–GPR27-a3 from the rest. This incompatibility is caused by parallel substitutions at the 456^th site on the branch from the GPR27-a1 node to the GPR27-a3 node and the branch from the GPR27-a4 node to the GPR27-a2 node. It may also be caused by parallel substitutions at the 306^th site on the branch from the GPR27-a3 node to the GPR27-a2 node and the branch from the GPR27-a1 node to the GPR27-a4 node. In Fig. 1A, four sites (207, 217, 337 and 717) divide haplotypes in the Suruga Bay population from those in the Akita population. In contrast, at the 273^rd site, all Suruga Bay haplotypes had guanine (G) except GPR4-s7 (adenine, A), while all Akita haplotypes had A except GPR4-a5 and GPR4-a7 (G). We can explain this pattern by postulating that after the divergence between the Suruga Bay (G) and Akita (A) haplotypes, backward or parallel substitutions occurred in GPR4-s7 (G to A) and in GPR4-a5 and GPR4-a7 (A to G). The 234^th and 540^th sites contained similar parallel and/or backward substitutions. At the 234^th site, all Suruga Bay haplotypes had thymine (T) except GPR4-s1, GPR4-s3 and GPR4-s4 (cytosine, C), while all Akita haplotypes had C except GPR4-a6 (T). At the 540^th site, all Suruga Bay haplotypes had A except GPR4-s4 (G), while all Akita haplotypes had G except GPR4-a3 (A).

Regardless of the exact mechanisms involved, all four phylogenetic networks (Figs. 1–2) indicated divergence between the two populations. Note that in the GPR63 phylogenetic network (Fig. 1C), haplotypes from the Akita population can be divided into two groups: ((GPR63-a1–a5) and (GPR63-a6 and GPR63-a7)). This suggests that there are subpopulation structures or some natural selection pressure acting on GPR63 in the Akita population, which calls for further investigation.

Furthermore, our results revealed inconsistency in the level of nucleotide diversity between the mtDNA and the three nuclear gene markers. While the nucleotide diversity of the mtDNA from the Suruga Bay population was approximately five times that of the Akita population (for which we obtained only two haplotypes; Fig. 2, Table 1), the differences in the diversity of the nuclear genes were much smaller. It is not clear why we obtained only two haplotypes of COX1 for the Akita population, but further investigation, using larger sample sizes, may answer this question.

Morphological differences between the two populations

McMillan and Wisner (2004) reported that the new species E. moki and E. walkeri could be distinguished from E. atami based on their different numbers of teeth. Eptatretus moki has 3/2 fused cusps, 6–8 AUC, 7–9 PUC and 38–42 cusps altogether (Table 2). Eptatretus walkeri has 3/2 fused cusps, 6–8 AUC, 6–9 PUC and 38–44 cusps altogether. In contrast, E. atami has 3/3 fused cusps, 9–10 AUC, 8–10 PUC and 47–52 cusps altogether (McMillan and Wisner, 2004). Although the number of fused cusps in our Suruga Bay sample (3/3) matched that previously reported for E. atami, that in our Akita sample (3/2) did not. In this respect, the Akita population was similar to E. walkeri and E. moki. Similarly to the genetic data, therefore, these morphological data suggest that the brown hagfish in the Sea of Japan, off the northwest coast of Honshu, is a distinct species from E. atami.

Table 2. Numbers of cusps in several hagfish species

	McMillan and Wisner (2004)				This study		Iwata (1993)		Nakabo and Kai (2013)
	E. burgeri	E. atami	E. moki	E. walkeri	Suruga Bay	Akita	E. burgeri	E. atami	E. burgeri	E. atami
Gill pouches	6	6	6	6	6	6	6a	6	6a	6
Fused cusps	3/2	3/3	3/2	3/2	3/3	3/2	3/2	3/2–3	3/2	3/2–3
AUC	6–8	9–10	6–8	6–8	8–9	6–8	8–9	10–13	8–11	8–13
PUC	7–9	8–10	7–9	6–9	5–9	7–9	9	10–13	9–11	9–13
Total cusps	35–42	47–52	38–42	38–44	42–47	37–42	–	–	–	–

^a  Occasionally 7.

Divergence of the two populations

The COX1 phylogenetic tree (Fig. 4) clearly indicates that the genetic differences between the two populations are equivalent to species differences. The significant positive value of Tajima’s (1989) D for the overall COX1 data, and the positive (although not statistically significant) Tajima’s D for the overall data for GPR4, GPR27 and GPR63 (Table 1), also support the divergence of the two populations.

Jansson et al. (1995) identified two sympatric Paramyxine (Eptatretus) species from Misaki by studying protein electrophoresis. They provisionally named the species Paramyxine (Eptatretus) sp. A and Paramyxine (Eptatretus) sp. B, which distinguished them from Paramyxine (Eptatretus) atami. It is not clear whether these two species are E. walkeri and E. moki, because the authors collected only tissue samples as experimental materials, and did not conduct morphological examinations (Jansson et al., 1995). Nevertheless, their discovery indicates that the brown hagfish complex is not monospecific, which is consistent with our results.

The phylogenetic tree (Fig. 4) indicates that the Akita population is identical to E. sp. ‘Korea’. The origin of E. sp. ‘Korea’ is described as only “South Korea” (Fernholm et al., 2013); since South Korea and the northwest coast of Honshu are geographically close, it is feasible that our Akita sample and E. sp. ‘Korea’ came from the same population.

Since our Akita sample was collected from within the distribution of E. walkeri, but not from within the distribution of E. moki (McMillan and Wisner, 2004), it may be appropriate to treat the Akita population as E. walkeri. However, the holotype of E. walkeri was recorded off the east coast of Honshu, and it is possible that there is genetic differentiation between the populations from the east and west coasts. Alternatively, it is possible that E. walkeri also inhabits the region between the two recorded populations, along the north coast, so that there is actually only one contiguous population. Either way, it is too early to designate the Akita population as E. walkeri, and we have tentatively named it E. sp. Akita. Based on our results, we propose “northern brown hagfish” as the common name, and “Kita kuro-nutaunagi” (kita = northern, kuro = black, nutaunagi = hagfish) as the standard Japanese name. Further investigation, involving samples from off Choshi to the north coast of Honshu, may resolve the question of whether the northern brown hagfish (E. sp. Akita) is, in fact, E. walkeri.

Recently, concern has emerged regarding the overexploitation of hagfish. Sound taxonomic knowledge and detailed examination of the diversity of hagfish species at the genetic level is thus needed for the fisheries industry to be able to plan sustainable fishing.

ACKNOWLEDGMENTS

We are grateful to Mr. Hisashi Hasegawa of the fishing boat Tyoukane-maru for capturing the brown hagfish samples in Suruga Bay, and to Mr. Choyu Sawaki of the fishing boat Mishima-maru for capturing the northern brown hagfish samples off Akita.

REFERENCES

•  April,  J.,  Mayden,  R. L.,  Hanner,  R. H., and  Bernatchez,  L. (2011) Genetic calibration of species diversity among North America’s freshwater fishes. Proc. Natl. Acad. Sci. USA 108, 10602–10607.
•  Bandelt,  H. J. (1994) Phylogenetic networks. Verh. Naturwiss. Ver. Hamburg 34, 51–71.
•  Chen,  Y. W.,  Chang,  H. W., and  Mok,  H. K. (2005) Phylogenetic position of Eptatretus chinesis (Myxinidae: Myxiniformes) inferred by 16S rRNA gene sequence and morphology. Zool. Stud. 44, 111–118.
•  Dean,  B. (1904) Notes on Japanese myxinoids. A new Paramyxine and a new species Homea okinoseana. Reference also to their eggs. J. Coll. Sci. Imp. Univ. Tokyo 19, 1–23.
•  Ewing,  B.,  Hillier,  L.,  Wendl,  M. C., and  Green,  P. (1998) Base-calling of automated sequencer traces using phred. I. Accuracy assessment. Genome Res. 8, 175–185.
•  Felsenstein,  J. (1985) Confidence limits on phylogenies: An approach using the bootstrap. Evolution 39, 783–791.
•  Fernholm,  B.,  Norén,  M.,  Kullander,  S. O.,  Quattrini,  A. M.,  Zintzen,  V.,  Roberts,  C. D.,  Mok,  H. K., and  Kuo,  C. H. (2013) Hagfish phylogeny and taxonomy, with description of the new genus Rubicundus (Craniata, Myxinidae). J. Zool. Syst. Evol. Res. 51, 296–307.
•  Gordon,  D.,  Abajian,  C., and  Green,  P. (1998) Consed: a graphical tool for sequence finishing. Genome Res. 8, 195–202.
•  Iwata,  A. (1993) Myxinidae. In: Fishes of Japan with pictorial keys of the species, first edition. (ed.: Nakabo, T.), pp. 83–84. Tokai University Press, Kanagawa. (in Japanese)
•  Jansson,  H.,  Wyöni,  P. I.,  Fernholm,  B.,  Bredwad,  M.,  Mierzykowaka,  A., and  Tegelström,  H. (1995) Genetic relationships among species of hagfish revealed by protein electrophoresis. J. Fish Biol. 47, 599–608.
•  Kitano,  T. (2012) Application of phylogenetic network. In: Post-Genome Biology of Primates. (eds.: Hirai, H., Imai, H., and Go, Y.), pp. 181–190. Springer, Tokyo.
•  Kumar,  S.,  Stecher,  G., and  Tamura,  K. (2016) MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874.
•  Kuo,  C. H.,  Huang,  H., and  Lee,  S. C. (2003) Phylogeny of hagfish based on the mitochondrial 16S rRNA gene. Mol. Phylogenet. Evol. 28, 448–457.
•  Kuo,  C. H.,  Lee,  S. C., and  Mok,  H. K. (2010) A new species of hagfish Eptatretus rubicundus (Myxinidae: Myxiniformes) from Taiwan, with reference to its phylogenetic position based on its mitochondrial DNA sequence. Zool. Stud. 49, 855–864.
•  Kuraku,  S., and  Kuratani,  S. (2006) Time scale for cyclostome evolution inferred with a phylogenetic diagnosis of hagfish and lamprey cDNA sequences. Zoolog. Sci. 23, 1053–1064.
•  Librado,  P., and  Rozas,  J. (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452.
•  Mabragaña,  E.,  Diaz de Astarloa,  J. M.,  Hanner,  R.,  Zhang,  J., and  González Castro,  M. (2011) DNA barcoding identifies argentine fishes from marine and brackish waters. PLoS One 6, e28655.
•  McMillan,  C. B., and  Wisner,  R. L. (2004) Review of the hagfishes (Myxinidae, Myxiniformes) of the northwestern Pacific Ocean, with descriptions of three species Eptatretus fernholmi, Paramyxine moki and P. walkeri. Zool. Stud. 43, 51–73.
•  Nakabo,  T., and  Kai,  Y. (2013) Myxini, In: Fishes of Japan with pictorial keys to the species, third edition. (ed.: Nakabo, T.), pp. 141–143. Tokai University Press, Kanagawa. (in Japanese)
•  Rasmussen,  A. S.,  Janke,  A., and  Arnason,  U. (1998) The mitochondrial DNA molecule of the hagfish (myxine glutinosa) and vertebrate phylogeny. J. Mol. Evol. 46, 382–388.
•  Saitou,  N., and  Nei,  M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, 406–425.
•  Saitou,  N., and  Yamamoto,  F. (1997) Evolution of primate ABO blood group genes and their homologous genes. Mol. Biol. Evol. 14, 399–411.
•  Steinke,  D.,  Zemlak,  T. S.,  Boutillier,  J. A., and  Hebert,  P. D. N. (2009) DNA barcoding fishes of the Canadian Pacific. Mar. Biol. 156, 2641–2647.
•  Stephens,  M., and  Scheet,  P. (2005) Accounting for decay of linkage disequilibrium in haplotype inference and missing data imputation. Am. J. Hum. Genet. 76, 449–462.
•  Stephens,  M.,  Smith,  N. J., and  Donnelly,  P. (2001) A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 68, 978–989.
•  Suzuki,  A.,  Komata,  H.,  Iwashita,  S.,  Seto,  S.,  Ikeya,  H.,  Tabata,  M., and  Kitano,  T. (2017) Evolution of the RH gene family in vertebrates revealed by brown hagfish (Eptatretus atami) genome sequences. Mol. Phylogenet. Evol. 107, 1–9.
•  Tajima,  F. (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.
•  Tamura,  K., and  Nei,  M. (1993) Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol. 10, 512–526.
•  Zintzen,  V.,  Roberts,  C. D.,  Shepherd,  L.,  Stewart,  A. L.,  Struthers,  C. D.,  Anderson,  M. J.,  McVeagh,  M.,  Noren,  M., and  Fernholm,  B. (2015) Review and phylogeny of the New Zealand hagfishes (Myxiniformes: Myxinidae), with a description of three new species. Zool. J. Linn. Soc. 174, 363–393.