In 2023, samples of three snow trout species—S. eurystomus (n = 23), S. fedtschenkoi (n = 12), and Schizothorax sp. (n = 3)—were collected from the Chirchik River (41°37.6553'N, 69°56.6954'E), Sentobsoy stream (40°36.5486'N, 66°40.4471'E), and Qoratog River, a left tributary of the Surkhandarya River, (38°21.5284'N, 68°4.1839'E) (Fig. 1a), respectively. Samples were collected using trammel and gill nets under standardized field protocols. Sampling and storing procedures adhered to the standard fisheries methods outlined by Kottelat and Freyhof (2007) and Jenkins et al. (2014). All the samples were identified based on morphological characteristics (Berg 1949; Turdakov 1963; Mitrofanov et al. 1988) and stored in 96% ethanol at -20 °C. Whole genome resequencing was performed on one individual from each of the three snow trout species sampled in this study. The Biomarker Technologies Corporation performed genomic DNA extraction, DNA library construction, and sequencing. Briefly, genomic DNA was extracted from ethanol-preserved fin clips, DNA library with an insert size of 300–500 bp was constructed for each species and sequenced on the Illumina NovaSeq platform with a 150 bp paired-end strategy. A total of 2 Gb Illumina paired-end short sequencing data were generated and used for mitochondrial genome assembly.

Figure 1. 

The sampling sites and mitochondrial genome maps of the three Schizothorax species. a. The sampling sites of three snow trout species. Amu Darya Basin, Syr Darya Basin, Tarim Basin, and Ili Basin, were shown in different colors; b–d. The mitochondrial genome map of S. eurystomus, S. fedtschenkoi, and Schizothorax sp, respectively.

The mitogenomes of three snow trout species were assembled using GetOrganelle v1.7.7.0 (Jin et al. 2020), applying 2 Gb of Illumina paired-end sequencing data. Multiple k-mer sizes (21, 35, 45, 65, 85, 105, 121, 135, and 145) were used to optimize assembly accuracy, with a minimum coverage threshold of 10× to ensure high-quality mitochondrial genome reconstruction. The resulting assemblies exhibited average base coverage of 195.9×, 126.1×, and 111.4× for S. eurystomus, S. fedtschenkoi, and Schizothorax sp., respectively. Consequently, the assembly was annotated using MitoAnnotator v3.95 (Zhu et al. 2023). The annotated tRNA genes were confirmed, and the secondary structure was predicted with the tRNAscan-SE search server v2.0 (Chan et al. 2021). The software MEGA11 (Tamura et al. 2021) was used to analyze the nucleotide composition and relative synonymous codon usage (RSCU). To understand the base compositional bias of the mitochondrial genome, AT skew and GC skew were estimated using the following formulae: AT skew = (A-T)/(A+T) and GC skew = (G-C)/(G+C) (Perna and Kocher 1995).

To clarify the phylogenetic positions of our three snow trout species within the genus Schizothorax, we downloaded the complete mitogenomes of thirty snow trout species and six outgroup species of Percocypris pingi, Gymnodiptychus dybowskii, Onychostoma macrolepis, Luciobarbus brachycephalus, Hypsibarbus vernayi, and Pethia ticto from NCBI (Table 1). Thirteen PCG sequences from thirty-nine individuals representing thirty-seven species were aligned using MUSCLE v3.8.31 (Edgar 2004) with default parameters. Phylogenetic relationships were inferred using the methods of maximum likelihood (ML) and Bayesian inference (BI). The ML phylogenetic tree was constructed with IQ-TREE v2.2.6 (Minh et al. 2020) with the following parameters: “-m MFP+MERGE -B 1000”. Using these parameters, the best-fit substitution model models (including FreeRate heterogeneity models) and partition schemes were inferred via the built-in ModelFinder (Kalyaanamoorthy et al. 2017). “MFP” (ModelFinderPlus) allows for extended model selection followed by tree inference, while “-B 1000” ensures that 1000 bootstrap searches will be performed to infer the consensus trees. The BI analysis was performed in MrBayes v3.2.7a (Ronquist et al. 2012), with the data partitioned by the gene. Before BI analysis, the best-fit model (GTR+I+G) was selected by Akaike Information Criterion (AICc) with jModeltest v2.1.10 (Darriba et al. 2012). Then, BI analysis was carried out using four simultaneous Markov Chain Monte Carlo chains for 1,000,000 generations and sampled every 100 generations, using a burn-in of 25% of the generations. The average standard deviation of split frequencies was < 0.01, showing that the sampling of the posterior distribution was enough. Phylogenetic trees generated from both ML and BI methods were visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

