Genetic structure of two sympatric gudgeon fishes (Xenophysogobio boulengeri and X. nudicorpa) in the upper reaches of Yangtze River Basin

Ethics statement

All handling of X. boulengeri and X. nudicorpa specimens was conducted in strict accordance with Animal Experimental Ethical Inspection of Laboratory Animal Centre, Yangtze River Fisheries Research Institute, Chinese Academy of Fishery Sciences (ID Number: FRE0006).

Sample collection and experimental methods

A total of 227 X. boulengeri individuals from five sites, Jiangjin (JJ), YB, Shuifu (SF), Yongshan (YS), and Qianwei (QW), and 126 X. nudicorpa individuals from four sites, JJ, YB, QJ and PZH, were collected between 2011 and 2018 (Table 1; Fig. 1). All samples were identified based on morphological characteristics (Ding, 1994), and small fins were clipped and preserved in 99% ethanol for DNA extraction.

Genomic DNA was isolated by proteinase K digestion followed by a salt extraction method (Aljanabi & Martinez, 1997). Nine polymorphic microsatellite primers (Xb2, Xb3, Xb4, Xb5, Xb11, LT-C5, LT-D2, LT-D7, LT-D8) in X. boulengeri and nine polymorphic microsatellite primers (LT-C5, LT-C6, LT-C7, LT-D1, LT-D2, LT-D3, LT-D7, LT-D8, LT-D9) in X. nudicorpa were used, according to previously reported protocols (Cheng et al., 2012; Zeng et al., 2015). The final products were used for SSR analysis based on capillary electrophoresis fluorescence on the ABI 377 DNA Analyzer (Tianyihuiyuan Bio-Technology Co., Ltd, Wuhan, China), and the results were analyzed by the GeneMarker 1.5 software (Hulce et al., 2011). The Cyt b gene was amplified by polymerase chain reaction (PCR) using primers L14724: (5′-GACTTGAAAAACCACCGTTG-3′) and H15915 (5′-CTCCGATCTCCGGATTACAAG-3′) (Xiao, Zhang & Liu, 2001), and CR using MitDI-F (5′-CACCCYTRRCTCCCAAAGCYA-3′) and MitDI-R (5′-GGTGCGGRKACTT GCATGTRTAA-3′) (Shen et al., 2017a). The PCR amplification was carried out in 50 μl volumes containing 25 μl of mix (Qingke Biological Technology Co., Ltd Wuhan, China), 2 μl of template DNA, 2 μl of each primer (10 mM/l), and 19 μl of ultrapure water. The PCR amplification conditions were performed by first denaturation step at 94 °C for 3 min, followed by 35 cycles at 94 °C for 30 s, at 55 °C for 30 s, and at 72 °C for 30 s, plus a final extension at 72 °C for 8 min. PCR products were sequenced in both directions with the same primers as PCR in an ABI 3730XL sequencer (Tianyihuiyuan Bio-Technology Co., Ltd, Wuhan, China).

Latitude, longitude and altitude were measured by a hand-held sub meter GPS (GEO-XT6000; Nanjing Jun can Instrument Co., Ltd., Nanjing, China).

Statistical analysis of microsatellite data

