Genetic structure of Rhodeus sinensis reflects historical imprints of East Asian paleo-drainages

ABSTRACT

A complete reconstruction of the history regarding the distribution of a species requires a separate analysis of the events before and after colonizing the current distribution. In a freshwater species, the historical pattern of linkages between water systems, as an event before colonization, can be determined by examining the phylogenetic relationship between the populations from different drainages. Given the insular nature of freshwater systems, on the other hand, intra-population genetic diversity is largely related to the effective sizes of the first colonizing populations and the habitat capacities. This study analyzed the genetic diversity and structure of Rhodeus sinensis, a freshwater acheilognathid inhabiting the southern and western parts of the Korean Peninsula using two mitochondrial and twenty microsatellite loci. Korean R. sinensis appears to have originated from mainland China along the paleo-Yellow river across the vast continental shelf of the Yellow Sea and first colonized the west coast of the Korean Peninsula, based on phylogenetic data obtained from mitochondrial loci. These west coast populations have remained relatively small, and populations that migrated south appear to have thrived in the southwestern part of the peninsula, where the populations thereby displayed relatively high genetic diversity. Southern populations appear to have been formed by the movement of relatively small numbers following geological fluctuations, which likely have left traces of genetic drift throughout history, as also reflected in our results. The population of this species on the southeast coast of the peninsula appears to have been established by anthropogenic translocation.

KEYWORDS:

• Rhodeus sinensis
• Acheilognathidae
• Korean Peninsula
• mitochondrial DNA
• microsatellites

Introduction

The distribution of a species is a by-product of biogeographical dispersal in response to past climate change and adaptation to newly colonized habitats (Thomas 2010). A complete reconstruction of the history associated with the distribution of a species requires a separate analysis of the events before and after colonizing the current distribution. The genetic data from different classes should be used to accomplish this. By studying the phylogenetic relationship among populations of a particular species and closely related species, it is possible to reveal the origin and dispersal pathways of the species (Avise et al. 1987). Quantifying genetic diversity of populations or lineages provides information on the historical demographic pattern since the events of colonization (Austerlitz et al. 1997; Avise 2000; Crawford and Whitney 2010).

In a freshwater system, climate change is not simply a unique historical issue to be considered as a trigger for many species to spread out in search of new niches. Decreases in global temperature have been accompanied by lower sea levels, greatly expanding the total area of the freshwater system (Lambeck and Chappell 2001; Clark and Mix 2002; Lambeck et al. 2014). Consequently, the areas of shallow seas have often formed huge lakes or rivers, which may have served as pathways for the migration of freshwater fish species by connecting rivers currently isolated by the sea (Unmack et al. 2013; Zúñiga-Vega et al. 2014). The historical connections between drainages that do not currently exist can be determined by examining the phylogenetic relationship between the populations from different water systems. Even in the absence of historical inter-drainage confluence formation, gene flow between spatially separated water systems can occur through events, such as watershed capture due to geological fluctuations (Lima et al. 2017; Jeon et al. 2018) or anthropogenic introduction (Kwan et al. 2014; Baltazar-Soares et al. 2020). Given the insular nature of freshwater systems, the level of population genetic diversity within drainages is largely related to the effective sizes of the first colonizing populations and the environmental factors, such as habitat capacity (Pavlova et al. 2017; Martinez et al. 2018; Bernos et al. 2023). Watershed capture or anthropogenic introduction only allowed relatively fewer individuals to migrate than inter-drainage confluence formations through sea level fluctuations, likely generating a greater potential for a genetic bottleneck in newly formed populations (McCall and Fluker 2021).

The Yellow Sea is a gulf area located on the eastern edge of the Eurasian continent, surrounded by the Korean Peninsula and mainland China. Because the large continental shelf extending from the lands occupies most of this sea, the water depth is low (<100 m) throughout most of its area (Yang and Liu 2007), and the seafloor topography has been exposed frequently to sea level fluctuations in response to the historical repetition of warm and cold climates during the Pleistocene (Wageman et al. 1970; Park 1992; Yi and Saito 2004; Kong et al. 2006). The exposed seafloor topography of the continental shelf may have served as a link connecting rivers in China, such as the paleo-Yellow River, which flowed east through the middle of China, with rivers flowing to the west coast of the Korean Peninsula (Song et al. 2010; Jeon and Suk 2014; Jang et al. 2017; Kwan et al. 2018). Considering the species composition of freshwater fish, various species migrated from the Eurasian continent through the Yellow Sea and colonized the rivers flowing to the west coast of the Korean Peninsula (Won et al. 2020; Jeon et al. 2022; Hu et al. 2023). In addition, the Yellow Sea frequently or intermittently formed confluences allowing gene flow between rivers flowing to the west coast of the peninsula (Jeon and Suk 2014; Won et al. 2020; Jeon et al. 2022). Because the continental shelf is relatively small on the southern coast (Park et al. 1996), it is unlikely that the rivers flowing south had less chance of gene flow through sea level fluctuations.

