In total, 140 individuals were collected from 12 localities in 10 river systems in the Seto Inland Sea and Ise Bay regions from 2015 to 2020 (Fig. 1, Suppl. material 1), covering the whole geographic range of A. t. tabira. The unsampled sites in this study are the Nagara, Kiso, and Kushida river systems investigated by previous studies (Kitamura et al. 2012; Umemura et al. 2012), and rivers that flow into the Sea of Japan, which are probably extinct. For the populations in the Nagara, Kiso, and Kushida River, we cited base sequence data from previous studies (Kitamura et al. 2012; Umemura et al. 2012, Table 1). Additionally, a captive population from the Gifu World Freshwater Aquarium was included in this study. This captive population is believed to have originated from a population collected and cultured by a citizen from the Nagara River system in the Ise Bay area, which was donated to the Gifu World Freshwater Aquarium in 2004 (Koki Ikeya, personal communication). The mtDNA lineage of this population is described as native to the Ise Bay region by Mukai (2019). However, it is not indicated whether it belongs to the Nobi Plain group I or II. Individual bitterling were collected using hand nets and fishing methods. Each specimen was subjected to caudal fin clipping, and the remaining specimens were fixed in 10% formalin and preserved in 70% ethanol. For habitat conservation, the number of individuals collected was limited to 20 or fewer. The fin clips were preserved in 99% ethanol and stored at −20 °C. The specimens were registered with the Gifu and Mie Prefectural Museum along with collection site information (GPM-Z-22109, MIE-Fi3500, 3506–3512, 3519–3534, 3536–3540, 3543–3579, 3581–3604, 3606–3624, 3626, 3630–3633, 3718, 4272).

Figure 1.

Sampling localities of Acheilognathus tabira tabira. The asterisks indicate localities used by Kitamura et al. (2012), Umemura et al. (2012), and Kitamura and Uchiyama (2020). Circles show mtDNA lineages (see Fig. 2). Gray areas represent the estimated natural distribution area of the five subspecies of A. tabira based on Kitamura et al. (2012). This altitude map was used with permission from the Geospatial Information Authority of Japan (https://maps.gsi.go.jp/) and digital national and information of Ministry of Land, Infrastructure, Transport and Tourism (https://nlftp.mlit.go.jp).

Total genomic DNA was extracted from a portion of each caudal fin using the Kaneka Easy DNA extraction kit version 2 (Kaneka, Hyogo, Japan) or the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). Total genomic DNA was used to amplify DNA fragments using polymerase chain reaction (PCR). For PCR, the following forward primer was used: L14690-Cb-AH, 5'-GGT CAT AAT TCT TGC TCG GA-3' (Kitanishi et al. 2016), and for the reverse primers H15913- Thr-AH, 5'-CCG ATC TTC GGA TTA CAA GAC CG-3' (Kitanishi et al. 2016), or Cytb-Rev, 5'-GAT CTT CGG ATT ACA AGA CC-3' (Hashiguchi et al. 2006). The sequencing protocol was previously described by Ito et al. (2020). All sequences were deposited in the DNA Data Bank of Japan (DDBJ), European Nucleotide Archive (EMBL), and GenBank databases under the accession numbers LC775317–LC775352.

Multiple alignments of nucleotide sequences were performed using MUSCLE (Edgar 2004). For the phylogenetic analysis, nucleotide sequence data for A. t. tabira were obtained from DDBJ, EMBL, and GenBank (AB620138, AB620141, AB620150, AB620159, AB759881–AB759890, and LC578851, Kitamura et al. 2012; Umemura et al. 2012; Ito et al. 2021; Table 1). In addition, we used Acheilognathus tabira jordani Arai, Fujikawa & Nagata, 2007, a sister group of A. t. tabira, as outgroups (AB620149, and AB620156, Kitamura et al. 2012). Phylogenetic analyses were conducted using maximum likelihood (ML) and Bayesian inference (BI) methods. Search for the best evolutionary model each partition and ML analyses were performed using IQ-TREE 2.2.2.6 (Minh et al. 2020), with the TIM2 + F + I, K2P + I, and F81 + F models for the first, second, and third codon positions as selected by ModelFinder (Chernomor et al. 2016; Kalyaanamoorthy et al. 2017) based on the Bayesian information criterion (BIC). The reliability of each internal branch was evaluated using Shimodaira-Hasegawa-like approximate likelihood ratio test (SH-aLRT, Guindon et al. 2010) and ultrafast bootstrap (UFBoot, Hoang et al. 2018) and with 1000 replicates. In this study, according to the manual, support values of 80% or higher for SH-aLRT and 95% or higher for UFBoot were considered as high support values. BI analyses were performed using MrBayes v3.2.7a (Ronquist et al. 2012), with the HKY + I, K80 + I + G4, F81 models for the first, second, and third codon positions as selected by ModelTest-NG (Darriba et al. 2017) based on the BIC. The MCMC analyses were performed using the following settings: ngen = 10000000, sample freq = 100, and burnin = 25000. Statistical parsimony networks were constructed using TCS v1.2.1 (Clement et al. 2000). In this study, the Bayesian posterior probability (PP) of 80% or higher was considered to indicate high support values.