Genetic diversity indexes of microsatellite loci were calculated and genetic diversity parameters, including the number of alleles, number of effective alleles, observed heterozygosity (Ho) and expected heterozygosity (He) were all detected by POPGENE version 1.32 (Yeh & Boyle, 1997). Duncan’s multiple comparison procedure was used to compare the means for the genetic diversity values between Jinsha River populations (PZH and QJ) and Yangtze River populations (JJ and YB) in X. nudicorpa by SPSS 17.0. The analysis of molecular variance (AMOVA) and pairwise F_ST based microsatellite data were assessed in Arlequin version 3.1.1 (Excoffier, Laval & Schneider, 2005). Benjamini Hochberg was applied to correct the P-values of F_ST whenever multiple tests were performed. Due to the long distance and the geographic isolation of two dams between Jinsha River populations and Yangtze River populations, we divided X. nudicorpa populations into two groups: Group 1 represented the Jinsha River populations (PZH and QJ), and Group 2 represented the Yangtze River populations (JJ and YB). The global F_ST was also calculated between Group 1 and Group 2. Correlations between geographical distances (ln km) and genetic distances (F_ST/(1-F_ST)) based on microsatellite data were tested in GENALEX version 6.5 (Peakall & Smouse, 2012). Geographic distances following the river network were measured in Google Earth (© 2018 Google) (Lisle, 2006). POPGENE was also performed to detect deviations from the Hardy–Weinberg equilibrium (HWE). The inbreeding coefficient (F_is) and allelic richness were obtained by FSTAT version 2.9.3 (Goudet, 1995) and the null alleles were checked by the software Micro-Checker version 2.2.3 (Oosterhout et al., 2004).

