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

RESEARCH ARTICLE The Himalayan uplift and evolution of aquatic biodiversity across Asia: Snowtrout (Cyprininae: Schizothorax) as a test case Binod Regmi1¤, Marlis R. Douglas1, Karma Wangchuk1,2, Zachery D. Zbinden1, David R. Edds3, Singye Tshering2, Michael E. DouglasID1* 1 Department of Biological Sciences, University of Arkansas, Fayetteville, Arkansas, United States of America, 2 National Research & Development Centre for Riverine and Lake Fisheries, Ministry of Agriculture & Forests, Royal Government of Bhutan, Thimphu, Bhutan, 3 Department of Biological Sciences, Emporia State University, Emporia, Kansas, United States of America ¤ Current address: National Institute for Arthritis, Musculoskeletal & Skin Diseases (NIAMS), National Institute of Health, Bethesda, Maryland, United States of America * med1@uark.edu Abstract Global biodiversity hotspots are often remote, tectonically active areas undergoing climatic fluctuations, such as the Himalaya Mountains and neighboring Qinghai-Tibetan Plateau (QTP). They provide biogeographic templates upon which endemic biodiversity can be mapped to infer diversification scenarios. Yet, this process can be somewhat opaque for the Himalaya, given substantial data gaps separating eastern and western regions. To help clar- ify, we evaluated phylogeographic and phylogenetic hypotheses for a widespread fish (Snowtrout: Cyprininae; Schizothorax) by sequencing 1,140 base pair of mtDNA cyto- chrome-b (cytb) from Central Himalaya samples (Nepal: N = 53; Bhutan: N = 19), aug-

mented with 68 GenBank sequences (N = 60 Schizothorax/N = 8 outgroups). Genealogical relationships (N = 132) were analyzed via maximum likelihood (ML), Bayesian (BA), and haplotype network clustering, with clade divergence estimated via TIMETREE. Snowtrout seemingly originated in Central Asia, dispersed across the QTP, then into Bhutan via south- ward-flowing tributaries of the east-flowing Yarlung-Tsangpo River (YLTR). Headwaters of five large Asian rivers provided dispersal corridors from Central into eastern/southeastern Asia. South of the Himalaya, the YLTR transitions into the Brahmaputra River, facilitating successive westward colonization of Himalayan drainages first in Bhutan, then Nepal, fol- lowed by far-western drainages subsequently captured by the (now) westward-flowing Indus River. Two distinct Bhutanese phylogenetic groups were recovered: Bhutan-1 (with three subclades) seemingly represents southward dispersal from the QTP; Bhutan-2 appar- ently illustrates northward colonization from the Lower Brahmaputra. The close phyloge- netic/phylogeographic relationships between the Indus River (Pakistan) and western tributaries of the Upper Ganges (India/Nepal) potentially implicate an historic, now disjunct connection. Greater species-divergences occurred across rather than within-basins, sug- gesting vicariance as a driver. The Himalaya is a component of the Earth’s largest glacial reservoir (i.e., the “third-pole”) separate from the Arctic/Antarctic. Its unique aquatic PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 1 / 27

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Regmi B, Douglas MR, Wangchuk K, Zbinden ZD, Edds DR, Tshering S, et al. (2023) The Himalayan uplift and evolution of aquatic biodiversity across Asia: Snowtrout (Cyprininae: Schizothorax) as a test case. PLoS ONE 18(10): e0289736. https://doi.org/10.1371/journal. pone.0289736 Editor: Tzen-Yuh Chiang, National Cheng Kung University, TAIWAN Received: September 5, 2022 Accepted: July 26, 2023 Published: October 24, 2023 Copyright: © 2023 Regmi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: Data available via Genbank repository: Haplotype GenBank Accession Number 58thun003 OR373078 58DAKP06 OR373079 58BERT01 OR373080 58khar001 OR373081 58pots002 OR373082 58puza001 OR373083 58thun008 OR373084 58danr004 OR373085 58NYAC01 OR373086 58haac001 OR373087 58DAKP011 OR373088 58DIKC03 OR373089 58KAME01 OR373090 58shen006 OR373091 58TOEB12 OR373092 58danr001

