Phylogeny and historical demography of endemic fishes in Lake Biwa: the ancient lake as a promoter of evolution and diversification of freshwater fishes in western Japan

Fish specimens

The mtDNA sequences were obtained from a total of 190 specimens of 44 species/subspecies, including 11 of the 16 endemic/semiendemic species of Lake Biwa and their related riverine species for phylogenetic analyses, and from a total of 2319 specimens of 42 species (including 14 endemic/semiendemic species) occurring in the lake for population genetic analyses (Table 1, Fig. 2, Tables S1 and S2). Furthermore, to reveal the relationships among intraspecific clades through phylogeographic analysis, we used 497 specimens of nine species from other regions, including some data from previous studies (Table S3). Five endemic species/subspecies were not used in the phylogenetic analysis because they or their relatives exhibit polyploidy (e.g., the cobitids Cobitis magnostriata and Cobitis minamorii oumiensis (see Saitoh 2000; Kitagawa et al. 2005); and the cyprinid Carassius buergeri grandoculis (see Mishina et al. 2014)), known mtDNA introgression (the gobiid Rhinogobius sp. BW (sensu Akihito et al. 2013; Yamasaki et al. 2015)), or the closest relative was uncertain (the cyprinid Ischikauia steenackeri; Takeuchi and Hosoya 2011).

mtDNA sequencing

Total genomic DNA was isolated from fin clips preserved in 100% ethanol using a Genomic DNA Purification kit (Promega, Tokyo, Japan). The mitochondrial nucleotide sequences of the 16S ribosomal RNA (16S rRNA), cytochrome c oxidase I (CO1), NADH dehydrogenase subunit 5 (ND5), cytochrome b (cyt b), and control region (CR) genes were chosen as molecular markers to infer phylogenetic relationships and divergence times because large amounts of data and information have been accumulated for these genes with respect to molecular evolution in fish (e.g., Briolay et al. 1998; Zardoya and Doadrio 1999; Watanabe and Takahashi 2010). These five different regions of mtDNA have different evolutionary rates (Meyer 1994; Avise 2000; Kuo et al. 2003); hence, we expected these genes to enable a comprehensive estimation of divergence times over relatively short or long time scales. In population genetic analyses, we used the cyt b gene (or the CR gene, see below) because of its higher substitution rate (see “Results”).

The primers described in Table S4 were used for polymerase chain reaction (PCR) amplification of the five mitochondrial regions. The PCR conditions consisted of 30 cycles of denaturation (94°C, 15 s), annealing (48°C to 55°C for each region, 15 s), and extension (72°C, 30 s) in a PC808 thermal cycler (Astec, Fukuoka, Japan). After purifying the PCR products by treatment with ExoSAP‐IT (USB Corp., Cleveland, OH, USA) at 37°C, they were sequenced using an ABI 3130xl DNA sequencer (Applied Biosystems, Foster City, CA, USA) with amplification primers using the BigDye Terminator Cycle Sequencing FS Ready Reaction kit ver. 3.1 (Applied Biosystems).