In populations of the Kinki-Sanyo region, we calculated genetic differentiation, estimated by genetic differentiation coefficient (Φ_st; Excoffier et al. 1992) using Arlequin version 3.5 (Excoffier and Lischer 2010). Critical significance levels for multiple testing were corrected using the sequential Bonferroni procedure (Rice 1989). In addition, we identified groups of populations using the spatial analysis of molecular variance version 2.0 (SAMOVA 2.0) (Dupanloup et al. 2002). To detect the number K of groups with the largest F_ct value, K was user-defined between 2 and 8, and 100 independent simulated annealing processes were performed in each run. However, collection sites with two or fewer individuals (the Kurashiki and Takahashi River systems) were excluded from the Φ_st and SAMOVA.

The divergence times of intraspecific lineage of Acheilognathus tabira tabira were estimated using BEAST ver. 2.7.6 (Bouckaert et al. 2019). We used sequences of A. t. jordani, A. cyanostigma Jordan & Fowler, 1903, and Opsariichthys platypus (Temminck & Schlegel, 1846) as outgroups. We adopted the estimated clock rate of cytb (0.76%, Zardoya and Doadrio 1999) that has been applied to the Acheilognathinae (Tominaga et al. 2020; Miyake et al. 2021). Divergence times were calibrated using the first appearance in the fossil record of Acheilognathinae from the early Miocene (ca. 20 Mya) found in Japan (Yasuno 1984). Constraints of the fossil were specified as a log-normal distribution, ranging from 16.5–23.0 Mya in the 95% range. Estimation was carried out using an optimised relaxed clock and applied the substitution model HKY+I selected by the BIC in ModelTest-NG. MCMC chains were run for 50000000 generations and sampled every 1000 generations, with the exclusion of the first 5000000 generations as burn-in. The convergence of MCMC was checked by calculating ESS values (> 200) using Tracer ver. 1.7.2 (Rambaut et al. 2018). TreeAnnotator ver. 2.7.6 in the BEAST package was used to obtain a maximum credibility tree with the annotation of average node ages and 95% highest posterior density (HPD) interval. The phylogenetic tree was visualized with FigTree ver. 1.4.4 (Rambaut 2018).

We sequenced 1069-bp mtDNA cytb gene nucleotide fragments from 140 individuals of A. t. tabira collected from 10 river systems and the captive population. As a result, 36 haplotypes were detected (T1–36), 5 of which [T1: AB759882 (Kiso River system), T3: AB759884 (Kiso R.), T4: AB620138, AB620159, and LC578851 (Lake Biwa, Kizu R., Kiso R., and northern Mie Prefecture), T23: AB759887 (Nagara R.), and T36: AB759888 (Nagara R.)] had been detected in previous studies (Kitamura et al. 2012; Umemura et al. 2012; Ito et al. 2021).

The topologies of the ML and BI phylogenetic trees were partially different (Fig. 2, Suppl. material 2). In the ML tree, Lineage I and III were supported as monophyletic with high values, respectively (SH-aLRT > 94%, UFBoot > 96%). The monophyletic clade of Lineage I and III was supported with high values (88.5%) by SH-aLRT, however with slightly lower value (88%) by UFBoot. In the BI tree, the monophyletic clade of Lineage I was supported by PP with a low value (0.54), while Lineage III was supported with a high value (0.99). The monophyletic clade of Lineage I and II was not supported, whereas the monophyletic clade of Lineage II and III was supported, albeit with a lower value (0.69). In addition, Lineage III was further subdivided into two sub-lineages: Sub-lineages III-i and ii in the ML and BI trees (SH-aLRT = 87.4%, UFBoot = 95%, and PP = 0.96, Fig. 2).