This study analyzed the genetic diversity and structure of Rhodeus sinensis, a freshwater bitterling (Acheilognathidae) that inhabits the southern and western parts of the Korean Peninsula, by the analysis of two mitochondrial and 20 microsatellite loci. Bitterlings are well known for their unusual reproductive ecology, in which females deposit the eggs inside freshwater mussels, typically Unionidae, and males spray the sperm into the inhalent current of the mussels, thereby fertilizing the eggs (Barnes et al. 1936). R. sinensis also exists in China (Akai and Arai 1998; Chang et al. 2014). Korean R. sinensis was known as a Korean endemic species under the name R. uyekii until recently. Nevertheless, it is now regarded as the synonym of Chinese R. sinensis, based on the results of recent studies (Akai and Arai 1998; Chang et al. 2014). Given its distribution, it can be assumed that Korean R. sinensis originated from mainland China, settled around the west coast of the Korean Peninsula, and later expanded its distribution to the southeast. R. sinensis is the species with the widest distribution among bitterlings on the Korean Peninsula, distributed throughout most rivers in the western and southern regions (Kim and Park 2002). On the other hand, the genetic structure of this species has not been reported at the population level to prove this assumption.

Our study was designed to achieve three objectives. First, our genetic analysis tested whether the continental shelf and historical repetition of warm and cold climates during the Pleistocene shaped the distribution and structure of the R. sinensis populations during the proliferating process of freshwater species from continental to distal. There had probably been direct connections allowing potential active gene flow between rivers flowing down the west coast, and R. sinensis may have dispersed along the continental shelf to the southwest part of the peninsula. Our study attempted to find genetic and phylogenetic traces related to these biogeographic processes. Second, populations distributed in rivers flowing down the southern coast, where the continental shelf is not well-developed, may have originated from relatively small levels of gene flow without ancient confluence. Our genetic analysis tested wheter these populations are more likely to exhibit low levels of genetic diversity and the signatures of genetic drift than the populations of west coastal rivers. Finally, our results were compared with the results of previous genetic studies of other Rhodeus species with similar geographic distributions to R. sinensis, R. notatus (Won et al. 2020), and R. pseudosericeus (Jeon and Suk 2014; Won et al. 2021), which partially corroborated the assumptions related to the continental shelf and sea level changes.

Materials and methods

Sample collection and DNA isolation

Two hundred and fifty-five Rhodeus sinensis individuals were collected from eight localities throughout the Korean Peninsula between 2012 and 2015 (Figure 1; Table S1). The samples for genetic analysis were prepared by cutting a fin-clip (2 × 2 mm) from the tip of the caudal fin that can be regenerated (Uemoto et al. 2020). All individuals were released to the site where they were collected. The fin-clips were transferred to the laboratory and stored in 95% ethanol at −20 °C until DNA isolation. The genomic DNA was isolated using a Wizard Genomic DNA purification kit (Promega). Sample collections were performed in strict accordance with the Inland Water Fisheries Act and Wildlife Protection and Management Act of the Republic of Korea.

Figure 1. The sampling locations of Rhodeus sinensis (see Table 1 for the details of the locations) and the spatiotemporal pattern of its ancestral populations. Black dash lines represent the lowest sea level (−125 m) at the last glacial maximum, and blue lines indicate the paleo-drainage systems, both of which were reconstructed based on the bathymetric information extracted from the digital elevation model implemented in ArcGIS (see Materials and Methods for the details). Orange dashed line represents the limit of ancestral distributions for the ancestral R. sinensis around 1 MYA (million years ago). Orange solid and purple lines indicate the inferred colonization pathways of ancestral populations and the inferred event(s) of anthropogenic introductions, assumed based on our results.

Mitochondrial sequencing and haplotype analysis

