Population Structure of Upstream Fat Minnow (Rhynchocypris oxycephalus) Based on Mitochondrial DNA Cytochrome b Analysis in Japan

Abstract

Population genetic analysis using mtDNA sequence analysis was carried out to clarify the population structure of Rhynchocypris oxycephalus in Japan. For each of 499 individuals from 71 localities, 506 bp of an mtDNA sequence from the 3ʹ end of the cytochrome b region was obtained and 187 haplotypes were identified. The haplotypes were differentiated into two large groups: Group I and Group II. Group I was divided into two subgroups (Subgroups A and B) and Group II was divided into five subgroups (Subgroups C, D, E, F, and G). The mean of the genetic heterogeneity value (F_ST) between each subgroup was 0.755, and the mean of the sequence divergence between each group was 5.5%. These results indicate a high level of genetic diversity and geographic differentiation.

Subgroup A was distributed from Tottori Prefecture and west to Kyushu on the Sea of Japan side. Subgroup B was limited to some rivers on the Pacific Ocean side of Shikoku. Subgroup C was distributed along the inland sea of the Kii Channel to the eastern Seto Inland Sea, mainly in Tokushima, Wakayama, and Hyogo Prefectures. Subgroup D was distributed only in the Kuzuryu River. Subgroup E was distributed in the area around Ise Bay and in central and western Shizuoka Prefecture. Subgroup F was collected from one river in Toyama Prefecture and two rivers in Aichi Prefecture. Subgroup G was distributed in the Ado River and several other rivers, and the Yahagi River was also likely a natural distribution region. The results suggest a relationship between the time of divergence between each subgroup and past geograpic events in Japan. The individuals collected in four localities (Tama River, Niizaki River, and two rivers in Yamagata Prefecture) are genetically closely related to the haplotypes in Lake Biwa, suggesting those fish may be alien species from Lake Biwa.

Translated Abstract

日本のタカハヤの集団構造を明らかにするため、mtDNA シーケンス分析による集団遺伝学的解析を行った。71標本群499個体について mtDNA cytochrome b 領域 3’末端側の506 bp の塩基配列を調べ、187種類のハプロタイプを決定した。その結果、日本のタカハヤ集団は、大きく2つのグループ（Group I と Group II）に分化し、グループ内でも Group I で2つ（Subgroup A, B）、Group II で5つ（Subgroup C～G）のサブグループを形成していた。グループ間の F_ST の平均値は0.755、各グループ間の塩基置換率の平均値は5.5%であり、遺伝的多様性および地理的分化レベルは高かった。Aグループは鳥取以西から九州に分布していた。Bグループは四国の太平洋側に、Cグループは徳島県・和歌山県・兵庫県を中心とした紀伊水道～東瀬戸内海の内海沿いに分布していた。Dグループの出現は九頭竜川のみであった。Eグループは、伊勢湾周辺域と静岡県中西部で出現した。Fグループの出現は、主に愛知県から富山県までにかけて南北であった。Gグループは琵琶湖で出現し、愛知県の矢作川は在来の可能性があった。各グループ間の分岐年代と日本における過去の地理的事象との関係性が示唆された。多摩川・新崎川と山形県に分布する本種は、琵琶湖流入河川のサンプルと遺伝的に近縁となり、琵琶湖からの侵入の可能性が高いと考えられた。

Introduction

Because the migratory dispersal of pure freshwater fishes is restricted to freshwater areas, genetic differentiation is likely to occur, and genetic population structures are formed at various geographic scales within each fish species^1-2). Information on genetic population structure is important not only for species conservation but also for considering conservation measures for local populations¹⁾. The Japanese archipelago has a complex topography with many mountain ranges and straits, which creates many geographical barriers for pure freshwater fishes, making it difficult for their distribution areas to expand. Under these conditions, differences are observed, with some species having a broad distribution and others having a local distribution. For example, the genetic differentiation of medaka ricefish species has occurred in various regions of the Japanese archipelago^3-4). Biwia zezera is locally distributed in the Nobi Plain, Lake Biwa-Yodo River system, Sanyo region, and northwestern Kyushu⁵⁾, and it is genetically differentiated into the Lake Biwa type, Gifu type, and Kyushu type⁶⁾. Topographic changes and climatic variations in the Japanese archipelago, as well as ecological differences among fish species, have affected the genetic population structure of each fish species. Genetic population structure has been studied in various fish species such as Odontobutis obscura⁷⁾, medaka ricefish⁴⁾, Cobitis biwae⁸⁾, Tanakia lanceolata⁹⁾, and Psudobagrus nudiceps¹⁰⁾. Comparison of these genetic population structures has provided knowledge on the distribution range formation process of Japanese freshwater fish species¹¹⁾. In recent years, differences in geographic and genetic differentiation have led to an increase in the number of new species described, including medaka ricefish (Oryzias sakaizumii, Oryzias latipes), Pseudogobio (Pseudogobio polystictus, Pseudogobio agathonectris), and Silurus tomodai^11-14).

