Human-mediated processes affecting distribution and genetic structure of Squalidus multimaculatus, a freshwater cyprinid with small spatial range

ABSTRACT Endemic species typically have a narrow niche breadth, and are likely more vulnerable to extinction than more widespread taxa. Squalidus multimaculatus is a small cyprinid endemic to the Korean Peninsula, and its reported geographical range was restricted to several small rivers located along the southeast coast. Several populations of S. multimaculatus have supposedly been subject to a variety of recent anthropogenic actions. Here, we analyzed the pattern of genetic diversity within and among populations of S. multimaculatus using nine microsatellite loci to quantify the relative contributions of human-mediated processes to the contemporary distribution and genetic structure. Overall, low levels of genetic diversity were exhibited in the populations of S. multimaculatus. Genetic differentiations among populations were not completely represented by their geographical proximity, likely resulting from the low intrapopulation genetic variability and anthropogenic transplants. The most conspicuous outcome of the anthropogenic activities was the introgression of alleles from a closely related species, S. gracilis majimae. Our study showed that anthropogenic transplanting, even with only a small number of individuals, can challenge our conservation goal to maintain the species integrity that has long been shaped in evolutionary processes.


Introduction
Evolutionary history and dispersal ability are key components in shaping natural distribution of a species, and each species has a unique geographic distribution (Goldberg and Lande 2007;Gaston 2009). While many species occupy wide geographic areas, other species are endemic, being restricted to a small range. Since endemic species typically have a narrow niche breadth, they are likely more prone to extinction than more widespread taxa (Gaston et al. 1997;Lamoreux et al. 2006). Thus, many of endemic species are primary targets of national conservation priorities (Lamoreux et al. 2006). The conservation of endemic species is particularly important in freshwater ecosystems, because endemism is exceptionally high (Revenga et al. 2005), and the global extinction rate of freshwater species is much higher than their terrestrial counterparts (Strayer and Dudgeon 2010;Collen et al. 2014). To develop systematic conservation strategies of endemic freshwater species, historical and contemporary factors that have shaped its current distribution and genetic structure must be determined.
Squalidus multimaculatus, the spotted barbel gudgeon, is a small cyprinid species endemic to the Korean Peninsula.
When this species was first reported, its geographical range was restricted to several eastern-flowing rivers located along the southeast coast (Hosoya and Jeon 1984; Figure  1). This species typically occupies a narrow range of habitats within tiny east coast rivers, raising concerns that the populations are relatively small and highly vulnerable to environmental disturbance. However, S. multimaculatus has never been officially listed as endangered in South Korea.
There are three other congeneric species occurring on the Korean Peninsula, S. gracilis majimae, S. chankaensis tsuchigae and S. japonicus coreanus (Froese and Pauly 2016). The distributions of these species are completely separated from that of S. multimaculatus because they are only found in rivers flowing west and south and their tributaries on southwest Korean Peninsula, while the natural populations of these three species have not been reported from any areas on the east coast ( Figure 1). The distributional pattern of these three congeneric species is likely related to the historical events of colonization from China through the Paleo-Yellow River (Nishimura 1974; see also Jeon and Suk 2014). Additionally considering that the natural populations of cyprinid species are very rarely found on the east coast, the distribution of S. multimaculatus is somewhat unusual. Based on the close phenotypic similarity between S. g. majimae and S. multimaculatus, it is likely that S. multimaculatus originated from divergence from S. g. majimae and colonization in the southern part of the east coast on the Korean Peninsula, although this speculation came before the detailed analysis.
However, the distributional range of S. multimaculatus has recently been expanded with reports of its presence in rivers further north of the known distribution (personal observations; Park et al. 2013). This rapid expansion may be a signature of anthropogenic transplant (Park et al. 2013). East coastal rivers on the Korean Peninsula are not physically connected with each other, and it is not likely that S. multimaculatus naturally migrate between rivers. Anthropogenic actions appeared to create hybridization problems in this species, because individuals showing a completely intermediate phenotype between S. multimaculatus and S. g. majimae have been found in the Jangan River, the southernmost population in the distribution of S. multimaculatus, from our recent field investigations. Genetic analyses are thus required for this species throughout its geographic range.
In this study, nine microsatellite loci were used to analyze the pattern of genetic diversity within and among populations of S. multimaculatus. Our ultimate objective was to quantify the relative contributions of human-mediated processes to the contemporary distribution and genetic structure of this species.