Two mitochondrial loci, cytochrome b (cyt b) and cytochrome oxidase subunit I (COI), were amplified using previously reported primers (Brito et al. 1997; Perdices et al. 2004; Ivanova et al. 2007). Amplification was performed in a 25 μL reaction volume containing 1 μL of 30–50 ng of genomic DNA, 2.5 μL of reaction buffer, 0.25 μM of each dNTP, 0.3 μM of each primer, 1 unit of Taq polymerase (Solgent) and 19.93 μL of water. Amplification was carried out using GenePro (Bioer) under the following thermal cycling profile: initial denaturation at 94 °C for 3 min, 30 cycles of denaturation at 95 °C for 1 min, annealing at 56 °C for 1 min, extension at 72 °C for 1 min, and final extension at 72 °C for 10 min. The PCR products were loaded on 1.5% agarose gels to visualize the products for the checking of quality. Successful products were purified using Wizard Genomic DNA purification kit (Promega) and sequenced by Genotech (Daejeon, South Korea) on an ABI3730XL (Applied Biosystems) with BigDye Terminator 3.1 Cycle Sequencing Kit (Applied Biosystems).

All mitochondrial sequences were checked using BLAST searches and edited using GENEIOUS v 9.1.8 (Kearse et al. 2012). Alignments were performed using ClustalX (Thompson et al. 1994) implemented in GENEIOUS. All haplotypes characterized were sorted using DnaSP v 5.1 (Librado and Rozas 2009). All novel haplotype sequences found in this study were deposited in the NCBI GenBank (cyt b: OR496837 – OR496873, COI: OR485643 – OR485661).

Phylogenetic structure

The genetic diversity indices (Nei 1987) of the two mitochondrial loci, including the number of haplotypes (h), haplotype diversity (h_d), and nucleotide diversity (π), were calculated using DnaSP. The unrooted haplotype network was generated for each locus and the sequences generated by concatenating the two loci (hereafter ‘cyt b + COI’) based on the connection limit above 0.95 probability using TCS 1.2 (Clement et al. 2000). The cyt b haplotypes were used to reconstruct phylogenetic trees based on Bayesian inference (BI) algorithms performed using MrBayes v 3.1.4 (Ronquist et al. 2012). GTR + I + G was determined as the best-fitting substitution model according to Akaike Information Criterion (AIC; Akaike 1974) using jModelTest V 2 (Darriba et al. 2012). Four independent runs of Markov chain Monte Carlo (MCMC) were performed under 40 million generations, with sampling every 1,000th step. The convergence between runs was determined based on whether the effective sample sizes (ESSs) reached more than 200 using TRACER v 1.7 (Rambaut et al. 2018). The node confidence was provided with the Bayesian posterior probabilities. The nodes with low values (<50% or 0.5) for both algorithms collapsed. FigTree v 1.4.2 was used to visualize the topology of a consensus tree.

The divergence time between the lineages was estimated using BEAST v 2.3.0 (Bouckaert et al. 2014) on cyt b. Two independent runs were performed based on 40 million MCMC generations with sampling every 1,000 generations under GTR + I + G and a strict clock model with known substitution rates in cyprinids (1% per million years; Zardoya and Doadrio 1999). The molecular clock was calibrated using the fossil record data of Acheilognathidae (Yasuno 1984) for the basal node in the tree and Tanakia lanceolata (Tomoda 1977) for node clustering Rhodeus and Tanakia. The final consensus tree with maximum clade credibilities was generated, after discarding the first 25% of trees as a burn-in using TreeAnnotator v 2.3.0 (Bouckaert et al. 2014). The node credibilities of the consensus tree were evaluated using the posterior probabilities. TRACER was used to determine if the ESSs exceeded 200 and the MCMC parameters of runs converged on the same stationary point with a unimodal distribution. The geographic regions where the common ancestors of R. sinensis populations and the sister species were distributed were estimated using the dispersal-extinction-cladogenesis (DEC) model, which was implemented in RASP (Reconstruct Ancestral State in Phylogenies; Yu et al. 2015). Nine areas, including China and eight Korean drainages collected in this study, were coded to calculate the likelihood of each area being the distributional region of the common ancestor corresponding to the nodes of the time-calibrated consensus tree produced by BEAST analysis. This analysis was conducted under the JC model with the default parameters.

Demographic history

The mitochondrial loci were used to find the best approximation regarding the pattern of the demographic history for each R. sinensis population based on the three following approaches. First, Fu’s F_S (Fu 1997) and Tajima’s D (Tajima 1989) statistics were calculated under the coalescent model using ARLEQUIN 3.5 (Excoffier and Lischer 2010) using each mitochondrial locus and ‘cyt b + COI’. Second, mismatch distribution analysis (MMD; Rogers and Harpending 1992) was performed using ARLEQUIN by examining whether the observed distribution of pairwise differences between haplotype sequences deviated from the predicted pattern using ‘cyt b + COI’ to identify the statistically significant signature of historical population expansion. Finally, extended Bayesian skyline plot (EBSP) analysis (Drummond et al. 2012) was used to describe the historical changes in the ESSs under the HKY model and strict-clock model with a constant mutation rate (1% per million years; Zardoya and Doadrio 1999) using ‘cyt b + COI’. EBSP analysis was performed for 40 million generations with sampling every 1,000 steps; the first 25% were discarded as burn-in. TRACER was used to assess the convergence between runs by examining whether ESSs exceeded 200.

