Molecular systematics of the Triplophysa robusta (Cobitoidea) complex: Extensive gene flow in a depauperate lineage

Accepted Manuscript Molecular systematics of the Triplophysa robusta (Cobitoidea) complex: extensive gene flow in a depauperate lineage Chenguang Feng, Weiwei Zhou, Yongtao Tang, Yun Gao, Jinmin Chen, Chao Tong, Sijia Liu, Kunyuan Wanghe, Kai Zhao PII: DOI: Reference:

S1055-7903(18)30536-0 https://doi.org/10.1016/j.ympev.2018.12.009 YMPEV 6365

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

17 August 2018 4 December 2018 7 December 2018

Please cite this article as: Feng, C., Zhou, W., Tang, Y., Gao, Y., Chen, J., Tong, C., Liu, S., Wanghe, K., Zhao, K., Molecular systematics of the Triplophysa robusta (Cobitoidea) complex: extensive gene flow in a depauperate lineage, Molecular Phylogenetics and Evolution (2018), doi: https://doi.org/10.1016/j.ympev.2018.12.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Target Journal: Molecular Phylogenetics and Evolution

Title: Molecular systematics of the Triplophysa robusta (Cobitoidea) complex: extensive gene flow in a depauperate lineage

Chenguang Fenga,b, Weiwei Zhouc, Yongtao Tanga,b, Yun Gaoc, Jinmin Chenb,c, Chao Tonga,b, Sijia Liua,b, Kunyuan Wanghea,b, Kai Zhaoa,*

a

Key Laboratory of Adaptation and Evolution of Plateau Biota, and Laboratory of Plateau Fish

Evolutionary and Functional Genomics, and Qinghai Key Laboratory of Animal Ecological Genomics, Northwest Institute of Plateau Biology, Chinese Academy of Sciences, Xining 810008, Qinghai, China b

University of Chinese Academy of Sciences, Beijing 100049, China

c

State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese

Academy of Sciences, Kunming 650223, Yunnan, China *

Correspondence:

Kai Zhao Northwest Institute of Plateau Biology, Chinese Academy of Sciences No. 23 Xinning Lu, Xining, Qinghai 810008, China E-mail: [email protected]

ABSTRACT

Gene flow between populations assumed to be isolated frequently leads to incorrect inferences of evolutionary history. Understanding gene flow and its causes has long been a key topic in evolutionary biology. In this study, we explored the evolutionary history of the Triplophysa robusta complex, using a combination of multilocus analyses and coalescent simulation. Our multilocus approach detected conspicuous mitonuclear discordances in the T. robusta complex. Mitochondrial results showed reticular clades, whereas the nuclear results corresponded with the morphological data. Coalescent simulation indicated that gene flow was the source of these discordances. Molecular clock analysis combined with geological processes suggest that intense geological upheavals have shaped a complicated evolutionary history for the T. robusta complex since the late Miocene, causing extensive gene flow which has distorted the molecular systematics of the T. robusta complex. We suggest that frequent gene flow may restrict speciation in the T. robusta complex, leading to such a depauperate lineage. Based on this comprehensive understanding, we provide our proposals for taxonomic revision of the T. robusta complex.

Keywords: Incomplete lineage sorting, Introgression, Phylogeny, Qinghai-Tibetan Plateau, Tibetan loach.

1. Introduction

Geological and climatic events generally cause range shifts of species, providing opportunities for gene flow between formerly diverged populations or even different species (Hewitt, 2000; Pereira et al., 2016; Qu et al., 2012; Yan et al., 2013). Gene flow during secondary contact confounds population structures, obscures species boundaries, and distorts phylogenies, which can lead to erroneous inferences of evolutionary history (Folk et al., 2017; Leavitt et al., 2017; Liu et al., 2010; McLean et al., 2016; Smith et al., 2017; Sonsthagen et al., 2016). Hence, understanding gene flow and its causes is a key topic in evolutionary biology in order to correct for erroneous inferences.

The northeastern part of the Qinghai-Tibetan Plateau (QTP) is an important region for evolution of the plateau’s biota (Qu et al., 2010; Zhang et al., 2005; Zhao et al., 2011), and has been one of the most active regions of tectonic activity in the QTP since the late Miocene (c. 8 Ma; Li et al., 2015). Studies on plants (Zhang et al., 2005), amphibians (Zhou et al., 2012), and birds (Qu et al., 2010) endemic to this region have provided empirical evidence of gene flow between seemingly isolated populations. However, our understanding of the prevalence of gene flow, as well as its primary role across taxa in this region, remains lacking. This is particularly true for fishes.

Various riverine systems have developed in the northeastern QTP region,

including the Hexi River system, the upper reaches of the Huanghe River, the Weihe River and the Jialing River (Fig. 1). The Hexi River system comprises the Shiyang, Heihe, and Shule rivers, listed from east to west. Stratigraphic studies suggest that geological upheaval caused by the uplift of the QTP led to complicated relationships among rivers in this region, with the possibility of multiple exchanges occurring between once-separated rivers (Fang et al., 2005; Feng, 1981; Guo et al., 2006; Li et al., 1996, 2001). Such exchanges have also been documented in the historical record (Zhou, 2010). Additionally, evidence from the ichthyofauna in this region also indicates species diffusion across rivers (Feng et al., 2017a). Frequent river diversions would have provided multiple opportunities for formerly diverged populations to come into contact, resulting in gene flow across species without effective reproductive isolation. As a result, there is likely to be extensive gene flow across the endemic fish populations.

The Triplophysa robusta complex (Cobitoidea: Nemacheilidae), or TRC, is a lineage with a handful species but their interrelationships are controversial, which presents an interesting case study of the issue. The TRC comprises five species: T. robusta, T. hsutschouensis, T. minxianensis, T. pappenheimi, and T. siluroides. All live only in the northeastern part of the QTP (Fig. 1; Wu and Wu, 1992; Zhu, 1989). Triplophysa robusta, T. hsutschouensis, and T. minxianensis can be easily distinguished from T. pappenheimi and T. siluroides morphologically, where caudal peduncle depth is nearly uniform towards the caudal fin in the first three species and