Sampling
The individuals used in this study were collected in accordance with the Inland Water Fisheries Act and Wildlife Protection and Management Act of Republic of Korea. This study was approved by the Yeungnam University Institutional Animal Care and Use Committee (Protocol # 2015013).
S. multimaculatus (N = 247) was sampled from 11 different drainages (populations) between February and March of 2013 using a kick-net (10 mm stretch mesh; Table S1 in supplementary material; Figure 1). S. g. majimae (N = 20) was also collected from a site in the Nakdong River ( Figure 1; Table S1 in supplementary material) to examine its introgression into S. multimaculatus populations. We rapidly counted the number of scales above the lateral-line (hereafter SAL), a morphological characteristic known to be different between S. multimaculatus and S. g. majimae (Hosoya and Jeon 1984), for all samples collected. A fin-clip (1 × 1 mm) was removed from the tip of each tail fin for genetic analyses, and all individuals collected were  Table S1 in supplementary material). released back to the original site. The tissue samples were placed in ∼1.5 ml of 95% ethanol until processing.

Microsatellite genotyping
Genomic DNA extraction was conducted using the Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer's protocols. Nine microsatellite loci were used for genotyping analyses (Table S2 in supplementary material). The forward primers of the loci were fluorescently labeled with FAM, HEX or NED (Applied Biosystems, Foster City, CA, USA). A GenePro Thermocycler (Bioer, Hangzhou, P. R. China) was used to amplify the microsatellites by subjecting the samples to the following conditions: denaturation at 94°C for 10 min, 35 cycles of 94°C for 30 sec, 55-60°C for 30 sec and 72°C for 30 sec, and final extension at 72°C for 10 min. PCR products were commercially run with an internal standard (ROX-500; Applied Biosystems) on an ABI3130 Genetic Analyzer. The sizes of the microsatellites were scored using Gene-Mapper 4.0 (Applied Biosystems).

Microsatellite analyses
The presence of possible genotyping errors, resulting from null alleles, large allele dropout and stutter peaks, was tested using MicroChecker 2.2.3 (Van Oosterhout et al. 2004). The level of diversity in each locus and each population was quantified as the mean number of alleles (A), observed (H O ) and expected heterozygosity (H E ) using Fstat 2.9.3.2 (Goudet 1995) and Arlequin 3.5 (Excoffier and Lischer 2010). For a proper comparison of interpopulation diversity, correction was done for sample sizes to quantify the allelic richness (A R ; Fstat). Departure of genotypic proportions from null hypothesis under Hardy-Weinberg equilibrium (HWE) was tested for each locus-population combination based on the exact test following Markov chain parameters with 100 batches and 1000 iterations per batch as implemented in Genepop 4.2 (Raymond and Rousset 1995). F IS was used to test the significance regarding the deviation from HWE in each population and locus.
The proportion of family relationships (full-, halfsibling and parent-offspring) was estimated for every possible pair of individuals within each population using ML-Relate (Kalinowski et al. 2006) and Colony 2.0 (Jones and Wang 2010). Bottleneck 1.2.02 (Piry et al. 1999) was applied to assess the likelihood of a recent decline in population sizes by Wilcoxon testing of the existence of significant excess in expected heterozygosity under mutation-drift equilibrium relative to HWE based on three different microsatellite mutation models, IAM (Infinite Allele Model), SMM (Stepwise Mutations Model) and TPM (Two-Phase Model with a setting of SMM = 70%). An additional examination was conducted based on a mode-shift away from the typical L-shaped allele distribution (using Bottleneck). Historical decline in population sizes was identified based on comparison of M-ratio (the mean ratio of the number of alleles to the range in allele size; Garza and Williamson 2001) calculated using Arlequin with the population specific M c (critical value of M ) threshold range generated in Critical_M (Garza and Williamson 2001). The M c threshold range was quantified using two different values of θ (4Neµ; 1-10) with a constant mutation rate of 5 × 10 −4 , the default (3.5) value of Δg (mean size of non-one-step changes) and P s of .9.
Three different approaches were applied to assess the level of the population genetic structure. First, genetic differentiation between populations was quantified based on pairwise-F ST and -R ST using Arlequin. These pairwise values were tested for the significance by Fisher's exact tests following 10 4 permutations in Arlequin. Second, bi-dimensional factorial correspondence analysis (FCA) was conducted on Genetix 4.05.2 (Belkhir et al. 2000) to visualize the genetic relationships among all genotyped individuals and all populations. The final approach was the model-based Bayesian structure analysis to examine the existence of distinct genetic clusters using Structure 2.3.4 (Pritchard et al. 2000). In this approach, posterior probabilities were calculated, while assuming HWE and minimizing linkage disequilibrium. Ten independent runs were performed under the admixture model, and each run was performed under 4 × 10 5 iterations while discarding the first 10 5 iterations as burn-in. The optimal value of K was determined based on the ΔK approach in Structure Harvester 0.6.93 (Evanno et al. 2005). The optimal K was also confirmed with highest similarity coefficient (H') calculated on Clumpak (Kopelman et al. 2015).