Rhynchocypris oxycephalus is a pure freshwater fish of the family Cyprinidae, which is mostly distributed over southeastern Japan from Shizuoka Prefecture facing the Pacific Ocean and from Toyama Prefecture facing the Sea of Japan. The habitation area of Rhynchocypris oxycephalus is a upper and middle reaches of rivers and mountain lakes and marshes. That distribution area is similar to the habitation area of Oncorhynchus masou ishikawae and is further upstream of the habitation area of the Candidia temminckii¹⁵⁾. As a pure freshwater fish with a wide distribution west of Shizuoka and Toyama Prefectures, Rhynchocypris oxycephalus is useful as an indicator organism for understanding the distribution range formation process of Japanese freshwater fish species that inhabit the upper reaches of rivers in western Japan. Understanding the genetic population structure of this species should also provide insight into release guidelines for commercially valuable fish species in the upper reaches of rivers. For example, it is difficult to determine the original genetic population structure of commercially valuable upstream fish species such as Salvelinus leucomaenis and Oncorhynchus masou ishikawae because their distribution areas have been disturbed by artificial releases for recreational fishing^16-23).Thus, genetic management of native populations of these fish species is already difficult. On the other hand, freshwater fish with low commercial value, such as Rhynchocypris oxycephalus, are less likely to be released intentionally. If we can understand the geographical differentiation of these fish species, it should be possible to infer the original distribution areas of Salvelinus leucomaenis and Oncorhynchus masou ishikawae, which are losing their original genetic population structure, and to provide guidelines for future releases. Although there have been reports on the morphological and genetic differences between Rhynchocypris oxycephalus and a closely related species, Rhynchocypris lagowskii, which is distributed mainly in eastern Japan, and on ecological differences in coexisting rivers^24-25), there has been no comprehensive genetic analysis of the domestic distribution area. Therefore, the objective of this study was to understand the genetic population structure of Rhynchocypris oxycephalus based on the mtDNA cytochrome b gene region¹¹⁾, which is commonly used in phylogeographical analyses of freshwater fishes, and to gain insight into the genetic diversity and distribution range formation process.

Materials and Methods

A total of 499 individuals of Rhynchocypris oxycephalus were collected from 71 sampling localities in 68 rivers by fishing or landing net from August 2008 to November 2012 (Fig. 1). After collecting the fish, the fins were cut and stored in pure ethanol.

Fig. 1. Sampling locations of the upstream fat minnow Rhynchocypris oxycephalus in Japan. Locality numbers correspond to those in Table 1.

DNA was extracted from the fins using the QuickGene SP kit DNA tissue (Kurabo Industries Ltd.). The extracted DNA was then subjected to PCR of the cytochrome b domain of mtDNA, using TaKaRa r Taq or Ex Taq (Takara Bio, Inc.) as the reaction enzyme. The primers referred to as L14391 (5´ATGGCAAGCCTACGAAAAAC 3´) and H15551 (5´GATTACAAGACCGATGCTTT 3´) were used for the PCR²⁶⁾. TaKaRa PCR Thermal Cycler MP was used, after performing 94℃ preheating for 1 min, a temperature cycle of 94℃ for 15 s (heat denaturation process), 50℃ for 15 s (annealing process), 72℃ for 30 s (extension process) repeated for 30 cycles, and finally an extension process at 72℃ for 7 min.

The PCR product was cleaned to remove excess primer using ExoSAP-IT (USB Corporation). The sequencing reactions were performed with H15551 primer using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystem). TaKaRa PCR Thermal Cycler MP was used for the sequencing reaction. DNA sequencing was performed using ABI3130 Genetic Analyzer Avant (Applied Biosystems). The obtained DNA sequences were aligned using the Clustal W method²⁷⁾. Determining the haplotypes, calculating the base substitution rate based on Kimura’s two-parameter method (K2P)²⁸⁾, and estimating the neighbor-joining tree²⁹⁾ were done using MEGA 5.1 software.

The cytochrome b sequence data of Rhynchocypris lagowskii (AB198969) was quoted from DDBJ as an outgroup at the time the genealogical tree was created. The reliability of each node of a phylogenetic tree is estimated as bootstrap probability by 1,000 repeated calculations³⁰⁾. Based on obtained haplotype sequence data, haplotype diversity (h), nucleotide diversity (π), and the genetic divergence coefficient F_ST value were estimated, and mismatch distribution among individuals of each sampling point were analyzed using Arlequin ver. 3.5.1.2³¹⁾. The false discovery rate (FDR) was controlled at q < 0.05 based on the Benjamini Hochberg method. The mismatch distribution represented a frequency distribution of the base sequence differences between individuals. The observed distribution was compared with the simulated one using the sudden expansion model with 1,000 bootstrap iterations. The significance of difference was evaluated based on the sum of squared deviation (SSD)³²⁾. To estimate the population interconnection by migration in each locality, we used the analysis of the spatial expansion³³⁾. The haplotype network was depicted using TCS ver. 1.21³⁴⁾. For presumption of branch age, the molecular clock (cytochrome b gene region: 0.0152 nucleotide substitution per site per one million years) of the Cyprinidae fishes in Europe was applied³⁵⁾. The DNA sequences determined in this study have been deposited in the DNA Data Bank of Japan (DDBJ) under accession numbers LC154940-LC154946.

Results

Sequence variation and haplotypes

For each of the 499 individuals from 71 localities, 506 bp of an mtDNA sequence from the 3ʹ end of the cytochrome b region was obtained. In this study, 187 haplotypes (Hap 1 to 187) were identified. One to seven haplotypes were observed in each locality (Table 1, Table 2).

Table 1. Sampling localities, number of individuals (N), sampling day, and haplotypes detected in each locality

Table 2. Mitochondrial sequence variation from 506 bp of the cytochrome b region in R. oxycephalus. Dots indicate the same state as in the first haplotype

Table 2. Continued (Hap91 to Hap187; Nucleotide position 1 to 182)

Table 2. Continued (Hap1 to Hap90; Nucleotide position 183 to 506)

Table 2. Continued (Hap91 to Hap187; Nucleotide position 183 to 506)

Nucleotide sequence divergences in each locality