tapered in the latter two. Triplophysa minxianensis is the only species with complete scale coverage and is restricted to the Taohe River (a tributary of the Huanghe River) and the upper reaches of the Weihe River. Triplophysa robusta is either scaleless or shows scale coverage only on the caudal peduncle. It is endemic to the Huanghe River and the upper reaches of the Jialing River. Meanwhile, T. hsutschouensis is scaleless and lives only in the Hexi River system. Since T. robusta and T. hsutschouensis show nearly the same superficial and anatomical morphology, Wu and Wu (1992) argued that the two should be synonymous. The other two sympatric species T. pappenheimi and T. siluroides are endemic to the Huanghe River. At the adult stage, T. siluroides has ring stripes and a comparatively elongated caudal peduncle relative to T. pappenheimi, but the two can be difficult to distinguish during juvenile stages. Nonetheless, there is no dispute about the species boundary between them. More than 100 non-cave dwelling species of Triplophysa have been reported, including the TRC (Froese and Pauly, 2017). Both our preliminary experiment and He et al. (2006) suggest that the TRC is a monophyletic clade and sister to all other non-cave congeners. Compared to its sister clade, the TRC comprises much fewer species. The great disparity in species richness between sibling clades indicates that the TRC represents a depauperate lineage. Unexpectedly, our preliminary experiment with mitochondrial data depicted reticulate clades in the TRC.

These counterintuitive phylogenetic relationships may stem from different causes, such as hybridization, incomplete lineage sorting (ILS), or imperfect taxonomy

(Fontenot et al., 2011; Leavitt et al., 2017; Takahata, 1995). Hybridization and ILS generally lead to mitonuclear discordance (Leavitt et al., 2017; Slatkin and Maddison, 1989; Smith et al., 2017). In our investigation, we employed multilocus analyses to infer the phylogenetic relationships and population structures of the TRC and detect potential mitonuclear discordances. Although hybridization and ILS generally leave very similar phylogenetic signatures, they present different coalescence forms in terms of species divergence (Joly et al., 2009). Thus, we used a coalescent simulation approach involving posterior predictive checks (Joly, 2012) to assess hybridization and ILS. When a hybridization scenario is accepted, we calibrate the phylogeny and date the nodes to assess the influence of geological events on the extent genetic structure of Triplophysa. Finally, based on those comprehensive means, we re-examined the putative cryptic species and re-evaluate the current taxonomy of Triplophysa.

2. Materials and methods

2.1 Samples and laboratory procedures