Interspecific hybridization
We attempted to select loci showing a clear difference between S. multimaculatus and S. g. majimae. Using the selected loci, we identified S. multimaculatus populations exhibiting alleles primarily found in S. g. majimae individuals. In the identified populations, NewHybrids v.1.1 beta (Anderson and Thompson 2002) was applied to allocate hybrids into each of hybrid categories (S. multimaculatus, S. g. majimae, F 1 , F 2 and backcross) based on a Bayesian posterior probability > 0.45. To examine the phenotypic expression of hybridization, we compared SAL measured for all individuals collected.

Microsatellite diversity
No significant evidence of genotyping errors was detected at nine microsatellite loci genotyped from 267 individuals of S. multimaculatus and S. g. majimae by MicroChecker. Low-to-moderate levels of genetic polymorphism were observed at all loci examined in S. multimaculatus, with the total number of alleles per locus ranging from 4 (LC27 and SQC4) to 17 (SQC60) and H O and H E ranged from 0.089 (LC27) to 0.502 (SQC60) and from 0.085 (LC27) to 0.586 (SQC60), respectively (Table 1). No significant deviation from HWE was observed at the locus level following Bonferroni correction (N = 9; α = 0.0056; Table 1). Population GM of S. g. majimae exhibited higher genetic diversity estimates than those observed from S. multimaculatus populations (Table 2). One population, CS, showed a significant departure from HWE following Bonferroni correction (N = 12; α = 0.0041; Table 2), presumably because of significant heterozygote deficiency in some loci.
As predicted from the genetic diversity data, S. multimaculatus populations exhibited relatively higher simulated proportions of family relationship than S. g. majimae (Table 3). Wilcoxon's tests for heterozygosity excess showed a strong signature of recent decline in population PN under all three mutation models when considering significance at α = 0.05, although no evidence was detected following Bonferroni correction (N = 11, α = 0.0045; Table 3). Two populations, MP and PN, likely experienced a recent significant bottleneck when tested based on the deviation from the normal L-shape of allele class distribution (Table 3). Although a strong signature of historical decline was detected in population PN, M ratios in other populations were similar to or higher than their M c threshold range (Figure 2).

