Phylogeography of Triplophysa stenura (Nemacheilidae): Responded to the Mid-Pleistocene Climate Transition in the Qinghai-Tibetan Plateau

© 2020 Academia Sinica, Taiwan Open Access Phylogeography of Triplophysa stenura (Nemacheilidae): Responded to the MidPleistocene Climate Transition in the QinghaiTibetan Plateau Jiaxin Hu1,2, Mingdian Liu3,*, and Dekui He1,* 1The Key Laboratory of Aquatic Biodiversity and Conservation, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan 430072, China. *Correspondence: E-mail: [email protected] (He) E-mail: [email protected] (Hu) 2University of Chinese Academy of Sciences, Beijing, China 3Yangtze River Fisheries Research Institute of Chinese Academy of Fishery Science, Scientific Observing and Experimental Station of Fishery and Environment in the Upper and Middle Reaches of the Yangtze River, Ministry of Agriculture, Wuhan, Hubei 430223, China. *Correspondence: E-mail: [email protected] (Liu) Received 18 July 2020 / Accepted 22 October 2020 / Published 30 November 2020 Communicated by Hin-Kiu Mok Pleistocene climatic oscillations have greatly influenced the evolutionary histories and distribution patterns of most extant species. In this study, the phylogeographic patterns and evolutionary history of Triplophysa stenura were investigated. A total of 359 individuals from 19 locations covering the species’ distribution range were collected, and two mitochondrial genes (COI and Cyt b) and the recombination activating protein 1 gene (Rag1) were analyzed. Two deeply divergent lineages, corresponding geographically to the northern and southern Tanggula Mountain, were observed, and shallow divergences were found within the southern and central Qinghai-Tibetan plateau (QTP). The estimated divergence time between the northern and southern Tanggula lineages was around 1.07 Mya. Within the southern Tanggula Mountain lineages, the Nu and Subansiri rivers populations were split about 0.74 Mya, and the southwestern and central QTP populations were divided with the southeastern QTP ones approximately 0.67 Mya. The divergence times of the lineages were matched with two major glaciations in QTP (the Xixiabangma Glaciation, 1.17–0.80 Mya and Kunlun Glaciation, 0.72–0.62 Mya). Together with demographic history analyses, our results highlighted that major glacial cycles during the mid-Pleistocene played a major role in sculpting the distribution pattern of T. stenura, and led to the gene homogenization crossing the drainage populations in the southwestern and central QTP. Key words: The Qinghai-Tibetan plateau, Phylogeography, Climate oscillation, Pleistocene glacial cycles, Triplophysa stenura. BACKGROUND The Qinghai-Tibetan Plateau (QTP), the highest and largest plateau in the world, is the greatest glaciated track outside the polar region in earth, with a total glacial area of ~100,000 km2 (Lehmkuhl and Owen 2005; Yao et al. 2012). The glaciation throughout eastern QTP is strongly controlled by precipitation changes related to oscillation in the South Asian monsoonal system combined with cooling events (Lehmkuhl and Owen 2005). However, the influence of the South Asian monsoon is weakened significantly across eastern mountain belt, leading to a sharp decrease in precipitation from east to west. Thus, the eastern Citation: Hu J, Liu M, He D. 2020. Phylogeography of Triplophysa stenura (Nemacheilidae): responded to the mid-Pleistocene climate transition in the Qinghai-Tibetan Plateau. Zool Stud 59:67. doi:10.6620/ZS.2020.59-67. Zoological Studies 59:67 (2020) doi:10.6620/ZS.2020.59-67 1 © 2020 Academia Sinica, Taiwan mountainous margins of the plateau are much more favorable for glacier development compared to the central and western parts (Dong et al. 2010). The climatic perturbations during the Quaternary have had a profound impact on the distribution pattern and have been ascribed a dominant role in shaping the distribution and genetic attribute of species in the north hemisphere (Hewitt 2000). Species went extinct in situ, survived in refugia or dispersed to new locations during periods of unsuitable climate and then underwent range expansions and colonized new habitations during periods of favorable conditions. Such repeated changes must have left imprints in the evolutionary history of these species (Hewitt 2000 2004; Yang et al. 2008; Wang et al. 2010; He et al. 2016). Phylogeographic analysis, especially the estimation of divergence time, provides tremendous insight into the patterns and underlying causes of the historical diversification of lineages and is a powerful tool for reconstructing the histories of living species (Funk and Wagner 1995; Avise et al. 2000; Wiens and Donoghue 2004; Lemey et al. 2009; Dincă et al. 2011). Recently, phylogeographic analyses showed that distribution and genetic structure of various species in QTP and its adjacent areas were affected by the climatic oscillations or geological barriers (Xia et al. 2018; Ye et al. 2018; Rana et al. 2019). However, the degrees to which geographic barriers affected for the population structure were depended on several factors such as dispersal ability and habitat preference (Qu et al. 2010; Fan et al. 2012; Xu et al. 2014). The highly migratory species (e.g., birds) were likely to be less affect by geographic barriers compared to less motile species (Zhao et al. 2019). Unlike that of terrestrial animals, the dispersal of freshwater fishes depended on direct connections between drainage basins (Roach et al. 2009; De Macedo‐Soares et al. 2010), and separation among drainage basins caused via the terrestrial environment probably led to the intra-species subdivision of freshwater fishes. However, glacial meltwater during the interglacial period triggered the erosion and incision of river systems which probably prompted species and populations exchanges, while expanded ice during the glacial period may have promoted vicariance due to the form of geographic barriers (e.g., glacier), both of which profoundly influenced the distribution and their genetic structure of cold-adapted fish species. Hence, the freshwater fish are likely serve as good candidates for phylogeographic studies in QTP, allowing examination of how broad-scale historical climatic events have influenced their distribution and genetic structure. Earth’s climate underwent a fundamental change between c. 1.25 and 0.70 Mya, marked by a shift of previously dominant 41 kya cycles to lately dominant 100 kya cycles. At the same time, the amplitude of glacial-interglacial climate variations became larger (Elderfield et al. 2012). This fundamental shift, known as the mid-Pleistocene climate transition (MPT), is one of an enigma in the Quaternary climate evolution. This climatic transition also arose in QTP, when the plateau entered cryosphere, and caused the appearance of the oldest glaciation (the Xixiabangma Glaciation, 1.17–0.8 Mya) and maximum glaciation (Kunlun Glaciation, 0.72–0.62 Mya) in QTP during the Quaternary. The Tibetan stone loaches, genus Triplophysa, is one of three principal components of ichthyofauna in QTP, and widely distribute in QTP and its surrounding areas (He et al. 2020). Triplophysa stenura (Nemacheilidae) is a small loach endemic to QTP and adjacent areas (Zhu 1989). It resides in the swift-flowing streams throughout the exorheic rivers in southeast and west QTP and endorheic drainages in north QTP (e.g., tributaries of Kongyu Tso, Tarok Tso and Selin Tso) as well as the Himalayas (e.g., tributaries of Paiku Tso and Dochen Tso). Due to its adaptation to high altitude and widely distribution, T. stenura is tolerant of the complex QTP environment, which makes it an ideal model organism for verifying the effect of QTP’s climate change on the genetic structure. However, prior to this study, genetic surveys of T. stenura was confined to the east QTP and small samples (Ren et al. 2018). In this study, we analyzed the phylogeographic structure and demographic history of T. stenura in QTP and its surrounding areas used three gene segments (COI, Cyt b and Rag1). We collected 359 individuals from 19 locations across the five river basins which almost covered the range of the species (Fig. 1). Thus, the evolutionary history of the species could thoroughly reflect the historic climate events in this region. The aims of this study are to (i) investigate the genetic diversity, phylogeographic structure and demographic history of T. stenura, and (ii) explore the underlying causes of the phylogeographical pattern concerning the mid-Pleistocene climatic oscillation in QTP. MATERIALS AND METHODS Ethics statement The methods involving animals in this study were conducted in accordance with the Laboratory Animal Management Principles of China. All experimental protocols were approved by the Ethics Committee of the Institute of Hydrobiology, Chinese Academy of Sciences. page 2 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan Sampling Specimens were sampled via cast nets, hand nets and shrimp coops (Table 1). Specimens were stored in 8% formaldehyde solution, and were identified according Zhu (1989). Fin clips were obtained and persevered in 95% ethanol (-20℃). Total DNA extraction and sequencing Total genomic DNA was extracted using the standard phenol-chloroform method (Maniatis et al. 1982). Fragments of two mitochondrial genes (cytochrome b, Cyt b; cytochrome oxidase subunit 1, COI) and one nuclear gene (recombination activating protein 1, Rag1) were amplified and sequenced based on He and Chen (2006), except that annealing temperature was 57℃ for Rag1. The corresponding primer sets are given in table S1. Molecular diversity and phylogenetic analyses Sequences were checked with chromatograms using visually in BioEdit 7.0 (Hall et al. 2011) after alignment with Clustal X 2.0 (Larkin et al. 2007). The haplotypes were identified using DnaSP 6.0 (Rozas et al. 2017) and have been deposited in the GenBank library under the Accession Nos MT681283-681640. Haplotype diversity (h), nucleotide diversity (π) were estimated in DnaSP 6.0 (Rozas et al. 2017). The TCS network (Clement et al. 2000) was constructed using PopART based on the haplotypes of Cyt b (Leigh and Bryant 2015). A spatial analysis of molecular variance was conducted using the spatial analysis of molecular variation (SAMOVA) software with the mitochondrial gene fragments (COI + Cyt b). The SAMOVA takes into consideration geographic locations of sampling, and maximizes the proportion of genetic variance among K groups of populations (Dupanloup et al. 2002). As such, this analysis was useful for statistically differentiating between historically isolated groups in the network. SAMOVA was run with several different group structures to determine the maximum value for FCT (the maximized proportion of total genetic variance due to differences between groups). Phylogenetic relationships were reconstructed under maximum likelihood (ML) and Bayesian inference (BI) in CIPRES (Miller et al. 2010) using concatenated sequence data and two other Triplophysa species, T. brevicauda and T. stewarti, were chosen as outgroups (GenBank accession nos. MW086881MW086892). We used PartitionFinder 2.1.1 to identify the appropriate model and partition scheme under the Akaike information criterion (AIC) (Lanfear et al. 2017), which divided our data into two partitions (subset1 = mtDNA; subset2 = Rag1). Both ML and BI were conducted with two Fig. 1. Map of the study area and geographical distribution of sampling sites in this study. The sites are numbered according to table 1. page 3 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan partitions scheme and GTR + G + I substitute model. For ML analysis, a majority consensus tree was obtained with 1000 pseudoreplicates in RAxML 8.2.4 (Stamatakis 2014). For BI, two independent runs with four Metropolis-coupled Monte-Carlo Markov Chains (MCMC) were simultaneously conducted with 3 million generations, and 1/1000 sample frequency in MrBayes 3.2 (Ronquist et al. 2012). The first trees with a standard deviation of split frequencies above 0.01 were discarded as the burn-in phase. Nodal support was assessed by calculating the mean posterior probability (BP) values of each node in the resulting consensus tree after burnin. Phylogenetic tree visualization was performed in FigTree 1.4 (Rambaut 2012). Divergence time and demographic history analyses The divergence times (Common Ancestor heights) with 95% highest posterior density (HPD) among lineages were tentatively estimated based the mitochondrial gene fragments using BEAST 2.4.7 (Bouckaert et al. 2014). To our knowledge, there is no well-dated calibration substitute rate for Tibetan stone loaches. The substitute rates of mtDNA protein coding genes for nemacheilian loaches were usually adopted the evolutionary rate of the European cobitids (0.34% to 0.42% substitution rate per million years) (Zhang et al. 2017) or cyprinids (1.05% substitution rate per million years) (Ren et al. 2018). However, the substitute rates may be differences depending on that of substitution models, taxa and molecular markers; thus they may not all be exactly comparable (He and Chen 2007). Here, we first recovered the time-calibrated phylogeny of Triplophysa and other loaches in China: The substitution rate was estimated by three calibration points for Triplophysa and other loaches using three calibration points: (1) the crown group of Cobitidae and Balitoridae were dated to 27.82–100.5 Mya (Chen et al. 2015); (2) the crown group of Botiidae was prior to 5.33–27.82 Mya; and (3) the crown group of genus Triplophysa was corrected to 7.25–27.82 Mya (Prokofiev 2007). We then adopted the lineage substitution rate of T. stenura (0.61% per million years), this substitution Table 1. Collection localities, drainages, sample size (N), numbers of haplotypes (H), and haplotype diversity (h) and nucleotide diversity (π) of the sampled populations of T. stenura in the Qinghai-Tibetan Plateau. Tajima’s D values and Fu’s Fs values shows in six different drainages. The single asterisk (*) for p ≤ 0.05; double asterisks (**) for p ≤ 0.01 Drainage Localities N H h πLongitude Latitude Altitude Tajima’ DFu’Fs upper Yangtze 40 29 0.98 ± 0.00008 0.0036 ± 0.00535 -1.2 -13.0** 1Budongquan 9 8 0.97 ± 0.00409 0.0033 ± 0.00411 93.97 35.51 4597 2 Tuotuo River 31 26 0.98 ± 0.00012 0.0038 ± 0.00539 92.43 34.21 4542 Nu River 47 39 0.97 ± 0.00021 0.0026 ± 0.00443 -1.2* -21.7** 3 Tsona 16 15 0.99 ± 0.00008 0.0024 ± 0.00347 90.87 32.03 4730 4 Nagchuka 14 14 1.00 ± 0.00073 0.0021 ± 0.00238 91.55 31.5 4532 5 Dza Chu 17 11 0.84 ± 0.00785 0.0012 ± 0.0021 89.84 32.13 4610 Yarlung Tsangpo 125 89 0.98 ± 0.00003 0.0056 ± 0.00870 -1.9** -33.9** 6 Lhasa 58 41 0.97 ± 0.00008 0.0030 ± 0.00422 91.15 29.64 3670 7 Tsetang 19 16 0.98 ± 0.00050 0.0036 ± 0.00434 92.68 29.19 3538 8 Saga 15 14 0.99 ± 0.00079 0.0044 ± 0.00515 85.15 29.32 4474 9 Zhongba 33 30 0.99 ± 0.00011 0.0024 ± 0.00550 83.98 29.76 4570 Brahmaputra 53 35 0.97 ± 0.00016 0.0024 ± 0.00439 -1.54 -0.09** 11 Dochen 22 21 0.99 ± 0.00023 0.0035 ± 0.00427 89.16 28.03 4384 12 Yadong 31 15 0.92 ± 0.00069 0.0010 ± 0.00150 89.03 27.39 4111 Subansiri River 33 13 0.72 ± 0.00633 0.0005 ± 0.00128 -2.1** -6.9** 10 Lhontse 33 13 0.72 ± 0.00633 0.0005 ± 0.00128 92.34 28.41 3940 upper Ganges 26 17 0.96 ± 0.00047 0.0015 ± 0.00282 -0.53 -7.3** 13 Sar 4 3 0.83 ± 0.04948 0.00040 ± 0.0004 87.87 28.17 4275 14 Gangkha 8 6 0.93 ± 0.00711 0.0015 ± 0.00200 86.61 28.58 4300 15 Phung Chu 12 12 1.00 ± 0.00160 0.0019 ± 0.00277 86.45 28.63 4362 16 Paiku Tso 2 2 1.00 ± 0.25000 0.0008 ± 0.00080 85.59 28.77 4592 Endorheic drainages 35 23 0.95 ± 0.00044 0.0018 ± 0.00310 -1.84** -5.5** 17 Kongyu Tso 4 4 1.00 ± 0.03125 0.0025 ± 0.00239 82.3 30.61 4924 18 Khinglo 28 16 0.93 ± 0.00091 0.0009 ± 0.00133 82.21 31.31 4922 19 Tarok Tso 3 3 1.00 ± 0.07407 0.0008 ± 0.00080 83.98 31.05 4572 page 4 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan rate was slow that of Schizothoracins in QTP, but rapid that of European cobitids (Doadrio and Perdices 2005). The substitution model was set GTR + G + I and a strict clock was used. The simulation was run for 60 million generations and sampling every 1,000 generations, and convergence was evaluated in Tracer 1.7.1 (Rambaut et al. 2018). The time tree was gained after 25% burnin in TreeAnnotator v 2.4.7. Time tree visualization was performed in FigTree. Two statistical tests, Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997), were calculated by two mitochondrial genes (COI, Cyt b) fragments with 10,000 pseudoreplicates in DnaSP 6.0 (Rozas et al. 2017). The significance of both values was calculated from 10000 reted as signatures of population expansion. In addition, we performed a mismatch distribution analysis in DnaSP to investigate population expansion signatures. A multimodal distribution of differences between haplotypes is usually found in samples drawn from populations at demographic equilibrium, whereas the distribution is usually unimodal in populations having undergone recent demographic expansion (Excoffier 2004). To estimate past population sizes, Bayesian Skyline Plot (BSP) was used with a strict molecular clock and 20 million generations in BEAST2.4.7. The substitution model and substitution rate were identical with divergence time analysis. ESS and trace plot were checked ensuring the Skyline parameters for acceptable effective sample sizes values (i.e., > 100). The trend plot of fluctuation of population size were reconstructed using 500 points in Tracer 1.7.1. RESULTS Molecular diversity From 359 individuals, 51 haplotypes of COI (666 bp), 147 haplotypes of Cyt b (1140 bp), 61 haplotypes of Rag1 (701 bp), 168 haplotypes of the mitochondrial genes (1806 bp), and 241 haplotypes of the concatenated gene (2507 bp) were identified. Estimates of nucleotide diversity and haplotype diversity for each sampled site and drainage basins are given in table 1. Haplotype diversity averaged over all populations was 0.99 with a range from 0.72 (0.72 ± 0.00633; Lhontse) to 0.99 (0.99 ± 0.00023; Dochen). The highest nucleotide diversity and the second-highest nucleotide diversity were observed in the Yarlung Tsangpo River (h = 0.98 ± 0.00003; π = 0.00323 ± 0.00731), which might have been caused by the large number of samples from this region. In spatial analysis (SAMOVA), the highest FCT (FCT = 77.9) was obtained when populations were grouped into six main clusters (Table 2), the biggest of which comprised nine populations including Khinglo, Saga, Zhongba, Tsetang, Yadong, Tarok Tso, Lhasa, Dochen, and Kongyu Tso. The samples from Yangtze were formed a cluster whereas the two different clusters comprised the three populations from Nu River basin (group 3 and group 4; Table 2). Phylogenetic relationships and divergence time Both ML and BI trees showed a consistency in topology (Fig. 2). Therefore, we only present the ML tree. We found substantial genetic structure within T. stenura, consisting of at least two significantly supported clades (BP = 1, Bootstrap > 90; clades 1–3, clade 4; Fig. 2). The clades 1–3 and clade 4 showed deeply divergent and separated by the intervening Tanggula Mountain. The clade 4 (Bootstrap = 100, BP = 1) geographically was restricted to the upper Yangtze River in the northern of the Tanggula Mountain. In the southern of the Tanggula Mountain, the clade 3 (Bootstrap = 94, BP = 1) composed for haplotypes limited to the Subansiri River (Lhontse) and clade 2 (Bootstrap = 96, BP = 1) consisted of Nu river systems. In the southwestern and central QTP, populations were distributed in two large rivers and three endorheic drainages (Ganges River and the Brahmaputra-Yarlung Tsangpo River; Kongyu Tso, Khinglo and Tarok Tso). The southwestern and central populations QTP were formed clade 1 and displayed a shallow divergence. Dated divergence time for the main mtDNA lineages are concentrated in the Pleistocene and shown with 95% HPD in figure 3. The oldest lineages arose at 2.05 Mya (95% HPD: 1.54–2.60 Mya) in the Yangtze River (clade 4), while clade 3 (Lhontse) diverged with the remaining clades at 1.07 Mya (95% CI: 0.80–1.38 Mya). Approximately 0.86 Mya (95% HPD: 0.65–1.08 Mya), clade 1 and clade 2 were separated with other southern Tanggula lineages. The estimated divergence times within clade 2 and clade 1 were tentatively dated back to 0.74 Mya (95% HPD: 0.52–0.93 Mya), 0.67 Mya (95% HPD: 0.50–0.85 Mya) respectively. The TCS analysis revealed four distinct separate networks (networks YR, NR, SR, SWT; Fig. 4) clearly corresponded to the four clades recovered by the phylogenetic analyses. Within the network YR, the haplotypes geographically were restricted to the upper Yangtze River, corresponding to phylogenetic clade 4. Haplotypes of population on the Nu River basin (network NR; clade 2) and Subansiri River (network SR; clade 3) formed two separate cluster. Network SWT (clade 1) was restricted among the southwestern and page 5 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan central QTP and showed several star-like topologies, indicative of recent expansions at least. Within network SWT, two regionally defined haplogroups (the Brahmaputra-Yarlung Tsangpo and endorheic drainage basin lineages; Ganges River lineages; Fig. 4) could be distinguished. Moreover, many mixed haplotypes were found to be mixed among the SWT haplogroups, suggesting recent gene flow occurred. Demographic history For all of T. stenura populations, Tajima’s D were negative, suggesting that these populations has experienced sudden population expansion, whereas the non-significant value were found in the populations of three different drainages (the Ganges River, Yangtze River and Brahmaputra River; p > 0.05). However, an extremely significant negative values of Fu’s Fs (p < 0.02) indicated tendency of sudden expansion for all of T. stenura populations. These contradiction might have resulted from different sensitiveness levels of statistics. Generally, Tajima’s D can be affected by various factors (e.g., bottleneck effect or founder effect) (Tajima 1989; Nei and Kumar 2000; Schmidt and Pool 2002), but Fu’s Fs is frequently considered as more sensitive to population expansion than Tajima’s D (Fu 1997; Chen et al. 2004). The unimodal distribution showed in the mismatch analysis for four clades indicated a high possibility of historical population expansion. Result of BSP revealed that three clades (clade 1; clade 2; clade 4; Fig. 5) populations showed increment in the effective population size starting approximately 120 thousand years ago (kya), corresponding to the Last Interglacial (Eemian Stage, 130–115 kya). The Subansiri River population (clade 3) had experienced population growth phases since approximately 20 kya. The limited gene flow in lowaltitude areas (Lhontse 3910 m; the lowest altitude among the localities we surveyed) may be the cause of the unsynchronized time of population expansion. DISCUSSION Phylogenetic analyses based on the concatenated dataset of three gene fragments (COI + Cyt b + Rag1) showed distinct phylogeographic patterns of Triplophysa stenura. The most obvious geographical pattern exhibited between northern and southern Tanggula Mountain with a deeply lineages divergence, Table 2. Results from a spatial analysis of molecular variance (SAMOVA) showing F values given different numbers of groupings. Sets of lineages that were combined within the groups are indicated Grouping Populations composition F value 4 groups 1. Budongquan + Tuotuo River FSC = 0.43 FST = 0.84 FCT = 0.72 2. Lhontse 3. Dza Chu + Nagchuka + Tsona 4. Phung Chu + Khinglo + Saga + Zhongba + Tsetang + Gangkha + Yadong + Tarok Tso + Lhasa + Dochen 5 groups 1. Budongquan + Tuotuo River FSC = 0.25 FST = 0.83 FCT = 0.77 2. Lhontse 3. Dza Chu + Nagchuka + Tsona 4. Sar + Phung Chu + Gangkha + Paiku Tso 5. Khinglo + Saga + Zhongba + Tsetang + Yadong + Tarok Tso + Lhasa + Dochen + Kongyu Tso 6 groups 1. Budongquan + Tuotuo River FSC = 0.23 FST = 0.83 FCT = 0.78 2. Lhontse 3. Dza Chu 4. Nagchuka + Tsona 5. Sar + Phung Chu + Gangkha + Paiku Tso 6. Khinglo + Saga + Zhongba + Tsetang + Yadong + Tarok Tso + Lhasa + Dochen + Kongyu Tso 7 groups 1. Budongquan + Tuotuo River FSC = 0.16 FST = 0.81 FCT = 0.77 2. Lhontse 3. Dza Chu 4. Nagchuka + Tsona 5. Yadong 6. Sar + Phung Chu + Gangkha + Paiku Tso 7. Khinglo + Saga + Zhongba + Tsetang + Tarok Tso + Lhasa + Dochen + Kongyu Tso page 6 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan while the southwestern and central QTP populations displayed a shallow divergence (Fig. 2). The previous phylogeographical study of T. stenura based on the two mitochondrial genes (Cyt b; D-loop) revealed that T. stenura is a single lineage with three well supported clades and the divergence occurring at 0.64 Mya (Ren et al. 2018), corresponding to the Kunlun-Huanghe Movement (1.1–0.6 Mya) (Cui et al. 1998). However, the presence of two different phylogeographic patterns in this study suggests that the evolutionary history of T. stenura is more complex than previously assumed. In our study, the southwestern and central QTP populations showed a shallow divergence which did not follow geographic patterns, indicating that other factors may have played a more dominant role than that of tectonic movement in shaping the evolutionary history and genetic divergence of the species. Interestingly, in most studies, tectonic movement was considered the important factor shaping the genetic structure of QTP species during the Pleistocene. For instance, Ding and Liao (2019) suggested that two major genetic differentiations of Cricetulus kamensis that occurred during the early Pleistocene were influenced by the Qing-Zang tectonic movements. Qu et al. (2010), Cun and Wang (2010), Li et al. (2017) and Wang et al. (2010) emphasized that the uplift of QTP and the subsequently induced climate fluctuations shaped the patterns of genetic diversity observed in extant QTP’s species during the Pleistocene. All these studies proposed that the interaction of historic climate changes and geological events shaped the evolutionary history and genetic divergence of the species in QTP. However, recent review suggested that QTP has been 4,000–5,000 m high since the mid-Eocene, and the Indian summer monsoon, South-east Asian summer monsoon, and Central Asian winter monsoon arose at Fig. 2. ML tree for T. stenura on QTP and in its adjacent drainages based on a combination of genes including two mitochondrial genes (COI, Cyt b) and recombination activating protein 1 gene (Rag1-701 bp). Clade credibility values of major lineages are given for nodes with bootstrap support for ML (above branch) and posterior probability for Bayesian inferences (below branch). Major clades referred to in the text are listed to the right. Different colors were assigned for each clade: Clade 1, green; Clade 2, red; Clade 3, black; Clade 4, blue; Outgroups, black. page 7 of 13 Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan different times and are unrelated to the uplift of QTP (Renner 2016). Therefore, the explanation that the current phylogeographic structure was influenced by the uplift of QTP during the Pleistocene does not seem plausible. Thus, factors other than distinct orogenic movements may have played the predominant role in shaping the evolutionary history and genetic divergence of the species. Based on the estimated divergence times obtained from the BEAST analyses, the Yangtze River is likely where T. stenura first arose, approximately at 2.05 Mya (95% HPD: 1.54–2.60 Mya), followed by the southwestward dispersal into the southern Tanggula Mountain (approximately 1.07 Mya). This date well coincides with a cold stage in the Xixiabangma Gaciation (1.17–0.8 Mya) which was the oldest glaciation took place in the late Early Pleistocene. Moreover, the divergences of clades 1 and clade 2 (clade Fig. 3. Bayesian estimates of divergence time for the lineages of T. stenura based on the mitochondrial genes (COI, Cyt b) data set. Dusty purple bars represent 95% highest posterior density for divergence estimates. The numbers on the nodes are million years ago (Mya). page 8 of 13Zoological Studies 59:67 (2020) © 2020 Academia Sinica, Taiwan 1, southwestern and central QTP; clade 2, Nu River basin, southern Tanggula Mountain) are approximately 0.68 Mya and 0.74 Mya, respectively. These two dates are well consistent with the largest glaciation (Kunlun Glaciation 0.78–0.62 Mya). During the most extensive glacier stage of Kunlun Glaciation on the plateau, the area covered by ice was found to be approximately 500,000 km2 (Zheng et al. 2002). Following the advancement of glaciers during these glaciations, some glacial refugia (e.g., river valley) might persist on QTP platform, where cold-adapted species (e.g., Triplohysa) could survive to the frigid climate and further promoted the intraspecific divergences. Taberlet and Cheddadi (2002) revealed that the refugia populations were often been accompanied by high levels of genetic diversity and unique haplotypes. In this study, unique haplotypes and high levels of genetic diversity were found in the Yangtze River, Subansiri River and Nu River. Our data implied that the lineages of T. stenura became restricted to at least three different regions that acted as refugia during the mid-Pleistocene glacial period. A logical explanation is that the private haplotypes existed in the glacial refugia of these regions over a long period and survived in situ. Our data indicated that T. stenura comprises four major clades (clade 1, clade 2, clade 3 and clade 4) and revealed different phylogeographic patterns in QTP. In contrast to northern and southern Tanggula Mountain populations, where populations were strictly nested by drainages (clade 2, 3 and 4 from Nu River, the Subansiri River and Yangtze River, respectively), exhibiting a marked geographical pattern, the southwestern and central QTP populations (clade 1) usually entangled together and lack of a clearly geographic pattern. These different phylogeographic patterns may be due to the asynchronous glacier development in different regions Fig. 4. TCS network generated by Popart based on cytochrome b haplotypes for the T. stenura in QTP. Numbers in the networks represent haplotype designations and the areas of the circles are proportional to haplotype frequency; Black dots represent missing intermediate haplotypes. page 9 of 13Zoological Studies 59:67 (2020)