Genetic structure analyses of populations identified using the microsatellite loci were conducted. We examined genetic relationships among populations in Structure version 2.3.4 (Pritchard, Stephens & Donnelly, 2000). Structure applies a Bayesian framework to identify the most likely number of cluster in the sample. The delta K method was used to determine the appropriate value of K (Evanno, Regnaut & Goudet, 2005). The lengths of the Markov chain Monte Carlo (MCMC) reps after burnin were set to 1,200,000 with a burn-in period of 200,000. The simulated K values ranged from 1 to 10, and 10 independent runs were used for each K with correlated allelic and without locprior (Wilson et al., 2015). The most likely K value was chosen according to the peak value of the mean log likelihood [Ln P(X/K)] and the Delta K statistic for a given K. The results were summarized on the online platform Structure Harvester (http://taylor0.biology.ucla.edu/struct_harvest/). Genetic relationships among populations were also examined by applying discriminant analysis of principal components (DAPC) (Jombart, Devillard & Balloux, 2010) using the adegenet 2.0.1 package in R 3.2.2 software (R Development Core Team, 2014). For microsatellite data, CONVERT version 1.3 (Glaubitz, 2004) was used to transform the input formats of the following programs: POPGENE, ARLEQUIN, and STRUCTURE.

Statistical analysis of mtDNA sequences

Sequences were aligned using the Clustal X program (Thompson et al., 1997) and pruned using the software MEGA version 6.06 (Tamura et al., 2013). Genetic diversity indexes of mtDNA sequences were calculated. The haplotype (h) and nucleotide (π) diversity, the number of haplotypes and the haplotype frequencies were obtained using the software Dnasp version 5.10 (Librado & Rozas, 2009). To estimate population differentiation, global F_ST and pairwise F_ST were calculated using Arlequin version 3.1.1 (Excoffier, Laval & Schneider, 2005). The global F_ST was also calculated between Group 1 and Group 2 based on Cyt b and CR datasets. To test for signatures of isolation by distance, correlations between geographical distances (ln km) and genetic distances (F_ST/(1-F_ST)) were tested using a MANTEL matrix correlation test in GENALEX version 6.5 (Peakall & Smouse, 2012).

Genetic structure analyses of populations identified using the mtDNA sequences were conducted. To visualize intraspecific genetic variation, median-joining networks of haplotypes were constructed in the software NETWORK version 5.0.0.3 (Bandelt, Forster & Röhl, 1999). Phylogenetic analyses were also conducted. Before reconstructing the phylogenetic trees, optimal X. boulengeri and X. nudicorpa DNA substitution models (TN+F+I; HKY+F+I and HKY+F+G4; HKY+F+I, respectively) were obtained for Cyt b and CR sequences using the Bayesian information criterion in IQ-TREE (Nguyen et al., 2015). X. boulengeri and X. nudicorpa are sister taxon, and they were used as outgroups each other. We downloaded X. boulengeri and X. nudicorpa mtDNA sequences from Genbank to perform phylogenetic analyses (AF375868, KM516103, KM373519, KM255691, NC-025300, KU314698). Bayesian inference was performed to reconstruct the phylogenetic tree among the Cyt b and CR haplotypes in X. boulengeri and X. nudicorpa in MRBAYES version 3.1.2 (Ronquist et al., 2012). Divergence times among the detected mitochondrial clades were evaluated in BEAST version 1.8.0 (Bouckaert et al., 2014) using an uncorrelated relaxed molecular clock Bayesian approach, in addition to a Yule prior approach and a random starting tree. The mean mutation rate was specified as a normal distribution, and estimates were calibrated using two age constraints. One constraint represented an upper bound of 5.3 Ma, the divergence time of Hemibarbus barbus (AB070241) and Hemibarbus labeo (DQ347953) or Gnathopogon elongates (AB218687) and Gobiocypris rarus (NC_018099) based on fossils from the Late Miocene (Zhou, 1990). The second time constraint represented a lower bound of 1.1 Ma, derived from the Kunlun-Yellow River tectonic movement (KP316067, NC_033403, AB239595), which occurred before this time (Cui et al., 1998; He & Chen, 2007). The MCMC chain was run for 1 × 10⁷ generations and was sampled every 10,000 generations. The first 25% were burn-in. TRACER version 1.5 (Goderniaux et al., 2009) was used to test the convergence of the chains to the stationary distribution, which was determined by an effective sample size (ESS) of more than 200. Trees were visualized in interactive tree of life (Ivica & Peer, 2016). The values of transition/transversion (Ts/Tv) were calculated in MEGA version 6.06 (Tamura et al., 2013). Automated barcode gap detection (ABGD) was used for species delimitation using the online version (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) (Eberle, Warnock & Ahrens, 2016). The first significant gap in the distribution of sequence distances beyond intraspecific sequence divergence can thus be used to infer operational taxonomic units that may be related to species. One critical parameter of the ABGD method is the prior maximum divergence of intraspecific diversity (P). Default settings were used for the prior range for maximum intraspecific divergence (0.001, 0.1), corrected distances (JC69) and minimum slope increase (X) of 1.5. Although there is still a lack of consensus of how to interpret discordant ABGD results (Kekkonen & Hebert, 2014), previous studies advocate using a P-value of ∼0.01, which has previously been demonstrated to recover species accurately (Eberle, Warnock & Ahrens, 2016; Puillandre et al., 2012).

Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu & Li, 1993) statistics were examined to infer the neutrality deviant. The mismatch distribution was also performed to evaluate the demographic history in populations of X. boulengeri and X. nudicorpa. Neutrality tests and Mismatch distribution were obtained in Dnasp. The Bayesian skyline plot method was used to infer the population dynamics of X. boulengeri and X. nudicorpa. After sequences were calculated in BEAST version 2.4.3 (Bouckaert et al., 2014), the chain convergence and the skyline plot graphic were visualized in Tracer version 1.5 (Goderniaux et al., 2009) with ESS of more than 200. In addition, Neutrality tests, Mismatch distribution and Bayesian skyline plot were also carried out separately for samples of X. nudicorpa between the Jinsha River and the Yangtze River because of a significant difference in genetic structure between these locations.

Population genetic diversity

A total of 271 X. boulengeri individuals from five sites and 164 X. nudicorpa individuals from four sites were screened across nine microsatellite loci, respectively. For X. boulengeri, the Ho and He per population varied from 0.766 (JJ) to 0.825 (SF) and from 0.852 (JJ) to 0.879 (SF), respectively, (Table 2). Moreover, the global level of genetic diversity for X. nudicorpa was moderately lower than for that of X. boulengeri. For X. nudicorpa, the Ho and He per population varied from 0.482 (QJ) to 0.613 (JJ) and from 0.469 (PZH) to 0.705 (YB), respectively, (Table 2). In sympatric populations (JJ and YB), the level of genetic diversity was similar for X. boulengeri and X. nudicorpa (Table 2).

