Ancient drainage networks mediated a large‐scale genetic introgression in the East Asian freshwater snails

Abstract Biogeography and genetic variation of freshwater organisms are influenced not only by current freshwater connections but also by past drainage networks. The Seto Inland Sea is a shallow enclosed sea in Japan, but geological evidence showed that a large freshwater drainage had intermittently appeared in this area between the late Pliocene and Pleistocene. Here, we demonstrated that this paleodrainage greatly affected the genetic variation of the East Asian freshwater snails, Semisulcospira spp. We found that the mtDNA haplotypes originated in the Lake Biwa endemic Semisulcospira species at the upstream side of the paleodrainage were frequently observed in the riverine Semisulcospira species at its downstream side. The genome‐wide DNA and morphological analyses consistently showed that there was no clear evidence of nuclear introgression between the Lake Biwa endemics and riverine species. These results suggest that the large paleodrainage had facilitated mitochondrial introgression and had broadly spread the introgressed mtDNA haplotypes to its downstream region around the Seto Inland Sea. Our study highlights the role of paleodrainages in shaping the genetic variation of freshwater organisms.


| INTRODUC TI ON
The pattern of biogeography and genetic variation in modern freshwater species should reflect not only the current landscapes but also the ancient geography (Hewitt, 2000;Petit et al., 2003). The historical change in drainage networks is a fundamental factor in shaping the distribution and genetic structure in obligate freshwater organisms, due to their limited dispersal opportunities among drainage systems. The drainages could temporally be extended by seawater regression and be contracted by transgression (de Bruyn et al., 2013;Thomaz, Malabarba, Bonatto, & Knowles, 2015;Thomaz, Malabarba, & Knowles, 2017;Unmack, Hammer, Adams, Johnson, & Dowling, 2013;Wu, Tsang, Chen, & Chu, 2016;Yan et al., 2013).
In addition, orogenesis can rearrange the drainage networks (Yan et al., 2013). These historical connections and dissociations of freshwater drainages should have changed the migration pattern of freshwater organisms and, therefore, intimately determine the distribution and genetic diversity of freshwater organisms.
The Seto Inland Sea is the largest enclosed sea in Japan and has historically experienced substantial changes in its topography and environment (Itihara, 1996;Kuwashiro, 1959). The region around the Seto Inland Sea is currently filled with seawater, but it had been a plain in the late Miocene (Itihara, 2001). Gentle block tilting during the late Pliocene and early Pleistocene produced a basin around the Seto Inland Sea area (the Setouchi Basin), which facilitated the | 8187 MIURA et Al.
formation of a large paleodrainage that is called the paleo-Setouchi drainage ( Figure 1a). The geological evidence exhibited that the paleodrainage included the ancient Lake Biwa on the upstream side and flowed west through the Setouchi Basin and Kyushu Island to the East China Sea (Itihara, 1996(Itihara, , 2001. The block tilting became active, and seawater had entered in this basin after the middle Pleistocene.
However, the successive regression and transgression by cyclical climatic changes during the Quaternary period produced temporal freshwater drainage in this basin several times because the Seto Inland Sea is a shallow sea with the current average depth of 31 m.
For example, a large part of the Setouchi Basin should be exposed and had converted into terrestrial and freshwater environments due to 140 m lowering of sea level during the Last Glacial Maximum about 20,000 years ago (Kuwashiro, 1959;Yashima, 1994). These geological events should rearrange the populations of freshwater organisms around the Setouchi Basin and affect their genetic diversity.
Specifically, since the river flowing from Lake Biwa currently end up at the Seto Inland Sea (Figure 1b), freshwater organisms in Lake Biwa hardly achieve gene flow at a broad spatial scale. However, when the large paleodrainage had formed around the Setouchi Basin, the river flowing from Lake Biwa was likely to be connected with other drainages around the basin, which may facilitate the connection of freshwater populations among the lake and rivers around the basin.
Semisulcospira is the freshwater snail genus broadly distributed in East Asia (Davis, 1969). Semisulcospira is viviparous with no planktonic larval stage and thus has limited dispersal ability.
Further, Semisulcospira is an obligate freshwater genus, which cannot traverse terrestrial or marine environments. Therefore, the Semisulcospira snails are ideal organisms to investigate the role of paleodrainages on the gene flow among the local populations.
Eighteen species of Semisulcospira are reported in Japan, while one species, Semisulcospira libertina, is the most common, being distributed in all regions in Japan (Kihira, Matsuda, & Uchiyama, 2003). In this study, we investigate mtDNA and genome DNA variations of Semisulcospira spp. to evaluate how and what extent the paleodrainage had affected the genetic variation of Semisulcospira spp. The entire configuration of the paleo-Setouchi drainage is still unclear, and the drainage might have connected to the river systems in the continental region (Itihara, 2001). We, therefore, collected the snails from broad geographical areas in Japan and Korea to evaluate the significance of the paleodrainage in shaping the genetic variation of the Semisulcospira snails.

| Sample collections
Approximately four hundred individuals of Semisulcospira spp. were collected at 82 sites in Japan and Korea (Figure 2 and Appendix S1).
The most common species, S. libertina, was collected from 67 sites.
Semisulcospira reiniana and S. kurodai were collected from four and one sites, respectively. The endemic species group in Lake Biwa (the Semisulcospira (Biwamelania) habei group and Semisulcospira (Biwamelania) decipiens group, see Miura, Urabe, Nishimura, Nakai, & Chiba, 2019) were obtained from 16 sites in Lake Biwa and its drainage. We identified these species based on adult shell features and the number and shape of embryos following the procedures of Davis (1969) and Watanabe and Nishino (1995). However, the taxonomy of Semisulcospira is currently confused. For example, although S. reiniana has a distinct shell morphology, it genetically belongs to a clade of S. libertina (Miura et al., 2019). Further, the species in the Lake Biwa endemic groups have very similar genome DNA sequences (Miura et al., 2019). We, therefore, use the species groups (the S. libertina group, S. (B.) habei group, and S. (B.) decipiens group), rather than morphospecies, as operational units in this study. 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, Torchin, Bermingham, Jacobs, and Hechinger (2012). Some of those samples were also used in the prior studies (Miura, Köhler, Lee, Li, & Foighil, 2013;Miura et al., 2019).

F I G U R E 1
The map of the Seto Inland Sea and Lake Biwa (a). The shaded area indicates the hypothetical location of the paleo-Setouchi drainage intermittently appeared between late Pliocene and Pleistocene. The area with horizontal lines indicates the Seto Inland Sea. The current drainage of Lake Biwa (b). The blue lines indicate the rivers flowing into Lake Biwa, and the red line indicates the river flowing from Lake Biwa.

| Mitochondrial DNA analyses
A total of 406 Semisulcospira spp. and related species were used for sequence analysis (Appendix S1). We analyzed mtDNA encoding the cytochrome oxidase c subunit I (COI) gene to identify major lineages and their phylogeographic patterns. The primer pair used in this study was COI-bf (5′-GGGGCTCCTGATATAGCTTTTCC-3′) (Miura, Torchin, Kuris, Hechinger, & Chiba, 2006) and COI-6 (5′-GGRTARTCNS WRTANCGNCGNGGYAT-3′) (Shimayama et al., 1990  The numbers in parentheses indicate the nodes used for the divergence time estimation (see Table 1). The upper map shows the distribution of the three subclades in mt-A clade (blue for mt-a1, red for mt-a2, and yellow for mt-a3), while the lower map shows the distribution of mt-B clade (gray). The detailed sampling locations for the Lake Biwa endemics are shown at the lower right. KX080149-KX080549 and LC500664-LC500669). We ensured the validity of sequences by inspecting the sequence chromatograms and by checking the quality scores of each base using Geneious 6.1.6 (Kearse et al., 2012). Sequences were aligned by ClustalW, implemented in Geneious for further phylogenetic analyses. Phylogenetic trees were constructed using the maximum likelihood (ML) algorithms. We searched for the best evolutionary models using MEGA 7.0.18 (Kumar, Stecher, & Tamura, 2016); the HKY + G + I model was selected. The ML analysis was also conducted using MEGA 7.0.18.

| Genome DNA analyses
Because of the unusual mitochondrial evolution in Semisulcospira (Köhler, 2017) and related snails (Whelan & Strong, 2016), other independent measurements are necessary to evaluate the phylogenetic relationship of Semisulcospira spp. We, therefore, analyzed genome-wide DNA variations in Semisulcospira spp. We used 39 S. libertina genome DNA samples; these samples were the subset of the samples used in the mtDNA analysis (Appendix S2). We analyzed the genome-wide DNA variations using a double digest restriction site-associated DNA library (ddRAD) sequencing tech- Raw sequence reads were processed using pyRAD 3.0.66 (Eaton, 2014). We integrated these sequence reads with the genome Sequences were demultiplexed 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 (Rognes, Flouri, Nichols, Quince, & Mahé, 2016) 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 across-sample 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 a clustering of paralogs (Hohenlohe, Amish, Catchen, Allendorf, & Luikart, 2011). Clusters shared among fewer than 50 individuals were excluded (the minimum taxon coverage = 50), and the remaining clusters were treated as ddRAD loci.
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 (GTR + G). Node robustness was assessed using bootstrapping and 100 replicates. J. plicifera, J. silicula, Koreoleptoxis sp., and K. tegulata were selected as an outgroup. We further estimated coancestry matrix to characterize genetic differentiation and admixture pattern, using the SNPs extracted from ddRAD loci. The coancestry values were calculated in fineRADstructure (Malinsky, Trucchi, Lawson, & Falush, 2018).

| Estimation of divergence time
To determine when the gene flow (mitochondrial introgression) occurred among the populations in the currently separated drainage systems, we estimated divergence time using mitochondrial DNA variations and a molecular clock. The oldest fossil species that is endemic to Lake Biwa is recorded from the Ueno and Iga Formation in the Kobiwako Group and Kameyama Formation in the Tokai Group (3.9-3.2 million years ago [MYA] [Matsuoka, 1985[Matsuoka, , 2001). This fossil species, Semisulcospira (Biwamelania) praemultigranosa, is the prospective earliest stem lineage of the Lake Biwa endemic snails (Matsuoka, 1985;Miura et al., 2019). Therefore, we set the split age of mt-B at 3.9 MYA, with the assumption of the Lake Biwa origin of mt-B. We selected eleven pairs of S. libertina and the Lake Biwa endemics, which are suspicious for mitochondrial introgression (see the numbered nodes in Figure 2) and estimated the age of their divergence. Because the relative rate test rejected equal evolutionary rate throughout the tree (p < .01), the divergence time was estimated based on uncorrelated relaxed lognormal clock implemented in BEAST v. 2.5.1 (Bouckaert et al., 2014). An exponential prior distribution was used for the split of mt-B. The Yule speciation model was used as a tree prior. The same evolutionary model, as the prior phylogenetic analysis, was used in the analysis. The analysis was run for 100 million generations, sampled every thousand steps, and the first thousand samples were discarded as burn-in. To check for convergence and to visualize the results, we used TRACER v. 1.6 and FIGTREE v. 1.4.2 (Drummond & Rambaut, 2007).

| Morphometrics
To evaluate whether introgressive hybridization affected the mor-

| RE SULTS
Based on 812 bp of the COI gene, 191 haplotypes were identified.
There were two largely distinct clades in Semisulcospira in Japan and Korea ( Figure 2). Mt-A was numerically predominant and was broadly distributed in Japan and Korea. Our dataset showed that this clade involved three subclades (mt-a1, mt-a2, and mt-a3). Mt-a1 mainly occurred in northeastern Japan, while mt-a2 occurred in southwestern Japan. In addition, mt-a3 was exclusively distributed in Korea ( Figure 2). The geographical distributions of these subclades were mostly separated, though mt-a1 and mt-a2 co-occurred in a few locations near the boundary of their distribution ranges (Figure 2).
On the other hand, mt-B involved highly diverged lineages and was distributed only around the Seto Inland Sea area and a few locations in Korea (Figure 2). Importantly, most of the endemic snails in Lake Biwa have mitochondrial DNA haplotypes that belong to mt-B, except for a few haplotypes that belong to mt-a2.  Figure 3). Further, this individual had a unique mtDNA haplotype that was located at the base of mt-A (Figures 2 and 4). Since we lack an adequate understanding of the evolutionary origin of this individual, we do not explore it in detail in this study.
The discordances between the mitochondrial and genome DNA phylogenies were observed (Figure 4a). The mt-B haplotypes were dominant in Lake Biwa, but these haplotypes also appeared in the  (Table 1).
The morphometric analyses demonstrated that there were no significant morphological differences between S. libertina with mt-A haplotypes and that with mt-B haplotypes (MANOVA, p = .96, Table 2) (see also Figure 5). However, the morphologies of these S. libertina groups were significantly different from those of the S.  Table 2).

| D ISCUSS I ON
We found that the clustering pattern based on mitochondrial DNA variations was largely concordant to that based on the genome DNA variations. However, there were apparent exceptions. The mitochondrial and genome DNA discordance was often observed in the populations around the Seto Inland Sea. The mtDNA haplotypes which were likely to be originated within the Lake Biwa endemics were broadly spread over the region around the Seto Inland Sea, while these snails had the genome DNA loci clustering with the local riverine species. These results suggest that the paleo-Setouchi drainage mediated mitochondrial introgression between the Lake Biwa endemics and riverine species and facilitated the spread of the introgressed haplotypes to the downstream area of the paleodrainage. Our results demonstrated that ancient connections of populations through the paleodrainage strongly affected the geographical genetic variations of the freshwater snails. Lee et al. (2007) first investigated the genetic variation of Semisulcospira spp. in Korea and found that there was a strange mitochondrial clade that contained phylogenetically divergent haplotypes. They suggested that this "enigmatic" divergent mitochondrial clade was originated due to ancestral polymorphisms or was captured through introgressive hybridization. Miura et al. (2013) demonstrated that the divergent mitochondrial clade in Korea historically migrated from Japan, potentially through hybridization between Korean and Japanese snails during Pleistocene Glacial Maxima when Korean Peninsula had been connected to the Japanese archipelago. Köhler (2016) showed that Lake Biwa is a potential reservoir for the divergent clade. We found that this divergent mitochondrial clade (= mt-B) was distributed exclusively around Lake Biwa and the Seto Inland Sea area in Japan and at the southern part of the Korean Peninsula ( Figure 2). Importantly, most of the divergent haplotypes  Table 2). These results further provided supportive evidence of the mitochondrial introgression between S. libertina and the Lake Biwa endemics, without the introgression of nuclear genes that control morphological characters of Semisulcospira spp.
Geological studies indicate that the large paleodrainage, which presumably involved several lakes and rivers, had appeared around the Setouchi Basin in the past (Itihara, 1996(Itihara, , 2001) (see   also Figure 1). This paleodrainage was first formed about 3 -1 MYA. After this period, the Setouchi Basin had been covered by seawater (Itihara, 1996). However, successive regression and transgression due to glaciations between 1 MYA and 20,000 years ago (until the Last Glacial Maximum) intermittently restored a temporal paleodrainage in this area (Kuwashiro, 1959;Yashima, 1994). This paleodrainage had strongly affected the migration patterns and phylogeography of freshwater fishes in Japan (Takahashi et al., 2020;Takehana, Nagai, Matsuda, Tsuchiya, & Sakaizumi, 2003;Watanabe et al., 2010). We assumed that the Lake Biwa endemics extended their distribution range through the paleodrainage and transferred their mtDNA haplotypes to local S. libertina populations through mitochondrial introgression. Matsuoka and Kitabayashi (2001) found several Pliocene fossils in the subgenus Biwamelania (the endemic subgenus name for the Lake Biwa Semisulcospira snails) at Oita Prefecture in Kyushu, which is located at the western edge of the Setouchi Basin, supporting the temporal range expansion of the Lake Biwa endemics. We evaluated the time of the range expansion by estimating the time for the mitochondrial introgression between S. libertina and the Lake Biwa endemics. Our results exhibited that most of the estimated times ranged between 1.0 and 0.2 MYA (Table 1), suggesting the temporal paleodrainage formed during the glacial regression could facilitate mitochondrial introgression and spread mt-B haplotypes around the Setouchi Basin. We also found that there was mitochondrial introgression in the reverse direction.
The mt-a2 haplotypes of S. libertina were introgressed to the Lake Biwa endemics about 1 MYA (the node [8] in Figure 2 and Table 1).
Although some of the estimated ages were much older (1.9-1.4 MYA, Table 1), we consider that these time estimations can be biased to older times because of the insufficient sampling of closer sisters, or the paleodrainage formed before the glacial period might also have contributed to the introgressive hybridization between S.
libertina and the Lake Biwa endemics.
The distribution of mt-B haplotypes overlapped with that of mt-a2 haplotypes ( Figure 2). However, although mt-a2 haplotypes were distributed in the Ryukyu islands, mt-B haplotypes were absent from these islands. Geological studies indicate that Ryukyu and Kyushu islands had often connected in the past, but they have been completely separated after the formation of the deep channel about 1.7 MYA (Kitamura & Kimoto, 2004). The majority of the divergence time estimates indicated that mt-B haplotypes moved outside Lake Biwa after the formation of the deep channel (Table 1) (Kawamura, 1998;Taruno, 2010). The estimated times for mitochondrial introgression between S. libertina and the Lake Biwa endemics were roughly consistent with the ages of the land bridge formation (Table 1). Our results indicated that the paleo-Setouchi drainage had connected to Korean drainage systems through the land bridge formed during the glacial regressions (see also Miura et al., 2013). The historical connection of Kyushu and the Korean Peninsula during the glacial period was also supported by phylogeographic studies of fishes (Tominaga, Nakajima, & Watanabe, 2016), insects (Park, Kitade, Schwarz, Kim, & Kim, 2006), and plants (Okaura, Quang, Ubukata, & Harada, 2007).
The high level of mtDNA variations observed in the Semisulcospira snails had been elusive, and several hypotheses were proposed to explain this enigmatic diversity (Köhler, 2016;Lee et al., 2007;Miura et al., 2013 The detailed hydrological network of the paleo-Setouchi drainage is still ambiguous despite geological research efforts (Itihara, 1996(Itihara, , 2001Kuwashiro, 1959;Yashima, 1994). Itihara (1996) suggested that the paleo-Setouchi drainage had potentially connected to the huge drainages in China, such as the Yellow River and Yangtze River, at the continental shelf exposed in the past. Although it is not easy to evaluate this hypothesis since these paleodrainage connections are currently hidden under the sea, we consider that this hypothesis will still be testable using Semisulcospira as a biological marker. If we can find mt-B haplotypes from the Semisulcospira individuals in the Chinese drainages, and they exhibit a sister relationship with the haplotypes of the Lake Biwa endemics, it will indicate the past connection between the paleo-Setouchi drainage and the Chinese drainages. Therefore, future studies on genetic variations in the Chinese Semisulcospira spp. will be fruitful, which will reveal not only the evolutionary history of the Semisulcospira snails but also the geographical evolution of the East Asian region. Our study demonstrates that the immobile freshwater species such as Semisulcospira spp. preserve the trace of ancient drainage networks within their DNA and thus can provide valuable information to reconstruct the ancient geography of focal regions.

ACK N OWLED G M ENTS
We thank M. Hayashi, Y. Sasaki, N. Takahashi

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest associated with this manuscript.