biodiversity must be defined and conserved through broad, trans-national collaborations. Our study provides an initial baseline for this process. Introduction The Himalaya is the most extensive and recently evolved mountain system on Earth (length = 2,400km; width = 240km; elevation = 75–8,800m), with global significance under- scored by large-scale lithospheric, cryospheric, and atmospheric interactions [1]. These have not only driven global climate, but also defined the cultural and biological endemism of the region [2]. The Himalaya is topographically separated into three components: Western Hima- laya: ~880 km wide region framed by Indus and Kali Gandaki rivers (southwestern Tibet Autonomous Region of China to western Nepal, respectively); Central Himalaya: ~800 km wide area contained between the Kali Gandaki and Tista rivers (from western Nepal to Sikkim, India, respectively) and; Eastern Himalaya: ~720 km stretch delimited by the Tista and Brah- maputra rivers (from Sikkim to Arunachal Pradesh in India) (https://www.pmfias.com/ himalayas-regional-divisions-punjab-himalayas-assam-himalayas-western-himalayas-central- himalayas-eastern-himalayas/). Massive, tectonically derived mountain chains such as the Himalaya are fundamental in establishing global biodiversity gradients, but strong signals of vicariance and local adaptation are generally noted within terrestrial rather than aquatic sys- tems [3]. Here we evaluate how orogeny (the deformation and folding of Earth’s crust by lat-

eral compression) has contributed to the potential diversification of freshwater fishes broadly across Asia by quantifying phylogeographic and phylogenetic patterns withing a an endemic high-elevation fish genus (Schizothorax). We use these data to infer a potential evolutionary hypothesis for the genus, and offer the resulting perspective as a surrogate for other aquatic biodiversity components in the region. Orogeny and biodiversity of the study region The Himalaya formed as a series of parallel ranges bordered by rivers to the west (Indus) and east (Brahmaputra). It encompasses most of northern Pakistan, northern India (Kashmir to the west, Assam to the east), as well as Nepal, Bhutan, and parts of China, with the extensive Qinghai-Tibetan Plateau (QTP) to the north and the alluvial plains of India and Bangladesh to the south (Fig 1A). It evolved over a 30 Ma span, as one component of the extensive impact between the Indian craton and the Eurasian landmass [4]. This uplift established the Himalaya as a topographic entity and profoundly influenced the connectivity of its regional rivers [5, 6]. In late Oligocene–early Miocene (26–19 Ma), the Yar- lung-Tsangpo River (= YLTR; average elevation 4,000m) formed on the southern border of the QTP as an east-west series of lakes with south-flowing tributaries (Fig 11 in [7]). These lakes eventually coalesced into a westward flowing river (early–mid Miocene; 19–15 Ma) that even- tually reversed its flow in response to a continuing plateau uplift. The YLTR was captured on

the eastern QTP by the Brahmaputra River in early-mid Miocene (i.e., ~18–15 Ma) [6, 8]. The uplift of the QTP was a signature event that seemingly allowed aquatic biodiversity to disperse broadly from west to east, with subsequent access to south-flowing Himalayan drainages even- tually truncated by continued uplift [9]. Three stratigraphic zones subsequently emerged, as manifested by the contemporary spe- cies richness of terrestrial biodiversity [4]. Although this endemism has been hypothesized as driven by tectonic activity, requisite data are either sporadic or absent [10]. Freshwater fish divergence seemingly occurred concomitantly, with tributaries increasingly isolated as the PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 2 / 27 OR373093 58puza006 OR373094 58tink001 OR373095 58karo003 OR373096. Funding: Author(s) received no specific funding for this work. Competing interests: Authors have declared that no competing interests exist.

QTP elevated more than 3,000m in late-Miocene [11]. The uplift also drove strong selection for high-elevation specialists that could subsist under more severe conditions [12]. However, data remain sparse for high elevation Himalayan fishes, an aspect inherently magnified by the lack of large-scale monitoring programs [13]. Such data-deficiencies, tightly linked with the geomorphic and climatic history of the region, underscore the strong potential for cryptic spe- cies yet to be identified and formally described. Schizothorax and its outgroups Snowtrout (Fig 1B) offers a unique opportunity to deciphering the geologic processes that drove fish diversification and distribution in the Himalaya. The genus belongs to the largest and most diverse freshwater fish family (Cyprinidae: 1,700+ valid species, with 11 recognized subfamilies) [14], and comprises 28 recognized species across lakes and rivers of the QTP, the endorheic Tarim River basin, as well as the snow-fed tributaries of the Himalaya/ sub-Hima- laya (Fig 1A) [15–17]. Several hypotheses have been offered to explain its origin and dispersal. Researchers first suggested the Schizothoracinae initially dispersed into the QTP from southeastern Asia [18], and subsequently became isolated by orogeny and climate change. Others proposed that a Fig 1. (A) Simplified drainage map of Asia, showing major basins from which Schizothorax were obtained. Area in yellow indicates the Tibetan Plateau, referred to as Qinghai-Tibetan Plateau (QTP) in China. The Tarim River, an

endorheic basin north of the QTP in Central Asia, is not depicted. Three major rivers, the Irrawaddy, Salween, and Mekong, drain the QTP into Southeastern Asia (SEA), whereas the Yangtze River flows east into the China Sea (Republished from Fig 1 [100] under a CC BY license, with permission of the Council on Foreign Relations, original copyright 2016). (B) Schizothorax sp. cf richardsonii (i.e., potentially S. richardsonii) captured 02 December 2016 at confluence of Bibigang Chhu and Mangde Chuu, Bhutan (Total length = 540 mm). https://doi.org/10.1371/journal.pone.0289736.g001 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 3 / 27

