Population genetics and sympatric divergence of the freshwater gudgeon, Gobiobotia filifer, in the Yangtze River inferred from mitochondrial DNA

Abstract The ecosystem and Pleistocene glaciations play important roles in population demography. The freshwater gudgeon, Gobiobotia filifer, is an endemic benthic fish in the Yangtze River and is a good model for ecological and evolutionary studies. This study aimed to decode the population structure of G. filifer in the Yangtze River and reveal whether divergence occurred before or after population radiation. A total of 292 specimens from eight locations in the upper and middle reaches of the Yangtze River were collected from 2014 to 2016 and analyzed via mitochondrial DNA Cyt b gene sequencing. A moderately high level of genetic diversity was found without structures among the population. However, phylogenetic and network topology showed two distinct haplotype groups, and each group contained a similar proportion of individuals from all sampled sites. This suggested the existence of two genetically divergent source populations in G. filifer. We deduced that a secondary contact of distinct glacial refugia was the main factor creating sympatric populations of G. filifer, and climate improvement promoted population expansion and colonization.

The subfamily Gobiobotinae (Cyprinidae, Cypriniformes) is a group of small freshwater gudgeon distributed in the rivers of East Asia in Korea and China (Chen, 1998). Only two genera, Gobiobotia and Xenophysogobio, with 17 species belong to this subfamily, and all the fishes are benthic. He (1991) examined the skeletons and classified Xenophysogobio as a primitive species and Gobiobotia filifer as a specialized species. Wang, He, and Chen (2002) used mitochondrial DNA (mtDNA) data to construct the phylogenetic relationship that supported this idea. They also inferred that Gobiobotinae originated in the upper reaches of the Yangtze River. At present, four fishes, including X. boulengeri, X. nudicorpa, G. abbreviata, and G. filifer, live in the upper reaches of the Yangtze River. However, except for these four fishes, sympatric habitation among Gobiobotinae fishes is rare, indicating that they could have quickly adapted and specialized to new ecological environments.
Gobiobotia filifer is endemic to the Yangtze River basin (Chen, 1998), and the limits of its native range extend from the Yibin City to downstream, including tributary rivers such as Min River, Chishui River, Hanjiang River, and Xiangjiang River ( Figure 1). This fish produces drifting eggs in spawning reason (Liu et al., 2018;Tian et al., 2017) and feeds on benthic organisms. The Gezhouba Dam was completed in 1988 in Yichang City, and then, the Three Gorges Dam with 175 m of water deep was completed in 2006, which form physical and ecological barriers to genetic exchange among populations. Recent surveys found that the abundance of G. filifer in upper Yangtze River was relatively high comparing with that before the construction of the Three Gorges Dam (Fan, Ba, & Duan, 2012;Xiong, Liu, Duan, Liu, & Chen, 2015;Yang, Xin, Ma, Kong, & Liu, 2010). However, it was less common in section and tributary rivers below the dams now (Fan et al., 2012).
Own to the specialized position in phylogeny and relatively broad distribution, G. filifer is a good mode for studying the local adaptation and the effect of dam on fish. For the last decades, mitochondrial DNA has been proven to be an important tool in population history, biogeography, genetic structure, and species delimitation (Hebert, Penton, Burns, Janzen, & Hallwachs, 2004;Moore, 1995). In this study, we used cytochrome b (Cyt b) to examine population genetic diversity of G. filifer. We aimed to (a) reveal whether genetic divergence occurred within the population and (b) determine whether divergence occurred before or after population radiation.

| Sampling and DNA extraction
A total of 292 individuals of G. filifer were collected from eight locations from 2014 to 2016; four locations were in the upper reaches of the Yangtze River and four in the middle reaches (Table 1; Figure 1).
A small fin sample was clipped and preserved in 95% ethanol for each specimen.
Genomic DNA was extracted using an easy-DNA Kit (Omega) following the manufacturer's instructions.