Population structure
Genetic divergences among populations of S. multimaculatus were quite strong with high levels of average F ST and R ST per locus that ranged from 0.135 (MISA08) to 0.716 (SQC4) and from 0.133 (MISA08) to 0.715 (SQC4), respectively (Table 1). Most pairwise-F ST and -R ST values were statistically significant (Table 4) Table 4). A single population of Table 1. Summary of nine microsatellite loci used in this study to quantify the genetic diversity for 11 populations of S. multimaculatus.    Table 4). S. multimaculatus populations were most likely partitioned into five groups, MP-PN-CS-YO-DJ, BR, GG, HS-HY-TH and JA on the FCA scaling plot (Figure 3). JA also showed a strong affinity to S. g. majimae on the plot. Structure Harvestor proposed the most reliable number of genetic clusters as K = 3 (Figure 4(a)), which was also confirmed by the highest similarity coefficient (H ′ = 0.997) at K = 3 in Clumpak. Upon Structure analysis, two major clusters were represented by MP-PN-CS-YO-DJ and BR-HS-HY-TH, respectively (Figure 4(b)). Population GG showed an admixed nature of the two major clusters, and JA formed a separate cluster with S. g. mijimae (Figure 4(b)).

Interspecific hybridization
In our allele frequency data, allele 131 in LC27 and allele 169 in MFW1 were observed from all S. g. majimae individuals analyzed (Table S3 in supplementary material).
Although these alleles were also found in several S. multimaculatus populations, the frequencies were not greater than 0.05, except for YO and JA (Table S3 in supplementary material; Figure 5). The frequency of allele 169 was 0.125 in population YO (Table S3 in  According to the results from NewHybrids analyses, most of individuals (24 out of 26) in population JA were identified as mixed (more than two categories with more than 0.05) between S. multimaculatus and S. g. majimae (Figure 6(b)). Using a threshold of posterior probability > 0.45, a large proportion of individuals (N = 16; 61.5%) was assigned to F 2 genotype class, whereas a single individual and seven individuals were assigned to S. multimaculatus and S. g. majimae, respectively ( Figure 6(b)). Using the threshold, all of the individuals in population YO were assigned to S. multimaculatus (Figure 6(b)), suggesting the occurrence of unidirectional   introgression with a small number of individuals. In both population YO and JA, no individual was not classified as F 1 or backcross. JA individuals showed an average SAL intermediate to S. multimaculatus and S. g. majimae ( Figure S1 in supplementary material).

Discussion
Overall, low levels of genetic diversity were exhibited in the populations of S. multimaculatus. This was as expected considering that populations of narrowly distributed species typically show low genetic variation (Frankham 1996;Cole 2003). The low genetic diversity in S. multimaculatus was specifically obvious when compared with S. g. majimae, since no S. multimaculatus population was higher in any diversity estimates than population GM. Two categories of evolutionary scenarios could be put forward to explain the low-genetic diversity observed in S. multimaculatus. First, colonization of S. multimaculatus or its ancestor to the east coast has likely involved a limited number of individuals, leading to a loss of genetic diversity in the newly founded population. This assumption was possible, considering the contempory distribution of natural cyprinids and Squalidus species on the Korean Peninsula, as mentioned earlier. In order to verify this assumption, insight should be provided into how the founding individuals were able to disperse and colonize in the geographically isolated freshwater systems in the east coast. Biogeographic exchange between the southwestern part and the east coast on the Korean Peninsula  has been postulated in a few previous studies (Kim et al. 1989;Bae and Suk 2015). It is conceivable that eastward migration of the founders might happen probably during the Late Quaternary deformation events through the formation of confluences between tributaries in the Nakdong River and east coastal rivers (Kim et al. 1989;Bae and Suk 2015). In future studies, more extensive sampling of S. g. majimae throughout the Korean Peninsula is required to examine whether the populations in the Nakdong River show higher genetic affinity to S. multimaculatus than populations from other drainages.
Second, the physical characteristics of the east coast rivers have likely contributed to the contemporary level of genetic diversity in S. multimaculatus populations, because they are relatively small and spatially isolated from each other, despite their geographical proximity to each other. Ample habitats might not be provided for the populations to increase the sizes, which could lead to a loss of genetic variability. However, it is not likely that narrow habitat ranges have induced significant historical declines in population sizes, since M ratios of most populations were similar to or higher than the critical M values. Southern east coastal rivers might have often formed confluence with one another, having historically prevented the populations from declines (see Bae and Suk 2015). In such respect, this second scenario is somewhat less likely. For whatever reason, our genetic diversity data have important implications for the conservation of this species. From a population genetic perspective, our results suggest that S. multimaculatus populations become increasingly subject to genetic drift, likely resulting in the loss of their evolutionary potential. Genetic differentiations between populations were not completely represented by their geographical proximity. For example, population DJ located in the south did not cluster with the populations in geographical proximity, probably indicating that historical dispersal and geographic isolation may not be the only factors shaping the contemporary pattern of genetic structuring. A previous study speculated that S. multimaculatus individuals were transplanted to the Myeongpa River (population MP), which was previously devoid of S. multimaculatus (Park et al. 2013). This speculation can be partially supported by the results of our bottleneck tests, since a significant shifted mode from the typical L-shape of allele class distribution in population MP was detected, indicating that a small number of individuals were transplanted. Although it is not certain that population PN originated from artificial translocation, our bottleneck tests also showed the strong possibility of recent decline (founder effect) of this population. Population YO showed a much greater level of genetic variability than other populations in this cluster, probably due to the introgression of S. g. majimae genes at least partially, as discussed below.
Population BR, GG and JA exhibited an extraordinary genetic distinctness from other populations in different types of analyses. Population BR showed high levels of pairwise-F ST and -R ST , and clustered with geographically distant populations, HS-TH-HY in the Structure analysis. Population GG also showed high pairwise-F ST and -R ST from other populations, despite a partial overlap with other populations in the FCA plot. Upon Structure analysis, population GG represented an admixture of two major clusters, HS-HY-TH and MP-PN-YO-CS-DJ. Population BR and GG exhibited nearly complete monomorphisms in several loci with certain alleles that were infrequent in many other S. multimaculatus populations. To summarize, a low genetic diversity or small effective population size of BR and GG could increase the possibility of genetic drift in some loci, eventually leading to fixation or loss of any given alleles.
Genetic distinctness of JA from other populations is the effects of hybridization with S. g. majimae. Two diagnostic alleles of S. g. majimae were frequently observed in this population. Ongoing genetic exchange between these two species is not conceivable, since no F 1 hybrids were detected in our analyses. No significant deviation from HWE was found and a majority of individuals were assigned to advanced generations of hybrid within this population, indicating that reproductive barriers may be weak between these two species. In future studies, S. multimaculatus and S. g. majimae could be regarded as model species for investigate of the evolution of reproductive isolation in the early stages of speciation.
Diagnostic alleles of S. g. majimae were also observed in several S. multimaculatus populations, although the frequency was quite low relative to that of population JA. Among them, population YO showed a much higher frequency compared to all others. This unidirectional introgression of S. g. majimae alleles appears to be the outcome of anthropogenic transplants. However, further studies should be conducted to examine whether this introgression derived from the original S. g. majimae populations (i.e. the Nakdong River) or population JA. Considering that only two diagnostic alleles were used in our study, the overall levels of introgression might be underestimated (Johane et al. 2011). The hybridization in population JA could not simply be explained by unidirectional introgression from the alleles of S. g. majimae into the gene pool of S. multimaculatus, because the frequency of S. g. majimae alleles was slightly higher than that of S. multimaculatus in the diagnostic loci. Although direct evidence cannot be provided from our data, it is clear that S. g. majimae was not simply introduced into the population of S. multimaculatus in this river by artificial introduction. Instead, this river might be occupied by S. g. majimae before the secondary contact, or both species might be introduced into this river. Further investigation should be performed to resolve how these two species regained contact in this area.
Taken together, our results showed that the genetic structure among populations of S. multimaculatus has been shaped by the interplay of historical and humanmediated processes, and provided a lesson that anthropogenic transplanting, even with only a small number of individuals, can challenge our conservation goal to maintain the species integrity.