We obtained 994 bp and 967 bp sequences (after being aligned) for Cyt b in 227 individuals of X. boulengeri and 106 individuals of X. nudicorpa, respectively. For both species, the overall Ts/Tv ratios were 41.014 and 339.736, respectively. A total of 92 X. boulengeri haplotype and 37 X. nudicorpa haplotype sequences were deposited in GenBank (accession numbers MK001561–MK001652, MK001696–MK001732). Overall, haplotype and nucleotide diversities were high for X. boulengeri (h = 0.963, π = 0.00405), but lower for X. nudicorpa (h = 0.718, π = 0.00348) (Table 2). Within the same population (JJ and YB), overall genetic diversities were similar for both species.

Control region sequences amplified from 225 X. boulengeri individuals and 126 X. nudicorpa individuals were distributed into 43 (784 bp fragments) and 23 haplotypes (854 bp fragments), respectively, (Table 2). These sequences were deposited in GenBank (accession numbers MK001653–MK001695, MK001733–MK001755). The CR sequences had 39 and 18 variable sites for X. boulengeri and X. nudicorpa, respectively. The overall Ts/Tv ratio was 7.908 in X. boulengeri and 1.273 in X. nudicorpa. The global haplotype and nucleotide diversities were similar for both X. boulengeri (h = 0.791, π = 0.00201) and X. nudicorpa (h = 0.752, π = 0.00304) (Table 2). Compared to the mtDNA Cyt b marker, the values of haplotype and nucleotide diversities based on the mtDNA CR were lower in samples from all populations among X. boulengeri and X. nudicorpa (excluding PZH and QJ populations) (Table 2).

In terms of single species, five X. boulengeri populations appeared similar in genetic diversity. However, there was a difference in genetic diversity among X. nudicorpa populations. Duncan’s test shown that the genetic diversity in Jinsha river populations were significantly lower than that in Yangtze River populations (Table 2; Table S1).

For SSR dataset, 16 out of 50 HWE tests (32%) displayed a departure in X. boulengeri (Table S2). For X. nudicorpa, five out of 30 (17%) deviated from HWE (Table S3). Estimates of inbreeding coefficients showed that the values were positive in all populations of X. boulengeri (0.061–0.102) and JJ and YB populations of X. nudicorpa (0.121 and 0.174), except the values were negative in QJ and PZH populations of X. nudicorpa (Table 2).

Population genetic structure

For mtDNA and SSR dataset, AMOVA analysis showed that weak signal of geographic structure was detected in X. boulengeri (SSR, F_ST = 0.0003, P = 0.71359; Cyt b, F_ST = 0.0047, P = 0.22092; CR, F_ST = 0.0352, P = 0.00196) (Table 3). High geographic structure was detected in comparison to X. boulengeri. A substantial proportion of genetic variation was related to differences among groups (SSR, 8.85%; Cyt b, 47.45%; CR, 44.03%), with highly significant F_ST values (SSR, F_ST = 0.0612, P < 0.0001; Cyt b, F_ST = 0.5291, P < 0.0001; CR, F_ST = 0.5792, P < 0.0001) (Wright, 1978) (Table 3).

Among geographic populations within species, genetic differentiation was estimated by pairwise F_ST. The pairwise F_ST were lower than 0.05 and significant in few cases for X. boulengeri (Wright, 1978) (Table 4), which indicated weak signal gene differentiation among the five X. boulengeri populations. In contrast, a modest divergence was observed by pairwise F_ST results between JJ and PZH in X. nudicorpa. For X. nudicorpa, the pairwise F_ST values revealed significant differentiation between the Jinsha River (PZH and QJ) and Yangtze River (JJ and YB) populations (Table 5). For mtDNA dataset, the values of pairwise F_ST are larger than that of SSR dataset.