| PCR and sequencing
The cytochrome b (Cyt b) gene was amplified using polymerase chain reaction (PCR) with primers L14742 and H15915 (Xiao, Zhang, & Liu, 2001). Each 50 μl PCR contained 1-10 ng of template DNA, 5 μl of 10 × PCR buffer, 1 μl of dNTP mix (10 mM), 10 pmol of each primer, and 2 U of rTaq polymerase (TaKaRa). PCRs were conducted in a F I G U R E 1 Map showing the eight sampling locations (black dots), sample sizes, abbreviated as in Table 1, with the relative frequency of the two HapGroups (HapGroup 1 and HapGroup 2) thermal cycler (T100; Bio-Rad) with the following program: one cycle of denaturation at 95°C for 4 min; 30 cycles at 94°C for 30 s, 56°C for 30 s, and 72°C for 90 s; and one final cycle at 72°C for 10 min.
Successful PCR products were separated by electrophoresis on a 1.0% agarose gel and purified using the Agarose Gel Purification Kit (Qiagen).
Several samples with distinct phylogenetic clades were also used to test the existence of cryptic species using COI gene. The primers and PCR information were as Ward, Zemlak, Innes, Last, and Hebert (2005).
Purified products then were sequenced with an ABI PRISM 3730 sequencer. Sequencing primers were the same as those used for PCR amplification. All unique sequences have been deposited in GenBank.
A median-joining haplotype network was constructed in Network v4.6 (http://www.fluxus-engin eering.com/). Phylogenetic analysis of haplotypes was conducted using the neighbor-joining (NJ) method and Bayesian inference (BI). The NJ analysis with the Kimura 2-parameter distance method was carried out in MEGA 7.0 (Kumar, Stecher, & Tamura, 2016). The tree nodes and branch lengths were statistically tested using the bootstrap method with 1,000 replicates and an interior branch test, respectively. BI analysis was carried out using MrBayes v3.1.2 (Ronquist & Huelsenbeck, 2003). Best-fit models for the Bayesian analysis were inferred by hierarchical likelihood ratio tests on the IQ-TREE Web server (http:// iqtree.cibiv.univie.ac.at/). Pairwise F ST values and analysis of molecular variance (AMOVA) were used to assess the population configuration in Arlequin v3.5 (Excoffier & Lischer, 2010). The pairwise F ST values between sites were calculated and assessed for significance by comparison with 10,000 permutations of data.
Migration patterns were estimated using the coalescent-based program MIGRATE-n v3.6.11, which estimates migration rates between two groups of populations-groups with four populations in upper reaches (UPG) and four populations in middle reaches (MPG).
Estimation of parameters in MIGRATE was done using Bayesian approach (Beerli, 2006).
Historical demographic expansions were investigated with neutral tests, such as Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997), and pairwise nucleotide mismatch distributions; all of these tests were implemented in DnaSP 6.0 (Rozas et al., 2017 Harpending's raggedness index (Harpending, 1994;Rogers & Harpending, 1992). To test the significance of the observed data from the simulated data with sudden expansion, 100 bootstrap replicates were applied. If the sudden expansion model was not rejected, the equation τ = 2µt was used to estimate the approximate expansion time, where τ is the date of the growth or decline measured in units of mutational time calculated in Arlequin v3.5 (Excoffier & Lischer, 2010), and u is the mutation rate per sequence and per generation. The u was calculated from u = 2μk, where μ is the substitution rate per site, and k is the number of nucleotides of the analyzed fragment (1,009 bp in this study).

| RE SULTS
A Cyt b dataset of 1,009 bp in size was obtained, and a total of 57  (Table S1).
Both of the phylogenetic topologies of 66 G. filifer mtDNA haplotypes generated with NJ and Bayesian inference were similar, which are presented with two distinct haplotype groups, HapGroup 1 and HapGroup 2 with high confidence values ( Figure S1). The sample distributions in the groups were not consistent with the geographical locations, and the two populations could each be found at all locations (Table S2).
Haplotype median-joining network supported the phylogenetic tree result that all of the mtDNA haplotypes fell into two haplotype groups ( Figure 2). All haplotypes from each group were separated by at least 11 mutation steps, whereas neighbor haplotypes differed by a maximum of four (HapGroup 1) and two (HapGroup 2)

HapGroup1
HapGroup 2 YYB Eighty-five COI sequences with 570 bp in size were also obtained and identified to 11 unique haplotypes. Genetic distance among haplotypes ranges from 0.002 to 0.022 (Table S4). The NJ tree constructed from the haplotypes was shown similar topology as that of Cyt b ( Figure S2).

| Genetic diversity and structure
MtDNA Cyt b gene sequences have been used to study population genetics of some fishes in the Yangtze River, for example, Zacco platypus (Perdices, Cunha, & Coelho, 2004), Leiocassis longirostris (Xiao, Xia, & Ma, 2012;Yang, Xiao, Yu, & Xu, 2012), Leptobotia elongata (Tian, Duan, Wang, Liu, & Chen, 2013), Siniperca chuatsi (Tian et al., 2015), and Saurogobio dabryi (Li, Tang, Yu, & Liu, 2016). Compared to these fishes, the genetic diversity of G. filifer is currently high. Although the F ST test did not detect genetic differentiation among the sampled populations, the pairwise analysis (Table 3) revealed moderate divergence between the YHJ population and the others (except YYB, with the minimum value; Wright, 1978). However, the river distance from YHJ to YYB, which is about 195 km, is longer than that to YJJ (about 105 km) and there are not flow barriers between these populations. It was irrational and could be a sampling error-for sampling, time from YHJ was a week, however one to several months from other sampling sites. The XXJ population also displayed moderate differentiation from the four populations in the upper reaches of the Yangtze River (Table 3)
Second, sympatric intraspecific divergences in mtDNA could be shaped by long-term isolation coupled to food niche and/or reproductive separation (Wimberger, 1994). Such genetic structuring is typically detected with phenotypic differences. Trophic and genetically separate sympatric populations have been reported in salmonid fishes such as Arctic char, brown trout, and whitefish inhabiting the postglacial lakes in the Northern Hemisphere, which are landlocked populations (e.g., Ferguson & Mason, 1981;Gowell, Quinn, & Taylor, 2012;May-McNally, Quinn, Woods, & Taylor, 2015;Power, Power, Reist, & Bajno, 2009;Praebel et al., 2013;Siwertsson et al., 2013). The food content of G. filifer has been examined, and most of its food sources are benthic organisms, such as mosquito larvae, Limnoperna fortunei, and aquatic insects (Wu et al., 2008). No distinct differences in food niche or morphological characteristics were found in the G. filifer population.
Third, distinct lineages in the same site can also be mixed populations from geographical subpopulations through invasion or introduction. Correlations between admixture genetic lineages and invasion events have been observed in sculpins Cottus spp. (Nolte, Freyhof, Stemshorn, & Tautz, 2005), guppies Poecilia reticulata (Lindholm et al., 2005), and zander Sander lucioperca (Eschbach et al., 2014). However, that might not be the case for G. filifer. Lineages were usually detected in invaded habitats, and the proportions of population size in each lineage were varied with the scales and plasticity of the introduced population (Eschbach et al., 2014;Nolte et al., 2005) (Ruskey & Taylor, 2016), and Salmo trutta in two tiny subarctic Swedish lakes (Andersson et al., 2017) and in the Loch Maree catchment, Scotland (Jacobs, Hughes, Robinson, Adams, & Elmer, 2018). This may be the most likely hypothesis for G. filifer. Global homogeneity implied that the two distinct populations came in contact postglacially at one position and then colonized novel environments synchronously. Population expansions (Table 5) occurred in the last glacial period (10,000-70,000 years ago), and temperature and rainfall began to increase at that time in east China (Song, Yu, & Zhu, 1998). Many fishes in the Yangtze River have been reported to experience population bottleneck followed by expansion, such as Leiocassis longirostris (Yang, Xiao, et al., 2012), Squalidus argentatus (Yang, Tang, et al., 2012), Gymnocypris dobulai (Chan, Li, Hu, Liu, & Xu, 2016), and Parabramis pekinensis (Chen et al., 2016). The asymmetric proportion of HapGroup 1 to HapGroup 2 populations could be interpreted as different founder sizes and panmixia after recontact of distinct populations, but nuclear maker is needed to support this hypothesis.

| Concluding remarks
To our knowledge, this is the first report about sympatric genetically populations of fish in the Yangtze River. These fishes are not a cryptic species, but instead represent a secondary contact of distinct glacial refugia. This study highlights that historical and ecological factors play an important role in population patterns for fishes in the Yangtze River. G. filifer is likely not the only example of a sympatric population in the Yangtze River.

ACK N OWLED G M ENTS
We thank the following persons for collecting samples: Huan Dai, Shaoyi Shen, Yan Pu, and Hao Yang. We also thank Zhe Wang for contributing to sequence assembly and for statistical assistance. This work was supported by the Natural Science Foundation of China (Grant No. 31602161) and project for the protection of species from the Ministry of Agriculture and Rural Affairs, the People's Republic of China.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All DNA sequences can be accessed via GenBank accessions: MK050193-MK050258 and MK834299-MK834309.