Figure 2.

Maximum likelihood (ML) tree of the 1069-bp cytochrome b gene sequences of Acheilognathus tabira tabira individuals from the Seto Inland Sea and Ise Bay regions. Numbers at nodes indicate Shimodaira-Hasegawa-like approximate likelihood ratio test values (left), ultrafast bootstrap values (middle) in the ML tree, and Bayesian posterior probabilities (right) in Bayesian inference tree. Each value is indicated when it exceeds 80%, 95%, and 0.80. Numbers in parentheses indicate the number of specimens. The parentheses after each Lineage name indicate the natural distribution area. The statistical parsimony network of A. t. tabira is shown to the left of the tree. Pie charts of Lineage I indicate the relative frequencies of haplotypes of the five groups defined by SAMOVA.

In Lineage I, Haplotypes T1–12 and T35 were mainly detected in Lake Biwa and the Yodo River system (Loc. 1–4), four of which (T3, T4, T11, and T12) were detected in the Yoshino River system (Loc. 12). Haplotypes L13–20 were detected only in the Muko River system (Loc. 5). Haplotypes T21 and T22 were only detected in the Kako River system (Loc. 6). Haplotype T23 was detected in the Kako, Yoshii Asahi, Sasagase, and Takahashi River systems (Loc. 6–9, and 11). Haplotypes T24–34 were detected in the Yoshii Asahi, Sasagase, Kurashiki, and Takahashi River systems (Loc. 7–11). Haplotype T36 in Lineage II was detected only in a captive Gifu World Freshwater Aquarium population collected from an unknown river system in Gifu Prefecture (Loc. 13).

In the statistical parsimony network, A. t. tabira exhibits a bottleneck pattern (Fig. 2). Lineage I was represented by a star-like pattern centered on the ancestral haplotype T4 (Fig. 2).

The results of the pairwise Φ_st among the local populations are shown in Table 2. The Φ_st values of Yodo1 and Muko populations showed significant genetic differentiation between populations, excluding Yodo3 and Biwa populations, respectively (Φ_st 0.235–392, 0.269–0.415; P < 0.05). The Φ_st value of the Yoshino population also showed significant genetic differentiation between Kako, Yoshii, Asahi, and Sasagase populations (Φ_st 0.191–0.301: P < 0.05). In addition, the Φ_st values of the Asahi and Sasagase populations showed significant genetic differentiation between Biwa, Yodo2, and Yodo3 populations (Φ_st 0.269–0.392: P < 0.05), and Kako population showed significant genetic differences from Yodo3 population (Φ_st 0.281–0.339, P < 0.05).

The results of the population group estimation using SAMOVA are shown in Table 3 and Suppl. material 3. The highest F_ct value (0.32166; P < 0.001) was obtained when the 10 populations were divided into K = 5 groups: Group A (Lake Biwa-Yodo River and Yoshino River: Loc. 1, 3, 4, and 12), Group B (tributary A of Yodo River: Loc. 2), Group C (Muko R. Loc. 5), Group D (Kako R.: Loc. 6) and Group E (Yoshii, Asahi, and Sasagase Rs: Loc. 7–9), and explained 32.17% of the variation among groups (P < 0.001), -0.3% of the variation among populations within groups (P > 0.1), and 68.13% of the variation within populations (P < 0.001).

We showed the divergence times of the three lineages of A. t. tabira in Fig. 3. In this tree, the exclusivity of Lineage I was supported (1.00), similar to the ML and BI trees described above, whereas the exclusivity of Lineages II and III was supported (0.98), similar to the topology of the BI tree. The estimated divergence time between Lineage I and II+III was approximately 1.53 Mya (95% HPD, 0.71–2.64 Mya; Fig. 3, node 1), while that between Lineage II and III was approximately 0.91 Mya (0.34–1.68 Mya; node 2). The time of the most recent common ancestor (tMRCA) of Lineage I was estimated approximately 0.96 Mya (0.44–1.66 Mya; Fig. 3, node 3).

Figure 3.