The net nucleotide sequence divergences in each locality ranged from 0.0 to 3.5%, and the overall mean of the net nucleotide sequence divergence was 0.4%. The largest value of the K2P distance was 3.5% in the Kaifu River (locality No.3). The lowest value was 0.0% in the None River (No.2), the Naka River (No.4), the Katsuura River (No.5), the Gokase River (No.17), the Beppu River (No.19), the Koishiwara River in the Chikugo River basin system (No.22), the Matsuura River (No.29), the Takano River in the Yodo River basin system (No.49), the Koza River (No.54), the Sakanai River (No.57), and the Azu River (No.69) (Table 3). The net nucleotide sequence divergence between each locality was 0.0 to 8.5%, and the mean of net nucleotide sequence divergence between each locality was 4.1%. The largest value was 8.5% between the Monobe River (No.1) and the Kurobe River (No.63). The lowest value was 0.0% between the Naka River (No.4) and the Katsuura River (No.5).

Table 3. K2P distances within each locality of R. oxycephalus

Genetic diversity (haplotype diversity and nucleotide diversity) of each locality

The genetic diversity (haplotype diversity and nucleotide diversity) of each locality is shown in Table 4. The largest value of haplotype diversity was 1.0000 ± 0.2722 in the Kumano River (locality No.55). The largest value of nucleotide diversity was 0.0333 ± 0.0194 in the Kaifu River (No.3). In nine rivers [the None River (No.2, Hap110), the Naka River (No.4, Hap169), the Katsuura River (No.5, Hap169), the Beppu River (No.19, Hap40), the Matsuura River (No.29, Hap42), the Takano River in the Yodo River basin system (No.49, Hap163), the Koza River (No.54, Hap178), the Sakanai River (No.57, Hap143), and the Azu River (No.69, Hap41)], the haplotype diversity and nucleotide diversity were 0.0000±0.0000.

Table 4. Haplotype and nucleotide diversities in 71 population samples of R. oxycephalus

Phylogenetic trees and haplotype networks

The dendrogram estimated from the genetic distance between each haplotype is shown in Fig. 2, and the haplotype network of each haplotype is shown in Fig. 3. Rhynchocypris oxycephalus of Japan was found to comprise two major groups (I and II). Group I was divided into two subgroups (Subgroups A and B) and Group II was divided into five subgroups (Subgroups C, D, E, F, and G). The mean net nucleotide sequence divergence was 6.4% between I and II. Subgroup A consisted of 95 haplotypes (Hap1-95), and Subgroup B consisted of 14 haplotypes (Hap96-109). Individuals in Subgroup A were collected from the Kyushu region, Chugoku region, and the West Shikoku region. Individuals in Subgroup B were collected from the East Shikoku region. The population in the area northeast of the Kinki region constituted five subgroups: Subgroup C (31 haplotypes, Hap157-187), Subgroup D (2 haplotypes, Hap155-156), Subgroup E (18 haplotypes, Hap137-154), Subgroup F (12 haplotypes, Hap125-136), and Subgroup G (15 haplotypes, Hap110-124).

Fig. 2. Phylogenetic relationships among haplotypes based on the cytochrome b sequence (506 bp): neighbor-joining (NJ) tree based on sequence divergence by Kimura (1980). Bootstrap probabilities with 1,000 replications are shown for each cluster. Some geographical groups and subgroups were formed based on bootstrap values of more than 46.

Fig. 3. aplotype network tree of 187 haplotypes. Numbered circles are observed haplotypes and black circles (●) are virtual haplotypes in this study. Each group circled for neighbor-joining (NJ) tree.

The net nucleotide sequence divergence in each subgroup and between each subgroup

The net nucleotide sequence divergence in each subgroup and between each subgroup was estimated using Kimura’s two-parameter method²⁸⁾. The largest value was 1.8% in Subgroup A, the lowest value was 0.8% in Subgroup D, and the average value of all subgroups was 1.3%. Among the values between each subgroup, the largest was 8.1% between Subgroups B and F, the lowest was 1.6% between Subgroups F and G, and the average between each subgroup was 5.5% (Table 5).

Table 5. K2P distances among 7 subgroups of R. oxycephalus

Genetic diversity (haplotype diversity and nucleotide diversity) of each subgroup

The genetic diversity (haplotype diversity and nucleotide diversity) of each subgroup is shown in Table 6. The highest value of haplotype diversity was 0.9755 ± 0.0036 (Subgroup A), the lowest value was 0.2500 ± 0.1802 (Subgroup D), and the average value of all subgroups was 0.8108. The highest value of nucleotide diversity was 0.0149 ± 0.0077 (Subgroup A), and the lowest value was 0.0020 ± 0.0017 (Subgroup D), while the average value of all subgroups was 0.0091. The pairwise F_ST values between each subgroup exhibited significant differences (p < 0.05; Benjamini Hochberg corrected FDR q < 0.05) in all pairwise combinations. The average pairwise F_ST value between subgroups was 0.755.

Table 6. Haplotype and nucleotide diversities in 7 subgroups of R. oxycephalus

The mismatch distribution in each subgroup

The mismatch distribution was analyzed for the cytochrome b region in each subgroup. Four subgroups (A, B, C, and F) presented unimodal curves. On the other hand, two subgroups (E and G) presented multimodal curves (Fig. 4). Compared with the simulated expected distribution based on the population expansion model, all p-values for each mismatch were p > 0.05 and were not significantly different.

Fig. 4. Mismatch distribution performed for cytochrome b region on each subgroup. Observed (bars) and simulated distributions (squares) of pairwise sequence differences under the spatial expansion model.

Distribution of each subgroup