The sequences of the 5′‐halves of the 16S rRNA gene (~1200 bp), CO1 gene (~600 bp), and ND5 gene (~1000 bp); the nearly complete cyt b gene (~1100 bp; or 580–750 bp for population genetic analyses); and the nearly complete (~900 bp) or 5′‐half of the CR gene (450 bp for silurids) were deposited in DDBJ/EMBL/GenBank (accession nos. LC097321–098734; Tables S1–S3). The haplotype frequencies for populations are available from GEDIMAP (http://gedimap.zool.kyoto-u.ac.jp/; Watanabe et al. 2010a) under accession numbers P2000–2195. For population genetic analyses, the sequence data of some species were drawn from Harada et al. (2002), Tominaga et al. (2009, 2015), Watanabe et al. (2010b), Komiya et al. (2011), Kakioka et al. (2013), Tabata and Watanabe (2013), and Mishina et al. (2014).

Phylogenetic reconstruction and estimation of divergence time

For phylogenetic reconstruction and estimation of the divergence times of endemic fishes in Lake Biwa, we used the mtDNA sequence data of 11 endemic species/subspecies, their closely related species, and outgroups, including some data from previous studies. Datasets were created for (1) Salmonidae, (2) Cyprininae (Cyprinidae), (3) Gobioninae (Cyprinidae), (4) Oxygastrinae (Cyprinidae), (5) Siluridae, (6) Gobiidae, and (7) Cottidae (Table S1).

The sequences of the five mitochondrial regions were aligned using the Muscle (Edgar 2004) and Mafft (Katoh et al. 2002) software. These alignments were automatically trimmed with trimAl (Capella‐Gutiérrez et al. 2009) to remove poorly aligned regions with gaps in more than 20% of the sequences or with a similarity score lower than 0.001, unless this removed more than 40% of the columns in the original alignment before the sequences were concatenated. For each dataset, the Bayesian approach was applied to estimate the phylogenetic trees, divergence times, and molecular evolutionary rates for the lineages of the endemic forms, with the fittest models selected by the Bayesian information criterion (BIC) in jModelTest v2.1.1 (Table S5; Darriba et al. 2012) and Yule (speciation) tree prior using BEAST v1.6.2 (Drummond and Rambaut 2007). We adopted the random local clock model, which assumes one or more independent rates on different branches (Drummond and Suchard 2010). To estimate the age of divergence between lineages, we imposed some a priori constrains on the node ages, mostly from geological events (such as upheavals of mountains) that divided regional populations of freshwater fish (Table 2). No calibration point was available for Gymnogobius isaza because of the lack of fossil data and the unclear geographic population structure of its closely related species. Thus, the previously estimated evolutionary rate of the G. isaza cyt b gene (Tabata and Watanabe 2013) was adopted (Table 2).

All other model parameters used default priors. For each Markov chain Monte Carlo (MCMC) analysis, we performed two independent runs of 100 million generations, sampling every 1000th generation and removing the initial 10% of samples as burn‐in. The convergence of the chains to the stationary distribution was confirmed by inspection of MCMC samples using Tracer v1.5 (Rambaut and Drummond 2009). These two runs were combined, and the consensus tree was calculated using TreeAnnotator v1.6.2 in the BEAST package and visualized using FigTree v1.3.1 (Rambaut 2009).

Population genetic and historical demographic analyses

For population genetic analyses, we used the cyt b sequences of a total of 2319 specimens from 42 species (11–249 specimens/species) occurring in Lake Biwa, including 11 endemic forms (Table S2). Because Silurus lithophilus exhibited no genetic diversity in the cyt b gene, we used the CR gene, which of silurids showed a faster evolutionary rate and genetic polymorphism (see “Results”), for analyses of S. lithophilus and Silurus asotus. We treated Sarcocheilichthys biwaensis and S. variegates microoculus as the species complex “Sarcocheilichthys spp.” in the population analysis because no genetic difference was detected between them in the mtDNA (cyt b) and microsatellite markers (Komiya et al. 2011).

The genetic diversity of endemic fish was evaluated based on haplotype diversity (h) and nucleotide diversity (π) using the Arlequin 3.5 software (Excoffier and Lischer 2010). The relationships among mtDNA haplotypes of each species were estimated as a statistical parsimony network using the TCS 1.2.1 software (Clement et al. 2000). When we found distinctive multiple intraspecific mitochondrial clades in Lake Biwa, we treated those clades both together and separately in the following analyses, and added the data of specimens from other regions to examine the origins of the clades (Table S3).

To evaluate whether sequences had evolved under strict neutrality, Fu's F _S (Fu 1997) was calculated using Arlequin. Fu (1997) noticed that the 2% percentile of the distribution corresponded to the 5% cutoff value. Therefore, when the statistic of Fu's F _S became negative at the significance level of 2%, population expansion or selection was suggested. To examine the possible occurrence of a historical demographic expansion of each species, we performed mismatch distribution analysis (Rogers and Harpending 1992), which was implemented in Arlequin. The fitness of the observed population data to a model of sudden population expansion (Rogers and Harpending 1992) was tested by calculating the sum of the squared deviation (SSD) and raggedness (rg) statistics (Harpending 1994). For each statistic, the sudden expansion model was rejected when P < 0.05 (Excoffier and Lischer 2010). The parameter τ is related to the time of population expansion by the formula τ = 2Tu, where T measures time in years when u is the mutation rate per nucleotide per year, multiplied by the sequence length (Slatkin and Hudson 1991). We estimated the time of population expansion with the rate (u) obtained from the Bayesian phylogenetic analysis based on the random local clock model mentioned above or the rate reported in a previous study (Zardoya and Doadrio 1999; Kawamura et al. 2014; see Table S5).

The demographic history of the fishes in Lake Biwa, including endemics, was also examined by Bayesian skyline plot (BSP) analysis (Drummond et al. 2005), implemented in BEAST. We performed runs with a MCMC chain length of 50 million generations, sampling every 1000th generation and removing the initial 10% of samples as burn‐in. The substitution model used was the model selected by BIC in jModelTest, and the time to expansion was estimated using the mutation rate obtained in the above phylogenetic analysis. The BSP result with the stepwise (constant) model was summarized using Tracer. When we found two or more distinct intraspecific mitochondrial clades in the TCS analysis, we carried out mismatch distribution and BSP analyses for each intraspecific clade separately.

Divergence time estimation based on mtDNA genes

The evolutionary rates (nucleotide substitutions per million years per lineage) of the respective genes at the endemic species lineages were estimated to be 0.18–0.41% (16S rRNA), 0.37–1.00% (CO1), 0.44–1.82% (ND5), 0.53–1.89% (cyt b), and 0.27–1.57% (CR) (Table S5). In most species, the cyt b or ND5 genes showed the highest average molecular evolution rate, whereas the 16S rRNA gene exhibited the lowest rate. Exceptionally, the rate of CR was the highest in the Silurus species (Table S5).

The degree of genetic divergence between the endemic species in Lake Biwa and their closely related species varied highly between species. Some species showed relatively large genetic distances (uncorrected P): 3.9% (between Opsariichthys uncirostris and Opsariichthys bidens) to 6.5% (between Silurus species). Other species, such as Squalidus spp. and Sarcocheilichthys spp., exhibited small genetic distances (<1%) (Table 3). These reflected the wide range of estimated divergence times of the endemic lineages (Table 3; Figs. 3 and 4).

Dividing the endemic species based on a divergence time of ~0.4 Myr (0.3–0.5 Myr), when the present Lake Biwa started to develop, the estimation showed that seven forms diverged prior to the development of the present lake, whereas the remaining four diverged more recently. Of the former group, the divergence time of Silurus biwaensis from the common ancestor of S. asotus and S. lithophilus was the oldest (13 Myr), and these species formed a monophyletic group that includes Eurasian continental species. The divergence time of S. lithophilus from S. asotus in western Japan was the second oldest (9 Myr), but the former did not show a sister relationship to the latter; however, it did show such a relationship to the cryptic clade of S. asotus occurring in the eastern area (1.24 Myr, calibration point C1). The divergence times of Opsariichthys uncirostris (from O. bidens), Gnathopogon caerulescens (from the west clade of G. elongatus), Carassius cuvieri (from Carassius species in China and Japan), and Gymnogobius isaza (from the common ancestor of G. urotaenia and G. petschiliensis) were estimated at ~2–3 Myr The divergence time of Cottus reinii lacustrine type (from the Ise Bay area population of Cottus reinii amphidromous type) was estimated at 1.3 Myr.

On the other hand, the other four endemic forms, Squalidus biwae biwae, Sarcocheilichthys biwaensis, S. variegates microoculus, and Oncorhynchus masou subsp. “Lake Biwa trout,” were estimated to have diverged from their closely related species more recently (Table 3, Figs. 3 and 4). Except for the Lake Biwa trout, the other three forms showed no distinct genetic divergence in their mtDNA from their closely related species.

Population genetic and historical demographic analyses

The endemic species showed a wide range of the haplotype diversity (h; 0.000–0.9358, average ± SD = 0.724 ± 0.244 in cyt b or CR; see “Materials and Methods”). Among 14 endemic fishes, 10 had h values higher than 0.7. By contrast, Ischikauia steenackeri had no genetic variation in its mtDNA cyt b gene (Table 1). The haplotype networks of the respective fishes in Lake Biwa, including endemics, showed various patterns including star‐like, bush‐like, and dumbbell‐like patterns (Figs. 5 and S2). The cyt b haplotypes of Sarcocheilichthys spp. and Gymnogobius isaza were clearly divided into three and two clades, respectively. Hereafter, the clades of Sarcocheilichthys spp. are treated as “clade A” (main clade), “clade B,” and “clade C;” and those of Gymnogobius isaza are treated as “clade A” (main clade) and “clade B” (Figs. 5 and S2).

Fu's F _S suggested that the null hypothesis of neutrality was rejected in nine of 14 endemic forms (i.e., Oncorhynchus masou subsp., Cobitis minamorii oumiensis, Opsariichthys uncirostris, Squalidus biwae biwae, Gnathopogon caerulescens, Sarcocheilichthys spp., Rhinogobius sp. BW, Gymnogobius isaza, and Cottus reinii (lacustrine type)) and three intraspecific clades (i.e., clade A of Sarcocheilichthys spp. and clades A and B of G. isaza) (Table 4). In the mismatch distribution analysis, these endemic species, with the exception of Sarcocheilichthys spp. and G. isaza, fit the model of sudden expansion and exhibited unimodal distributions (Table 4, Figs. 5 and S2). The mismatch distributions for Sarcocheilichthys spp., Carassius buergeri grandoculis, C. cuvieri, and G. isaza were multimodal (Figs. 5 and S2). Using the evolutionary rate estimated by the above phylogenetic analysis, the population expansion times of the respective endemic forms and intraspecific clades were estimated at 0.03–0.22 Ma (average ± SD = 0.13 ± 0.06; Table 4, Fig. 6) (see Table S5).

On the other hand, in 12 of 27 nonendemic forms, Fu's F _S suggested that the null hypothesis of neutrality was rejected (P < 0.02). The mismatch distribution of those 12 forms fit the model of sudden expansion and, with the exception of Nipponocypris sieboldii, which showed a clearly bimodal distribution, all exhibited a unimodal distribution (Table 4, Figs. 5 and S2). The population expansion times of these 11 nonendemic forms were estimated at 0.01–0.35 Ma (average ± SD = 0.14 ± 0.09; Table 4, Fig. 6). For the other 15 nonendemic forms, the null hypothesis of neutrality was not rejected by Fu's F _S (P > 0.02) or the model of sudden expansion was not supported by the mismatch distribution analysis (P < 0.05).

The BSP analysis also revealed the population expansion pattern in most of the endemic forms, which began ca. 0.4 Ma or later (Figs. 6 and S2). However, some differences in the pattern were observed among the forms; that is, the forms were roughly classified into those with clear rapid expansion, those with calm expansion, and those with intermediate expansion. Clear rapid expansion was typically shown in clade A of Gymnogobius isaza and clade A of Sarcocheilichthys spp., whereas the expansion pattern was relatively calm in Squalidus biwae biwae (Figs. 5 and S2A). The species in which the null hypothesis of neutrality was not rejected by Fu's F _S exhibited an almost stable population size in BSP (e.g., Cobitis magnostriata, Carassius cuvieri, and Silurus lithophilus; Fig. 5B). The times to the most recent common ancestor (tMRCAs) of most endemic forms and five intraspecific clades (i.e., three clades of Sarcocheilichthys spp. and two clades of G. isaza) were generally shorter than 0.4 Myr (average ± SD = 0.15 ± 0.09, range 0.03–0.32, n = 13; Table 4, Figs. 6 and S2). On the other hand, the tMRCAs of the species that included two or more distinct intraspecific clades (i.e., C. buergeri grandoculis, C. cuvieri, Sarcocheilichthys spp., and G. isaza) were older than 0.4 Myr (average ± SD = 0.88 ± 0.67, range 0.40–1.88, n = 4; Table 4, Figs. 6 and S2).

The BSP analysis of 11 nonendemics also suggested population expansions that began c. 0.4 Ma or later (Figs. 6 and S2). The tMRCAs of these nonendemic fishes, except for Tanakia lanceolata, were mostly shorter than 0.4 Myr (Table 4, Figs. 6 and S2). The nine nonendemic forms that did not show the population expansion pattern included distinct intraspecific clades (e.g., Tribolodon hakonensis, Silurus asotus), and the tMRCAs of those species, as with some endemic forms, were longer than 0.4 Myr (Table 4, Figs. 6 and S2).

In several species groups, haplotypes common or very close to those found in Lake Biwa were detected from other riverine populations in western Japan (Fig. S3), suggesting multiple colonizations into Lake Biwa by some intraspecific lineages. For example, both clade A and B haplotypes of Sarcocheilichthys spp. were obtained from the Lake Biwa, Kinki, and Chugoku regions, and clade C haplotypes were obtained from the Lake Biwa and Kyushu regions. Common or very close haplotypes of Carassius buergeri grandoculis (Lake Biwa endemics) were shared by C. buergeri subspp. populations from the Kinki and Kyushu regions of western Japan. Lake Biwa populations of seven nonendemic species (i.e., Tanakia limbata, G. elongatus, Pseudogobio esocinus, Rhynchocypris oxycephalus jouyi, Tribolodon hakonensis, Nipponocypris sieboldii, and Silurus asotus) had haplotypes of two or more clades that were shared with the Chugoku, Kinki, and Tokai regions of western Japan.

Divergence times and origins of endemic fish in Lake Biwa

Based on a Bayesian approach with geological calibrations, the evolutionary rates of each mitochondrial gene were estimated to be ~0.3–1.5% Myr⁻¹ lineage⁻¹. No discrepancies were found between our rates and those previously reported for each region (e.g., Avise 2000; Kuo et al. 2003; Watanabe and Takahashi 2010). Ho et al. (2005) pointed out an apparent acceleration of the divergence rate over a short time period (<1–2 Myr or ~1% substitution). However, except for a few cases, the divergences of endemic forms from their closest relatives exhibited greater genetic distances (>1% substitution in the uncorrected P‐distance) than would occur with such a threshold. Additionally, the estimated rates were within in the range of known molecular clocks of teleost mitochondrial genes (~0.1–2.5% Myr⁻¹ lineage⁻¹; e.g., Burridge et al. 2008; Watanabe and Takahashi 2010). Therefore, it is reasonable to infer divergence times using these estimated evolutionary rates. Nevertheless, it should be noted that the estimations could be biased toward older divergence times when these long‐term evolutionary rates are used in the analysis of a recent historical demography (Ho et al. 2005) (see the next section).

Bayesian phylogenetic analysis using five mitochondrial genes suggested that seven of the 11 examined endemic species diverged from their closest species long before the origin of the present Lake Biwa environment (~0.4 Ma) (hereafter the “old‐endemic” group), whereas four other species diverged ca. 0.4 Ma (or did not clearly diverge from their closest relatives; hereafter the “new‐endemic” group). Previous studies presumed that endemic fish that had been inferred to be adapted to the present Lake Biwa environment (the rocky substrate, extensive deep pelagic area, etc.) differentiated from their ancestors after the establishment of the particular environments of the lake (e.g., Nakajima 1987; Takahashi 1989). However, many of the endemic species were inferred to have diverged from their closest species long before 0.4 Ma.

It has been inferred that the new‐endemics diverged recently from or maintained gene flow with their closest relatives. No remarkable genetic differentiation was found among Sarcocheilichthys biwaensis, S. variegates microoculus (endemic to Lake Biwa), and S. variegates variegates (distributed widely in rivers in western Japan but not in Lake Biwa) in their mtDNA (Komiya et al. 2011; 2013). The head shape divergence and the rapid speciation in this group in Lake Biwa were inferred to be associated with feeding habits corresponding to substrate types in the lake (i.e., rocky vs. pebbly areas; Komiya et al. 2011; 2013). Squalidus biwae biwae (endemic to Lake Biwa) and its closest form Squalidus biwae tsuchigae (distributed in the southern basin of Lake Biwa, the rivers flowing into the lake, and around the Seto Inland Sea and Ise Bay regions) exhibit steep but continuous genetic and morphological variations along the environmental gradient around the lake (Kakioka 2013). This could be explained by adaptations in body shape to pelagic life for efficient swimming performance, along with some gene flow. On the other hand, it has been inferred that the origin of Oncorhynchus masou subsp. “Lake Biwa trout” diverged a little earlier than those of the above groups (0.52 Ma, 95%HPD: 0.29–0.76 Ma), and a sister relationship between this trout and all other masu trout subspecies (i.e., Oncorhynchus masou masou, Oncorhynchus masou ishikawae, and Oncorhynchus masou formosanum) has been amply supported. The masu trout group basically has an anadromous life history. Our results support that the Lake Biwa trout was differentiated from the ancestral anadromous or landlocked form by its use of the lake instead of the sea after the pelagic zone of the present Lake Biwa began to form (Oohara and Okazaki 1996; Fujioka 2009).

The old‐endemic group was inferred to have originated from lineages that had already diverged from any other existing lineages much earlier than the present Lake Biwa stage (i.e., prior to 0.4 Ma). It has been suggested that those endemic forms exhibit apparent adaptations to the lake environment in terms of their morphologies, behaviors, and life histories (Tomoda 1978; Takahashi 1989). This suggests that those endemic forms originated from lineages different from those of the present riverine closest relatives and that they colonized and adapted to the present lake following the extinction of their direct riverine ancestral populations. This pattern differs from the rapid evolution following adaptation to new environments observed in the cichlids of Lake Victoria (Seehausen 2006; Turner 2007) or the sticklebacks and salmonids (Coregonus and Salvelinus) of the postglacial lakes of North America and northern Europe (Schluter 1995; Bernatchez et al. 1999; Gíslason et al. 1999).

The older origins of Gymnogobius isaza, Gnathopogon caerulescens, and Carassius cuvieri have been suggested by previous molecular phylogenetic studies based on shorter mtDNA sequences (Harada et al. 2002; Takada et al. 2010; Kakioka et al. 2013, 2015; Tabata and Watanabe 2013). The present study supports these results and demonstrates that Silurus lithophilus, which presumably has traits that are adaptive in rocky areas (e.g., anteriorly positioned eyes and yellowish brown body color), also has an older origin, and its sister group is a cryptic lineage that has been included in Silurus asotus (see Table 3, Fig. 3). Cottus reinii (lacustrine type) (a cottoid, landlocked and presumably adapted to the bottom of the deep area) and Opsariichthys. uncirostris (a cyprinid, regarded as a “relict species” because of the absence of extant relative species in Japan) were also inferred to have diverged from their closest forms before the establishment of the present Lake Biwa. Silurus biwaensis was confirmed to have distant relationships with the two other silurids and is considered a relict. Although the endemic goby Rhinogobius sp. BW was excluded from our analysis because of suspected mtDNA introgression, Yamasaki et al. (2015) showed that the species is also included in the old‐endemic group (divergence ca. 0.9–1.1 Ma from its closest relative) using multiple nuclear genes.

The divergence times of the old‐endemic lineages from their closest lineages ranged from about 1.3 Myr (C. reinii lacustrine type–a lineage of C. reinii amphidromous type) to 3.0 Myr (O. uncirostris–O. bidens), which are mainly in the early Pleistocene (Fig. 4). One possible explanation is the Pleistocene climate changes accompanying global cooling and drying in the glacial periods (Provan and Bennett 2008), which could affect the divergence and speciation of the endemic species found in the present Lake Biwa. Although western Japan has never been covered with ice sheets, such climatic changes caused remarkable changes in the freshwater environment and may have led to the extinction of species (Lindberg 1972; Clark et al. 1999; Provan and Bennett 2008). In addition, such environmental changes, together with drastic crustal movements (e.g., uplifting of mountains), would have promoted the population isolation and differentiation as well as the adaptation to novel environments of freshwater fish.

The possible important roles played by Paleo‐Lake Biwa in the divergence and adaptation of endemic species should be taken into consideration. Paleo‐Lake Biwa is the historical series of lake or marsh environments that existed from c. 4 Ma (the Pliocene) to 0.4 Ma, prior to the present Lake Biwa (Yokoyama 1984; Kawabe 1989, 1994). After the Pliocene, the topographical and limnological features of Paleo‐Lake Biwa (and Lake Tokai, located east of Paleo‐Lake Biwa) successively changed: the paleo‐lake experienced (1) shallow‐and‐large stages (“Lake Ohyamada,” >3.2 Ma; “Lake Ayama,” 3.0–2.7 Ma); (2) a deep‐and‐large stage (“Lake Koka,” 2.7–2.5 Ma), similar to the northern basin of the present Lake Biwa; (3) a second shallow‐and‐large stage (“Lake Gamo,” 2.5–1.8 Ma); (4) a swamp‐like stage (1.8–1.0 Ma); and (5) a final shallow‐and‐large stage (“Lake Katata,” 1.0–0.5 Ma) (Kawabe 1989, 1994). Past lakes and river systems have contributed to the present species richness and distribution in some cases. The foundation of the diversity of present cichlids, generated in Paleo‐Lake Makgadikgadi in Africa (Joyce et al. 2005), and the distribution of Barbus in the Balkan Peninsula (Marková et al. 2010) have been cited as examples. Although the confidence intervals of estimated divergence times are wide, the divergences of the old‐endemic lineages were likely concentrated between 2.0 and 3.0 Ma, that is, during the “Lake Koka” stage. It might be that Lake Koka, which existed as a deep‐and‐large lake, induced the divergence with adaptation of the ancestral populations of the extant species to the limnetic environment (Harada et al. 2002). During the subsequent desiccated periods, the alleles that had adapted to the paleo‐lake might have been retained as standing genetic variations in the populations that survived in rivers or marshes, as was the case with the cichlids in Lake Malawi and Lake Victoria (Elmer et al. 2009; Loh et al. 2013).

The mtDNA loci consist of a single linkage group; hence, they might cause a biased estimation of phylogenetic relationships and divergence times due to discrepancies between mtDNA and nuclear DNA trees and/or between gene trees and species trees (Avise 2000; Edwards and Beerli 2000; Walstrom et al. 2012). Our estimations, therefore, must be tested with a multilocus nuclear DNA dataset. Furthermore, studies on genetic backgrounds and the origins of adaptive genetic variations to lake environments using evolutionary genomic approaches are necessary (Gilad et al. 2009; Schoville et al. 2012; Bernardi 2013).

Historical demography and establishment of the ichthyofauna in the present lake

The demographic history inferred based on mtDNA diversity suggests that most of the present fish species in Lake Biwa experienced some historical bottleneck or founder effect after the present Lake Biwa environment began to develop. That is, the mismatch distribution and BSP analyses revealed expansion patterns of the populations (more precisely, of the mitochondrial effective population sizes) after 0.3 Ma (predominantly after 0.15 Ma) in more than half the analyzed fish species in Lake Biwa, including the old‐ and new‐endemics. In addition, the tMRCAs of almost all of the examined fish populations were estimated to be <0.4 Myr. These figures might be biased toward earlier dates because the molecular evolutionary rates used were calibrated by older events and hence they might be too small to apply to recent population histories (Ho et al. 2005). Even in such cases, our results regarding the expansion in recent geological periods are still supported. These demographic patterns suggest drastic changes and the establishment of the ichthyofauna during the present Lake Biwa stage.

The new‐endemic species in Lake Biwa are inferred to have expanded their populations following divergences from their closest species (or subspecies). The old‐endemics, on the other hand, had already diverged from their closest relatives long before the present Lake Biwa stage (prior to 0.4 Ma) and had colonized present Lake Biwa. The latter species might have acquired their presently observed phenotypic characteristics by adapting to the lacustrine environment after the colonization of their ancestral populations. That is, there would have been a time lag between the phylogenetic divergences of the endemic lineages and the acquisition (or reacquisition) of their present adaptive traits. Such a time lag might be explained by rapid phenotypic acquisition (or reacquisition) through strong selective pressure when the fish entered the novel limnetic environments.

Alternatively, the reduction and rapid expansion of the mtDNA effective population size might not reflect the dynamics of actual population size but instead might have resulted from a selective sweep followed by the recovery of genetic diversity in mtDNA. The mitochondrion is an important organelle that bears genes related to aerobic respiration and energy production that are involved in exercise and metabolism. Therefore, the mtDNA genome (and nuclear genes related to mitochondrial functions) can be a target of natural selection (e.g., Doiron et al. 2002; Hurst and Jiggins 2005; Ballard and Melvin 2010; Silva et al. 2014); i.e., natural selection might have favored improvements in exercise and metabolic capacity, when fish colonized the large pelagic area or the hypoxic, low‐temperature bottom area in the present Lake Biwa. Selection for genetic variation of a single gene in a mitochondrial genome could cause a selective sweep of an entire mitochondrial genome due to the lack of recombination; this could manifest as a pattern similar to a population bottleneck. However, selective sweep alone would not be able to explain the entire expansion pattern because the pattern was observed not only in endemic species but also in other species around the lake. The historical population demography of Lake Biwa fish should be reexamined using multiple, highly variable nuclear DNA loci to differentiate selective sweep caused by selection on a single locus.

The roles of Lake Biwa in the formation of ichthyofauna in western Japan

Lake Biwa is presumed to have functioned not only as the center of adaptive evolution and speciation but also as an important refugium (or “reservoir”) for freshwater fish in western Japan since the middle Pleistocene because of its huge capacity and historical stability. In some endemic and non‐endemic species, common or very close haplotypes were shared between Lake Biwa and rivers/lakes in other regions of western Japan. This suggests some historical gene flow between populations in these regions as well as secondary contact following multiple colonizations in Lake Biwa (Fig. S3). This situation differs from the cases in other typical ancient lakes, such as the African Great Lakes and Lake Baikal, with much higher rates of endemic species (>80% vs. <25% in Lake Biwa). The divergent mtDNA lineages within the endemic population of Gymnogobius isaza also support the secondary contact in Lake Biwa between past divergent populations (Tabata and Watanabe 2013; present study). Global climate changes (e.g. Clark et al. 1999; Provan and Bennett 2008) and drastic changes in inland water environments in the Pleistocene (Bernatchez and Wilson 1998; Yokoyama et al. 2000) would have resulted in the reduction or extinction of many freshwater fish populations in western Japan (Watanabe and Takahashi 2010). Such extinctions outside Lake Biwa would be evidenced by some relict species in the lake, such as Opsariichthys uncirostris and Ischikauia steenackeri, whose distributions have been expanded through artificial introduction (Sato et al. 2010; Watanabe 2012). The significantly lower genetic diversity in riverine populations compared with the Lake Biwa population of Pseudobagrus nudiceps may also support the severe riverine conditions in the past (Watanabe and Nishida 2003). The role of Lake Biwa as a reservoir of species and genetic diversity and, probably, as a source of dispersal, would have contributed to the formation and maintenance of ichthyofauna in western Japan.

The faunal history of Lake Biwa involving the absorption and accumulation of divergent lineages within species may have also contributed to the evolvability of Lake Biwa populations by increasing and preserving genetic diversity. Unlike those in the high‐latitude lakes of Europe and North America that were severely affected by glaciations (Avise 2000; Hewitt 2000, 2001), the ichthyofauna in the temperate zone of western Japan may have continuously developed since the Pliocene (e.g., Hewitt 2004; Watanabe and Takahashi 2010). The lake would have functioned as a biogeographic reservoir not only in the present Lake Biwa but also in the Paleo‐Lake Biwa stages. Standing genetic variations accumulated in the long history of Lake Biwa may have promoted rapid evolutionary responses in changing lake environments (Renaut et al. 2011; Kakioka et al. 2013, 2015). The new insights on the origin of the endemic ichthyofauna of Lake Biwa obtained through this study will be an important starting point for progressive studies of adaptive evolution in this ancient lake.