primitive barbin clade was isolated by tectonism on the QTP [19], thus serving as the progeni- tor for contemporary Schizothorax. The subfamily Cyprininae was recently revised [20], based on a combination of molecular markers [i.e., cyt-b (N = 791); mitogenomes (N = 85); nuclear RAG-1 sequences (N = 97)], with Schizothorax being retained within the monophyletic tribe Schizothoracini, which evolved either pre- or post-Barbini, depending upon the molecular marker employed. Six major biogeographic groups of Snowtrout have been proposed [21]: A Central Asiatic clade, followed by five drainage-specific sister-clades: YLTR/ Brahmaputra; Indus; Irrawaddy; Mekong/Salween; and Yangtze. However, this perspective remains topographically incomplete as data for major Himalayan drainages were unavailable. We extended these results by focus- ing on two previously unrepresented regions: Bhutan (eastern Himalaya) and Nepal (central Himalaya). In tandem, they link the far eastern Himalaya with its western terminus, thus allowing for a more complete evaluation of tectonism, orogeny, and the potential diversifica- tion of aquatic biodiversity in the region. By filling these biogeographic gaps, Snowtrout can be more appropriately recognize as a biodiversity component extending from Central Asia, across the QTP into eastern, southern, and western Asia. This, in turn, provides a formal tem- plate of regional biodiversity from which potential impacts of Quaternary climate oscillations can be more appropriately interpreted [22]. Contemporary significance

The QTP (and surrounding mountains) are termed the “third pole” because the region con- tains the greatest amount of global snow and ice beyond the Arctic/Antarctic. It also represents the headwater source for 10 large Asian rivers (Fig 1A), which collectively provide ecosystem services for >1.5 billion people [23]. Glaciers in the eastern Himalaya, for example, represent an estimated 14.5% of the global total, a quarter of which has been lost since 1970 [24], due largely to a regional warming 2x greater than the global average [2]. These cumulative impacts are recorded as increasing precipitation, decreasing glacial surface area, diminishing snow cover, expanding glacial lakes, and a loss of permafrost [25]. These biophysical occurrences have substantially impacted the distribution of Schizothorax, provoking a migration into higher-elevation streams with a concomitant contraction of trail- ing-edge habitats [26]. This redistribution has also facilitated contact with invasive aquatic predators [27], and represents an essential element of ecological resilience for Schizothorax (e.g., its adaptive capacity). This shift further adds to the taxonomic uncertainties within the genus, as well as emphasizing the strong potential for unrecognized, cryptic diversity. A comprehensive understanding of endemic biodiversity is a necessary requirement for holistic, scalable approaches that mitigate the impacts of regional climate-change. These must be derived from standardized metrics that quantify interspecific variation across political

boundaries. Such an approach can provide templates for cooperative regional conservation [25] and help define the ecological resilience of native biota, against which climate impacts are defined and incorporated. Herein, our approach is to evaluate a biologically-complex genus of freshwater fish as a focal point for the development of such components. Materials and methods Ethics statement All methods were performed in accordance with relevant guidelines and regulations. Nepali samples were collected as one component of a Fulbright Award to DRE (see Acknowledg- ments), whereas Bhutanese collecting permits were in conjunction with the National Research and Development Centre for Riverine & Lake Fisheries (NRDCR&LF), Ministry of Agriculture PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 4 / 27