The divergence times of Schizothorax fishes were estimated using MCMCTree within the PAML v4.10 package (Yang 2007), employing the 13 mitochondrial protein-coding genes and ML tree topology as the input files. We applied the HKY85 substitution model with the overall substitution rate (rgene_gamma) and rate-drift parameter (sigma2_gamma), which were set as G (1, 2) and G (1, 10), respectively. The MCMC run was first executed for 20,000 iterations as burn-in, followed by an additional 100,000 generations with a sample frequency of 10. We used four calibration points: > 0.09 and < 0.20 million years ago (Mya) for the ancestor of S. davidi and S. grahami; > 2.30 and < 3.23 Mya for the ancestor of S. yunnanensis paoshanensis and S. gongshanensis; > 0.07 and < 0.17 Mya for the ancestor of S. wangchiachii and S. oconnori; and > 8.13 and < 11.43 Mya for the common ancestor of Percocypris pingi alongside all Schizothorax species. These calibration points were derived from the study by Li and Guo (2020). We checked convergence by running all analyses twice and using Tracer v1.7.2 (Rambaut et al. 2018) to monitor effective sample sizes (ESS). ESS values for all parameters over 200 indicated sufficient sampling and good convergence.

The partial fragments of 652 bp from the 5’ end of the mitochondrial COI gene were retrieved from the mitochondrial genomes of three newly sequenced species, alongside four previously published mitochondrial genomes of Schizothorax eurystomus, S. fedtschenkoi, Schizothorax sp. (Sheraliev and Peng 2021), Schizothorax argentatus (Wang et al. 2022) and S. pseudoaksaiensis (Bao et al. 2016). These COI barcode sequences were then aligned, and the pairwise genetic distances between or within species were calculated using the Kimura 2-Parameter (K2P) model (Kimura 1980) implemented in MEGA11 software.