The mixing ratio of each subgroup in each locality is shown in Fig. 5. Subgroup A was distributed from Tottori Prefecture and west to Kyushu on the Sea of Japan side. The island groups (localities No.69-71), such as Tsushima Island and the Goto Islands, belong to this subgroup. Subgroup B was localized to some rivers on the Pacific Ocean side of Shikoku. Subgroup C was distributed along the inland sea of the Kii Channel to the eastern Seto Inland Sea, mainly in Tokushima, Wakayama, and Hyogo Prefectures. Subgroup D was distributed only in the Kuzuryu River (No.64). Subgroup E was distributed in the area around Ise Bay, such as the Sakanai River (No.57), the Kumozu River (No.58), the Toyo River (No.60), the Abe River (No.62), and the Ooi River (No.61). Subgroup F was collected from one river in Toyama Prefecture (No.63) and two rivers in Aichi Prefecture (No.59 and 60). Subgroup G was distributed in the Ado River (No.48), the Yahagi River (No.59), the Niizaki River (No.66), the Tama River (No.65), the Aka River (No.67), and the Ira River (No.68).

Fig. 5. The mixing ratio of each subgroup in each locality. Locality numbers correspond to those in Table 1.

Discussion

Genetic diversity of Rhynchocypris oxycephalus

The pairwise F_ST values between each subgroup exhibited significant differences (p < 0.05; Benjamini Hochberg corrected FDR q < 0.05) in all pairwise combinations. The average pairwise F_ST value between subgroups was 0.755. Generally, “values of F_ST above 0.25 indicate very great genetic differentiation”³⁶⁾ in wild populations. The high F_ST values in this study indicated that Rhynchocypris oxycephalus had genetic differentiation in each locality and maintained high genetic divergence in total populations.

The average value (5.5%) of the base substitution rate between different subgroups was close to 5.1%, which is the value between Japanese rose bitterling Rhodeus ocellatus kurumeus and Chinese rose bitterling (between a subspecies)³⁷⁾. The largest value (8.1% between Subgroups B and F) among different subgroups in this data was close to the value (8.0%) between each species of striated spined loach kind groups in the Kyushu district and the Korean Peninsula³⁸⁾, which, accordingly, suggests that the level of genetic diversity and geographic differentiation within Rhynchocypris oxycephalus was high.

Four subgroups (A, B, C, and F) presented a unimodal curve with a mismatch distribution, suggesting that three subgroups (A, C, and F) have maintained high genetic diversity because the shapes of the mismatch distribution histograms present a stable maintenance type³⁹⁾. On the other hand, Subgroup B presented a unimodal curve with a peak at three base substitutions, representing a mismatch distribution of simultaneous diffusion type. These results suggested that the area of distribution of Subgroup B began to expand in relatively recent years³⁹⁾.

Subgroup D had two haplotypes (Hap155 and Hap156). The four bases of substitution had formed between Hap155 and Hap156. These unique haplotypes were collected only in the Kuzuryu River. These haplotypes varied in 18 bases from the nearest haplotype (Hap 163 in Subgroup C), indicating that this population has decreased based on bottleneck effects in the past. However, only eight individuals were analyzed from the group. More individuals from the Kuzuryu River and outskirts are needed to investigate the characteristics of the area in detail.

Subgroup E presented bimodal curves and two gentle peaks. The simulated curve has one peak and a large dispersion of nucleotide substitution distances with a unimodal function type (p = 0.34), suggesting that this pattern is formed with a stable population maintained for long periods of time³⁹⁾.

Subgroup G presented bimodal curves and two peaks around one base and five bases. However, it was not significantly different from a unimodal function type with a peak around one base (p = 0.59), suggesting that the distribution area has expanded in relatively recent years. The type appears to be intermediate between the stable maintenance type and the simultaneous diffused type.

Most subgroups had a few peaks and a large dispersion of nucleotide substitution distances with a unimodal function type. We believe that Rhynchocypris oxycephalus has maintained a stable population in each group.

Genetic population structure and formation process of the population

・Diverging between Group I (west group) and Group II (east group)

Our phylogenetic analysis (NJ tree) suggests Rhynchocypris oxycephalus consists of two distinctly diverged groups. Group I (west group) includes the region of Shikoku except a part of East Shikoku (around Tokushima) and of Chugoku and Kyushu, and Group II (east group) includes the areas of East Shikoku, Kinki, and Chubu. The average nucleotide sequence divergence between Group I and Group II was estimated as 6.4%. The estimated divergence time of Group I from Group II was calculated to be about 2.1 million years using the molecular clock of Zardoya and Doadrio³⁵⁾. The Seto Inland Sea did not exist at the time (1.8-0.78 million years ago). The islands of Honshu, Shikoku, and Kyushu were connected and existed as one big island. In the first half of the Pleistocene era (about 1.8-0.78 million years ago), the shoreline was extended by the incursion of the Bungo Channel and Kii Channel, forming the present Seto Inland Sea. This estimated branch age (about 2.1 million years ago) between the west group (Group I) and the east group (Group II) was consistent with the time in the first half of the Pleistocene era⁴⁰⁾. The specialization (about 2.1 million years ago) between the two groups was possibly caused by the old river dividing to the west and the east, which had flowed to the west from old Lake Biwa to the second Seto Inland sedimentation zone.

・Diverging in the Kinki and Chubu regions