No correlation were detected between geographical distances (ln km) and genetic distances (F_ST/(1−F_ST)) for all the three markers dataset in X. boulengeri. For X. nudicorpa, a significant correlation was detected in the SSR (r = 0.898, P = 0.020) and Cyt b dataset (r = 0.839, P = 0.040), but not in CR (r = 0.906, P = 0.160) (Fig. 2; Table S4 and S5).

For SSR dataset, a structure analysis was also used to detect population structure. For X. boulengeri and X. nudicorpa, the Delta K values were highest (8.02 and 667.01, respectively) when K = 2 (Fig. S1). The structure analysis indicated greater differences among both the Jinsha and Yangtze River populations in X. nudicorpa, but not in X. boulengeri (Figs. 3A and 3B), supporting the results of the global differentiation. The DAPC clusters showed obvious clusters between the Jinsha and Yangtze River populations of X. nudicorpa (Fig. S2B), but not in X. boulengeri (Fig. S2A).

For mtDNA dataset, the haplotype network was performed to detected population structure. The network presented a star shape, dominated by four common haplotypes for the Cyt b dataset and two common haplotypes for the CR dataset, with no geographic aggregation for haplotypes in X. boulengeri (Figs. 4A and 4B). The Cyt b haplotypes of X. boulengeri were divided into two clades (Clade I and Clade II) in network topology. There were five mutational steps between clades (Fig. 4A). Based on CR data, clades could not be observed in the haplotypes network of X. boulengeri (Fig. 4B). The AMOVA analysis also corroborated the divergence between Clade I and Clade II with a high and significant F_ST value (F_ST = 0.4428, P < 0.0001) (Table 3).

For X. nudicorpa, the results of haplotype network showed two clades based on Cyt b (Clade A and Clade B) and CR (Clade C and Clade D) respectively. Clade A and Clade C encompassed individuals from four populations with the most frequent haplotype usually being present in the Jinsha River populations (PZH and QJ), whereas Clade B and Clade D individuals were found only in JJ and YB (Figs. 4C and 4D). For PZH population, 41 out of 44 (93%) and 47 out of 55 (85%) individuals were located in Hnb-1 (Fig. 4C) and Hnc-1 (Fig. 4D), respectively.

The Bayesian inference phylogenetic trees supported the monophyly of Clade I, Clade A and Clade C (Figs. S3A, S3C and S3D), which were consistent with the haplotype networks. Divergence time between Clade I and Clade II was 5.0 Ma from the Cyt b sequence in X. boulengeri (Fig. 5A). For X. nudicorpa, the Cyt b and CR datasets showed that divergence time was not consistent in clades, 5.8 Ma between Clade A and Clade B from Cyt b sequence and 2.4 Ma between Clade C and Clade D from the CR sequence (Figs. 5C and 5D).

Further species delimitation analysis has been performed using ABGD based on the Cyt b and CR datasets, which suggested one species based on recursive partitioning over a range of prior values for maximum intraspecific divergence in X. boulengeri (Puillandre et al., 2012), as well as X. nudicorpa, when P-value is closest to 0.01 (Table S6). For all Cyt b or CR sequences of X. boulengeri and X. nudicorpa, ABGD analysis suggested a total of two species (Table S6), which indicated that no new species has been formed.

Past demographic history

Neutrality tests were used to detect population expansions using Cyt b and CR sequences. For X. boulengeri and X. nudicorpa, Tajima’s D and Fu’s Fs were significant with negative values (except Tajima’s D for the CR dataset in X. nudicorpa) (Table 3), which indicated demographic expansion events. Due to genetic divergence, neutrality tests (Tajima’s D and Fu’s Fs) for the Jinsha River and Yangtze River X. nudicorpa populations were calculated separately. The values of Tajima’s D and Fu’s Fs were significant and negative in Yangtze River populations (except Tajima’s D for the CR dataset), but not significant in Jinsha River populations.