Triplophysa siluroides and T. pappenheimi are local protected species and rare in the wild (http://www.gansu.gov.cn/art/2007/8/28/art_4843_216013.html). Under the supervision of the fisheries departments of the Chinese government, we sampled the entire range of the TRC (Fig. 1). The Taohe River population of T. minxianensis, which was denoted by T. minxianensis (Taohe), was only detected in one site with 10

individuals. In total, 302 samples of the TRC were used in the analyses (Table S1). In addition, 17 non-TRC congeners were used as outgroups based on previous studies (He et al., 2006; Wang et al., 2016; Table S1).

Total genomic DNA was extracted from fin or muscle using the standard 3-step phenol-chloroform method (Sambrook et al., 1989). Sequences from two mitochondrial genes (cytochrome b, Cyt b; 16S ribosomal RNA, 16S) and three widely used single-copy nuclear genes (rhodopsin, RH1; interphotoreceptor retinoid binding protein, IRBP; myosin heavy polypeptide 6, myh6) were amplified and sequenced. We chose widely used single-copy nuclear genes in our analyses because paralogs may also create reticulate topologies (Schulte et al., 2015; Smith et al., 2017). Primer pairs and PCR conditions are given in Table S2.

Sequences were assembled in SeqMan (DNASTAR, Madison, WI, USA), and were aligned using CLUSTAL W (Thompson et al., 1994) in MEGA v6.0 (Tamura et al., 2013) with default parameters. As per Song et al. (2016), we used the PHI test in Splitstree4 (Huson and Bryant, 2006) and the GARD method (Pond et al., 2006) to search the signatures of recombination in the nuclear sequences. We used the PHASE algorithm (Stephens et al., 2001) in DnaSP v5.1 (Librado and Rozas, 2009) with default parameters to determine the gametic phase of heterozygous nuclear genes. Results with posterior probabilities greater than 85% were retained. Identical haplotypes for mtDNA and alleles of phased nuDNA were collapsed using DnaSP

v5.1 (Librado and Rozas, 2009). All sequences were deposited in GenBank (Table S1).

2.2 Phylogenetic analyses

Cyt b and 16S were concatenated as one mitochondrial locus due to the linked constraint. Population genetic information, such as the number of polymorphic sites (S), the number of haplotypes (H), haplotype diversity (Hd), and nucleotide diversity (π) were computed in DnaSP v5.1 (Librado and Rozas, 2009). We inferred the phylogenetic relationships of the TRC using maximum-likelihood (ML) and Bayesian inference (BI) methods. Gene trees were estimated separately for the mitochondrial loci and concatenated nuclear loci. The optimal partitioning strategy was identified by PARTITIONFINDER 1.1.1 (Lanfear et al., 2012) based on the Bayesian information criterion. ML trees were constructed in RAxML v8.2.9 (Stamatakis, 2014) with the GTRGAMMA model under the optimal partitioning scheme. Statistical support for hypothesized clades was assessed by performing 1000 bootstrap replicates for the ML trees. For the BI analyses, the posterior distributions were obtained by Markov Chain Monte Carlo (MCMC) analysis with one cold chain and three heated chains in MrBayes v3.2.2 (Ronquist et al., 2012) using the same partitioning strategy above. Samples of trees and parameters were drawn every 5000 steps from a total of 50,000,000 MCMC generations, with the first 25 % samples discarded as ‘burn-in’. Convergence was judged according to the average standard deviation of split

frequencies (< 0.01) and the effective sample size (ESS) values (> 200). Both ML and BI analyses were implemented in the CIPRES Science Gateway (http://www.phylo.org/index.php). We performed the approximately unbiased (AU) test (Shimodaira, 2002) using the CONSEL package (Shimodaira and Hasegawa, 2001) to assess prominent conflict between mitochondrial and nuclear topologies.

2.3 Population genetic structure analyses

We used unrooted networks to visualize relationships among haplotypes of mitochondrial locus. The network was constructed using the parsimony software TCS v1.21 (Clement et al., 2000). We used the default 95 % to set the parsimoniously plausible branch connections.

To show population structure in nuDNA data, we conducted model-based clustering analyses in STRUCTURE v2.3 (Pritchard et al., 2000) using the genotype matrix of the three nuclear loci. The number of clusters (K) was set from 2 to 10, each with 10 iterations. In total, 1.2 2

106 MCMC repeats were retained after a burn-in of

105. The best K was identified by ΔK in Structure Harvester (Earl, 2012).

Following this, we used CLUMPP v1.1.2 (Jakobsson and Rosenberg, 2007) to find the optimal alignments of replicate analyses of each K. Finally, DISTRUCT v1.1 (Rosenberg, 2004) was used to generate a bar diagram of population structure.

2.4 Bayesian species delimitation

Analyses based on mitochondrial and nuclear datasets reveal discordant results (AU test, P < 0.001; see Section 3.2). Before testing potential causes of mitonuclear discordance, we need to evaluate the species delimitation in the TRC. An initial species assignment was carried out based on the geographical populations of known species (detailed in Fig. 1). Because samples of the Weihe River population of T. minxianensis were split into two distinct clades in the mitochondrial gene trees (see Section 3.2), these two clades were treated as two candidate species (labeled as ‘Weihe1’ and ‘Weihe2’). Samples of the TRC were finally assigned to 10 candidate species. To test the validity of these candidate species, we employed Bayesian Phylogenetics and Phylogeography (BP&P) v3.3 (Yang, 2015; Yang and Rannala, 2014), which performs well with small datasets (Olave et al., 2014).

Conspicuous mitonuclear discordances were detected (see Section 3.2), with the nuclear results more in line with the morphological data. To avoid errors due to potential gene flow (Zhang et al., 2014), we used only nuclear datasets to infer the guide tree and conduct the BP&P analyses. The guide tree was first constructed using *

BEAST (Heled and Drummond, 2010). To reduce the computational burden, we

randomly resampled each group with 5 samples. We set the substitution parameters for each locus according to the PARTITIONFINDER 1.1.1 (Lanfear et al., 2012) and used the lognormal relaxed clock model. We chose the Yule process as the species tree prior and the piecewise linear with constant root as the population size prior. We ran 200 million generations with sampling frequency every 2000 generations. The first 20%

of the samples were discarded as burn-in, and the convergence of posterior distributions was evaluated using Tracer 1.5 (Rambaut and Drummond, 2007). A maximum-credibility tree was created in TreeAnnotator v1.8.0 (Drummond et al., 2012).

This maximum-credibility tree was used as a guide in the BP&P analysis. Following Leaché and Fujita (2010), we used varying priors for ancestral population size (θ) and the root age (τ0): θ ~ G (1, 10) and τ0 ~ G (1, 10); θ ~ G (2, 2000) and τ0 ~ G (2, 2000); and θ ~ G (1, 10) and τ0 ~ G (2, 2000). The rjMCMC algorithm 0 was implemented, and fine-turn parameter (ε) was set to 10. For each run, 500,000 samples (a sampling interval of 5) were taken after the first 10,000 generations were discarded as burn-in. We ran these twice for each analysis to confirm the consistency between them.

2.5 Testing the scenario of hybridization

To test the causes of mitonuclear discordances (hybridization and/or ILS), we conducted the posterior predictive checks (Joly et al., 2009) in JML 1.3.0 (Joly, 2012). This method is based on the theory that we can reject ILS and accept hybridization as the explanation for the data when the observed minimum distance between two lineages is smaller than 95 % of the expected distance under a no-hybridization scenario (ILS; e.g. *BEAST model).

Under the species delimitation achieved in the BP&P (Section 3.3), we ran *BEAST analysis again using both mitochondrial and nuclear datasets. We chose the piecewise constant as the population size prior and ran 100 million generations with sampling frequency every 2000 generations. The output file containing 50,000 trees from the *BEAST analysis was used as the input file for JML. In the JML analysis, we discarded the first 20 % of trees and ran simulations under the remaining trees using the mitochondrial sequences. The heredity scalar was set to 0.5. Locusrate used the mean value for sequences parameters from *BEAST log file as per the JML manual. Seq-gen command used the model parameters from jModelTest 0.1.1 (Posada, 2008). Finally, 1000 simulations were carried out with the thinning set to 40.

2.6 Time estimation of the most recent common ancestor

Results of the JML analysis supported hybridization as the causes of mitonuclear discordances (see Section 3.4). To understand the background of gene flow, we further conducted the time estimation for the TRC. Conspicuous mitonuclear discordances implied distinct evolutionary histories between mitochondrial and nuclear loci. Hence, we estimated the time of the most recent common ancestor (TMRCA) of the TRC using mitochondrial and nuclear datasets separately.

We employed two validated time-calibration priors to infer the time (Feng et al., 2017b). According to the normal prior, the TMRCA of T. scleroptera in Lake Qinghai and Yellow River was 0.15 Ma, and the TMRCA of T. strauchii in Ili River and

Junggar River was 2.92 Ma. The time estimations were implemented in BEAST v1.8.0 (Drummond et al., 2012) using the lognormal relaxed clock model and the Constant size tree prior. For the mitochondrial dataset, we carried out two independent runs for 500 million generations with sampling frequency every 50,000 generations. For the nuclear dataset, we conducted three independent runs for 200 million generations with sampling frequency every 2000 generations. The first 10 % of these samples were discarded as burn-in and the remaining trees were combined in LogCombiner v1.8.0 (Drummond et al., 2012). The convergence evaluation and the maximum-credibility tree were conducted as described in Section 2.4.

3. Results

3.1 Sequence information

In total, 4345 bp sequence datasets were obtained and aligned from 319 individuals. These comprised five gene fragments: 1140 bp from Cyt b, 1102 bp from 16S, 712 bp from RH1, 734 bp from myh6 and 657 bp from IRBP. Cyt b and 16S were concatenated as one mitochondrial locus. For the TRC, mtDNA had 397 variable sites and 209 resolved haplotypes. A signature of recombination was not detected in the nuclear sequences. RH1, myh6, and IRBP had 18, 22, and 30 variable sites, and resolved 24, 33, and 44 haplotypes, respectively. Details were listed in Table 1.

3.2 Phylogeny and population structure

Both mitochondrial and nuclear analyses suggested the monophyly of the TRC (PP = 1, BS = 98 and 94; Figs. 2 and 3a). The conservative AU test indicated a significant conflict between mitochondrial and nuclear topologies (P < 0.001).

The mitochondrial gene tree recovered four major clades (A – D) in the TRC (Fig. 2). Specifically, clade A (PP = 1, BS = 99) comprised T. pappenheimi, T. siluroides, and most individuals of T. minxianensis (Weihe), and was rooted at the base. Clade B (PP = 1, BS = 100) was made up of only the rest of the individuals of T. minxianensis (Weihe). Individuals of T. minxianensis (Weihe) in clades A and B hereinafter are referred to as T. minxianensis (Weihe1) and T. minxianensis (Weihe2), respectively. Clade C (PP = 1, BS = 91) comprised T. robusta and T. minxianensis (Taohe), and depicted an extremely complex phylogenetic structure between species as well as populations. Clade D (PP = 1, BS = 100) was sister to clade C, and contained only T. hsutschouensis. In addition, T. hsutschouensis from the Shiyang river was separate from T. hsutschouensis from the Heihe and Shule rivers. The TCS networks depicted the detailed structure (Fig. 2).

Intriguingly, the nuclear gene tree suggested a vastly different phylogenetic structure (Fig. 3a). Triplophysa minxianensis (Taohe) was a monophyletic clade (PP = 1, BS = 88), and showed a distant relationship with T. robusta. Triplophysa minxianensis (Weihe1 and Weihe2) was recovered as monophyletic (PP = 0.94, BS = 82), but was not clustered with T. pappenheimi and T. siluroides. According to the

Structure analyses, K = 4 and 7 were two accepted clustering values (Fig. S1). When K = 7, the bar plots showed a subdivided genetic partition where K = 4 (Fig. 3b). The genetic distinctiveness between T. minxianensis (Taohe) and T. minxianensis (Weihe1-Weihe2), T. hsutschouensis (Shiyang) and T. hsutschouensis (Shule-Heihe), even T. robusta (Huanghe) and T. robusta (Jialing) were further supported. The genetic partitions in K = 7 were consistent with the nuclear phylogenetic structure.

In addition to these mitonuclear discordances, both mitochondrial and nuclear analyses failed to delineate the species boundary between T. pappenheimi and T. siluroides.

3.3 Species delimitation

Terminal branches among candidate species were well resolved in the guide tree (Fig. 4). The BP&P analysis suggested seven putative species by speciation probabilities ≥ 0.95 (Fig. 4), including T. hsutschouensis (Shule-Heihe), T. hsutschouensis (Shiyang), T. robusta (Huanghe), T. robusta (Jialing), T. minxianensis (Weihe1-Weihe2), T. minxianensis (Taohe) and T. pappenheimi-T. siluroides.

3.4 Statistical supports of hybridization

Results of the JML analysis demonstrated that hybridization, rather than ILS, best explained the mitonuclear discordances (Table 2). The observed genetic distances between T. robusta (Huanghe) and T. robusta (Jialing), T. robusta (Huanghe) and T.

minxianensis (Taohe), T. robusta (Jialing) and T. minxianensis (Taohe), T. minxianensis (Weihe) and T. pappenheimi-T. siluroides were less than 95 % of the corresponding simulated distances (P = 0.001~0.041).

3.5 Time estimation of the T. robusta complex

The estimated time with 95 % highest posterior density (HPD) was shown in Fig. 5 and Fig. S2. The TMRCA of the TRC was estimated to be around 6.09 Ma (95 % HPD: 3.70~9.26) in mitochondrial analyses, and around 5.17 Ma (95 % HPD: 2.00~9.67) in nuclear analyses. Additionally, mitochondrial and nuclear analyses indicated the time that T. hsutschouensis separated from its sibling species was about 2.15 Ma (95 % HPD: 1.20~3.33) and 2.46 Ma (95 % HPD: 1.44~3.80), respectively. The mitochondrial analyses also inferred a divergence time of about 0.54 Ma (95 % HPD: 0.25~0.85) between T. hsutschouensis (Shiyang) and T. hsutschouensis (Heihe-Shule).

4. Discussion

A comprehensive molecular evaluation of the TRC provided conspicuous mitonuclear discordances (Figs. 2 and 3), occurring mainly between the following four pairs: T. robusta (Huanghe) and T. robusta (Jialing), T. robusta (Huanghe) and T. minxianensis (Taohe), T. robusta (Jialing) and T. minxianensis (Taohe), T. minxianensis (Weihe) and T. pappenheimi-T. siluroides.

4.1 Gene flow explains mitonuclear discordances

The JML analysis suggests that a hybridization model explains all the mitonuclear discordances within these species pairs (Table 2). However, the JML approach may increase the risk of Type I error, thereby reducing the effectiveness of rejecting ILS, in cases where species trees are poorly resolved (Joly et al., 2009; McLean et al., 2016), such as in all species pairs except T. robusta (Huanghe) and T. robusta (Jialing). Generally, lineage sorting in mtDNA occurs about four times faster than in nuDNA (McLean et al., 2016; Palumbi et al., 2001). In our case, because our nuDNA results delineate clear species boundaries (Fig. 2 and 3), it is unlikely that nuDNA has completed lineage sorting while the more rapidly evolving mtDNA has not. Hence, we can confidently reject ILS and accept the hybridization scenario.

The hybridization scenario within different species pairs implies different types of introgression. In pairs T. minxianensis (Taohe) with T. robusta and T. minxianensis (Weihe) with T. pappenheimi-T. siluroides, the introgressions are apparently asymmetric. Mitochondrial genomes of all T. minxianensis (Taohe) are replaced by that of T. robusta, which likely represents recent or sustained mitochondrial introgression events from T. robusta to T. minxianensis (Taohe). In T. minxianensis (Weihe), a considerable part of the population captures their mitochondrial genomes from T. pappenheimi-T. siluroides via historical hybridization events. In contrast, the introgressions between T. robusta (Huanghe) and T. robusta (Jialing) appear to be

bidirectional. The genetic structure between T. robusta (Huanghe) and T. robusta (Jialing) is strongly intertwined (Fig. 2), which indicates secondary admixture. In summary, gene flow is prevalent and explains all mitonuclear discordances in the TRC.

4.2 Continual tectonic events have led to extensive gene flow in the T. robusta complex

The TMRCA of the TRC was at around 5.17 Ma (or c. 6.09 Ma; Fig. 5), coincident with an early phase of intense tectonic activities in the northeastern part of the QTP (Li et al., 2015). We therefore speculate that geological activities at around 5.17 Ma caused the initial divergence within the TRC. Subsequently, frequent geological activities led to further recurrent vicariance and coalescence among populations and species, resulting in gene flow.

Triplophysa minxianensis (Taohe) and T. robusta are sympatric in the Taohe River (Fig. 1). Gene flow between them indicates secondary contact after vicariance, which was presumably prompted by river capture events. Triplophysa minxianensis (Weihe) and T. pappenheimi-T. siluroides are allopatric. The mitochondrial introgression between them demonstrates past exchange between the Weihe River and the upper reaches of the Huanghe River. In addition, the intertwined genetic structure in the pair T. robusta (Huanghe) and T. robusta (Jialing) suggests the exchange between the Huanghe River and the Jialing River. Fossil evidence indicates that the

Taohe River, a tributary of the Huanghe River, was once connected to the Jialing River (c. 0.03 Ma; Guo et al., 2006). Additionally, the study region is a seismically active zone, with strong earthquakes recorded in historical literature (Guo et al., 2006). Many strong earthquakes, such as the Wudu earthquake during the Han Dynasty in 186 BC, have changed river patterns (Zhou, 2010). Such earthquakes would likely have occurred numerous times over millions of years. We therefore suggest that continual tectonic activities have led to a complicated evolutionary process and extensive gene flow in the TRC.

Similar phenomena have also been detected in plants (Zhang et al., 2005), amphibians (Zhou et al., 2012), and birds (Qu et al., 2010) endemic to this region. Interestingly, our study indicates comparatively more extensive gene flow. One potential explanation is that fish populations are more sensitive to tectonic activities than terrestrial populations, because slight tectonic activities may divide rivers without blocking the migration of terrestrial populations. Frequent gene flow is prone to hinder lineage sorting among populations or recently separated species, thereby prolonging or even reversing the speciation process (Webb et al., 2011; Leaché et al., 2014). We therefore speculate that frequent gene flow has restricted speciation in the TRC, resulting in a depauperate lineage.

Triplophysa hsutschouensis, endemic to the Hexi River system, is the only species in the TRC with a simple history. The uplift of the Qilian Mountains separated

the Hexi River system from the Huanghe River (Feng, 1981; Fang et al., 2005). We estimate the time T. hsutschouensis separated from its sibling species (Huanghe River) at about 2.46 Ma (or c. 2.15 Ma; Fig. 5), closely matching previous physiographic studies in the eastern Qilian Mountains (c. 2.4~1.8 Ma: Li et al., 1996; Li et al., 2001). Additionally, the divergence time between T. hsutschouensis (Shiyang) and T. hsutschouensis (Heihe-Shule) (0.54 Ma, 95 % HPD: 0.25~0.85) is coincident with the sympatric Gymnocypris chilianensis (0.19~0.56 Ma; Zhao et al., 2011). Zhao et al. (2011) pointed out that the G. chilianensis has experienced a sequential westward colonization from the Shiyang River to the Shule River, and exhibits a geographic pattern of decreasing genetic diversity. T. hsutschouensis demonstrates the same decreasing genetic diversity pattern (Hd and π; Table 1). These lines of evidence suggest that T. hsutschouensis stems from the Huanghe River and may have experienced the same biogeographic scenario as the G. chilianensis, with sequential westward colonization.

4.3 Taxonomic implications

Ideally, the taxonomy of a group should be consistent with its evolutionary history and confer stability (Burbrink et al., 2000; Chen et al., 2017). However, visible morphology and gene flow frequently mislead the taxonomy. Our findings suggest that a reevaluation of the TRC lineage is advisable.

Due to their extremely similar morphology, the validity of T. hsutschouensis and

T. robusta has been disputed (Wu and Wu, 1992; Zhu, 1989). Wu and Wu (1992) argued that T. hsutschouensis should be a synonym of T. robusta, and even suggested that T. minxianensis was a synonym of T. robusta. Our results suggest that neither T. hsutschouensis nor T. minxianensis is synonymous with T. robusta. T. minxianensis can be distinguished from T. hsutschouensis and T. robusta by the distribution of scales over the body (Zhu, 1989). Additionally, our results indicate that T. minxianensis (Taohe) and T. minxianensis (Weihe) are two distinct species (Figs. 3 and 4), with different evolutionary histories and distinct regional populations. T. minxianensis was first described from the Taohe River and was established based solely on its dense scales (Wang and Zhu, 1979). Therefore, we suggest retaining the name T. minxianensis for T. minxianensis (Taohe), and suggest that T. minxianensis (Weihe) should be redescribed and treated as a new species. Although less morphological divergence is observed between T. hsutschouensis and T. robusta, evolutionary history and biogeographic differentiation suggest that they should constitute separate species. Furthermore, BP&P analysis suggests possible subdivisions of each species: T. hsutschouensis (Shiyang), T. hsutschouensis (Heihe-Shule), T. robusta (Huanghe), and T. robusta (Jialing) (Fig. 4). However, the result of BP&P merely reflects the barriers between populations. To avoid the over-segmentation in species, we suggest a conservative approach in retaining T. hsutschouensis and T. robusta for the above described populations.

Finally, the validity of T. pappenheimi and T. siluroides is challenged by the

BP&P analysis (Fig. 4). This is presumably due to recent speciation or insufficient information in these loci. However, T. siluroides has intraspecific morphological variations (Zhang et al., 1998). No prior studies have examined morphological variation in T. siluroides. All samples of less than 100 mm body length collected in the field appear consistent with the T. pappenheimi phenotype, and it remains difficult to distinguish the two species even in the sub-adult life stage. We cannot preclude the possibility that T. pappenheimi is a life cycle stage of T. siluroides. Hence, we recommend further research should be conducted to provide a more detailed picture in delimiting species boundaries.

5. Conclusions

This study offers a novel examination of the demographic genetics of species belonging to the Triplophysa genus. Intense geological upheavals since the late Miocene in the northeastern part of the QTP have shaped complicated evolutionary processes for the TRC and caused extensive gene flow between seemingly isolated populations. This gene flow distorts the molecular systematics of the TRC, creating conspicuous mitonuclear discordances. Frequent gene flow may restrict speciation within the TRC, resulting in a depauperate lineage. Finally, our analyses suggest a revision for the taxonomy of the TRC is necessary.

Acknowledgements

We extend our sincere gratitude to the editor and two anonymous reviewers for their constructive comments that helped to improve the manuscript. We appreciate Prof. Jing Che for her special help, Prof. Simon Joly, Dr. Narayan Koju for their help in JML, Drs. Ruiying Zhang, Gang Song for their advisements in data handling. We gave our thanks to Prof. Jody Hey, Jared Strasburg for their help in gene flow, and Prof. Ziheng Yang for his crucial aids in BP&P. Prof. Zuogang Peng deserves special thanks for providing spelean loaches samples. We are also grateful to Drs. Guogang Li, Renyi Zhang, Baolin Zhang, and Yongliu Yan for their kindly help and suggestions. This work is supported by grants from the National Natural Science Foundation of China [31870365 and 31572258] and the China Biodiversity Observation and Research Network.

References

Burbrink, F.T., Lawson, R., Slowinski, J.B., 2000. Mitochondrial DNA phylogeography of the polytypic north American rat snake (Elaphe obsoleta): a critique of the subspecies concept. Evolution 54, 2107–2118.

Chen, J.M., Zhou, W.W., Poyarkov, N.A., Stuart, B.L., Brown, R.M., Lathrop, A., Wang, Y.Y., Yuan, Z.Y., Jiang, K., Hou, M., Chen, H.M., Suwannapoom, C., Nguyen, S.N., Duong, T.V., Papenfuss, T.J., Murphy, R.W., Zhang, Y.P., Che, J., 2017. A novel multilocus phylogenetic estimation reveals unrecognized diversity in Asian horned toads, genus Megophrys sensu lato (Anura: Megophryidae). Mol. Phylogenet. Evol. 106, 28–43.

Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659.

Drummond, A.J., Suchard, M.A., Xie, D., Rambaut, A., 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973.

Earl, D.A., 2012. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361.

Fang, X.M., Zhao, Z.J., Li, J.J., Yan, M.D., Pan, B.T., Song, C.H., Dai, S., 2005.

Magnetostratigraphy of the late Cenozoic Laojunmiao anticline in the northern Qilian Mountains and its implications for the northern Tibetan Plateau uplift. Sci. China Ser. D. 48, 1040–1051 (In Chinese).

Feng, C.G., Tong, C., Zhang, R.Y., Li, G.G., Wanghe, K.Y., Tang, Y.T., Zhang, C.F., Zhao, K., 2017a. Biodiversity and distribution patterns of Triplophysa species in the northeastern margin of the Tibetan Plateau. Biodivers. Sci. 25, 53–61.

Feng, C.G., Wu, Y.J., Tian, F., Tong, C., Tang, Y.T., Zhang, R.Y., Li, G.G., Zhao, K., 2017b. Elevational diversity gradients of Tibetan loaches: The relative roles of ecological and evolutionary processes. Ecol. Evol. 7, 9970–9977.

Feng, S.W., 1981. Evolution of Hexi water system in Gansu. J. Lanzhou Univ. 17, 125–129 (In Chinese).

Folk, R.A., Mandel, J.R., Freudenstein, J.V., 2017. Ancestral Gene Flow and Parallel Organellar Genome Capture Result in Extreme Phylogenomic Discord in a Lineage of Angiosperms. Syst. Biol. 66, 320–337.

Fontenot, B.E., Makowsky, R., Chippindale, P.T., 2011. Nuclear–mitochondrial discordance and gene flow in a recent radiation of toads. Mol. Phylogenet. Evol. 59, 66–80.

Froese, R., Pauly, D., 2017. FishBase. World Wide Web electronic publication.,

http://www.fishbase.org/ (Accessed 4 March 2017).

Guo, J., Han, W., Liang, S., Li, Z., 2006. Study and comparison on the terrace of Taohe River and Minjiang River in Minxian–Dangchang area of the Western Qinling. Geol. Surv. Res. 29, 271–278 (In Chinese).

He, D., Chen, Y., Chen, Y., 2006. The molecular phylogeny and biogeography of genus Triplophysa. Prog. Natural Sci. 16, 1395–1404 (In Chinese).

Heled, J., Drummond, A.J., 2010. Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580.

Hewitt, G.M., 2000. The genetic legacy of the Quaternary ice ages. Nature 405, 907–913.

Huson, D.H., Bryant, D., 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267.

Jakobsson, M., Rosenberg, N.A., 2007. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801–1806.

Joly, S., 2012. JML: testing hybridization from species trees. Mol. Ecol. Resour. 12, 179–184.

Joly, S., McLenachan, P.A., Lockhart, P.J., 2009. A statistical approach for distinguishing hybridization and incomplete lineage sorting. Am. Nat. 174, E54–E70.

Lanfear, R., Calcott, B., Ho, S.Y., Guindon, S., 2012. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695–1701.

Leaché, A.D., Fujita, M.K., 2010. Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). P. Roy. Soc. Lond. B Bio. 277, 3071–3077.

Leaché, A.D., Harris, R.B., Rannala, B., Yang, Z., 2014. The influence of gene flow on species tree estimation: a simulation study. Syst. Biol. 63, 17–30.

Leavitt, D.H., Marion, A.B., Hollingsworth, B.D., Reeder, T.W., 2017. Multilocus phylogeny of alligator lizards (Elgaria, Anguidae): Testing mtDNA introgression as the source of discordant molecular phylogenetic hypotheses. Mol. Phylogenet. Evol. 110, 104–121.

Li, J., Fang, X., Ma, H., Zhu, J., Pan, B., Chen, H., 1996. Geomorphic evolution in upper reaches of Yellow River and the uplift of Qinghai-Xizang Plateau in late Cenozoic. Sci. China 26, 316–322 (In Chinese).

Li, J., Fang, X., Pan, B., Zhao, Z., Song, Y., 2001. Late Cenozoic intensive uplift of Qinghai-Xizang Plateau and its impacts on environments in surrounding area. Quatern.

Sci. 21, 381–391 (In Chinese).

Li, J., Zhou, S., Zhao, Z., Zhang, J., 2015. The Qingzang movement: the major uplift of the Qinghai-Tibetan plateau. Sci. China Earth Sci. 58, 2113–2122 (In Chinese).

Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452.

Liu, K., Wang, F., Chen, W., Tu, L., Min, M.S., Bi, K., Fu, J., 2010. Rampant historical mitochondrial genome introgression between two species of green pond frogs, Pelophylax nigromaculatus and P. plancyi. BMC Evol. Biol. 10, 201.

McLean, B.S., Jackson, D.J., Cook, J.A., 2016. Rapid divergence and gene flow at high latitudes shape the history of Holarctic ground squirrels (Urocitellus). Mol. Phylogenet. Evol. 102, 174–188.

Olave, M., Solà, E., Knowles, L.L., 2014. Upstream analyses create problems with DNA-based species delimitation. Syst. Biol. 63, 263–271.

Palumbi, S.R., Cipriano, F., Hare, M.P., 2001. Predicting nuclear gene coalescence from mitochondrial data: the three-times rule. Evolution 55, 859–868.

Pereira, R.J., Martínez-Solano, I., Buckley, D., 2016. Hybridization during altitudinal range shifts: nuclear introgression leads to extensive cyto-nuclear discordance in the fire salamander. Mol. Ecol. 25, 1551–1565.

Pond, S.L.K., Posada, D., Gravenor, M.B., Woelk, C.H., Frost, S.D., 2006. Automated phylogenetic detection of recombination using a genetic algorithm. Mol. Biol. Evol. 23, 1891–1901.

Posada, D., 2008. jModelTest: phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256.

Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945–959.

Qu, Y., Lei, F., Zhang, R., Lu, X., 2010. Comparative phylogeography of five avian species: implications for Pleistocene evolutionary history in the Qinghai-Tibetan plateau. Mol. Ecol. 19, 338–351.

Qu, Y.H., Zhang, R.Y., Quan, Q., Song, G., Li, S.H., Lei, F.M., 2012. Incomplete lineage sorting or secondary admixture: disentangling historical divergence from recent gene flow in the Vinous-throated parrotbill (Paradoxornis webbianus). Mol. Ecol. 21, 6117–6133.

Rambaut, A., Drummond, A., 2007. TRACER version 1.5. Available: http. beast. bio. ed. ac. uk/Tracer.

Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D.L., Darling, A., Höhna, S., Larget, B., Liu, L., Suchard, M.A., Huelsenbeck, J.P., 2012. MrBayes 3.2: efficient

Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542.

Rosenberg, N.A., 2004. DISTRUCT: a program for the graphical display of population structure. Mol. Ecol. Notes 4, 137–138.

Sambrook, J., Fritsch, E.F., Maniatis, T., 1989. Molecular cloning. Cold spring harbor laboratory press, New York.

Schulte, L.J., Clark, J.L., Novak, S.J., Jeffries, S.K., Smith, J.F., 2015. Speciation within Columnea section Angustiflora (Gesneriaceae): Islands, pollinators and climate. Mol. Phylogenet. Evol. 84, 125–144.

Shimodaira, H., 2002. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508.

Shimodaira, H., Hasegawa, M., 2001. CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics 17, 1246–1247.

Slatkin, M., Maddison, W.P., 1989. A cladistic measure of gene flow inferred from the phylogenies of alleles. Genetics 123, 603–613.

Smith, J.F., Clark, J.L., Amaya-Márquez, M., Marín-Gómez, O.H., 2017. Resolving incongruence: Species of hybrid origin in Columnea (Gesneriaceae). Mol. Phylogenet. Evol. 106, 228–240.

Song, G., Zhang, R., Qu, Y., Wang, Z., Dong, L., Kristin, A., Alström, P., Ericson, P.G., Lambert, D.M., Fjeldså, J., 2016. A zoogeographical boundary between the Palaearctic and Sino-Japanese realms documented by consistent north/south phylogeographical divergences in three woodland birds in eastern China. J. Biogeogr. 43, 2099–2112.

Sonsthagen, S.A., Wilson, R.E., Chesser, R.T., Pons, J.-M., Crochet, P.-A., Driskell, A., Dove, C., 2016. Recurrent hybridization and recent origin obscure phylogenetic relationships within the ‘white-headed’gull (Larus sp.) complex. Mol. Phylogenet. Evol. 103, 41–54.

Stamatakis, A., 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313.

Stephens, M., Smith, N.J., Donnelly, P., 2001. A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 68, 978–989.

Takahata, N., 1995. A genetic perspective on the origin and history of humans. Annu. Rev. Ecol. Syst. 26, 343–372.

Tamura, K., Stecher, G., Peterson, D., Filipski, A., Kumar, S., 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729.

Thompson, J.D., Higgins, D.G., Gibson, T.J., 1994. CLUSTAL W: improving the

sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680.

Wang, X.T., Zhu, S.Q., 1979. On a New Species of the Genus Nemachilus in Gansu Province, China. J. Lanzhou Univ. 4, 129–132 (In Chinese).

Wang, Y., Shen, Y.J., Feng, C.G., Zhao, K., Song, Z.B., Zhang, Y.G., Yang, L.D., He, S.P., 2016. Mitogenomic perspectives on the origin of Tibetan loaches and their adaptation to high altitude. Sci. Rep. 6, 1–10.

Webb, W.C., Marzluff, J.M., Omland, K.E., 2011. Random interbreeding between cryptic lineages of the common raven: evidence for speciation in reverse. Mol. Ecol. 20, 2390–2402.

Wu, Y., Wu, C., 1992. The fishes of the Qinghai-Xizang plateau. Sichuan Publishing House of Science & Technology, Chengdu (In Chinese).

Yan, F., Zhou, W.W., Zhao, H.T., Yuan, Z.Y., Wang, Y.Y., Jiang, K., Jin, J.Q., Murphy, R.W., Che, J., Zhang, Y.P., 2013. Geological events play a larger role than Pleistocene climatic fluctuations in driving the genetic structure of Quasipaa boulengeri (Anura: Dicroglossidae). Mol. Ecol. 22, 1120–1133.

Yang, Z., 2015. The BPP program for species tree estimation and species delimitation.

Curr. Zool. 61, 854–865.

Yang, Z., Rannala, B., 2014. Unguided species delimitation using DNA sequence data from multiple loci. Mol. Biol. Evol. 31, 3125–3135.

Zhang, C., Rannala, B., Yang, Z., 2014. Bayesian species delimitation can be robust to guide-tree inference errors. Syst. Biol. 63, 993-1004.

Zhang, D., Zhang, X., Dai, Z., Dei, P., Lu, J., 1998. Morphological change of Triplophysa siluroides in Ningxia. J. Ningxia Agric. Coll. 19, 24–26 (In Chinese).

Zhang, Q., Chiang, T.Y., George, M., Liu, J.Q., Abbott, R.J., 2005. Phylogeography of the Qinghai-Tibetan Plateau endemic Juniperus przewalskii (Cupressaceae) inferred from chloroplast DNA sequence variation. Mol. Ecol. 14, 3513–3524.

Zhao, K., Duan, Z., Peng, Z., Gan, X., Zhang, R., He, S., Zhao, X., 2011. Phylogeography of the endemic Gymnocypris chilianensis (Cyprinidae): sequential westward colonization followed by allopatric evolution in response to cyclical Pleistocene glaciations on the Tibetan Plateau. Mol. Phylogenet. Evol. 59, 303–310.

Zhou, H.W., 2010. The Relationship between the Wudu Earthquake in the Early Han Dynasty and Changes in the Upper Reaches of the Han River Water System. Hist. Res. 4, 49–69 (In Chinese).

Zhou, W.W., Wen, Y., Fu, J.Z., Xu, Y.B., Jin, J.Q., Ding, L., Min, M.S., Che, J., Zhang,

Y.P., 2012. Speciation in the Rana chensinensis species complex and its relationship to the uplift of the Qinghai–Tibetan Plateau. Mol. Ecol. 21, 960–973.

Zhu, S., 1989. The loaches of the subfamily Nemacheilinae in China (Cypriniformes: Cobitidae). Jiangsu Science and Technology Publishing House, Nanjing (In Chinese).

Supplementary material

Supplementary 1 Tables S1 and S2

Supplementary 2 Figs. S1 and S2

Table Table 1 Genetic statistics in the Triplophysa robusta complex

mtDNA

nuDNA

Cyt b

16S

Combined

RH1

myh6

IRBP

Length(bp)

1140

1102

2242

712

734

657

T. hsutschouensis (Shiyang)

S H Hd π

26 22 0.9586 0.0020

14 17 0.8805 0.0015

40 28 0.9954 0.0017

3 5 0.5559 0.0010

8 12 0.8461 0.0027

6 8 0.7441 0.0019

T. hsutschouensis (Heihe)

S H Hd π

13 12 0.8460 0.0018

7 8 0.5126 0.0005

20 17 0.9149 0.0012

2 3 0.1282 0.0002

5 6 0.3948 0.0017

1 2 0.1351 0.0002

T. hsutschouensis (Shule)

S H Hd π

3 4 0.3333 0.0003

10 11 0.7965 0.0010

13 13 0.8398 0.0006

0 1 0 0

4 2 0.4059 0.0022

1 2 0.0455 0.0001

T. pappenheimi

S H Hd π

18 13 0.9476 0.0032

11 9 0.8000 0.0021

29 16 0.9667 0.0026

4 6 0.6609 0.0013

2 3 0.2321 0.0003

6 10 0.8734 0.0038

T. siluroides

S H Hd π

7 6 0.7426 0.0014

2 3 0.2279 0.0002

9 7 0.7647 0.0008

2 3 0.5989 0.0010

0 1 0 0

9 8 0.7379 0.0040

T. robusta (Huanghe)

S H Hd π

64 45 0.9808 0.0055

32 30 0.8889 0.0024

96 54 0.9933 0.0039

4 6 0.2528 0.0005

9 10 0.6597 0.0013

9 14 0.8841 0.0032

T. robusta (Jialing)

S H Hd

40 28 0.9430

21 19 0.7982

61 38 0.9774

3 5 0.6927

9 9 0.6651

8 11 0.8050

π

0.0049

0.0024

0.0037

0.0014

0.0019

0.0036

T. minxianensis (Weihe)

S H Hd π

157 21 0.9005 0.0239

48 23 0.8673 0.0066

205 32 0.9524 0.0154

0 1 0 0

1 2 0.0204 0.0000

5 6 0.1000 0.0002

T. minxianensis (Taohe)

S H Hd π

15 6 0.8667 0.0056

7 3 0.6000 0.0025

22 6 0.8667 0.0035

3 4 0.6474 0.0013

3 2 0.1000 0.0004

1 2 0.4848 0.0007

S indicates the number of polymorphic sites, H indicates the number of haplotypes, Hd indicates haplotype diversity, and π indicates nucleotide diversity.

Table 2 Results of posterior predictive checks in JML analysis. Species comparisons listed are under the P = 0.05 significance level.

Species comparison

Minimum genetic distance

P

0.0004

0.041

0

0.001

T. robusta (Jialing) vs. T. minxianensis (Taohe)

0.0027

0.004

T. minxianensis (Weihe) vs. T. siluroides-T. pappenheimi

0.0031

0.001

T. robusta (Huanghe) vs. T. robusta (Jialing) T. robusta (Huanghe) vs. T. minxianensis (Taohe)

Figure captions

Fig. 1 Distribution area and sampling sites of five species of the Triplophysa robusta complex. Light blue patches and lines depict the drainage system. Colored dashed lines correspond to the approximate range of populations of the five species. Point shapes on the map indicate different species; colors indicate distinct populations. Site numbers are detailed in Table S1. Fish illustrations are derived from Wu and Wu (1992) and Zhu (1989).

Fig. 2 Bayesian inference tree and TCS networks of the Triplophysa robusta complex based on mitochondrial haplotypes. Numbers on branches indicate posterior probabilities (PP) obtained in the BI analyses followed by bootstrap supports (BS) obtained in the ML analyses. Colored bars at the tips of the tree represent the corresponding populations (species). Numbered circles in TCS networks indicate mitochondrial haplotypes in the phylogenetic tree. Circle sizes roughly represent the number of individuals. Black dots refer to missing haplotypes.

Fig. 3 Phylogram and population structures of the Triplophysa robusta complex based on nuDNA. (a) Bayesian inference tree. Numbers on nodes indicate posterior probabilities (PP) followed by bootstrap supports (BS). Colored branches indicate corresponding populations (species). (b) STRUCTURE plots of K = 4 and 7. Different colors represent pre-set clusters (K). Grouping under the figures corresponds to the initial species assignment in BP&P analysis. “Weihe1” and “Weihe2” indicate T.

minxianensis (Weihe) in clade A and clade B in Fig. 2, respectively.

Fig. 4 Bayesian species delimitation results for the Triplophysa robusta complex. The guide tree is recovered by *BEAST analysis assuming 10 species. Posterior probabilities of the tree are shown near the nodes. Numbers in the gray box indicate speciation probabilities under different combinations of priors in BP&P analyses. Numerical values of less than 0.95 are set to red. “Weihe1” and “Weihe2” indicate T. minxianensis in clade A and clade B in Fig. 2, respectively.

Fig. 5 Chronogram of the Triplophysa robusta complex derived from BEAST analysis. The shaded bars on the nodes represent 95 % highest posterior density of the divergence time estimated from mtDNA (purple) and nuDNA (pink). Numbers above the branches give mean values of divergence time from nuDNA; numbers below give values from mtDNA. Red dots indicate time-calibration markers.

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Highlights

1.

Evolutionary history of the Triplophysa robusta complex (TRC) was studied;

2.

Extensive gene flow between isolated populations of the TRC was detected;

3.

Intense geological upheavals may have shaped a complicated history for the TRC since the late Miocene;

4.

Taxonomy of the TRC was reevaluated.

Graphical Abstract