Geographical differentiation of the Euchiloglanis fish complex (Teleostei: Siluriformes) in the Hengduan Mountain Region, China: Phylogeographic evidence of altered drainage patterns

Abstract The uplift of the Tibetan Plateau caused significant ecogeographical changes that had a major impact on the exchange and isolation of regional fauna and flora. Furthermore, Pleistocene glacial oscillations were linked to temporal large‐scale landmass and drainage system reconfigurations near the Hengduan Mountain Region and might have facilitated speciation and promoted biodiversity in southwestern China. However, strong biotic evidence supporting this role is lacking. Here, we use the Euchiloglanis fish species complex as a model to demonstrate the compound effects of the Tibetan Plateau uplift and Pleistocene glacial oscillations on species formation in this region. The genetic structure and geographical differentiation of the Euchiloglanis complex in four river systems within the Hengduan Mountain Region were deduced using the cytochrome b (cyt b) gene and 10 microsatellite loci from 360 to 192 individuals, respectively. The results indicated that the populations were divided into four independently evolving lineages, in which the populations from the Qingyi River and Jinsha River formed two sub‐lineages. Phylogenetic relationships were structured by geographical isolation, especially near drainage systems. Divergence time estimation analyses showed that the Euchiloglanis complex diverged from its sister clade Pareuchiloglanis sinensis at around 1.3 Million years ago (Ma). Within the Euchiloglanis complex, the divergence time between the Dadu–Yalong and Jinsha–Qingyi River populations occurred at 1.0 Ma. This divergence time was in concordance with recent geological events, including the Kun‐Huang Movement (1.2–0.6 Ma) and the lag time (<2.0 Ma) of river incision in the Hengduan Mountain Region. Population expansion signals were detected from mismatched distribution analyses, and the expansion times were concurrent with Pleistocene glacier fluctuations. Therefore, current phylogeographic patterns of the Euchiloglanis fish complex in the Hengduan Mountain Region were influenced by the uplift event of the Tibetan Plateau and were subsequently altered by paleo‐river transitions during the late Pleistocene glacial oscillations.