A mismatch distribution and Bayesian skyline plot indicated recent demographic expansion events in X. boulengeri and X. nudicorpa based on the Cyt b and CR datasets. For X. boulengeri, a smooth unimodal mismatch distribution was observed, which was compatible with a single expansion (Rogers & Harpending, 1992). X. nudicorpa showed a bimodal distribution pattern. Jinsha River populations had narrow mismatch distributions that differed from those of the Yangtze River populations, although the mismatch distribution graphs were both unimodal (Fig. S4). Mismatch distributions were very conservative for detecting population growth (Ramos-Onsins & Rozas, 2002). A unimodal mismatch distribution, consisting of only two types of pairwise differences in Jinsha River populations of X. nudicorpa, failed to indicate a recent population expansion. Growth in population size according for Bayesian skyline plot suggested demographic expansion among populations in X. boulengeri and X. nudicorpa (Figs. 6A and 6B). However, we found Jinsha River populations remained at a stable population size (Fig. 6C), while Yangtze River populations had increased (Fig. 6D).

Taking into account a substitution rate of 0.65% per million years in the Cyt b gene of mtDNA (Macey et al., 1998), this episode of demographic expansion of X. boulengeri and X. nudicorpa began separately approximately 20–200 Ka (Fig. 6A) and 75–135 Ka (Fig. 6B), respectively. The Jinsha River X. nudicorpa populations had a stable population size between approximately 0–0.5 Ka (Fig. 6C), and the Yangtze River X. nudicorpa populations had an increasing population size between approximately 125–175 Ka (Fig. 6D).

For the CR dataset, demographic expansions in Yangtze River X. boulengeri and X. nudicorpa populations were consistent with the Cyt b dataset. Taking into account a substitution rate of 3.6% per million years in the CR of mtDNA (Donaldson & Wilson, 1999), these episodes of demographic expansion began separately in X. boulengeri and X. nudicorpa approximately 1–20 Ka (Fig. 6E) and 1–10 Ka (Fig. 6F), respectively. The Jinsha River X. nudicorpa populations remained at a stable size between approximately 0–0.04 Ka (Fig. 6G), and the Yangtze River X. nudicorpa populations had an increasing population size between approximately 2–15 Ka (Fig. 6H).

Population genetic diversity

Genetic diversity forms the basis for species to respond to natural selection (Vrijenhoek, 1994). In the case of mtDNA and SSR markers, the genetic diversity (h and π) of X. nudicorpa was similar to X. boulengeri in JJ and YB populations (Table 2), indicating similar natural selection pressure in sympatric habitats for both species. The upper reaches of the Yangtze River, as the origin center of Gobiobotinae fish, is a sympatric habitat for these two species. In sympatric habitats, similar natural selection pressures can result in similar genetic variation, which can contribute to similar levels of genetic diversity.

For X. boulengeri, the level of genetic diversity was similar among populations, which was consistent with the short distances and absence of any long-term geographical barriers among these five populations. However, different levels of genetic diversity were detected between the Jinsha and Yangtze River X. nudicorpa populations. The low level of genetic diversity in Jinsha River populations (PZH and QJ) might have been caused by the founder effect or genetic drift (Franklin & Frankham, 1998). In addition, X. boulengeri was only rarely present in the QJ-PZH section of the Jinsha River. In contrast, the abundance of X. nudicorpa in the Jinsha River was higher than that in the Yangtze River (Gao et al., 2013). We surmise that X. nudicorpa, particularly individuals with the most frequent mtDNA haplotype, is better adapted to the local environment. However, low genetic diversity would reduce the capacity to cope with environmental change in the Jinsha River populations of X. nudicorpa (Franklin & Frankham, 1998).

