The Himalayan uplift and evolution of aquatic biodiversity across Asia: Snowtrout (Cyprininae: Schizothorax) as a test case

We derived a network representing the 132 Schizothorax haplotypes using PopART (POPu- lation Analysis with Reticulate Trees, v.1.7; https://popart.maths.otago.ac.nz/download/ [58]), which employs the TCS metric ([59], as implemented in [60]). To promote visualization, each haplotype was colored per its basin of origin (S3 Table), with colors corresponding to groups depicted in phylogenetic trees (Figs 3 and 4; black circles reflect hypothetical haplotypes not represented in our dataset. Each line indicates a single mutation connecting two haplotypes, but larger numbers of mutational steps are summarized within a black circle). Results Sequence data Our phylogenetic analyses incorporated 140 cytochrome-b (cyt-b) haplotypes, each spanning 1,140 bp. Of these, 643 sites were monomorphic (= 56%), 432 parsimony-informative (38%), and 65 (6%) represented singletons. Based on the Bayesian Information Criterion (BIC), the best fitting model of sequence evolution was TN+F+I+G4 (i.e., unequal evolutionary rates among purines/pyrimidines, empirical base frequencies, proportion of invariant sites, and a 4-rate category discrete Gamma model). Schizothorax across Asia and the Himalaya Our Bayesian approach without a priori assumptions detected evidence for 13 discrete, geo- graphically designated groups, (Fig 3). A ’ghost lineage’ emerged as an additional group but was subsequently disregarded as an artifact of isolation-by-distance (per [61]). Groups were designated according to geographic locations of constituent samples. The AMOVA indicated

sequence variation was primarily apportioned among the 13 groups, accommodating 78% of the total observed variance. Additionally, all were found to be significantly different one from another (Table 1). Lineages within regions also differed significantly, as did haplotypes within lineages (Fig 4). Probabilities for all comparisons were based on 1,000 permutations, with sig- nificance gauged by Bonferroni-correction (p<0.005 for 13 regions; p<0.0017 for 29 lineages). The visualization of phylogenetic relationship among geographic groups and constituents using maximum likelihood (ML) and Bayesian (BA) approaches resulted in largely consistent topologies (ML tree: Figs 3 and 4; Bayesian tree: S1 Fig). Statistical tests of topology subse- quently excluded the BA tree (p<0.05; Table 2). Given this, our interpretations primarily focused on the ML topology, but also highlighted differences between the two. Most nodes in the ML trees were strongly supported by SH-ALRT values and ultrafast boot- strap approximations [e.g., composite (Fig 3) and extended (Fig 4)]. Nodes were likewise sup- ported in the BA tree, with the majority displaying a posterior probability ranging from 0.9– 1.0 (S1 Fig). Two surprising results emerged in both analyses: Bhutanese snowtrout were consistently non-monophyletic, nesting instead within two spatially and temporally disparate clades (i.e., Central QTP/ YLTR versus Koshi basin of Nepal). This suggests multiple, independent coloni- zation of high-elevation Bhutanese drainages. In contrast, Nepali snowtrout clustered in a dis-

tinct east-to-west pattern, diverging in a sequential manner consistent with the broad-scale topographies of major tributaries, and with sub-clades being approximately congruent from a taxonomic standpoint. The specifics of group membership are provided in S3 Table. Groups recovered using different genetic clustering approaches corresponded generally with geographic regions (N = 13), as depicted in the ML tree (Fig 3). These are: (1) Central Asia; (2) Central QTP; (3) YLTR West-A; (4) YLTR West-B; (5) Bhutan-1A; (6) Bhutan-1B; (7) YLTR East; (8) Bhutan-1C; (9) Eastern QTP/Southeast Asia-Clade; (10) Bhutan-2 (via Brahmaputra River); (11) Koshi basin (eastern Nepal); (12) Gandaki drainage (central Nepal); PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 10 / 27