Divergence time estimation by Bayesian inference tree of the 1069-bp cytochrome b gene sequences of Acheilognathus tabira tabira and outgroups. The blue rectangular bars on the nodes indicate the 95% highest probability density. Bayesian posterior probabilities are indicated at nodes, with values exceeding 0.90 shown. The node marked with an asterisk indicates the calibration point based on fossil record for the Acheilognathinae. Nodes with circled numbers are referenced in the text.

We estimated the phylogenetic tree of A. t. tabira based on the sequence of the cytochrome b region of the mtDNA, and in the samples used in the present study, three lineages (Lineages I, II, and III) were identified primarily based on the ML tree. The results were similar to those of a previous study (Umemura et al. 2012), where Lineage I to III were referred to as the Kinki-Sanyo Group, Nobi Plain Group I, and Nobi Plain Group II, respectively. The populations newly analyzed in the present study (Muko, Kako, Sasagase, Takahashi, Kurashiki, and Asahikawa River systems) were all included in Lineage I, and the captive Gifu World Freshwater Aquarium population collected from Gifu Prefecture was included in Lineage II. In a previous study, a haplotype of Lineage II was detected only in an individual collected from a tributary of the Nagara River system in the previous study (Umemura et al. 2012). Lineages I (non-native lineage) and III (native lineage) were identified in this tributary; therefore, it was unclear whether Lineage II was native to the Ise Bay region or non-native to the Seto Inland Sea region (Umemura et al. 2012). In the present study, the haplotype belonging to Lineage II was not identified in any water system in the Seto Inland Sea region; therefore, the natural distribution area of Lineage II was considered to be the Ise Bay region.

The populations of many freshwater fishes [e.g., Sarcocheilichthys variegatus variegatus (Temminck & Schlegel, 1846) and Opsariichthys platypus] in the Ise Bay region are thought to have been divided from the populations of the Seto Inland Sea region by the uplift of the Suzuka Mountains approximately one million years ago (Mya) (Watanabe et al. 2017). The genetic lineages of populations in many freshwater fishes in each region of Ise Bay and the Seto Inland Sea are generally exclusive (e.g., S. v. variegatus, Komiya et al. 2014; O. platypus, Kitanishi et al. 2016; Pseudogobio esocinus (Temminck & Schlegel, 1846) and P. agathonectris Tominaga & Kawase, 2019, Tominaga et al. 2016; Tanakia lanceolata (Temminck & Schlegel, 1846), Tominaga et al. 2020). However, in the case of A. t. tabira, the exclusivity of the two native lineages in the Ise Bay region (lineages II and III) was not supported in the ML tree and was supported with low values in the BI tree. In the ML tree, these two lineages were shown as paraphyletic groups, and each branch was highly supported. Furthermore, the results of the statistical parsimony network also indicated that lineages II and III were closely related to different haplotypes of Lineage I. On the other hand, the BI tree does not have the same topology as the ML tree. Divergence time estimates suggest that the supported topology is similar to that of the BI tree, with Lineages II and III diverging from Lineage I approximately 1.53 Mya, and Lineages II and III diverging approximately 0.91 Mya. However, due to the varying exclusivity of Lineages II and III across different phylogenetic trees, accepting this divergence time estimation as it is would be risky. And more, the original distribution range of Lineage II and III is unclear, as A. t. tabira in the Ise Bay region has already become extinct in many river systems (Mukai 2019). If DNA analysis of specimens collected decades ago and stored in museums becomes possible, it will be possible to verify this issue. To consider the divergence order among the three lineages, further studies using longer sequences, such as mitogenomes, including sequence data from historical specimens, are necessary to re-estimate divergence times.

In Lineage I, which was detected only in the Seto Inland Sea region, genetic differentiation has not been recognized in previous studies because of the small number of sampling sites and individuals (Kitamura et al. 2012). In the present study, we greatly increased the number of sampling sites and individuals; as a result, five genetic groups (A–E) were distinguished within Lineage I using SAMOVA. These five groups comprised adjacent river systems, indicating that the genetic population of A. t. tabira within the Seto Inland Sea region was genetically differentiated into narrow regions. The river systems flowing into the Seto Inland Sea are thought to have connected as a single paleo-river system during the glacial periods of the Pleistocene (ca. 0.01–2.5 Mya) and were isolated during the interglacial periods (Kuwashiro 1959; Ota et al. 2004). Additionally, the uplift of the mountains areas surrounding the Seto Inland Sea is thought to have become active since the Pleistocene (Ota et al. 2004), making the migration of freshwater fish between river systems difficult during this period. The tMRCA of Lineage I is estimated to be 0.96 Mya (95% HPD, 0.44–1.66 Mya), overlapping with the periods of connection and isolation of the paleo-river systems and the uplift of mountains around the Seto Inland Sea. Therefore, the isolation factors between the regional groups are suggested to be related to the paleo-river systems and active uplift of mountains that occurred during the Pleistocene. Genetic differentiation in the same period due to similar factors in the Seto Inland Sea region has also been suggested in P. esocinus (Tominaga et al. 2016).

