Recent lake expansion triggered the adaptive radiation of freshwater snails in the ancient Lake Biwa

Abstract Lake expansion that leads to the formation of new habitats has potential to drive intralacustrine diversification. The ancient Lake Biwa in central Japan has historically experienced substantial changes in the lake size, and it provides a useful system for evaluating the role of lake‐size fluctuations in the diversification of endemic fauna. Here, we used genome‐wide DNA analyses and reconstructed the diversification history of the endemic freshwater snails belonging to the subgenus Biwamelania with respect to the geological history of Lake Biwa. We found that two genetically distinct snail lineages independently colonized Lake Biwa and they concurrently and rapidly radiated into 15 extant Biwamelania species. A combination of paleontological evidence and molecular dating technique demonstrated that the radiation of Biwamelania was tightly linked to the latest enlargement of the lake about 0.4 million years ago and suggested that increased ecological opportunity associated with the lake expansion drove the rapid adaptive radiation. We propose that the Biwamelania snails in Lake Biwa offer a promising new system for understanding the association between the geological history of the lake and rapid intralacustrine diversification.

Lake-size changes are one of the major drivers of extinction and speciation in lake fauna (Sturmbauer and Meyer 1992;Sturmbauer 1998;Kornfield and Smith 2000;Sturmbauer et al. 2001;Sturmbauer et al. 2011). A shallowing and reduction of lake area will elevate the extinction rate because of intense competition in the shrinking environments. In contrast, a deepening and enlargement of lake will promote population expansion and increase the speciation rate because of colonization to new habitats, reduced competition, and population subdivisions (Sturmbauer 1998). For example, the cichlids in the Great East African Lakes have shown that intralacustrine radiations occur upon lake refilling after severe drought (Owen et al. 1990;Sturmbauer and Meyer 1992;Sturmbauer 1998;Kornfield and Smith 2000;Sturmbauer et al. 2001;Genner et al. 2010;Sturmbauer et al. 2011). Severe drought could have extinguished many cichlid species while survived cichlids could have encountered several empty niches upon the lake expansion, which are considered to partly facilitate the rapid adaptive radiation in the African cichlids (Sturmbauer et al. 2011).
The cichlids in the Great East African Lakes have caused a major development in models for adaptive radiation (Albertson et al. 1999;Sturmbauer et al. 2001;Kocher 2004;Seehausen 2006;Sturmbauer et al. 2011). While these species flocks exhibit intralacustrine radiations and will certainly help in obtaining more information on evolutionary dynamics of species diversification, the intralacustrine radiations of endemic invertebrates such as mollusks may also have a great potential in contributing to our knowledge on species diversification (von Rintelen et al. 2004;Albrecht et al. 2006;Albrecht et al. 2008;Glaubrecht 2008;Schultheiß et al. 2009;Van Bocxlaer and Hunt 2013;Van Bocxlaer 2017). Because mollusks are less mobile, they are likely to become geographically isolated to a higher degree. Further, because of their calcareous shells, mollusks have a high fossilization potential. These behavioral and structural characters of mollusks enable us to effectively access to the information on both recent and past species histories.
The subgenus Biwamelania (Mollusca: Caenogastropoda: Semisulcospiridae) is composed of 15 extant species and 10 fossil species (Matsuoka 1987;Nishino and Watanabe 2000;Matsuoka and Miura 2018) and it is the most diverse endemic group in Lake Biwa (Watanabe and Nishino 1995;Nishino and Watanabe 2000;Kihira et al. 2003). Lake Biwa is an ancient lake located in central Honshu in Japan (Fig. 1A). The lake has experienced substantial lake-size changes throughout its history (Satoguchi 2012) (see Fig. 1B). First, the relatively shallow and small lake was formed south of the current Lake Biwa about four million years ago, and the lake depth and area substantially increased between 3.2 and 2.6 million years ago. After this deep and large lake period, the lake eventually moved north and almost dried up (but remained as swamps) until 1.8 million years ago. The lake basin moved to the location near the current Lake Biwa and was refilled about 1.2 million years ago. The enlargement of Lake Biwa to its present volume was estimated to have occurred approximately 0.4 million years ago (Satoguchi 2012). These lake-size changes would have repeatedly affected lacustrine habitat conditions and should have influenced the distribution and diversity of the endemic snails (Matsuoka 1987).
The fossil records have demonstrated that Biwamelania species experienced extinction and speciation events associated with the lake-size changes (Fig. 1C). The oldest Biwamelania species, Semisulcospira (Biwamelania) praemultigranosa, first appeared at Paleo-lake Tokai about 3.9 million years ago and it colonized Lake Biwa about 3.6 million years ago (Matsuoka, 1987(Matsuoka, , 2001; however, it became extinct during the deepening event of the lake about 3.2 million years ago. Semisulcospira (Biwamelania) reticulataformis MS and Semisulcospira (Biwamelania) nojirina MS subsequently emerged in the deep and large lake, but they had been disappeared before or during shallowing event of the lake about 2.6 million years ago. The shallow swamps between 2.6 and 1.8 million years ago harbored Semisulcospira (Biwamelania) gamoensis MS and Semisulcospira (Biwamelania) tagaensis MS, and these species were replaced by six Biwamelania species when the lake was refilled about 1.2 million years ago. Of the six species, five species were extirpated before the latest expansion of the lake about 0.4 million years ago (Matsuoka 1987;Matsuoka and Miura 2018).
Fifteen Biwamelania species are present in the current Lake Biwa, but most of them were not found from fossil records, except for Semisulcospira (Biwamelania) habei (Matsuoka 1987). This suggests that the most of extant Biwamelania species may have rapidly radiated after the latest expansion of the lake (Nishino and Watanabe 2000). However, this scenario should be carefully tested because the fossil record is generally less complete at lower taxonomic levels (Benton 1995). Although it has become possible to reconstruct the evolutionary history of a focal group on the basis of DNA variations in the extant species, previous molecular studies have demonstrated that the molecular phylogenies of the genus Semisulcospira on the basis of a single or few genes did not accurately reflect their evolutionary relationship because of the retention of ancestral polymorphisms or introgression (Lee et al. 2007;Miura et al. 2013;Köhler 2016). A previous study performed using several allozyme loci also failed to resolve their phylogenetic relationships because of insufficient variations in these loci (Kamiya et al. 2011). Only a method that assays variations in a large number of unlinked loci is likely to minimize the risk of inaccurate inference of phylogenetic relationships (Albertson et al. 1999;Carstens and Knowles 2007). In this study, we used a genome-wide DNA analyses based on a double digest restriction site associated DNA library sequencing technique to infer the robust phylogenetic relationship of Biwamelania and related riverine species. We tested the hypothesis of recent radiation of the Biwamelania snails associated with the latest expansion of the lake and we further evaluated what ecological factors facilitated the radiation of the Biwamelania snails. Finally, we highlighted how this system can serve as a useful system to provide insight into what processes drive adaptive radiation in ancient lakes.

SAMPLE COLLECTIONS
Fifty-four individuals of 14 endemic species in Lake Biwa were obtained from 15 sites in Lake Biwa and its drainage (see Table S1 for details, see also Fig. S1 for their shell morphology). We used snorkeling and SCUBA to collect the samples from deeper habitats and obtained Semisulcospira (Biwamelania) multigranosa and one specimen of Semisulcospira (Biwamelania) reticulata from a fisherman who dredged near Okishima Island in the lake. We could not find Semisulcospira (Biwamelania) ourense de-spite our sampling efforts. Semisulcospira (Biwamelania) ourense is an extremely rare species, and only a handful of individuals were recorded after its description (Watanabe and Nishino 1995;Kihira et al. 2003). Most of the Biwamelania species were collected from their type localities, although the type localities of Semisulcospira (Biwamelania) decipiens, Semisulcospira (Biwamelania) niponica, and Semisulcospira (Biwamelania) multigranosa are recorded just as "Lake Biwa". In addition to endemic species in Lake Biwa, 42 individuals of the riverine species, Semisulcospira libertina, Semisulcospira reiniana, and Semisulcospira kurodai were collected from 29, four, and one sites, respectively. We identified these species based on adult shell features, and number and shape of embryos following the procedures of Davis (1969) and Watanabe and Nishino (1995). Snails were either fixed in 95% ethanol, stored at -30°C or both for molecular analyses. We isolated genome DNA using a modified CTAB procedure described by Miura et al. (2012). Some of these specimens were used in the former mtDNA study (Miura et al. 2013).

OF RAD LOCI
We used 104 genome DNA samples (Table S1) for a double digest restriction site associated DNA library (ddRAD) sequencing technique, as described by Peterson et al. (2012) with a slight modification. Briefly, the extracted DNA from each individual was further purified using a nucleospin gDNA clean-up kit (Macherey-Nagel) with addition of RNase A. Approximately 30 ng of DNA were digested using two restriction enzymes (EcoRI and MspI). P1 and P2 adapters from Peterson et al. (2012) were ligated to the DNA fragments of each individual. The ligated samples were multiplexed and purified using a nucleospin gDNA clean-up kit. An E-gel size select agarose gel (Invitrogen, CA) was used to collect 300-350 bp DNA fragments. We amplified the DNA fragments in eight single PCR reactions. The PCR products were combined and cleaned using E-gel size select agarose gel and the nucleospin gDNA clean-up kit. The constructed DNA library was sent to Genome Quebec Innovation Center and sequenced using Illumina Hiseq 2000 single-end sequencing, yielding maximum read lengths of 100 bp. Raw sequence reads were processed using pyRAD 3.0.66 (Eaton 2014). Sequences were de-multiplexed using their sample-specific barcode without allowing any mismatches. The restriction site and barcode were removed from each sequence. A nucleotide base with a FASTQ quality score less than 20 was replaced with N. Sequences having more than 5% Ns were discarded. Sequences within each sample were clustered using VSEARCH (https://github.com/torognes/vsearch) with an 85% similarity threshold, following the pyRAD SE ddRAD tutorial. Within-sample clusters with fewer than 10 sequences were excluded to ensure accurate base calls. Consensus sequences were created based on the clusters with consideration of the error rate and heterozygosity. Consensus sequences from all samples were clustered using the same similarity threshold that was applied in the within-sample clustering. The resulting acrosssample clusters were aligned with MUSCLE (Edgar 2004). Any clusters having more than 5% shared polymorphic sites were discarded, because a shared heterozygous site across many samples likely represents clustering of paralogs (Hohenlohe et al. 2011). Clusters shared among fewer than 50 individuals were excluded, and the remaining clusters were treated as ddRAD loci (Table S1).

INFERENCES
The ddRAD sequences were concatenated into a single sequence alignment by an output function in pyRAD. Phylogenetic analysis was conducted by maximum likelihood (ML) algorism, using RAxML 8.0.20 (Stamatakis 2014) with general time reversible and gamma model. Node robustness was assessed using bootstrapping and 100 replicates. Other semisulcospirid species, Semisulcospira extensa, Parajuga sp., Juga silicula, and Juga plicifera were selected as outgroups, based on Lee et al. (2007), Strong andKöhler (2009), andKöhler (2016). As concatenation of genome DNA sequence data can be problematic because of spuriously high bootstrap supports for incorrect partitions (Gadagkar et al. 2005), we also inferred species trees based on multispecies coalescent model. Gene tree for each locus was estimated by RAxML with general time reversible and gamma model using the MAGNET pipeline (http://github.com/justincbagley/MAGNET). We then used ASTRAL III (Mirarab and Warnow 2015), which estimates a species tree that agrees with the largest number of quartet trees within a set of unrooted gene trees.
We estimated divergence time of the Biwamelania species using BPP 3.3 (Yang 2015). This analysis assumes no recombination among loci, neutral clock-like evolution with JC69 mutation model, and no migration among species. Because the Bayesian analyses using the multispecies coalescent are computationally expensive, we used the subset of individuals and ddRAD loci shared among all those individuals to reduce the size of the dataset. We used species tree topology estimated by ASTRAL III and estimated relative divergence time (τ) and population size (θ) at each node. The prior for τ was Gamma (2, 250) and θ is Gamma (2, 1000). We converted τ value to actual time (t) using fossil calibration. The oldest fossil Biwamelania species, S. (B.) praemultigranosa, is recorded from the Ueno and Iga Formation of the Kobiwako Group and Kameyama Formation of the Tokai Group (3.9-3.2 million years ago) (Matsuoka, 1987(Matsuoka, , 2001. Semisulcospira (Biwamelania) praemultigranosa has an elongated conical shell outline and a small number (two or three) of basal cords on the body whorl, which correspond to diagnostic characters of the subgenus Biwamelania (Matsuoka 1985). This fossil species is the prospective earliest stem lineage of Biwamelania. Therefore, the age of the node representing the split of the basal group of Biwamelania [S. (B.) habei group] from the other Japanese Semisulcospira was set for 3.9 million years ago, which represents the first appearance of the subgenus Biwamelania (Matsuoka 2001) (see also Fig. 1). We used TRACER v. 1.6 and FIGTREE v. 1.4.2 (Drummond and Rambaut 2007) to check for convergence and to visualize the results. We also estimated net diversification rate of the Biwamelania and related riverine species using BAMM 2.5.0 (Rabosky 2014), based on the tree obtained by the BPP analysis. MCMC simulations in BAMM were run for 100 million steps, sampling parameters every 10,000 steps, and the MCMC results were analyzed within the R package BAMMtools .
The demographic histories of the Biwamelania species were reconstructed using Extended Bayesian Skyline Plot (EBSP) analysis (Heled and Drummond 2008), implemented in BEAST2 (Bouckaert et al. 2014). For each Biwamelania species, we randomly selected 200 loci that are shared by all individuals within the species. We used the same evolutionary model (JC69) as in the above BPP analysis and used the mutation rate estimated by the BPP analysis. The analyses were run for 1.5 billion generations, sampled every 50,000 steps and first 10-20% of samples were discarded as burnin. We used TRACER v. 1.6 to check convergence, and up to three independent runs were combined if the runs were not converged. We reduced the number of loci to 100 for Semisulcospira (Biwamelania) nakasekoae and Semisulcospira (Biwamelania) takeshimensis because their parameters did not converge in a reasonable computation time. The estimated relative population size was plotted against time for each Biwamelania species using R v3.3 (https://www.r-project.org).

ECOLOGICAL ASSESSMENTS
To evaluate habitat usage patterns of the Biwamelania species, we classified substrate types of the sampling locations into the six categories: mud, sandy mud, sand, sandy gravel, pebble, and rock. Sandy mud is a mixture of mud and sand, and sandy gravel is a mixture of sand and gravel. We roughly identified the substrate types by eye at the field based on dominant particle size. We also recorded the depth range at the sampling locations. We then integrated our dataset and the habitat information of each Biwamelania species reported in the study of Watanabe and Nishino (1995), which contain about a hundred sampling points in Lake Biwa. Note that the numbers of the sampling locations for some species are limited as these species are distributed in restricted locations (see Fig. S2). We used chi-square test to evaluate difference in habitat types among Biwamelania species. We used mid-range as a representative depth at each sampling location and evaluated the difference in habitat depth among Biwamelania species. We compared the distribution of Biwamelania species along the lake depth using a general linear model. These statistical tests were performed using JMP V. 9.0 (SAS Institute).

Results
We obtained 0.4 to 15.3 million reads for each individual after de-multiplexing the Illumina Hiseq raw dataset (Table S1). The average number of clusters greater than nine sequences was 15,271.1; the average coverage achieved per individual per loci was 55.7 (Table S1). There were 6,285 loci with a total alignment length of 597,943 bp and the dataset had 31% of missing data.
The ML tree based on the concatenated ddRAD sequences is shown in Fig. 2A. There were five genetically well-supported clades detected in Semisulcospira spp. All of these clades were supported by the highest bootstrap value (100%). Semisulcospira libertina L1 was distributed at the east side of Japan while S. libertina L2 was distributed at the west side of Japan. Semisulcospira reiniana was included in S. libertina L2. Semisulcospira libertina L3 was exclusively observed at the sites in Korea. The clear monophyly and a high level of divergence among the clades L1-3 suggest that these clades represent distinct biological species. However, we do not assign scientific names to these clades because the taxonomic revision is not the goal of this study. There were two largely separated clades in the Lake Biwa endemic species (Fig. 2A). One clade at basal position included Semisulcospira (Biwamelania) habei, Semisulcospira (Biwamelania) dilatata, Semisulcospira (Biwamelania) rugosa, Semisulcospira (Biwamelania) fuscata, Semisulcospira (Biwamelania) niponica, and Semisulcospira (Biwamelania) reticulata; this clade also included non-endemic species, Semisulcospira kurodai. The other clade included Semisulcospira (Biwamelania) decipiens, Semisulcospira (Biwamelania) nakasekoae, Semisulcospira (Biwamelania) fluvialis, Semisulcospira (Biwamelania) arenicola, Semisulcospira (Biwamelania) takeshimensis, Semisulcospira (Biwamelania) shiraishiensis, Semisulcospira (Biwamelania) multigranosa, and Semisulcospira (Biwamelania) morii. We followed the reports by Kamiya et al. (2011) by calling the former clade as the S. (B.) habei group and the later clade as the S. (B.) decipiens group. The species tree based on multispecies coalescent models was similar to that of the concatenated ML tree (Fig. 2B). However, the species tree showed that the phylogenetic positions of the species within the S. (B.) decipiens groups were often uncertain. The ddRAD loci shared among all selected individuals (in total of 403 loci) were used for the divergence time estimates. The estimated divergence times for major clades and the Lake Biwa endemics are shown in Fig. 3A. The estimated mutation rate of Biwamelania was 1.06 ×10 −9 substitution per year, which is comparable to the mutation rate in other organisms in diverse taxonomic group (Lynch 2010). The BAMM analysis demonstrated that the net diversification rate was increased after the latest enlargement of Lake Biwa (Fig. 3B). The EBSP analyses further exhibited that the most of Biwamelania species have experienced population expansions during and after the latest enlargement of the lake (Fig. 4).

Discussion
Molecular phylogeny based on the ddRAD sequencing technique showed that there are two largely distinct groups in the subgenus Biwamelania ( Fig. 2A, B) (L1-L3). These two endemic species groups and three riverine clades should be taxonomically revised in future study. Molecular phylogenies provide a framework for studying endemic radiations. Our tree rejected the hypothesis that the Biwamelania species flock originated from a single colonization event ( Fig. 2A, B). The tree topology is most plausibly interpreted by assuming two independent colonizations of Biwamelania in Lake Biwa. The first colonization by the S. (B.) habei group occurred at the initial stage of Lake Biwa, and the second colonization by the S. (B.) decipiens group occurred about 1.90 million years ago (1.28-2.49 million years ago: 95% HPD; Fig. 3A) when Lake Biwa was a group of shallow swamps (Satoguchi 2012). Multiple colonizations have often been reported in other endemic radiations in ancient lakes. For example, the cichlids in Lake Malawi were colonized in the lake at least two times (Joyce et al. 2011). Further, Tylomelania snails in ancient lakes in Sulawesi colonized four times from surrounding rivers (von Rintelen et al. 2004). These colonization events often facilitated diversifications of lineages in ancient lakes, while not all colonizations resulted in diversification scenario (e.g., Peart et al. 2014). Matsuoka and Miura (2018)  these five old fossil species are the ancestral species in the S. (B.) habei group. Detailed morphological examination of these fossil species will provide an opportunity to evaluate the validity of this hypothesis.
The fossil records and the geology of the lake suggested that lake-size change is a major factor for extinction and speciation (Fig. 1C). The latest extinction event in Lake Biwa occurred about 0.4 million years ago, when the lake basin was substantially enlarged and deepened as a result of a fault-block movement called Rokko Movements (Matsuoka 1987). Six Biwamelania species were present in the fossil record before this extinction event (Matsuoka 1987;Matsuoka and Miura 2018). However, all but one species were extirpated, and S. (B.) habei is the only species that appears both in the past and current lake (Fig. 1C). Therefore, the fossil evidence suggests that more than 10 species of extant Biwamelania snails rapidly radiated following the latest expansion of the lake (Nishino and Watanabe 2000). If this fossil-based inference is correct, we should expect to find a shallow, bushy phylogeny among Biwamelania species. Consistent with this expectation, we found that the phylogenetic relationships within the S. (B.) habei and S. (S.) decipiens groups were characterized by short branches (Fig. 2). The molecular dating analysis demonstrated that the divergence event of the S. (B.) habei group began about 0.35 million years ago and the S. (S.) decipiens group began to diverge about 0.17 million years ago (Fig. 3A), demonstrating that the radiation events in Biwamelania occurred following the latest increase in lake size. The net diversification rate also increased after the enlargement of the lake (Fig. 3B), further supporting the recent radiation of the Biwamelania snails.
The demographic reconstruction of the Biwamelania species demonstrated that the population of the most of Biwamelania snail expanded after the latest enlargement of the lake about 0.4 million years ago (Fig. 4). Similar patterns were reported in several fish populations in Lake Biwa. The demographic inferences based on mtDNA variations showed that 22 fish species or clades in Lake Biwa expanded their populations after the latest enlargement of the lake (Tabata et al. 2016). These demographic patterns suggest that the latest enlargement of the lake provided new stable No. sites=41 No. sites=4 No. sites=4 No. sites=4 No. sites=23 No. sites=8 No. sites=3 No. sites=1 No. sites=3 No. sites=2 No. sites=18 No. sites=2  habitats for faunas in Lake Biwa and resulted in the concurrent demographic increases across the diverse faunas.
New ecological niches should become available with the expansion of Lake Biwa. Adaptive radiation can take place when founders enter a new environment with a number of discrete ecological niches (Gavrilets and Losos 2009), and thus, the exploitation and specialization to new habitats is likely to be an important process in rapid diversification of the Biwamelania snails. Consistent with this idea, the Biwamelania species radiated from two ancestral species after the latest lake expansion now use a variety of habitats in Lake Biwa. In the S.  (Fig. 5B). These habitat usage patterns can be independently evolved between two groups. The exploitation of new habitats could result in the reduction of resource competition and are suggestive of a significant role for ecological factors in rapid diversification in Biwamelania snails.
The studies on species radiation in the Great East African Lakes have shown that cichlids inhabit the littoral zone were allopatrically differentiated at a small geographical scale during the expansion of the lake, owing to extreme territoriality and lack of dispersal opportunities during any life stage (Sturmbauer et al. 2011). Indeed, many cichlids are narrow endemics present only in a single stretch of continuous habitat (Ribbink et al. 1983) and exhibited low level of gene flow among local populations (Genner et al. 2010). This pattern was also reported in a shallow water catfish in Lake Tanganyika (Peart et al. 2018). Similar to the case of the cichlids and catfish in African ancient lakes, the geographic distribution of nine Biwamelania species is confined to small regions in the lake (Watanabe and Nishino 1995, see  (B) nakasekoae are distributed at only a part in Uji River, which is a major outlet of Lake Biwa (Fig. S2). Biwamelania species are ovoviviparous snails with no planktonic stages, and thereby, the different coasts within the lake may have been sufficient to isolate their populations. These distribution patterns suggest that, in addition to ecological factors, spatial factors can also play an important role in the radiation of the Biwamelania snails in the lake.
Several other factors may facilitate radiation in the Biwamelania species. Variation in the radula morphology of S. (B.) decipiens, S. (B.) multigranosa, and S. (B.) reticulata has been observed (Watanabe 1970;Prozorova and Rasshepkina 2006). The radula morphology often reflects the trophic system, such as food resource usage pattern (Hawkins et al. 1989). It suggests that the trophic specialization may also, in part, contribute to the diversification in Biwamelania, as it observed in the radiation of Tylomelania snails at ancient lakes in Sulawesi (von Rintelen et al. 2004). In addition, the Biwamelania species exhibited a high level of karyotype variations among species (Burch 1968;Kobayashi 1986). Reproduction between species with different karyotypes can yield hybrids that are heterozygous for chromosomal rearrangements, and these hybrids typically have reduced fertility because of error during the first meiotic division (Forejt 1996). Therefore, these karyotype variations may also account for the evolution of rapid reproductive isolation among the Biwamelania species.
We demonstrated that two distinct lineages of the Biwamelania snails were concurrently radiated during the latest expansion of Lake Biwa. Our results exhibited the potential of the Biwamelania snails in Lake Biwa to serve as a useful system for determining the evolutionary factors in speciation and adaptive radiation. Further development of evolutionary models for Biwamelania snails with ecological, paleontological, and karyotypic perspectives should provide insight into the relative importance of each evolutionary factor on radiation in Biwamelania snails.