Microsatellite genotyping

All 255 individuals used in the mitochondrial analysis were genotyped using 20 previously reported microsatellite primer sets (Supplementary Table 2). The 5’ end of each forward primer was labeled fluorescently with either 6-FAM, VIC, NED, or PET (Applied Biosystems). PCR amplification was performed in a mixture of the same composition as that used for mitochondrial analysis under the following thermal cycling profile: initial denaturation at 95 °C for 5 min, 35 cycles of denaturation at 95 °C for 30 sec, annealing at 54–58 °C for 30 sec, extension at 72 °C for 45 sec, and a final extension at 72 °C for 10 min. The fluorescently labeled PCR products were genotyped on an ABI 3730xl Genetic Analyzer by Macrogen (Seoul, South Korea) and scored using GeneMapper 3.7 (Applied Biosystems). The allele size was determined using Peak Scanner 1.0 (Applied Biosystems).

Intrapopulation diversity

The average number of alleles per locus (A), allelic richness (A_R), heterozygosity observed (H_O) and expected (H_E), and fixation index (F_IS) were quantified using ARLEQUIN and GenAlEx v 6.503 (Peakall and Smouse 2012). The likelihood of a recent bottleneck was tested for each population using Bottleneck v 1.273 (Piry et al. 1999) by checking whether significant excess of expected heterozygosity quantified under the mutation-drift equilibrium relative to expected heterozygosity quantified under Hardy-Weinberg equilibrium existed based on the Wilcoxon sign-rank test under three different microsatellite mutation models, infinite alleles model (IAM), stepwise mutations model (SMM), and two-phase model (TPM) with a setting of 70% SMM and 30% IAM. Bottleneck was also used to determine if there was a marked deviation from the typical L-shaped distribution of the allelic frequencies that often occurs when the population has recently undergone a significant size reduction. The likelihood of a historical population decline was checked based on whether the M-ratio, the mean ratio of the number of alleles to the range in the allele size, was significantly lower than the commonly used threshold level, 0.68 (Garza and Williamson 2001). The critical significance values of all statistics used were adjusted for multiple comparisons based on the Bonferroni procedure.

Population structure

Three different analyses were applied to detect the genetic structure among the populations of R. sinensis. First, the pairwise-F_ST and -R_ST were quantified among the populations using ARLEQUIN. Second, STRUCTURE v 2.3.4 (Pritchard et al. 2000) was used to examine the pattern of population structuring under the Bayesian framework. The most likely numbers of genetically different clusters (K) of R. sinensis populations were predicted by ΔK statistics (Evanno et al. 2005) implemented in STRUCTURE HARVESTOR v0.6.94 (Earl and vonHoldt 2012). STRUCTURE analysis was performed for each K (2–8), with 20 independent MCMC runs, each consisting of 2 × 10⁶ generations after discarding 10% as burn-in. Finally, principal component analysis (PCA) was performed using GenAlEx to visualize the relationship among populations at the individual level in 2-dimensional space. The critical significance values of all statistics used were adjusted for multiple comparisons based on the Bonferroni procedure.

One population, HS, is located on the east coast of the Korean Peninsula. Because Acheilognathidae has not been reported on the east coast of this peninsula until relatively recently, it is necessary to determine if it is a population formed by anthropogenic introduction (Kim and An 2010). NewHybrids v 1.1 (Anderson and Thompson 2002) was thus used to identify the origin(s) and the level of admixture in this population (due to the potential origin from multiple populations). The parameters were set to MCMC of 10 generations after discarding the first 10% as burn-in.

Reconstruction of paleo-drainages