Interestingly, the genetic diversity of the Cyt b gene was higher than that of the CR in these two species, except for X. nudicorpa populations in the Jinsha River (Table 2). The ratio of Ts/Tv was also lower in the CR than in the coding Cyt b. This phenomenon is not unique to X. boulengeri and X. nudicorpa, although most species have higher rates of mutation in CR than in Cyt b. Low rates of evolution in the CR have been reported in salmonids (Shedlock et al., 1992), Melanotaenia (Zhu et al., 1994) and Saurogobio (Zhang, 2001).

Population genetic structure

In sympatric habitats, the genetic structure of X. boulengeri and X. nudicorpa followed similar patterns, with weak signal of geographic population structure. Extrinsic factors, such as historical vicariant events and the same environmental features, can promote the development of similar patterns of population structure (McMillen-Jackson & Bert, 2003) observed for X. boulengeri and X. nudicorpa. Moreover, similar spawning patterns and migration behaviors for both species likely also contributed to the similar population structures.

For X. boulengeri, genetic homogeneity among populations indicated high levels of gene flow (Nm) (Table S7), which might be related to the habit of spawning drifting eggs (Liu, Zhou & Zhou, 2012). In addition, similar environment features, such as low altitude and a subtropical climate (except the YS which had a temperate climate) (Table 1), may have some influence on the genetic structure in X. boulengeri (Zhou et al., 2016). The genetic homogeneity among populations of X. boulengeri was congruent with genetic structure reports in other taxa from these regions (Liu et al., 2017; Shen et al., 2017a; Liu, Zhou & Zhou, 2012).

For X. nudicorpa, a significant correlation between genetic and geographic distance in Mantel’s test was observed using the Cyt b and SSR dataset, although the result was not supported by the CR datasets. Hence, population genetic differentiation cannot be explained by the isolation-by-distance model using the present data. The lack of strong evidence in the CR (Fig. 2) further suggests that the present pattern was shaped by a historical connection rather than ongoing gene flow (Hubert et al., 2007). We also cannot conclude that dams have affected the genetic structure of the two fish species. The reasons are that: (1) populations of X. boulengeri above and below the Xiangjiaba Dam did not show distinct divergence; and (2) accumulation of mutational steps between the two clades of X. nudicorpa, corresponding to either side of the dams, need a longer period of time than that lapsed since dam construction (2007 for the Xiluodu Dam and 2008 for the Xiangjiaba Dam). Thus, the current genetic patterns of the fish are likely to be derived from historical events. The Jinsha River was insulated in the middle Pliocene (Clark et al., 2004), which was congruent with the divergence time between the Jinsha and Yangtze River populations (5.8−2.4 Ma). Hence, geographic structure might be related to long-term isolation between the Jinsha and Yangtze River in the middle Pliocene. The climate and environmental heterogeneity between the Jinsha and Yangtze Rivers might also be the factors that lead to differences in genetic diversity and significant genetic structure between the fish populations of the two rivers (Zhou et al., 2016). The Yangtze River is characterized by a typical subtropical climate and the Jinsha River by a temperate climate (QJ) or three-dimensional climate based on the south subtropical zone (PZH) (Table 1).

Hardy–Weinberg Departure (HWD) was mostly due to excess of homozygosity. The micro-checker detected excess of homozygosity for most HWD tests (Table S8 and S9). The positive values of F_is, indicating heterozygote deficit, might be due to the Wahlund effect caused by the subpopulation structure (Dharmarajan, Beatty & Rhodes, 2013). If the sample includes a mixture of individuals from more than one breeding unit, then (on average) Ho will be less than He (Waples, 2015). The phylogenetic trees showed the genetic divergence within populations with low support rate based Cyt b in X. boulengeri. In addition, inbreeding behavior might also be one of the reasons for heterozygote deficit. Population size of X. boulengeri and X. nudicorpa were decreased due to environmental deterioration and intensive fishery exploitation. Moreover, the presence of null alleles may one cause of heterozygote deficit.