Group II was divided into five subgroups (C, D, E, F, and G) in the NJ tree. The Suzuka Mountains and the Kii Mountains exist at the border between Subgroup C and the other three subgroups (E, F, and G). The branch age between Subgroup C and those subgroups was calculated to be about 1.75 million years. The Kii Mountains underwent upheaval in the late Quaternary (about 1.8 million years ago)⁴¹⁾ and the Suzuka Mountains about 1.5-1.0 million years ago⁴²⁾. The genetic divergence between Subgroup C and the other three subgroups (E, F, and G) may have resulted from the formation of these mountains. There are some reports about morphological and genetical differences in freshwater fishes between the east and west sides of the Suzuka Mountains, for example, the medaka Oryzias latipes⁴³⁾, the three-spined stickleback Gasterosteus aculeatus⁴⁴⁾, and the dark chub Candidia temminckii⁴⁵⁾. Furthermore, it has been reported that this area was a distribution boundary of stumpy bagrid catfish Pseudobagrus ichikawai and forktail bullhead Psudobagrus nudiceps. The Suzuka Mountains may have been a physical barrier and resulted in the geographical isolation of many freshwater fishes.

・Diverging around the Noto Peninsula

The branch age between Subgroups D and C was estimated at about 1.6 million years. This branch age was a value near the second half of the Pliocene (about 3.6-2.58 million years ago), when the present Fukui Prefecture including Noto Peninsula connected with the Japanese Islands⁴⁰⁾. After the second half of the Pliocene, the area of Noto Peninsula dissociated from the Japanese Islands, and the detached island was formed in the first half of the Pleistocene era (about 2.0-0.8 million years ago)⁴⁶⁾. Subgroup D may have been an isolated group on the detached island, and that isolated group formed a new subgroup (D) by invading this place when the detached island was connected with the mainland. Also, the local population distributed on the west of Noto Peninsula has been reported as a unique population of medaka Oryzias latipes⁴⁾. The geological changes of the past in this area have been suggested as a common physical barrier for several freshwater fishes.

・Diverging in the Shikoku region

Group I was divided into two subgroups (A and B). Subgroup B was isolated geographically and only distributed southwards bordering the Shikoku Mountains. The branch age between Subgroup B and Subgroup A has been estimated to be about 0.98 million years, suggesting that the upheaval of the Shikoku Mountains took place by the start of new subduction movement of the Philippine Sea Plate in the first half of the Pleistocene era (about 1.0 million years ago)⁸⁾. This time was almost consistent with the branch age of 0.98 million years ago. Consequently, one subpopulation in the primitive group of Subgroup A may have been isolated on the south side of the Shikoku Mountains, and differentiated into Subgroup B and Subgroup A from the primitive group. The genetic population structures of the spined loach Cobitis biwae and medaka Oryzias latipes are consistent with this result^4,47). The phylogeographic barrier of the Shikoku Mountains has possibly influenced the formation of a unique population structure in the freshwater fish of south Shikoku including Kochi Prefecture. However, three subgroups (A, B, and C) existed in Shikoku, and Subgroups B and C existed in the eastern and western parts of Kochi⁴⁸⁾. It has been reported that high rank terrace formation had occurred in the Pleistocene era (about 0.8-0.14 million years ago) in southern Shikoku, and one big basin system (old Tosa River including the Niyodo River and the Monobe River) had formed there. After that, the maximum regression term came around the time of lower order terrace formation of the Pleistocene era (second half about 20,000-10,000 years ago). The old Niyodo River, the old Urado River (present-day Kagami River), and the old Monobe River were then formed⁴⁹⁾. Therefore, we believe that the three rivers (Niyodo River, Kagami River, and Monobe River) formed one group (Subgroup B).

・Diverging on Tsushima Island and the Goto Islands and in the Kyushu region

Tsushima Island and the Goto Islands were included in Subgroup A. It is thought that Tsushima Strait appeared in the second half of the Pleistocene era (about 10,000 years ago), and Tsushima Island and the Goto Islands were divided with Kyushu about 8,500 years ago⁵⁰⁾. Genetic difference was not seen in Kyushu, Tsushima Island, and the Goto Islands. It has been suggested that gene exchange between these islands and Kyushu (Subgroup A) occurred relatively recently.

・Diverging by landward migration of the Bungo and Kii Channels

Subgroups A and C were distributed along the inland sea from the Bungo Channel to Kii Channel including the western side of the Sea of Japan and the East China Sea in Japan. The sea level had fluctuated by about 100 m through repeated glacial epochs and the interglacial period in the second half of the Pleistocene era (about 1.8 million-10,000 years ago). Twenty thousand years ago, the sea level had fallen 100-140 m, and most of the Seto Island Sea had become land⁵¹⁾. In this area, it is assumed that the present-day Shiwaku Islands divided the Seto Inland Sea into two water systems in the west and the east⁵²⁾. The genetic divergence between Subgroups A and C accorded with that branch age. It has been reported that the populations in the east and west in the Seto Island Sea were genetically different, for example, Japanese freshwater goby Rhinogobius flumineus⁵³⁾, freshwater goby Odontobutis obscura⁷⁾, forktail bullhead Pseudobagrus nudiceps¹⁰⁾, and medaka Oryzias latipes⁵⁴⁾, etc. This finding indicated that the formation of the two water systems greatly affected the genetic differentiation of the east-west populations in the Seto Island Sea.

Origin of the domestic alien species

The natural distribution of Rhynchocypris oxycephalus is east of Shizuoka Prefecture and Toyama Prefecture¹⁵⁾. The populations of the Tama River or the rivers in Yamagata Prefecture have been non-native fish. However, the origin of Rhynchocypris oxycephalus in the Tama River and some rivers in the Yamagata populations is not clear. In this report, the haplotypes of one locality of the Tama River and two localities of Yamagata Prefecture were similar to the haplotypes of individuals collected in the Ado River locality, Lake Biwa. The haplotypes Hap110-115 from three localities of the Tama River and Yamagata Prefecture were different from Hap116-117 of the Ado River in 2-3 bases of those sequences. For example, it has been reported that genetic differentiation has produced Japanese freshwater goby in the river that flows into the west side of Lake Biwa and the river that flows into the east side⁵³⁾. If genetic differences also occur between rivers west and east of Lake Biwa for Rhynchocypris oxycephalus, then Rhynchocypris oxycephalus in the Tama River, Aka River, and Ira River may be a group other than those on the west side of Lake Biwa.