A unique genetic group (Group D) was identified in the Kako River system. However, Φ_st showed no significant genetic differentiation (P > 0.05) between the Kako River system and the other three river systems (Yoshii, Asahi, and Sasagase) included in Group E. SAMOVA results indicated that most of the genetic variation in this subspecies was within populations (68.13%) and that differentiation among groups was relatively small (32.17%). Genome-wide analysis of nuclear DNA may be useful for more detailed elucidation of the genetic population structure of A. t. tabira.

Populations collected from the Lake Biwa-Yodo and Yoshino River systems were included in Group A. In addition, four haplotypes detected in the Yoshino River system were similar to those in the Lake Biwa-Yodo River system. This study demonstrates that populations in the Seto Inland Sea region are genetically differentiated by localized areas. The reason for this is thought to be the same as with other species: the disappearance of the paleo-river system and isolation due to the uplift of mountains. Therefore, it is unlikely that the population in the Yoshino River system has the same haplotype as the population in the Lake Biwa–Yodo River system, which is across the Seto Inland Sea. In the Yoshino River system, non-native freshwater fishes [e.g., Acheilognathus cyanostigma and Acheilognathus rhombeus (Temminck & Schlegel, 1846)] were estimated to have been artificially introduced from Lake Biwa (Hosoya 2019; Miyake et al. 2021). Ministry of the Environment of Japan (2015) also indicated that A. t. tabira populations were artificially introduced into the Yoshino River system. Therefore, the finding that only haplotypes of the Yoshino River system are common to those of the Lake Biwa-Yodo River system supports the hypothesis that the Yoshino River system population was artificially introduced from the Lake Biwa-Yodo River system.

The captive population of the Gifu World Freshwater Aquarium was identified as Lineage II, which is thought to be native to the Ise Bay region. Non-native populations belonging to Lineage I have been artificially introduced into all habitats of native populations in the Ise Bay region (Umemura et al. 2012; Kitamura and Uchiyama 2020). Therefore, the captive population of the Gifu World Freshwater Aquarium may be a native population of the Ise Bay region that has not undergone genetic introgression. However, to confirm that the population is not genetically introgressed, examining the possibility of hybridization with Lineage I using nuclear DNA is necessary.

Conservation units need to focus on levels below species (Moritz 1994, Frankham et al. 2010). Evolutionarily significant units (ESUs) refer to phylogenetically unique intraspecific population groups, established by factors such as mitochondrial DNA monophyly (Moritz 1994). In the case of A. t. tabira, the exclusivity of Lineages I–III was supported; thus, we propose to designate them as ESUs.

Furthermore, Management Units (MUs) are established based on allele frequencies among populations (Moritz 1994, Frankham et al. 2010). In the case of A. t. tabira, Groups A–E found within Lineage I were significantly genetically differentiated by SAMOVA. Therefore, we propose that these five groups be conserved as MUs of A. t. tabira in the Seto Inland Sea region. However, it has been proposed that adaptive traits should also be taken into account in the establishment of ESUs and MUs (Crandall et al. 2000). In the future, in order to establish better conservation units, it is necessary to study adaptive traits among populations of A. t. tabira.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This study was supported by JSPS KAKENHI (22K14908).

Gen Ito: Conceptualization, Data curation, Formal Analysis, Funding acquisition, and Writing – original draft. Naoto Koyama, Ryota Noguchi, Ryoichi Tabata, Seigo Kawase, and Jyun-ichi Kitamura: Investigation, Resources, and Writing – review & editing. Yasunori Koya: Supervision, Funding acquisition, and Writing – review & editing.

Gen Ito https://orcid.org/0000-0002-9781-7206

All of the data that support the findings of this study are available in the main text or Supplementary Information.