& Forests (MoAF), Royal Government of Bhutan. The export of fin clips was authorized through a Material Transfer Agreement (MTA) provided by the Bhutan Agricultural and Food Regulatory Authority (BAFRA). Sampling protocols were approved by the University of Arkansas Institutional Animal Care and Use Committee (UA_IACUC_17064). Specifically, sampling was accomplished using standard fisheries methods, to include hand- held seines, dip nets, backpack shockers, as well as passive gear such as hoop/ trammel nets and minnow traps. These methods and devices are approved for use by the American Society of Ichthyologists and Herpetologists (ASIH; https://www.asih.org/resources) and the Ameri- can Fisheries Society (AFS; https://fisheries.org/books-journals/online-resources/). In addi- tion, educational resources [Association for the Advancement of Laboratory Animal Science (AALAS): https://www.aalas.org/] were employed to assist in the instruction of Bhutanese cadre. Small-bodied fishes (<10 cm total length) were removed from capture devices with <20 placed concurrently within a 19-liter bucket of well-aerated stream water. Direct handling of fish was minimized, and individuals were retrieved via dip net for identification and subse- quent release. Larger fishes were gently removed from capture devices and identified/ released at site of capture. Care was taken in all cases to ensure that water temperatures were identical between capture site and holding bucket. Once captured, fishes were continuously

monitored, and euthanasia immediately employed if sampling injuries were detected. Fishes were so recognized if they could not maintain equilibrium in the holding bucket, and/or dis- played erratic or impaired opercular movements (<20 respiratory cycles/minute). Euthana- sia was conducted via an overdose of clove oil (Eugenol), followed by fixation in 10% formalin. Clove oil satisfies specified criteria and is considered acceptable when used at high concen- trations (>400 mg/L), as it causes a reliable and rapid loss of consciousness and induces hyp- oxia (critical components for eliminating potential pain; American Veterinary Medical Association Guidelines on Euthanasia; https://www.avma.org/sites/default/files/2020-02/ Guidelines-on-Euthanasia-2020.pdf). A concentration of 50–150 mg/L was initially employed, then gradually increased to >400 mg/L through periodical additions (~50 mg/L). Complete cessation of opercular movement was used as an indicator of death, with the fish then being retained in solution for >10 additional min. Tissues for DNA analysis were non-lethally sampled by removing a small piece of anal fin (a non-propulsive appendage). The fish was first placed in a small bucket of stream water treated with clove oil, and subsequently removed via dipnet immediately following a loss of righting ability. A small piece of fin was removed with surgical scissors cleaned in 100% EtOH. The clip was transferred to a Whatman FTA© card for preservation and storage, and the fish revived in a small bucket of untreated stream water until righting ability was regained, at

which time it was released at the collecting site. Wounds caused by clipping of fins will nor- mally heal satisfactorily without application of antibiotics or topical analgesics, which are not recommended for use in wild fishes (per AFS resources). Sample sites: Drainages in Nepal and Bhutan Nepal is drained by three major east-to-west rivers (Koshi, Gandaki, and Karnali; Fig 2A) that collectively support eight tributaries, several of which originate on the QTP. In eastern Nepal the stronger effects of the monsoon are reflected by more eroded and abruptly-elevated moun- tainous slopes of the Koshi River. Central Nepal, dominated by the Gandaki Basin and its west- ern tributary (Kali Gandaki), serves as a distinct east-west vicariant break. Western Nepal is drained by the Karnali Basin and includes Rara Lake (2,990m elevation), a RAMSAR site with PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 5 / 27

Fig 2. Map depicting sampling locations for Schizothorax, with metadata for ’Site,’ ’River,’ ’Location,’ ’Lat-Long,’ ’Haplotype,’ and ’Group’ compiled in S1–S3 Tables, respectively. Red circles in map indicate sampling sites with samples identified by Group number in (A) Nepali as N = 1–9, and (B) Bhutan as B = 1–9. Maps were generated in ArcGIS 10.4.1 using publicly available data from HydroRIVERS and HydroBASIN (https://www.hydrosheds.org/products/hydrorivers). https://doi.org/10.1371/journal.pone.0289736.g002 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 6 / 27

a surface area of 10.6 km2. All three basins drain south into the east-flowing Ganges River which joins with the Brahmaputra River to terminate in the Bay of Bengal (Figs 1 and 2A). Bhutan has five major north-south river systems (Fig 2B): The Amo Chhu and Wang Chhu descend from the QTP, flow southeasterly through western Bhutan and then into West Bengal (India). The Punatsang Chhu (= Sunkosh) drains the Great Himalaya Range in west-central Bhutan, and flows into Assam (India). The Manas River and its major tributaries drain cen- tral-eastern Bhutan. The Dangme Chhu originates in Arunachal Pradesh (India) and flows southwesterly through Bhutan, whereas the Kuri Chhu (from the QTP), the Chamkhar Chhu, and the Mangde Chhu drain central Bhutan. The eastern-most river, the Nyera Ama Chhu, is the only drainage that originates entirely within Bhutan. All river systems eventually drain into the Brahmaputra River (northern India), which curves south and west through Bangladesh to merge first with the Ganges River before flowing into the Bay of Bengal (Figs 1, 2A and 2B). Amplification and sequencing of mtDNA samples Nepali samples (N = 53; S1 Table, Fig 2A) represent five Schizothorax species. Three [S. macro- phthalmus (N = 3), S. nepalensis (N = 4), and S. raraensis (N = 3); total N = 10] represent a landlocked species-flock within Rara Lake (Karnali River drainage; Fig 2A). Two other riverine species [S. richardsonii (N = 26), and S. progastus (N = 17)] were sampled from (west-to-east): Karnali River (N = 12/N = 3, respectively); Gandaki River (N = 11/9); Koshi River (N = 5/3).

Nepali samples were supplemented with those from Bhutan (N = 19; S2 Table, Fig 2B), col- lected from four west-to-east drainage basins: Wang Chhu (N = 2); Punatsang Chhu (N = 11); Mangde Chhu (N = 3); Dangme Chhu (N = 3). Genomic DNA was extracted using the Qiagen DNeasy kit, per manufacturer’s instruc- tions. We utilized published primers to amplify mtDNA cyt-b [28], with products enzymati- cally purified, sequenced using BIGDYE (ver.3.1) chemistry (Applied Biosystems Inc. [ABI], Forest City, CA), and analyzed on an ABI Prism 3700 Genetic Analyzer (W.M. Keck Center for Comparative and Functional Genomics, University of Illinois, Urbana/Champaign). In addition, we accessed 68 mtDNA cyt-b sequences representing 39 species (N = 60 Schizothorax plus eight outgroups [Acrossocheilus (N = 2); Pethia (N = 2); Puntius (N = 1); Tor (N = 2); Neo- lissochilus (N = 1); S3 Table] (National Institute of Health; https://www.ncbi.nlm.nih.gov/ genbank/). To avoid potential ambiguities with regards to evolutionary trajectories, we chose outgroups from select tribes identified as having originated prior to evolution of the Schi- zothoracini. These include (in sequence): Torini (i.e., Neolissochilus); Smilogastrini (i.e., Pun- tius, Pethia); and Acrossocheilini (i.e., Acrossocheilus, the most immediate outgroup) (S3 Table). Additional details on generating molecular data are provided [29]. Sequences were edited/aligned using SEQUENCHER v.5.3 (Gene Codes Inc., Ann Arbor, MI), with nucleotide composition and uncorrected p-distances derived among groups (MEGA X; [30]).

Polyploidy and choice of marker system The subfamily Cyprininae displays extensive polyploidy (>400 species; [20, 31]), and given the number of polyploid species, different ploidy levels, and the variety of independent polyploidi- zation events identified, it may represent the most complex vertebrate polyploid system cur- rently known [32]. Ploidy is an important consideration when molecular phylogenies are reconstructed in that it significantly impacts statistical approaches, results, and subsequent interpretations. For example, polyploidy complicates the detection of population structure, assignment tests, and principal component analyses [33]. These limitations are particularly problematic given the greater redundancy of nuclear genes in polyploid genomes. All may potentially diverge one from another over evolutionary time, such that results can differ PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 7 / 27

markedly, depending upon which ‘sub-genome’ is being accessed when polyploid data are phylogenetically evaluated [34]. Given this, phylogenetic analyses of polyploid species are often limited in practice to either mtDNA (given its haploid nature), or nuclear genes cloned into bacteria and subsequently iso- lated as single-copy (such as RAG-1) [20]. However, the latter are evolutionarily more conser- vative and, for example, cannot be utilized in a TIMETREE analysis [35] as herein. Yet polyploidization in Schizothorax also represents an evolutionary mechanism that may have facilitated its invasion of high-elevation habitat [31]. We thus restricted our study to mtDNA haplotypes and relied upon GENBANK to increase our sample number and geographic representation. Available numbers are often limited, due either to restricted sampling and/or logistic difficulties in accessing those regions where study species occurs. In this context, Schizothorax inhabits regions notoriously difficult to access, both topographically as well as socio-politically. In addition, GENBANK samples often represent output from ongoing investigations and their metadata thus lack precise GPS coordinates required for specific analyses. This, as well as the large variance in number of base pairs avail- able per sample, has constrained our total GENBANK sample number. Many of the above issue can also confound analyses and, to avoid such issues, we rigorously assessed GENBANK meta- data prior to incorporating sequences. Analyses to infer evolutionary history

We employed a four-step approach to infer potential evolutionary scenarios that could explain the biodiversity in Schizothorax across the Himalaya. This involved: (a) Evaluating congruence between taxonomic and genetic designation of samples; (b) Detecting cryptic biodiversity; (c) Characterizing regional and basin-specific genetic lineages; and (d) Quantifying phylogeo- graphic patterns. We employed a Bayesian clustering analysis to identify regional gene pools without a priori assumptions pertaining to location and/or taxonomic classification. We then used these groupings to statistically evaluate genetic variation within and among regions and basins (AMOVA). We generated phylogenetic hypotheses by deriving ML and Bayesian tree topologies, which were statistically and qualitatively compared for congruence (IQ-TREE). We also derived a haplotype network as a means of visualizing the number of stepwise mutations separating distinct sequences and groups. Fossil-calibrated methods were then employed to estimate divergence times for major groupings in the context of regional geologic events. After compiling a dataset of 132 Schizothorax samples using the criteria described above, we derived a dataset that included only variable sites (N = 365bp: DNASP v.6; [36]). We employed a Bayesian approach (GENELAND V.4.0.3: [37, 38]) to infer the number of distinct groups (K) and to designate group membership for each sample. Ten runs of 100,000 MCMC iterations and a thinning interval of 100 were used, with K varying between 1 and 15, but with-

out use of optional geo-referencing, given the absence of such coordinates for most GENBANK sequences. The estimated number of populations (K) was defined using the maximum poste- rior estimate. We conducted an analysis of molecular variance (AMOVA: ARLEQUIN v.3.5.1.2; [39]) to gauge how genetic variation was partitioned statistically among geographic regions (as informed by GENELAND). Here, we defined populations (= distinct gene pools) as being haplo- types within regions [40]. Probabilities were ascertained using 1,000 permutations, and statisti- cal significance was assigned using a Bonferroni-correction. We derived phylogenetic hypotheses using a maximum likelihood (ML) analysis (IQ-TREE v2; [41]), with the best-fitting substitution model being generated via MODELFINDER, selecting flexible rate heterogeneity (e.g., gamma-distributed) across sites [42]. Model space was hierar- chically explored by sequentially adding parameters (e.g., rate categories) until the model failed PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 8 / 27

to improve. Support for the resulting tree, inferred using the selected model, was explored using 1,000,000 ultrafast bootstrap approximations (UFBOOT2; [43]). We also assessed branch support using the SH-ALRT test (Shimodaira-Hasegawa approximate Likelihood Ratio Test, based on a null distribution derived from bootstrap replicates of varying sizes). Values �80% indicated strong support for clades, whereas those for UFBOOT were �95% (IQ-TREE v.2 Man- ual: http://www.iqtree.org). To employ a Bayesian approach, we first derived an input XML file using a relaxed molecu- lar clock with an uncorrelated lognormal rate distribution and a Yule process prior (graphical user interface (i.e., BEAUTY) within the phylogenetic software suite BEAST (v2.61; [44]). Each analysis was run for 100 million generations, with samples taken every 1,000th iteration, fol- lowing a burn-in of 20%. Effective sample size (ESS) was verified as >200 across all parameters with MCMC chain convergence visually inspected using TRACER [45]. Trees from three inde- pendent runs were then combined (via LOGCOMBINER) as a maximum clade credibility (MCC; TREEANNOTATOR [44]), with nodal support representing posterior probabilities. The resulting topology was visualized using FIGTREE (v.1.4.3; [46]). We statistically compared topologies from ML and BA analyses with respect to the original sequence alignment (IQ-TREE2 v2; [43]). We then examined the support displayed by each regarding RELL approximation (10,000 replicates), as well as: (a) Raw change in log-likeli-

hoods; (b) Bootstrap proportions [47]; (c) Kishino-Hasegawa test [48]; (d) Shimodaira-Hase- gawa test [49]; (e) Approximately Unbiased test [50]; and (f) Expected Likelihood weights [51]. Developing a TIMETREE For molecular dating, we employed a rapid, high-performance ML-based model (RELTIME, MEGA X; [33]) that incorporates a relative rate framework (RRF) within which the evolution- ary rates of sister lineages are contrasted against a minimum rate change, as diagnosed between evolutionary lineages and their respective descendants. The computational speed and accuracy of the model are appropriate for the analysis of large data sets (>100 ingroup sequences; [52], with procedural steps for RELTIME previously documented [53]). The evolu- tionary model for our ML-based phylogenetic tree was also implemented in RELTIME. We estimated divergence times based on the earliest fossil evidence for Schizothoracini (e.g., Paleoschizothorax qaidamensis; [54], as gauged in the Qaidam Basin [55]). This placed initial divergence at the onset of the Oligocene (33 Ma) (see also [56]). Phylogeographic analyses A haplotype network provides a visualization of genealogical relationships based on individual mutations separating haplotypes, and as such, provides insights into genetic diversity within populations or closely related species. It is beneficial in that it effectively separates those indi- viduals with similar haplotypes th generally depicted within a phylogenetic tree as unresolved polytomies (i.e., ’combs’). However, it differs from a strictly bifurcating phylogenetic visualiza-

tions in that each haplotype (depicted as a circle) can potentially be connected to multiple, rather than just a single, node. A second difference is that a comparative phylogenetic tree represents a platform from which genetic variation can be interpreted within the context of geologic history, as supported by statistical and/or Bayesian analyses [57]. In addition, a phylogenetic tree is rooted from ancestral to more derived, with geographically co-distributed species (and their shared lineage histories) transcending. In short, relationships can be well-depicted within a species using a haplotype network, but often much less so for distinct species, given the increase in genetic variability and the lack of an evolutionary trajectory. PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 9 / 27

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

Fig 3. A collapsed maximum likelihood (ML) tree reflecting phylogenetic relationships among biogeographic groups. The constituent samples comprising each biogeographic clade are provided in S3 Table under ’Species Name,’ ’Haplotype,’ ’Accession #,’ ’Location,’ and ’Group’ headings. Data were derived from sequence analysis of the cytochrome-b mitochondrial gene (1,140 bp, 140 haplotypes, clades collapsed). Red asterisk = SH-LRT (Shimodaira-Hasegawa approximate likelihood ratio test) �80% and UFBOOT (ultrafast bootstrap approximations) �95% (IQ-TREE v.2 Manual: http://www.iqtree.org). Numbers at colored triangles indicate the number of haplotypes PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 11 / 27

(13) a composite western Himalaya clade that included Rara Lake and the Karnali drainage (western Nepal), the upper Ganges (northwest India), and the Indus River Basin (Pakistan). A more detailed perspective is presented in Fig 4, with regional composites listed in S3 Table. Most noteworthy is the partition of the Eastern QTP/ Southeast Asia (E-QTP/ SEA) clade into three distinct geographic clusters: E-QTP/ SEA-A (Salween-Irrawaddy drainages); E-QTP/ SEA-B (upper Yangtze drainage); (9) E-QTP/SEA-C (Mekong-Salween drainages). In the BA tree (S1 Fig) the western Himalaya was additionally partitioned into two subclades: Indus/Upper Ganges versus Karnali/Rara Lake. In the ML tree, samples from Central Asia formed a distinct regional clade that clusters as sister to all other samples, including. S. argentatus, which is sister to all remaining species. These form two distinct composite clades. In the first, samples representing the Central QTP form a group, followed by the YLTR West-A and YLTR West-B, then Bhutan-1A separated from Bhutan-1B, and finally a group representing by YLTR East and Bhutan-1C. The second large ML cluster is subdivided into two distinct sister clades, each partitioned into subclades consistent with geographic regions. The first comprises the Eastern QTP/ Southeast Asia and is represented by three significantly different groups (labeled A, B, C). The second contains the Koshi River (eastern Nepal) as sister to the Bhutan-2 group (representing eight undescribed Bhutanese samples), followed by the Gandaki River as sister to a mixed

western Himalaya clade containing the Karnali River/ Rara Lake, the upper Ganges River, and the Indus River (Figs 3 and 4). Clustering of S. richardsonii and S. progastus in the Koshi, Kar- nali, Ganges, and Indus River drainages is noteworthy because individuals do not cluster sepa- rately by taxonomic designation within each drainage, but instead group as co-mingled components. Major differences between ML (Figs 3 and 4) and BA (S1 Fig) trees include: Position of the Central Asia and E-QTP/ SEA clades, relationship of YLTR and Bhutan-1 subgroups, and position of S. argentatus in relation to the Central Asia group. The regions are the same within the BA tree, but distributed across two large regional clusters, with Central Asia included within the first, rather than being distinct from it. YLTR West is subdivided into groups by the presence of Bhutan-1A, whereas Bhutan-1B is instead subsumed within it rather than being sister to YLTR West/Bhutan-1C. The E-QTP/ SEA-A clade is also further subdivided within the BA tree (a pattern also distinct in the haplotype network, Fig 5). And lastly, the Ganges and Indus rivers appear as distinct western Himalaya clades, whereas they are unresolved in the ML tree. Temporal divergences We employed RELTIME to derive major time-calibrated events, 13 of which are depicted in our ML-based tree (Fig 3; other dates annotated within discussion; See also S1 Appendix for the associations between divergence intervals and geomorphic/ climatic drivers in the Himalaya).

The average divergence time for Central Asia was 11.6 Ma (CI = 15.4–8.7). Two large sub- clades sister to Central Asia then separated at ~7.3 Ma (CI = 11.3–4.7). In the first, the Central QTP diverged at 6.2 Ma (10.2–3.7), followed by the YLTR West (4.2 Ma; CI = 7.4–2.4), then Bhutan-1A (3.3 Ma; CI = 6.7–1.5), with YLTR East and Bhutan-1B at 2.7 Ma (CI = 5.9–1.2). In the second large subclade, the Eastern QTP/ Southeast Asia region separated at 6.6 Ma (CI = 10.7–4.1), followed by a second large composite clade at 4.8 Ma (CI = 8.8–2.6). Here, the grouped within each. Numbers at nodes represent RELTIME TIMETREE mean divergence estimates (in Ma; confidence intervals in text; See also S1 Appendix). https://doi.org/10.1371/journal.pone.0289736.g003 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 12 / 27

Fig 4. Maximum likelihood (ML) tree reflecting phylogenetic relationships among biogeographic samples. The constituent samples representing each biogeographic clade are provided in S3 Table under ’Species Name,’ ’Haplotype,’ ’Accession #,’ ’Location,’ and ’Group’ headings.’ Data were derived from sequence analysis of the cytochrome-b mitochondrial gene (1,140 bp, 140 distinct haplotypes). Red asterisk = SH-ALRT (Shimodaira-Hasegawa approximate likelihood ratio test) � 80% and UFBOOT (ultrafast bootstrap approximations) � 95% (IQ-TREE v.2 Manual: http://www. PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 13 / 27

from other Western Himalaya haplotypes (>50 mutations). Our network exhibits two distinct patterns: (a) Clusters of similar haplotypes, such as the Karnali and Gandaki river groups, reflecting scant internal divergence normally characteristic of within/among population genetic diversity; and (b) Groups comprised of considerably divergent haplotypes, such as the E-QTP/ SEA assemblage (sensu lato, lower left). The latter is to be expected given the geographic spread of its constituent samples ([21], S2 Table). Discussion Synergy between tectonism and climate Orogeny and tectonism drive the diversification of aquatic ecosystems, and act synergistically with climate to impact both hydrology [62] and resident freshwater fishes [63]. For example, iqtree.org). Numbers only at specific nodes represent UFBOOT values. Letters within composite species-triangles at tips of the tree represent species, with R = S. richardsonii, P = S. progastus, and x = Both species, with numbers representing individuals contained. https://doi.org/10.1371/journal.pone.0289736.g004 Table 1. Results of an analysis of molecular variance (AMOVA) based on 132 haplotypes of snowtrout (Cyprini- dae: Schizothorax). Source represents: (1) Geographic regions (= Among regions); (2) Lineages within regions (= Within regions/among lineages); (3) Haplotypes within lineages (= Within lineages). Source = Grouping; D.F. = Degrees of freedom; S.S. = Sums of squares; Var. = Variance; % = Percent variance. Probability values for each Fst-

value (derived via 1,000 permutations) are significant at a Bonferroni-adjusted value of p<0.005. Source D.F. S.S. Var. % Among regions 10 3958.7 31.4 77.9 Within regions/among lineages 18 466.1 5.6 13.5 Within lineages 104 350.6 3.4 8.5 Total 132 4775.3 40.3 99.9 https://doi.org/10.1371/journal.pone.0289736.t001 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 14 / 27

historic drought over evolutionary time in western North America gradually shaped tectoni- cally-derived drainage systems, resulting in the erosion of genetic diversity, with populations of endemic fish repeatedly collapsing into refugia then expanding outward during more pluvial periods [64]. Episodic flooding, on the other hand, allowed antecedent streams to down-cut, with headwaters subsequently eroding, basins concomitantly expanding, with large-scale dispers- als being promoted [65]. Such alterations (i.e., captures, diversions, beheadings; [66]) not only extended fishes into adjacent basins [67] but also facilitated hybridization [68–70], a process that has (and continues to) confound both taxonomic resolution and management strategies. Flows within basins are often abridged by longitudinal barriers that are vicariant for larger-bodied species, yet serve as passive filters for fishes smaller and more limited in dis- persal capacity [71]. As a result, fish distributions represent a series of biogeographic ’islands’ within and among basins, each reflecting local expansions, contractions, and extinctions. This process provides a necessary template from which to gauge the response of fishes to similar tectonic, orogenic, and climatic impacts within dendritic riverscapes less well moni- tored (as herein). The Himalaya provides such a template as shaped by tectonic and climatic events. Its sharp elevational gradients derived via orogeny provided ample opportunities for diversification. For

example, many deeply incised valleys emerged along the border of the elevating QTP, and given their increasing gradients, subsequently transitioned into prominent geomorphic fea- tures. As a result, many unique physiographic characteristics emerged, particularly when com- pared with other large global river networks. For example, five (of 10) QTP rivers [i.e., Yellow, upper Yangtze (= Jinsha), upper Mekong (= Lancang), upper Salween (= Nujiang), and Yar- lung-Tsangpo (= upper Brahmaputra)] originate at >4,500m elevation, flow >1,000km to the border of the QTP, and subsequently drop >3,000m to Asian lowlands. Interestingly, the Mekong, Salween, and Yangtze flow in parallel for ~170km, and achieve their closest proximity at <70km (Fig 1B; [72]). In S1 Appendix, we provide references that link divergence intervals for our phylogeographic clades with geomorphic/ climatic drivers in the Himalaya. These events, in turn, have established the region as a global hotspot for biodiversity ende- mism [21, 73] with Snowtrout as a characteristic element of the aquatic fauna. A more specific example is found within the YLTR, where a deep canyon derived via orogeny effectively parti- tioned S. oconnori [74], promoting several distinct ESUs (Evolutionarily Significant Units; [75]). This vicariance similarly impacted S. waltoni (BA tree; [76]), but those components are evolutionarily less distinct and designated instead as MUs (Management Units; [75]). Studies of Himalayan freshwater most often focused on select regions, rivers, and/or moun-

tains [77–80]. Our approach broadens that perspective, in that we fill biogeographic gaps in the Himalaya and leverage existing data to more fully evaluate origin and subsequent dispersal of Snowtrout across Asia. Table 2. Statistical topology tests comparing maximum likelihood (ML) and Bayesian (BA) phylogenetic results (= Tree). logL = Log Likelihood value; deltaL = logL difference from maximal logL in the set; Bp-RELL = Bootstrap proportions using RELL method (weights sum to 1 across trees); p-KH = p-value of one-sided Kishino- Hasegawa test; p-SH = p-value of Shimodaira-Hasegawa test; c-ELW = Expected Likelihood Weight (weights sum to 1 across trees); p-AU = p-value of Approximately Unbiased test. A plus sign (+) next to a p-value denotes 95% confidence sets, whereas a minus sign (-) denotes a significant exclusion (i.e., the tree is rejected). All tests stemmed from 1,000 resamples using the RELL method. Tree logL deltaL bp-RELL p-KH p-SH c-ELW p-AU ML -9526.83 0.0 0.986 + 0.986 + 1.00 + 0.978 + 0.985 + BA -9582.76 55.93 0.014 - 0.013 - 0.013 - 0.014 - 0.015 - https://doi.org/10.1371/journal.pone.0289736.t002 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 15 / 27

Fig 5. TCS haplotype network reflecting stepwise mutational differences between samples and among geographic groups. The contents of each biogeographic clade are provided in S3 Table under ’Species Name,’ ’Haplotype,’ ’Accession #,’ ’Location,’ and ’Group’ headings. Data were derived from sequence analysis of the cytochrome-b mitochondrial gene (1,140 bp, 132 haplotypes). Geographic groupings (Fig 3, S3 Table) are labeled within rectangles, with group colors identical to those in Figs 3 and 4. Black circles represent hypothetical haplotypes to reflect stepwise mutations >14. See S1 Appendix for the associations between divergence intervals for clades and geomorphic/ climatic drivers in the Himalaya. https://doi.org/10.1371/journal.pone.0289736.g005 PLOS ONE Himalayan uplift and aquatic evolution PLOS ONE | https://doi.org/10.1371/journal.pone.0289736 October 24, 2023 16 / 27

Central Asia Five species from Central Asia form a distinct group, but their relationships within Schi- zothorax differs per ML and BA analyses, as does S. argentatus. In the ML tree (Figs 3 and 4), the Central Asia group clusters as a clade apart from all other Schizothorax, with S. argentatus positioned as sister to all remaining species. In contrast, the BA tree (S1 Fig) places the Eastern QTP/ Southeast Asia as a separate clade outside other species, with the Central Asia group clustering with S. argentatus as sister to the remaining species. These six species represent the same geographic region (Fig 3) and exhibit potential regional and/or taxonomic substructure, with S. pseudoaksaiensis and S. argentatus distributed in Lake Balkhash and its catchment basin (Kazakhstan, central Asia), as well as the Ili River (northwest China; [81]). Similarly, S. biddulphi [80] and S. eurystomus [82] occur in the Tarim River and tributaries (an endorheic basin in the Central Asiatic Desert of northwest China), whereas S. intermedius (represented by two distinct haplotypes) has a more extensive range through Iran, Uzbekistan, Kazakhstan, and Kyrgyzstan. Our TIMETREE estimation indicated the Central Asia group emerged at 11.6 Ma (CI = 14.3–6.2), an estimate comparable with previous results (11.4–10.5 Ma; [5]). Our estimate for S. argentatus (9.4 Ma; CI = 14.3–6.2) is also corroborated by earlier estimates (9.6–8.9 Ma; [5]). Based on the results of our ML tree, it seems reasonable to hypothesize Cen- tral Asia as the schizothoracine center of origin.