| INTRODUCTION
Climatic fluctuations and geological events in the Pleistocene induced an accelerated change in the genetic structure of species and populations (Yan et al., 2013;Yu, Chen, Tang, Li, & Liu, 2014). The advancement and retreat of ice sheets in the Pleistocene led to population divergence and the generation of new lineages, as well as shaping population demographics. Specifically, the activity (or distribution) range of organisms was limited to refuges during the glacial periods, and later dispersed to open habitats during the interglacial periods. This repeated process shaped the geographical distribution of populations and the genetic variation within species, which, in turn, stimulated adaptation and allopatric speciation (Hewitt, 2000(Hewitt, , 2004. Several studies have confirmed that species primarily distributed in the Tibetan Plateau region experienced population expansion after glacial retreat, suggesting that the eastern Tibetan Plateau might have been a refuge during the major Pleistocene glaciations (Qu & Lei, 2009;Qu, Lei, Zhang, & Lu, 2010). Gongga Mountain is located in the eastern region of the Tibetan Plateau, and it is a significant monsoonal maritime glacier center in the Hengduan Mountain Region. Quaternary glaciers remain widespread today, and glacial accumulation landforms are well preserved, due to the repeated glaciation of this region (Thomas, 1997). Glaciation cycles could drive the postglacial expansion of populations, thus shaping patterns in genetic variation (Li et al., 2009).
However, geographical events have also markedly affected genetic differentiation in this region (Yu et al., 2014). The uplift of mountain systems and the formation of river systems could lead to isolating events that result in limited gene flow between populations, which would consequently provide opportunities for genetic diversification and speciation due to genetic drift and natural selection (Che et al., 2010;Streelman & Danley, 2003).
Freshwater fishes that are strictly constrained by drainage systems could provide unique insights into the relationships between current species distributions and the historical evolution of the paleo-environment (Hewitt, 2004;Qi, Guo, Zhao, Yang, & Tangi, 2007). Historic basin connection events, which resulted from geological alterations, might have shaped the genetic structure of the fish population in the Hengduan Mountain Region (Durand, Templeton, Guinand, Imsiridou, & Bouvet, 1999). The main drainage trajectories have changed remarkably since the late Pliocene because the geomorphology has changed extensively (Clark et al., 2004). Researchers have advocated that the paleo-drainage configurations of the main continental East Asian rivers that drain the southeastern section of the Tibetan Plateau were noticeably different to current patterns (Clark et al., 2004;He & Chen, 2006;Qi et al., 2015). The rivers that currently drain the plateau margin were once attributed to a single paleo-Red River, which flowed southward and discharged into the South China Sea (Clark et al., 2004). River capture and reversal events related to the uplift of the Tibetan Plateau led to the subsequent reorganization of this river system into the current-day major river drainage systems. The evolution of distribution patterns of primary freshwater fishes responded to the complex paleo-geographical structure in the Tibetan Plateau, as well as to processes leading to their isolation or interconnection during the uplift event (Hurwood & Hughes, 1998;Montoya-Burgos, 2003;Zhang et al., 2015).
The Hengduan Mountain Region lies on the southeast edge of the Tibetan Plateau, and it is considered as an important biodiversity hotspot because of its unique geological history and complex topography (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000). The geomorphic evolution of this region resulted in the differentiation or isolation of many plant and animal populations (Fan et al., 2012). The Tibetan Plateau uplift resulted in major ecogeographical changes and hydrographic fluctuations; thus, species in the Hengduan Mountain Region represent appropriate models for examining the contributions of climate and geography to contemporary genetic diversification.
The Euchiloglanis fish species complex is composed of two species (E. kishinouyei and E. longibarbatus) that have high morphological similarity. This complex is part of a group of demersal freshwater catfish (Siluriformes: Sisoridae) that is distributed in the upstream region of the Yangtze River basin. The genetic divergence and population structure of the fishes are easily influenced by geographical events because of their weak mobility. Thus, the Euchiloglanis species complex is an ideal subject for investigating how paleo-drainage shifts affect speciation in connection with the historic uplift of the Tibetan Plateau.
However, morphology-based taxonomy or molecular-based phylogeny techniques have previously failed to recognize these fishes correctly (Guo, Zhang, & He, 2004;Zhou, Li, & Thomson, 2011). Consequently, the lack of a robust phylogenetic relationship for these fishes hindered detailed research on biogeography (Guo, He, & Zhang, 2007;Yu & He, 2012), and, hence, our understanding of the evolution of Euchiloglanis species in Southwestern China. In the current study, we considered the Euchiloglanis species distributed in the Hengduan Mountain Region as a species complex, and we attempted to determine the phylogeographical patterns of Euchiloglanis within this region. Based on the hypothesized links between paleo-drainage systems in southeastern Tibet (Clark et al., 2004), we speculated several important geographical separations that might have shaped the patterns of Euchiloglanis distribution in this region.
We analyzed the geographical differentiation of the Euchiloglanis complex in the Hengduan Mountain Region, using complete sequences of mitochondrial cytochrome b (cyt b) gene and 10 microsatellite markers.
Our goals were to: (1) infer the genetic structure and geographical differentiation of populations belonging to the Euchiloglanis species complex throughout the Hengduan Mountain Region; and (2) verify a vicariant speciation hypothesis (i.e. whereby the geographical range is split into discontinuous parts by the formation of a physical or biotic barrier to gene flow or dispersal) based on geological evidence of massive scale paleo-drainage shifts that are related to the uplift of the Tibetan Plateau.

| Sample collection
Samples were collected using fishhooks from 2012 to 2015. In all, 360 samples were gathered from 11 populations across four river systems, with 155 specimens being collected from the Dadu River, 117 from the Yalong River, 62 from the Jinsha River, and 26 from the Qingyi River (Table 1). A photograph of a Euchiloglanis fish specimen sampled from the Dadu River is shown in Figure 1. Detailed information about the sampling sites is presented in Figure 2. All individuals were used for mtDNA amplification, while 192 specimens from seven populations were selected for microsatellite genotyping. For microsatellite analyses, we chose specimens based on the range of the river. If the same river system contained more than two populations, we chose only two of them. Fin or muscle samples were pre-
PCR products were tested via electrophoresis through 1% agarose gels and were purified with the Qiagen Gel Extraction Kit (Qiagen).
Both strands of each product were sequenced using PCR primers.
Ten microsatellite loci (EK7, EK11, EK13, EK17, EK26, EK34, EK35, EK41, EK48, and EK66) were specifically developed for E. kishinouyei and were selected for genotyping analyses (Li, Wang, Zhao, Xie, & Peng, 2014). All loci were fluorescently labeled with FAM dye and amplified, as previously described ). An ABI 3730xl DNA Analyzer with ROX 500 was used as the internal size standard to determine the size of the PCR products. GENEMAPPER version 4.0 software (Applied Biosystems, USA) was used to score the allele designation.

| Phylogenetic analysis and population structure
The cyt b sequences of three Pareuchiloglanis sinensis individuals were amplified and used as outgroups because the species is the sister taxon of Euchiloglanis. Glyptosternon maculatum (DQ192471) was also chosen as an outgroup to avoid any bias. Before reconstructing the phylogenetic trees, an optimal DNA substitution model (GTR + I + G) was obtained based on model-averaged parameters using the Akaike representing burn-in. Posterior probabilities were obtained from the 50% majority rule consensus tree of the remaining topologies. PAUP* version 4.0b10 was used to perform NJ and MP analyses (Swofford, 2002). Nodal support for NJ phylogram was calculated using 1,000 bootstrap replicates. For the MP analysis, a heuristic search strategy was employed with the tree bisection and reconnection branchswapping algorithm, including the random addition of taxa and 1,000 replicates per search. Nodal support for MP trees was evaluated using 1,000 bootstrap replicates. To visualize intraspecific genetic variation within Euchiloglanis better, the haplotype median-joining network for cyt b was performed in NETWORK version 4.6.1 (Bandelt, Forster, & Rohl, 1999).
The genetic structure analyses of populations identified using the microsatellite loci were conducted using the Bayesian clustering analyses (Pritchard, Stephens, & Donnelly, 2000). Admixture models were chosen to assess possible clusters (K value). The lengths of the MCMC iterations were set to 50,000 with a burn-in period of 5,000.
The K value range was set to 1-7, and each K was replicated 20 times.
The most likely K value was chosen according to peak value of the mean log likelihood [Ln P(X/K)] and the Delta K statistic for a given K (Evanno, Regnaut, & Goudet, 2005).

| Divergence time estimation
Divergence times among the detected mitochondrial clades were evaluated in BEAST version 1.8.0 (Drummond, Suchard, Xie, & Rambaut, 2012), using an uncorrelated relaxed molecular clock Bayesian approach, following a lognormal distribution with the GTR + I + G substitution model proposed by JMODELTEST, 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 4 Ma, derived from the capture of Tsangpo by the Brahmaputra River, which occurred before this time (Clark et al., 2004). The second internal time constraint was the divergence of P. sinensis and E. davidi (Peng, Ho, Zhang, & He, 2006 generations and was sampled every 1,000 generations. The first 10% were burn-in. TRACER version 1.5 (Rambaut & Drummond, 2007) was used to test the convergence of the chains to the stationary distribution, which was determined by an effective size (ESS) of more than 200 (Rambaut & Drummond, 2007). Moreover, three analyses with different random seeds were conducted to verify convergence.
The corresponding tree files were merged with LOGCOMBINER1.8.0 (part of the BEAST package). TREEANNOTATOR version 1.8.0 was used to obtain a maximum credibility tree with the annotation of average node ages and the 95% highest posterior density (HPD) interval.
Phylogenetic tree visualization was performed in FIGTREE version 1.4 ).

| Historical demographic analyses
Demographic historical diversification in the population size of the Euchiloglanis complex in the Hengduan Mountain Region was explored using several approaches. Specifically, we completed neutrality tests, including Tajima's D (Tajima, 1989), Fu's Fs tests (Fu, 1997), and R 2 analyses (Ramos-Onsins & Rozas, 2002). Pairwise differences between haplotypes and mismatch distributions were evaluated for each clade using ARLEQUIN. Sum of squares deviations (SSD) and raggedness statistics (Rag) significance values were evaluated with 10,000 permutations (Harpending, 1994). Mismatch distribution and neutrality tests, except R 2 , were calculated in ARLEQUIN, and R 2 was performed in DNASP.
However, mismatch distribution and a neutrality test based on DNA data do not always catch historical signals because they depend only on the segregating sites and haplotype patterns (Fitzpatrick, Brasileiro, Haddad, & Zamudio, 2009 The results were consistent across runs, and a substitution rate of 2% was used in the Euchiloglanis cyt b region. Previously, the mtDNA substitution rate indicated that the speciation of a lacustrine fish species from its riverine ancestor (corrected based on mtDNA substitutions) was 0.02, which sufficiently pre-dated the formation of the lake where speciation likely happened (Ovenden, White, & Adams, 1993;Waters et al., 2007).

| Genetic diversity and population differentiation
The 1  F I G U R E 3 Scatter plots of genetic distance vs. geographical distance (km: kilometer) for pairwise population comparisons inferred from cyt b (a) and microsatellite data (b) with significant divergence being found between any two populations (Table S3). In addition, the Mantel test generated r values of 0.523 (p = .0034) and 0.467 (p = .0461) for mitochondrial and microsatellite data, respectively, when evaluating the genetic diversity and geographical distance in the Euchiloglanis populations (Figure 3).
The mean values of F IS and F IT were 0.0081 and 0.4632, respectively, based on microsatellite data (Table S6).

| Estimation of divergence times
The divergence time analyses indicated that the ingroup diverged  Tajima's D and Fu's Fs values associated with the Dadu River and Yalong River lineages were negative and highly significant (Table 3).   (Clark et al., 2004). In addition, the results showed high haplotype diversity (h = 0.9521 ± 0.0059) and nucleotide diversity (π = 0.01360 ± 0.00059) (h > 0.5 and π > 0.5%); therefore, the high differentiation between haplotypes might be ascribed to secondary contact between differentiated allopatric lineages and to the long evolutionary history of a large, steady population (Grant & Bowen, 1998). A haplotype network based on the cyt b data provided an enhanced visualization of intraspecific genetic variation within Euchiloglanis fishes. The results indicated that there were no shared haplotypes among tributaries or different reaches, suggesting that the Dadu and Yalong groups were completely isolated. Moreover, the Bayesian structure clustering analysis based on microsatellite data verified this pattern.

| Divergence times and historic demography
Geological research has suggested that the evolutionary drainage systems of the Tibetan Plateau are marked by significant changes in paleo-drainage patterns (Clark et al., 2004). The rivers that currently drain the plateau margin were historically a single paleo-Red River that flowed southward and discharged into the South China Sea. However, the river patterns have drastically changed because of nearby river capture and drainage direction reversal (Barbour, 1936;Lee, 1933). In the middle Pliocene, the Jinsha River was insulated, with this isolation stimulating the genetic diversification of its inhabitants at the genus level. However, these populations subsequently expanded to the Yunnan and Sichuan rivers during the uplift event of the Himalayan region. Clark et al. (2004) suggested that the current Dadu River is most likely the product of an ancient river capture with the Anning River. The current Dadu River has a short, sharp segment that runs transversely to the main mountain range, and a relatively large, low-gradient segment that flows parallel to or behind the main mountain range. Moreover, the high terraces of the Dadu River and low longitudinal river gradients on the Anning River potentially define a paleo-longitudinal profile of the paleo-Dadu/Anning River. Thus, the initially south flowing Anning River might have been captured by the high-gradient river that became present-day Dadu River (Barbour, 1936;Wang, 1998 (Zheng et al., 2002).
Apart from the EGP that occurred at 0.5-0.17 Ma, the last glacial period (LGP) occurred at 0.08-0.01 Ma, and the last glacial maximum (LGM) occurred at 0.021-0.017 Ma (Shi, 1998). During the EGP, ice coverage permanently existed at high elevations and middle areas of the Tibetan Plateau (Shi, 2002;Yang, Rost, Lehmkuhl, Zhenda, & Dodson, 2004). The mismatch analyses and neutrality tests detected significant signals of rapid expansion, with R 2 statistics being small and significant, suggesting recent demographic expansion in some parts