If Rhynchocypris oxycephalus in Yamagata Prefecture is a native species, then Rhynchocypris oxycephalus should be distributed along the Sea of Japan north of Fukui Prefecture. Also, if we assume that Rhynchocypris oxycephalus north of Fukui Prefecture became extinct and only the Rhynchocypris oxycephalus of Yamagata Prefecture survived, the Rhynchocypris oxycephalus of Yamagata Prefecture should be genetically very different from those in the Ado River sample group as a result of the natural expansion of their distribution area. This finding suggests that it is highly unlikely that the Yamagata Prefecture Rhynchocypris oxycephalus is native to the area. In addition, if the Tama River Rhynchocypris oxycephalus is a native species, it should be genetically closer to the geographically adjacent Abe and Oi Rivers than to the geographically distant Ado River. The distribution of other groups (Mie and Shizuoka groups and Hokuriku and Aichi groups) between the Ado River and Tama River sample groups also makes it highly likely that the Tama River Rhynchocypris oxycephalus is an invasive species. It will be necessary to investigate the genes of Rhynchocypris oxycephalus from the rivers of the Lake Biwa basin system to determine the origin of the Tama River, Aka River, and Ira River groups.

It has been reported that the Yahagi River was a natural distribution region of Rhynchocypris oxycephalus. Those localities had some haplotypes of Subgroup G including the Ado River (Lake Biwa). Those haplotypes had substitutions of 3-8 bases between the Ado River and the Yahagi River. Of the samples from the Yahagi River, 57% had some haplotypes in Subgroup G. The number of base exchanges (3-8 bases) between the Yahagi River and the Ado River populations is not small. Another question is why did the invader inhabit only the Yahagi River without inhabiting nearby rivers, if the fish in the Yahagi River were the invaders from Lake Biwa. In other words, the fish in the Yahagi River may not have invaded from Lake Biwa. It is thought that the old Lake Biwa was the Iga area, about 5.0 million years ago. The old lake had gradually moved north over 1.0 million years, and reached its present position about 1.0 million years ago⁵⁵⁾. The genetic similarity between the Yahagi River and Lake Biwa may have formed through the genetic exchange of each individual when those areas were previously the same place. Therefore, it is necessary to investigate the populations of Nobi Plain (east side of Lake Biwa).

On the other hand, the haplotype difference between the Niizaki River and the Ado river population was only 1-2 bases. That difference was the same as that of the Tama River and the Ira River (Yamagata Prefecture). The Ado River and the Niizaki River were geographically separated. Moreover, the populations that inhabited the area between the Ado River and the Niizaki River were identified in other subgroups (E and F). Therefore, the Niizaki River locality was geographically similar to the Tama River locality. In addition, the sample from the Tama River was genetically close to the sample from the Niizaki River. This genetic closeness suggests that the Niizaki River samples were invaders from Lake Biwa. To determine the eastern limit of the natural distribution region of Rhynchocypris oxycephalus on the Pacific Ocean side of Honshu, it will be necessary to investigate the area around Izu Peninsula.

Acknowledgments

We appreciate the help of Akira Ishida and Yusuke Ishida. In addition, without the English support and guidance of Winnie Naa Adjorkor Sowah, this paper would not have been possible. We also thank Takanori Ishikawa (Nihon University) and Natsuo Okabe (Yamagata Prefecture) for their help in collecting the specimens.

References