The concept of management unit (MU) is useful for identifying and prioritizing conservation units within a species. MUs are recognized as populations with significant divergence of allele frequencies at mitochondrial or nuclear loci, regardless of the phylogenetic distinctiveness of the alleles (Huang et al., 2007; Moritz, 1994). According to the definitions above, a single MU was identified within X. boulengeri, and two separate MUs within X. nudicorpa including Jinsha River Unit and Yangtze River Unit.

Past demographic history

The inferences about demographic history based on Cyt b and CR provide evidence of a recent population expansion in X. boulengeri and X. nudicorpa derived from the Yangtze River. However, the population expansion times obtained by Cyt b and the CR gene were not coincident. The interval might have been caused by the substitution rates of the markers we used. The average substitution rates used in the present study were those of vertebrate mtDNA (0.65% for Cyt b and 3.6% for CR). However, variation of the CR was lower than that of Cyt b in this study, against the common property of vertebrate mtDNA, and therefore the time deduced from the CR may be underestimated. For both X. boulengeri and the Yangtze River populations of X. nudicorpa, the beginnings of the expansion took place in a similar period in the Late Pleistocene. X. boulengeri has undergone a greater long-term expansion period than Yangtze River X. nudicorpa populations. Population expansion for both species has provided evidence that the main stream in the upper reaches of the Yangtze River could have been a refuge for X. boulengeri and X. nudicorpa in the Late Pleistocene. Similar population expansion events might contribute to the similar genetic structure pattern found in X. boulengeri and X. nudicorpa (Silva et al., 2018). Any expansion events might be related to past climatological changes and geological events (Huang et al., 2007). In the Late Pleistocene, a warm climate (Avise, Walker & Johns, 1998) and fluctuating water levels (Yang, 1986) provided population expansion chances for X. boulengeri and the Yangtze River populations of X. nudicorpa. The effects of a warming postglacial climate upon the demography of many species have been well documented in the upper reaches of the Yangtze River Basin, including Leptobotia microphthalma (Shen et al., 2017a), L. elongate (Liu, Zhou & Zhou, 2012), L. rubrilabris (Shen et al., 2017b) and Coreius heterodon (Yuan et al., 2010).

Conversely, the Bayesian skyline plot suggested that Jinsha River populations of X. nudicorpa were stable, and no recent population expansion has occurred. Stable populations should present a lower number of recently evolved mutations than expanding populations (Templeton, 2006). The low frequency of recent mutations, singletons or doubletons, in the Jinsha River populations of X. nudicorpa supports the hypothesis of a stable population, which is consistent with the neutrality tests and Bayesian skyline plot. Most likely, this genetic pattern reflects the distinct demographic history traits of X. nudicorpa between the Jinsha and Yangtze River populations, which can be related to differences in genetic diversity. In Yangtze River, the Pleistocene refuges might have contributed to the colonization for X. nudicorpa by creating opportunities for suitable establishment in stable areas during the climatic fluctuations of the Pleistocene. More individuals could survive and propagate to increase population size and the level of genetic diversity (Hubert et al., 2007).

Divergence time between Jinsha and Yangtze River populations could provide evidence for geographic events. In the middle Pliocene, the Jinsha River was insulated from the Yangtze River (Clark et al., 2004). During the Late Pliocene (≤3.4 Ma), the uplift of the eastern Qinghai-Tibetan Plateau resulted in river capture events, and drainage rearrangements occurred. The Jinsha River flowed into the Yangtze River, and the modern drainage basin morphology was formed (Clark et al., 2004). The two population genetic groups of X. nudicorpa corresponds to the geography of the past rather than modern drainage systems, and likely results from the insulated Jinsha River in the middle Pliocene, and associated river captures (Zhang, Comes & Sun, 2011). The divergence time between the Jinsha and Yangtze River populations falls into the middle Pliocene. This timing broadly agrees with both geological and molecular data, indicating that geographic isolation occurred between the Jinsha and Yangtze Rivers during the middle Pliocene. The similar genetic divergence has been reported in the Euchiloglanis fish complex (Li, Ludwig & Peng, 2017).