The paleo-drainage structure around the Korean Peninsula during the Last Glacial Maximum (LGM) was reconstructed from the bathymetric information extracted from the digital elevation model (DEM) using the hydrological implemented in ArcGIS10 (Figure 1). A 30-arc-second resolution (c. 1 km) of the General Bathymetric Chart of the Oceans (GEBCO, http://www.gebco.net/) was used as an input file. The sea level at LGM entered was 125 m, combining the results of previous studies (Elias et al. 1996). The watershed areas of paleo-drainages were delineated and extracted using the steps presented in a previous study (Thomaz et al. 2015). The sea level at LGM was determined using the Contour tool after masking the area around the Korean Peninsula using the Mask tool. The basin borders and flow direction were identified using the tools, Fill, Flow direction, and Basin. Paleo-drainages system for the masked areas was reconstructed using the Stream order tool.

Results

Mitochondrial analysis

Thirty-seven cyt b and 19 COI haplotypes were discovered from the analysis of 255 Rhodeus sinensis individuals. Among them, nine cyt b haplotypes were observed from at least two populations. In contrast, no COI haplotype was shared between populations, except for the population HS, which had the haplotypes of populations MG and ND. Population HS consisted of the cyt b haplotypes of populations HG, DC, MG, and ND. Haplotype network analysis showed that the haplotypes from the same population tended to be clustered closely together (Figure 2). From Bayesian phylogenetic analysis using cyt b data, Rhodeus species were likely divided into three major clades, R. sinensis, R. smithii group (R. notatus, R. atremius, R. fangi, and R. suigensis) and R. ocellatus complex (Figure 3). All Korean R. sinensis populations were resolved to form a single monophyly (Figure 3). On the other hand, Korean and Chinese (AB366522) R. sinensis did not form complete monophyly. The Korean populations even showed a closer phylogenetic affinity to R. amurensis, which raises the suspicion that R. sinensis may be a species complex (Figure 3). In the cyt b data, no distinct phylogenetic divisions between the west and south coastal populations were found, and the haplotypes of populations HG and DC were found in multiple clades (Figure 3).

Figure 2. (a) Ancestral distributions of Rhodeus sinensis populations reconstructed based on the two loci of mitochondrial haplotypes (‘cyt b + COI’); the phylogenetic tree for the ancestral area reconstruction was inferred using BEAST with the fossil calibration points for the major clades, and the outgroups were trimmed out to obtain the visual clearance. (b) Median-joining network analysis of R. sinensis haplotypes estimated using ‘cyt b + COI’; the numbers on the branches in the network represent the mutation steps.

Figure 3. Divergence time estimation of Rhodeus sinensis populations analyzed based on cyt b sequences. The blue shade represents the R. sinensis group, while the green shade indicates the R. smithii group that shows the sister relationship. R. ocellatus group are located at the basal placement in this phylogenetic tree.

Based on cyt b data, the common ancestor of R. sinensis and R. smithii group diverged from the R. ocellatus complex at approximately 17 MYA (Figure 3). The R. sinensis and R. smithii groups likely diverged at around 15 MYA (Figure 3). It was estimated that R. sinensis was geographically separated from Chinese R. sinensis and R. amurensis at around 10 and 3 MYA, respectively, and began to distribute in the freshwater systems on the Korean Peninsula (Figure 3). Considering such a phylogenetic pattern, R. sinensis likely colonized the Korean Peninsula for the first time around the Han River region at approximately 2 MYA. Subsequently, they took over the MG area and the southwest coastal corner by dispersing south (Figure 3). The BEAST results indicated that colonization at the population level occurred after 0.6 MYA (late Pleistocene; Figure 3). s-DIVA analysis showed that the geographical locations of the common ancestors of the major clades were found to be HG and DC, indicating that multiple dispersal routes may have formed from the population(s) located around the Han River and the south (Figure 2).

Fu's Fs and Tajima's D analysis yielded slightly different results of population demographic patterns depending on the mitochondrial loci. An analysis of the cyt b revealed the signatures of historical expansion in only two populations, MG and TJ (Supplementary Table 1). From COI analysis, the signatures of population expansion were observed in the DC, SJ, TJ, and ND populations (Supplementary Table 1). Giving high reliability for these estimates was difficult because the haplotype diversity of COI was low at the population level. Traces of historical expansion were found in populations MG, TJ, and ND when both mitochondrial loci were considered simultaneously (Supplementary Table 1). In mismatch distribution analysis, a unimodal distribution showing the historical expansion of the population was observed in populations MG, SJ, TJ, and YS (Supplementary Figure 1). Despite its multimodal pattern, population ND showed a distribution close to unimodal (Supplementary Figure 1). Populations HG and DC showed a perfect multimodal distribution pattern, confirming the existence of several isolates that may have been the beginning of a dispersal route to other regions, given that they are phylogenetically close to other populations (Supplementary Figure 1). Population HS showed a distinct multimodal distribution pattern, suggesting immigration originating from several different populations (Supplementary Figure 1). In the EBSP analysis, the signature of historical population size expansion was observed in populations MG, SJ, TJ, and YS (Supplementary Figure 2). As predicted from mismatch distribution analysis, a distinct increase pattern was not observed in the HG and DC populations (Figure S4). Populations ND and HS showed a sharp increase relatively recently, with no statistical signature of a historical increase (Supplementary Figure 2).

Microsatellite analysis

Extensive microsatellite polymorphism was detected at the twenty loci genotyped. For example, the observed number of alleles at a locus ranged from nine (RC741) to 60 (RY73), and H_E ranged from 0.168 (RC741) to 0.690 (RY73; Table S2). Significant departure from HWE was detected at two loci (RY73 and RY414; Supplementary Table 1) due to the heterozygote deficiency in several populations (data not shown). No significant deviation from HWE was observed at the population level (Table 1). Three populations located around the southwest tip of the Korean Peninsula, MG, TJ, and YS, exhibited a higher level of genetic diversity than the others (Table 1; Mann–Whitney U: P < 0.05 for all indices). In the heterozygosity excess test or mode-shift examination analyzed using BOTTLENECK, the signature of the contemporary bottleneck did not appear in any population (data not shown). On the other hand, an estimation of the M-ratio gave insight into the historical growth and decline of R. sinensis populations. Populations HG, DC, SJ, and HS appear to have undergone historical population size declines with values below 0.68. The three populations with high genetic diversity, MG, TJ, and YS, exhibited relatively high M-ratio values (Table 1).

Table 1. Genetic diversity indices estimated from the analysis of twenty microsatellites for the eight Rhodeus sinensis populations.

Pop	River	GPS coordinate	N	A	AR	HO	HE	M-ratio
HG	Han	36°01'23"/128°58'08"	33	4.00	2.832	0.278	0.296	0.66
DC	Daecheon	35°59'59"/129°14'58"	32	4.75	3.396	0.304	0.354	0.51
MG	Mangyeong	34°45'26"/126°59'27"	33	8.05	4.966	0.428	0.456	0.83
YS	Youngsan	34°39'16"/126°53'13"	33	6.50	4.444	0.527	0.547	0.79
TJ	Tamjin	34°57'22"/126°58'29"	33	9.10	5.557	0.539	0.551	0.83
SJ	Seomjin	35°51'15"/127°10'25"	33	4.10	2.826	0.334	0.337	0.61
ND	Nakdong	37°20'02"/127°40'37"	34	3.50	2.452	0.248	0.249	0.68
HS	Hyeongsan	36°23'28"/126°39'48"	24	4.20	3.594	0.498	0.529	0.50

A: the average number of alleles per locus; A_R: allelic richness; H_O: heterozygosity observed; H_E: heterozygosity expected; The information for drainage (River) and location of population (GPS coordinate) were also provided.

The pairwise-F_ST was 0.356–0.714 (mean = 0.531), and the pairwise-R_ST was 0.295–0.876 (mean = 0.674), indicating that the degree of genetic differentiation between populations was relatively high overall (Supplementary Table 3). The eight populations were allocated into six genetic clusters (best K = 6) based on the results of K statistics using STRUCTURE HARVESTER, indicating that each population has distinctive genetic characteristics (Supplementary Figure 4). From a more conservative point of view at K = 3, the STRUCTURE results were generally divided into three groups, HG-DC-MG, YS-TJ-SY, and ND-HS, in relation to the geographical location (Figure 4(a)). Populations YS-SJ and ND-HS were noted because they existed in the same cluster, even at different K values (Figure 4(a)). In PCA analysis, the first principal component (PC; 17.04%) showed similar results to the STRUCTURE analysis results at K = 3, but only the HG population was clearly separated in the second PC (13.73%; Supplementary Figure 3). In the present STRUCTURE and PCA results, population HS appeared to be the closest genetically to population ND, and shared some genetic traits with population MG, but not as much as ND (Figure 4(a); Supplementary Figure 3). In NewHybrids analysis, 11 out of 24 individuals had similar genetic characteristics to that of population ND individuals, and the remaining individuals exhibited a pattern of a genomic admixture (Figure 4(b)). On the other hand, none of the individuals from population HS possessed the genetic characteristics of individuals in the population MG (Figure 4(b)).

Figure 4. (a) Population structure of Rhodeus sinensis populations that was estimated using STRUCTURE analyses at K (cluster number) of 5–8 in descending order. (b) Contemporary genomic admixture analysis of population HS to discover the populational origin. Based on this analysis, population HS was at least partially originated from population ND. Each color bar represents posterior probability of individual ancestry, in which orange and green indicates the potential hybrid generations between populations (green: backcross, orange: F₂).

Discussion

In the robust phylogenetic data, the Korean Rhodeus sinensis populations were monophyletic, and the sister shared a recent common ancestry with R. amurensis occurring in the Russian Amur River. The clade to which R. sinensis belongs and the R. smithii group appeared to have differentiated in the middle Miocene (15 MYA) when the global temperature dropped sharply (Flower and Kennett 1994). A decline in sea levels accompanied the middle Miocene climate transition to such an extent that mainland China, the Korean Peninsula, and the Japanese Archipelago were connected (Noda 2004; Collins et al. 2018). A recent study suggested that the common ancestor of acheilognathid species may have originated in the Japanese Archipelago, based on the abundance of endemic bitterling species (Kawamura et al. 2014). Assuming that this hypothetical estimation is correct, a common ancestor of R. sinensis, which diverged from species of the R. smithii group abundant in the Japanese Archipelago in the middle Miocene, would have dispersed into the paleo-Yellow region, as supported by s-DIVA analysis. On the other hand, the biogeographic origin or dispersion of acheilognathid species distributed in East Asia is expected to become clear only when comprehensive sampling and phylogenetic studies based on substantial genetic information are performed.

Based on the phylogenetic data, the common ancestor of R. amurensis and Korean R. sinensis was distributed in the paleo-Amur and paleo-Yellow rivers. It is presumed that the Liao River in China flowing into the present Yellow Sea was, at least intermittently, connected to the paleo-Amur River until the late Pleistocene (Qiu et al. 2012; Xu et al. 2014). R. amurensis and Korean R. sinensis appear to have been differentiated from the Late Pliocene (3 MYA), based on the phylogenetic data. In the present study, Chinese R. sinensis diverged during Pliocene warming from the group of Korean R. sinensis and R. amurensis. It is difficult to draw clear conclusions from only the single Chinese R. sinensis used in this analysis. A clearer phylogenetic estimation will be possible if sampling and direct genetic information of Chinese R. sinensis are obtained from various regions.

s-DIVA analysis suggested that Korean R. sinensis was differentiated from R. amurensis in the Late Pliocene or Early Pleistocene and appeared to be distributed around populations HG (Han River) and DC (Daecheon River) for the first time on the Korean Peninsula. Considering that the haplotypes of populations HG and DC were distributed across several clades in phylogenetic analysis, and there was no graphic pattern of population growth by mismatch distribution analysis or EBSP analysis, these two populations were presumed to have consisted of small and isolated subpopulations. The reason that the populations in this area were fragmented without growing was probably due to environmental factors at the time of colonization, which can only be judged strictly through genetic comparisons with other freshwater fish species distributed widely on the Korean Peninsula. In the case of R. notatus, which does not inhabit the Daecheon River, but has a similar distribution to R. sinensis, the size of the Han River population was smaller than that of other western and southwestern coastal populations (Won et al. 2020). Some HG and DC haplotypes were clustered as phylogenetic sisters, indicating at least partial gene flow between these two populations until differentiation in the Late Pleistocene (ca. 0.5–0.6 MYA). Evidence of a biogeographical connection between the Han River and the Daecheon River was also found in a study of R. pseudosericeus, a closely related species to R. sinensis (Jeon and Suk 2014).

Based on the s-DIVA analysis and geographic location, it appears that some of the subpopulations of the population DC migrated to the southwestern corner of the Korean Peninsula and formed several populations. Population MG (Mangyeong River), which is geographically closest to population DC, appears to have formed first in this region, probably approximately 700,000 years ago. The mitochondrial and microsatellite data showed that the populations MG, TJ (Tamjin River), and YS (Yeongsan River) located near this area exhibited rapid growth and could retain high genetic diversity. This pattern of high genetic diversity in the southwestern corner populations of the Korean Peninsula was also found in a study of R. notatus (Won et al. 2020). When these populations colonized these rivers, it was presumed there were a wide variety of habitats that could accommodate large populations. Considering that these three populations formed distinctly different clades in the haplotype network analysis or phylogenetic trees, it was estimated that the level of gene flow between these rivers was not that active, but this corner area was connected to the paleo-Yellow confluence while dropping in sea levels during the Pleistocene.

Colonization in the rivers flowing to the south coast, where the continental shelf size is relatively small, occurred more recently than in the rivers flowing to the west coast, assuming no gene flow through a direct connection with paleo-yellow confluence. The historical demographic dynamics shown in the EBSP results indicated that the populations of the two rivers flowing down the south coast, populations SJ (Seomjin River) and ND (Nakdong River), grew rapidly at a significantly later time than that of the other populations. In these results, it is more a matter of where the populations SJ and ND came from rather than when they were formed. Population SJ showed the highest genetic affinity with population YS in mitochondrial and microsatellite data. Both the Tamjin and Yeongsan rivers are located in the southwestern corner of the Korean Peninsula, but given the estuary location, the Tamjin River is geographically closer to the Seomjin River. These results indicated that the population SJ did not colonize the present habitat through estuary coalescence. Because the upper areas of the Yeongsan River are adjacent to the upper areas of the Seomjin River, it is believed that capture accompanying the geological fluctuations in the headwater regions likely occurred (Kang 2005), which may have led to the migration of R. sinensis into the Seomjin River. According to the estimates of the time of divergence made in the present study, the populations SJ and YS were estimated to have differentiated at approximately 0.4 MYA.

Although it is possible to presume from a geographical point of view that the population of the Nakdong River (population ND) may have originated from the nearest adjacent population (Seomjin River), our results did not show any genetic or phylogenetic affinity between the populations SJ and ND. Based on the mitochondrial data, population TJ appears to be the origin of population ND. No empirical evidence has been reported to suggest that there may have been any geological events that would enable gene flow between these two geographically distant populations. Events that would have formed a confluence between the Nakdong River and Tamjin River regions can also be presumed, considering that the sea level off the southern coast of the Korean Peninsula was 50–60 m lower than at present, approximately 10,000 years ago (Lee et al. 2008). During this period of sea level descent, the Korean Peninsula was connected to the southernmost islands (e.g. Jeju Island), which can be supported by the floral similarity between Jeju Island and the southern coast of the Korean Peninsula (Chung et al. 2013). Because no empirical studies that make it possible to estimate biogeographic migrations in the southern coast of the Korean Peninsula have been conducted, more relevant data needs to be accumulated in the future.

Population HS was assumed to be formed by artificial introduction. The Hyeongsan River, where this population is located, is not likely to have been associated with paleo-Yellow confluences, either from its geographical location or from the surrounding mountainous terrain. Several species appearing with high frequency in this river flowing to the eastern coast of the Korean Peninsula (e.g. Squalidus multimaculatus, Iksookimia yongdokensis, and Liobagrus hyeongsanensis) are not found in the rivers flowing to the western and southern coasts (Kim and Park 2002). A previous study reported that anthropogenic inflow was frequently observed in the rivers flowing to the east coast of the Korean Peninsula (Lee et al. 2017). Microsatellite analysis showed that population HS showed a clear genetic affinity only with the geographically closest population ND. However, the mitochondrial haplotypes originating from population MG and ND were found frequently in population HS. Hence, population HS appears to have been introduced anthropogenically from the Nakdong River and other water systems, such as the Mangyeong River (population MG). Probably because relatively few individuals were introduced from the Mangyeong River, the genetic characteristics of microsatellites have been diluted over time, while the mitochondrial haplotype remained and was detected in the present study.

Based on these results, the management units for establishing a conservation strategy with Korean R. sinensis can be determined. Although R. sinensis inhabits the Korean Peninsula and mainland China, this study shows that the Korean and Chinese species are in a relatively distant evolutionary relationship to be considered the same species. This is why Korean R. sinensis should be recognized as a local native species, which can be the purpose for devising a conservation strategy, but extensive sampling and investigation in China are needed in the future. It is difficult to say that the current freshwater system on the Korean Peninsula is well preserved in terms of the ecological conditions that allow the natural populations of R. sinensis to maintain the historically formed size and genetic diversity. Korean R. sinensis is rare in many areas where large populations were naturally maintained 20 years ago (Choi 1973; Jang et al. 2008). Korean R. sinensis should be divided into three regions from the standpoint of natural population management: west coastal rivers (populations HG and DC), southwest coastal rivers (populations MG, TJ, and YS), and south coastal rivers (populations SJ and ND). The west coastal rivers can be considered the historical starting point of Korea R. sinensis, but the formation of large populations was complex, probably because of environmental factors, which likely led to relatively low population genetic diversity. For this reason, continuous monitoring of the populations and the surrounding ecosystems is needed to identify natural and anthropogenic factors that interfere with a population habitat or growth in this region. The southwest coastal rivers are a hotspot of Korean R. sinensis genetic diversity. Nevertheless, each population in this region has unique genetic characteristics. Populations in the south coastal rivers were formed by a small number of individuals during the dispersal process and are likely to have suffered historical bottlenecks.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The datasets generated during the current study are available from the corresponding author on reasonable request.

Additional information

Funding

This work was supported by National Research Foundation of Korea [grant number 2021R1A2C1014474]; National Institute of Biological Resources [grant number NIBR202005102].