The complete mitochondrial genomes of S. eurystomus, S. fedtschenkoi, and Schizothorax sp. were successfully assembled, revealing a consistent structure across these species. Each genome is a full DNA molecule, with lengths of 16,588 bp, 16,584 bp, and 16,592 bp (Fig. 1b–d), respectively. The genomes contain 37 genes: 13 PCGs, 22 tRNAs, and 2 rRNAs, and a non-coding control region (CR), following the typical vertebrate mitochondrial organization (Fig. 1b–d and Suppl. material 1: table S1). The nucleotide composition of the mitochondrial genome is nearly identical across species, with A (29.4–29.8%), T (25.3–25.6%), G (17.9–18.1%), and C (26.8–27.1%) (Suppl. material 1: table S2). The lengths of intergenic spacers varied from 1 to 34 bp, with notable long intergenic spacers between tRNA^Asn and tRNA^Cys (34 bp) and tRNA^Asp and COI (13 bp) (Suppl. material 1: table S1). The thirteen PCGs across the three species are highly conserved. The A+T content for the PCGs is similar, ranging from 54.2% to 54.6% (Suppl. material 1: table S2). Both the AT skew values and GC skew values for the PCGs of three species are nearly identical, with AT skew values of -0.007, -0.007, and 0, and GC skew values of -0.220, -0.218, and -0.211 for S. eurystomus, S. fedtschenkoi, and Schizothorax sp., respectively (Suppl. material 1: fig. S1, table S2). The start codon of most PCGs was ATG, but COI started with GTG (Suppl. material 1: fig. S2). The standard stop codon (TAA) was used by ND1, COI, ND4L, ND5, and ND6, and two incomplete codons (T and TA) were used by ND2, COII, ATPase6, COIII, ND3, ND4, and Cyt b (Suppl. material 1: table S1). The 13 protein-coding genes expressed 5,242, 5,018, and 5,111 amino acid triplets, respectively, excluding the stop codon (Suppl. material 1: table S3). These mitochondrial genomes also have one 12S rRNA and one 16S rRNA gene, located between tRNA^Phe and tRNA^Val, and between tRNA^Val and tRNA^Leu(UUR (Fig. 1b–d and Suppl. material 1: table S1), respectively. The size of the 12S rRNA were 956 bp, 956 bp, and 955 bp, and the size of the 16S rRNA were 1677 bp, 1678 bp, and 1679 bp in S. eurystomus, S. fedtschenkoi, and Schizothorax sp., respectively (Suppl. material 1: table S1). In Schizothorax sp., the A+T content of the rRNA genes was 54%, higher than in S. eurystomus (46.3%) and S. fedtschenkoi (46.2%) (Suppl. material 1: table S2). All three species showed a positive AT skew (0.270, 0.260, and 0.263, respectively) of the rRNA and a negative GC skew (Suppl. material 1: fig. S1, table S2), indicating higher A and C content than T and G content in the rRNAs (Suppl. material 1: table S2). The tRNA genes ranged in length from 67 (tRNA^Cys) to 76 bp (tRNA^Lys) (Suppl. material 1: table S1). Of these, 14 tRNA genes were located on the heavy (+) strand, and 8 tRNA genes were located on the light (–) strand (Fig. 1b–d and Suppl. material 1: table S1). All the tRNAs displayed the clover-leaf secondary structure (Suppl. material 1: fig. S3), excluding the tRNA^Ser(AGY), which cannot form a stable secondary structure because of the lack of the DHU (dihydrouridine) arm (Suppl. material 1: fig. S3). This loss was reported in many fishes, such as Alcichthys elongatus, Linichthys laticeps, and Sebasticus marmoratus (Jia et al. 2020; Patil et al. 2023; Zhang et al. 2023).

We used 13 mitochondrial PCGs of 31 Schizothorax species and six outgroup species to build the BI and ML trees to determine phylogenetic relationships for the three snow trout. The BI and ML trees were resolved into five distinct clusters (Fig. 2a), supporting the results of previous studies (He and Chen 2006; Yang et al. 2012) and confirming the consistency of the observed phylogenetic patterns. The phylogenetic tree showed snow trout from the Tien Shan Mountains group into two distinct clades, I and II. Clade I comprises a pair of sister species, S. fedtschenkoi and Schizothorax sp. (Fig. 2a), primarily distributed in the western Tien Shan Mountains, specifically in the right tributaries of the Amu Darya River. Clade II includes species from other regions of the Tien Shan Mountains (Fig. 2a).

Figure 2. 

Maximum likelihood and Bayesian inference phylogenetic tree of 33 Schizothorax fishes based on 13 protein-coding genes. a. Maximum likelihood phylogenetic tree. The numbers at the nodes indicate the bootstrap value for the maximum likelihood tree and the posterior probability for the Bayesian inference tree, respectively. Tip labels in bold indicate those for which whole genome sequencing was conducted in this study; b. Genetic distances calculated based on DNA barcode sequences. The bolded species names indicate the sequenced species or individual in this study, while the regular species names indicate the species or individual for which NCBI data was obtained. SE, S. eurystomus; SF, S. fedtschenkoi; Ssp, Schizothorax sp.; SE, S. eurystomus KY436758; SF, S. fedtschenkoi (not displayed in the tree); Ssp, Schizothorax sp. (not displayed in the tree); SP, S. pseudoaksaiensis KM079630; SA, S. argentatus SRR26304028.

Interestingly, our findings reveal that two individuals of S. eurystomus do not form a monophyletic group. One individual of S. eurystomus from the Tarim Basin (KY436758) and the individual of S. biddulphi (OR812523) from the same basin are sister taxa, jointly forming a clade with one individual of S. eurystomus sampled from the western Tien Shan Mountains. Although incomplete lineage sorting and hybridization might both account for this observation, we propose that interspecific hybridization between S. eurystomus and S. biddulphi is likely the main driver, given that reproductive isolation is not absolute among congeneric fishes (Campton and Utter 1985; Müller et al. 2010; Ma et al. 2018) and the two species share not only overlapping habitats but also breeding seasons (Guo et al. 2012). Further supporting our hypothesis, a study based on the COI gene revealed significant gene flow among four snow trout species (S. eurystomus, S. biddulphi, S. irregularis, and S. barbatus) in the major tributaries of the Tarim River (Haysa et al. 2021). Future research based on population genomics will provide deeper insights into the hybridization history of these fish species, enhancing our understanding of their evolutionary dynamics.

Notably, the phylogenetic analyses distinctly identified Schizothorax sp. as an independent species within the evolutionary tree, underscoring its unique taxonomic status. The K2P distance between Schizothorax sp. and other valid species ranges from 3.18% to 6.72%, while the distance between our sequenced Schizothorax sp. individual and that from previous study (Sheraliev and Peng 2021) is 0 (Fig. 2b). These results confirm its status as a valid independent species, meeting the established criterion of a genetic distance greater than a 2% interspecific distance and less than a 1% intraspecific distance, as outlined by Hubert and Hanner (2015). However, since the K2P genetic distance calculations did not account for variability, there may be issues with insufficient statistical robustness.

Future research should employ population genomics and large-scale sampling to deepen our understanding of the hybridization history and evolutionary dynamics of these snow trout species in the Tien Shan Mountains. Nuclear genome-level studies are critical to overcome the limitations and potential inaccuracies of phylogenetic inferences based solely on mitochondrial DNA, particularly in polyploid species with complex evolutionary histories. Moreover, population-level investigations using genomic data are essential for uncovering finer-scale hybridization patterns and offering more profound insights into their speciation processes and adaptive evolution.

Notably, the phylogeny based on the 13 protein-coding genes revealed that snow trout species from the Tien Shan Mountains formed a well-supported monophyletic clade (bootstrap value of 100 for ML tree and posterior probability of 1 for BI tree) (Fig. 2a). This clade originated around 7.12 Mya with 95% confidence interval (CI) of 8.46–5.70 Mya (Fig. 3). This clade subsequently diverged into two distinct clades (I and II) (Fig. 3) around 6.41 Mya (95% CI: 8.01–4.75 Mya). Clade I, which comprises Schizothorax sp. and S. fedtschenkoi, both endemic to the right tributaries of the Amu Darya River, diverged at 3.01 Mya (95% CI: 5.42–1.07 Mya) (Fig. 3). In clade II, S. pseudoaksaiensis from the Ili River, S. argentatus from the northern Tien Shan Mountains, and S. eurystomus from the Tarim Basin independently emerged at 3.91 Mya (95% CI: 5.97–1.78 Mya) and 1.31 Mya (95% CI: 2.29–0.45 Mya) (Fig. 3), respectively. Then, S. eurystomus diverged from S. biddulphi approximately 0.41 Mya (95% CI: 0.91–0.09 Mya) and then spread from the eastern Tien Shan Mountains to the western Tien Shan Mountains. The timing of these cladogenetic events coincides with the period of significant uplift of the Tien Shan Mountains that occurred 7.00–2.58 Mya (Sun et al. 2004). We infer that allopatric divergence of snow trout was likely driven by vicariance associated with the rapid orogeny of the Tien Shan Mountains from the Late Miocene to the Pliocene. Similar allopatric divergence patterns have been observed in other species, such as sand snake Psammophis turpanensis (Chen et al. 2021).

Figure 3. 

Time-calibrated phylogenetic tree of 33 Schizothorax species based on ML topology from 13 mitochondrial protein-coding genes. White circles with black borders represent calibrated nodes. The bolded species names indicate the sequenced species or individual in this study, while the regular species names indicate the species or individual for which NCBI data was obtained.

The uplift of the Tien Shan Mountain and its adjacent area had driven the diversification of these snow trout (Fig. 4). The extensive uplift of the Tien Mountain during the Late Miocene to the Pliocene separated the Amu Darya drainage and the Syr Darya drainage from the Tarim Basin, which divided the ancestor of these snow trout species into two groups, A1 and A2 (Figs 3, 4a, b). A1 is the ancestor of S. fedtschenkoi and Schizothorax sp., and A2 is the ancestor of S. pseudoaksaiensis, S. argentatus, S. eurystomus, and S. biddulphi. The orogeny of the Pamir-Alai mountains, western part of the Tien Shan Mountains, induced geographical isolation that led to the speciation within A1, resulting in two distinct species, S. fedtschenkoi and Schizothorax sp. (Fig. 4c). Concurrently, the uplift of the Tien Shan Mountains, led to the formation of the Ili River valley, sequentially triggered speciation processed within A2, leading to the emergence of S. pseudoaksaiensis (Fig. 4c), S. argentatus, and S. eurystomus (Fig. 4d). After that, S. eurystomus dispersed into the Syr Darya Basin along the Pamir-Tien Shan corridor, establishing the current biogeographic distribution pattern on both sides of the Tien Shan Mountains. Subsequent uplifts of the Tien Shan and the resulting changes in river systems facilitated the divergence of S. biddulphi from S. eurystomus, which is distributed in the Tarim Basin (Fig. 4e).

Figure 4. 

The biogeography history of snow trout in the Tien Shan Mountains. a. Current geographic distribution of these species. Occurrence data were compiled from the Global Biodiversity Information Facility (GBIF) and our field sampling sites; b. The barrier of the Pamirs-Tien Shan separated the ancestor of snow trout in the Tien Shan Mountains into two groups, A1 and A2, approximately 6.41 Mya. A1 is the ancestor of S. fedtschenkoi and Schizothorax sp., and A2 is the ancestor of S. pseudoaksaiensis, S. argentatus, S. eurystomus, and S. biddulphi. Thick white dashed lines indicate the barrier of the Pamirs-Tien Shan; c. The uplift of the northern Tien Shan Mountains around 3.91Mya facilitated the geographical isolation and subsequent speciation of S. pseudoaksaiensis from A2 along the Ili River. Additionally, the uplift of the Pamir-Alai mountains around 3.01 Mya prompted allopatric isolation, resulting in A1 diverging into two distinct species, S. fedtschenkoi and Schizothorax sp. Thin white dashed lines indicate different parts of the Tien Shan Mountains; d. Further uplift of the Tien Shan Mountains around 1.31 Mya and 0.41 Mya triggered allopatric speciation within A2, leading to the speciation of S. argentatus in the northern Tien Shan and S. eurystomus in the southern Tien Shan, respectively. Subsequently, S. eurystomus spread westward along the Pamir-Tien Shan corridor, resulting in the current geographic distribution pattern that crosses the Pamir-Tien Shan barrier. The white arrow indicates the dispersal event of S. eurystomus from east to west across the Pamir-Tien Shan barrier 0.41 Mya. e. The potential interspecific gene flow between S. biddulphi and the eastern population of S. eurystomus. The thin white dashed line with bidirectional arrows indicates the possible gene flow occurred between the two species.

The uplift of the Tien Shan Mountains has significantly influenced several species’ evolutionary trajectories and biogeographic patterns across this region. Similar to the observed allopatric divergences and origins in snow trout in this study, the uplift has influenced other taxa in comparable ways. For example, population differentiation in Meriones meridianus was driven by the uplift of the Tien Shan Mountains during the Middle Pleistocene (Wang et al. 2013), while amphipods, such as Gammarus lacustris, originated in the Tien Shan region and subsequently dispersed through alpine lakes to high-latitude and high-altitude areas across the Northern Hemisphere (Hou et al. 2022). Our study demonstrates that the uplift of the Tien Shan Mountains has played an essential role in shaping regional biodiversity through allopatric speciation, impacting both in the range’s western and eastern sections.