• 1)   Avise , J. C. (2000) Phylogeography: The history and formation of species, Harvard University Press, Cambridge, pp. 447.
• 2)   Watanabe,  K. and  M.  Nishida (2003) Freshwater fishes, in Conservation Genetics, H. Koike and M. Matsui (eds.), University of Tokyo Press, Tokyo, pp. 227-240. (In Japanese)
• 3)   Sakaizumi,  M. (1987) Biogeography of killifish using genetic methods. Biological Science: Genetics, 41(12): 17-22. (In Japanese)
• 4)  Takehana, Y., N. Nagai, M. Matsuda, et al. (5 co-authors) (2003) Geographic variation and diversity of the cytochrome b gene in Japanese wild populations of medaka, Oryzias latipes. Zoological Science, 20(10): 1279-1291.
• 5)   Hosoya,  K. (2002) Biwia zezera, in Freshwater Fishes of Japan Revised Edition, H. Kawanabe, N. Mizuno and K. Hosoya (eds.), Yama-Kei Publishers Co., Ltd., Tokyo, pp. 317. (In Japanese)
• 6)   Horikawa,  M.,  J.  Nakajima and  T.  Mukai (2007) Distribution of indigenous and non-indigenous mtDNA haplotypes of Biwia zezera (Cyprinidae) in northern Kyushu. Japanese Journal of Ichthyology, 54(2): 149-159. (In Japanese with English abstract)
• 7)   Sakai,  H.,  C.  Yamamoto and  A.  Iwata (1998) Genetic divergence, variation and zoogeography of a freshwater goby, Odontobutis obscura. Ichthyological Research, 45: 363-376.
• 8)  Kitagawa, T., M. Watanabe, E. Kitagawa, et al. (6 co-authors) (2003) Phylogeography and the maternal origin of the tetraploid form of the Japanese spined loach, Cobitis biwae, revealed by mitochondrial DNA analysis. Ichthyological Research, 50: 318-325.
• 9)   Itai,  K.,  N.  Nagata,  J.  Kitamura,  K.  Watanabe and  T.  Sota (2020) Phylogeography of the bitterling Tanakia lanceolata (Teleostei: Cyprinidae) in Japan inferred from mitochondrial cytochrome b gene sequences. Ichthyological Research, 67: 105-116.
• 10)   Watanabe,  K. and  M.  Nishida (2003) Genetic population structure of Japanese bagrid catfishes. Ichthyological Research, 50: 140-148.
• 11)  Watanabe, K., H. Takahashi, A. Kitamura, et al. (12 co-authors) (2006) Biogeographical history of Japanese freshwater fishes: Phylogeographic approaches and perspectives. Japanese Journal of Ichthyology, 53(1): 1-38. (In Japanese with English abstract)
• 12)   Asai,  T.,  H.  Senou and  K.  Hosoya (2011) Oryzias sakaizumii, a new ricefish from northern Japan (Teleostei: Adrianichthyidae). Ichthyological Exploration of Freshwaters, 22(4): 289-299.
• 13)   Tominaga,  K. and  S.  Kawase (2019) Two new species of Pseudogobio pike gudgeon (Cypriniformes: Cyprinidae: Gobioninae) from Japan, and redescription of P. esocinus (Temminck and Schlegel 1846). Ichthyological Research, 66: 488-508.
• 14)   Hibino,  Y. and  R.  Tabata (2018) Description of a new catfish, Silurus tomodai (Siluriformes: Siluridae) from central Japan. Zootaxa, 4459(3): 507-524.
• 15)   Itai,  T. (2002) Rhynchocypris oxycephalus, in Freshwater Fishes of Japan Revised Edition, H. Kawanabe, N. Mizuno and K. Hosoya (eds.), Yama-Kei Publishers Co., Ltd., Tokyo, pp. 274. (In Japanese)
• 16)  Yamamoto, S., K. Morita, S. Kitanomura, et al. (7 co-authors) (2004) Phylogeography of white-spotted charr (Salvelinus leucomaenis) inferred from mitochondrial DNA sequences. Zoological Science, 21: 229-240.
• 17)   Kubota,  H.,  T.  Doi,  S.  Yamamoto and  S.  Watanabe (2007) Genetic identification of native populations of fluvial white-spotted charr Salvelinus leucomaenis in the upper Tone River drainage. Fisheries Science, 73: 270-284.
• 18)   Kawamura,  K.,  M.  Kubota,  M.  Furukawa and  Y.  Harada (2007) The genetic structure of endangered indigenous populations of the amago salmon, Oncorhynchus masou ishikawae, in Japan. Conservation Genetics, 8: 1163-1176.
• 19)  Kikko, T., M. Kuwahara, K. Iguchi, et al. (7 co-authors) (2008) Mitochondrial DNA population structure of white-spotted charr (Salvelinus leucomaenis) in the Lake Biwa water system. Zoological Science, 25: 146-153.
• 20)   Kawamura,  K.,  M.  Furukawa,  M.  Kubota and  Y.  Harada (2012) Effects of stocking hatchery fish on the phenotype of indigenous populations in the amago salmon Oncorhynchus masou ishikawae in Japan. Journal of Fish Biology, 81: 94-109.
• 21)   Okabe,  M. and  A.  Komatsu (2012) Investigation of population structure and identification of a native population of amago salmon Oncorhynchus masou ishikawae in the Nahari River system, Kochi Prefecture, Japan. Aquaculture Science, 60(1): 89-97. (In Japanese with English abstract)
• 22)  Miyahara, H., H. Yamada, T. Sato, et al. (6 co-authors) (2012) Mitochondrial-nuclear discordance in the amago salmon, Oncorhynchus masou ishikawae, in the River Miya, Japan. Conservation Genetics, 13: 1343-1353.
• 23)  Iwatsuki, Y., H. Tanaka, S. Seki, et al. (5 co-authors) (2020) Biological characteristics of the six genetic group of complete cytochrome b genes (1141bp) on four subspecies of the Oncorhynchus masou complex and biogeographical genetic strains from northwestern pacific by Japanese explanation with notes (Iwatsuki et al., 2019). Nature of Kagoshima, 47: 5-16. (In Japanese)
• 24)   Fujita,  T. and  K.  Hosoya (2003). Biochemical and morphological comparison between two Japanese daces, Phoxinus lagowskii steindachneri and P. oxycephalus jouyi in the sympatric sites. Japanese Journal of Ichthyology, 50(1): 55-62. (In Japanese with English abstract)
• 25)   Higuchi,  F. and  K.  Watanabe (2005) Genetic diversity and hybridization in the cyprinid, Rhynchocypris lagowskii steindachneri, from Yokohama, central Honshu. Japanese Journal of Ichthyology, 52(1): 41-46. (In Japanese with English abstract)
• 26)  Sasaki, T., Y. P. Kartavtsev, S. N. Chiba, et al. (6 co-authors) (2007) Genetic divergence and phylogenetic independence of Far Eastern species in subfamily Leuciscinae (Pisces: Cyprinidae) inferred from mitochondrial DNA analyses. Genes & Genetic Systems, 82: 329-340.
• 27)   Thompson , J. D., D. G.  Higgins  and T. J.  Gibson  (1994) CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Research, 22: 4673-4680.
• 28)   Kimura,  M. (1980) A simple method for estimating evolutionary rates of base substitution through comparative studies of nucleotide sequences. Journal of Molecular Evolution, 16: 111-120.
• 29)   Saitou,  N. and  M.  Nei (1987) The neighbor-joining method: A new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4(4): 406-425.
• 30)   Felsenstein,  J. (1985) Confidence limits on phylogenies: An approach using the bootstrap. Evolution, 39(4): 783-791.
• 31)   Excoffier,  L. and H. E. L.  Lischer  (2010) Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10(3): 564-567.
• 32)   Schneider,  S. and  L.  Excoffier (1999) Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: Application to human mitochondrial DNA. Genetics, 152(3): 1079-1089.
• 33)   Excoffier,  L.,  G.  Laval and  S.  Schneider (2007) Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evolutionary Bioinformatics, 23(1): 47-50.
• 34)   Clement,  M.,  D.  Posada and K. A.  Crandall  (2000) TCS: A computer program to estimate gene genealogies. Molecular Ecology, 9: 1657-1659.
• 35)   Zardoya,  R. and  I.  Doadrio (1999) Molecular evidence on the evolutionary and biogeographical patterns of European cyprinids. Journal of Molecular Evolution, 49: 227-237.
• 36)   Hartl , D. L. and A. G.  Clark  (1997) Principles of population genetics, Third Edition, pp. 117-120, Sinauer, Sunderland, Massachusetts.
• 37)  Miyake, K., H. Tachida, Y. Oshima, et al. (7 co-authors) (2001) Genetic variation of the cytochrome b gene in the rosy bitterling, Rhodeus ocellatus (Cyprinidae) in Japan. Ichthyological Research, 48: 105-110.
• 38)  Kitagawa, T., S. R. Jeon, E. Kitagawa, et al. (6 co-authors) (2005) Genetic relationships among the Japanese and Korean striated spined loach complex (Cobitidae: Cobitis) and their phylogenetic positions. Ichthyological Research, 52: 111-122.
• 39)   Koike,  H. (2003) Intraspecific polymorphism and conservation genetics, in Conservation Genetics, H. Koike and M. Matsui (eds.), University of Tokyo Press, Tokyo, pp. 51-52. (In Japanese)
• 40)   Ichihara,  M. (2001) Osaka formation, old Seto inland river lake system. Urban Kubota, 39: 10-13. (In Japanese)
• 41)   Kitamura,  A. (2010) Formation of the Japanese Islands and paleoenvironment, in Natural History in Geography of Freshwater Fish Kinds Concerning Diversity and Differentiation, K. Watanabe and H. Takahashi (eds.), Hokkaido University Press, Hokkaido, pp. 16-17. (In Japanese)
• 42)   Kawabe,  T. (1989) Stratigraphy of the lower part of the Kobiwako Group around Ueno basin, Kinki district, Japan. Journal of Geosciences Osaka City University, 32: 39-52.
• 43)   Sakaizumi,  M.,  K.  Moriwaki and  N.  Egami (1983) Allozymic variation and regional differentiation in wild populations of the fish Oryzias latipes. Copeia, 2: 311-318.
• 44)   Taniguchi,  N.,  Y.  Honma and  K.  Kawamata (1990) Genetic differentiation of freshwater and anadromous threespine sticklebacks (Gasterosteus aculeatus) from Northern Japan. Japanese Journal of Ichthyology, 37(3): 230-238.
• 45)  Okazaki, T., M. Watanabe, K. Mizuguchi, et al. (4 co-authors) (1991) Genetic differentiation between two types of dark chub, Zacco temmincki, in Japan. Japanese Journal of Ichthyology, 38(2): 133-140.
• 46)   Gouhara,  Y. (1975) Japan in ice age Japan of quaternary era. Urban Kubota, 11: 4-11. (In Japanese)
• 47)   Kagami,  H.,  T.  Mitsushio and  H.  Takemasa (1991) Study on a Morasse sediment from the Shikoku Mountains: Joyama formation of red conglomerate facies in Kochi, Japan. Josai University annual study report (natural science), 15: 1-19. (In Japanese with English abstract)
• 48)   Okino,  T.,  Y.  Tanaka and  S.  Seki (2012) Population structure of the upstream fat minnow (Rhynchocypris oxycephalus) based on mitochondrial DNA cytochrome b analysis in Shikoku. Fish Genetics and Breeding Science, 41(2): 181-186. (In Japanese with English abstract)
• 49)   Mitsushio,  T. (2001) Environmental change of quaternary era in Kochi prefecture. Report of the Institute of the Kuroshio Sphere, Kochi University. Special series, 13: 3-15. (In Japanese)
• 50)   Oshima,  K. (1990) The history of straits around the Japanese Islands in the late-quaternary. The Quaternary Research, 29(3): 193-208. (In Japanese with English abstract)
• 51)   Nagai,  K. (1987) Guide of Ehime Prefecture earth science, K. Nagai (ed.), Corona Publishing Co., Ltd., Tokyo, pp. 18-19. (In Japanese)
• 52)   Kuwashiro,  I. (1959) Topography of the seabed in Inland Sea of Seto. Geographical Criticism, 32(1): 24-35. (In Japanese)
• 53)   Shimizu,  T.,  N.  Taniguchi and  N.  Mizuno (1993) An electrophoretic study of genetic differentiation of a Japanese freshwater goby, Rhinogobius flumineus. Japanese Journal of Ichthyology, 39(4): 329-343.
• 54)   Sakaizumi,  M. (2000) Lineage and intraspecific structure of medaka. Protein, Nucleic Acid and Enzyme, 45: 2909-2917. (In Japanese)
• 55)   Miyoshi,  T. and  K.  Ishibashi (2008) Neotectonics in and around the Kinki Triangle in relation to the geometry of the subducted Philippine Sea plate. The Quaternary Research, 47: 223-232. (In Japanese with English abstract)