Geometric morphometric analyses define riverine and lacustrine species flocks of Himalayan snowtrout (Cyprinidae: Schizothorax) in Nepal

Freshwater fishes in the river and lake systems in the Himalayas and Tibetan Plateau are morphologically diverged but the evolutionary relationship of putative subspecies separated in these freshwater systems has not been explored. Snowtrout (Schizothorax spp.) are minnows (Cyprinidae) broadly distributed in Asia. Body shapes of 3 Lake Rara (northwest Nepal) endemics (S. macrophthalmus, S. nepalensis, S. raraensis) and 2 widely distributed riverine species (S. progastus, S. richardsonii) across 3 drainages in Nepal (i.e. Karnali, Gandaki, and Koshi Rivers) were studied using geometric morphometry. Data were derived from museum voucher specimens/ tissues collected in 1984−1986 and 1996 (Lake Rara). Cartesian coordinates of 18 anatomical points (Type I landmarks) from 528 individuals were digitized; shape variation was then quantified with principal component analysis and visualized with thin-plate splines derived from a Procrustes analysis. Models of shape variation (i.e. taxonomy versus geography) were tested with a multivariate analysis of variance and a morphological distance matrix. Phylogeographic relationships were examined with a haplotype network (N = 115) derived from 1140 base pairs of the mitochondrial DNA cyto chrome b gene, and selected GenBank sequences (N = 5). Koshi River snowtrout diverged morphologically from conspecifics, consistent with the phylogeographic data. In contrast, Gandaki and Karnali River snowtrout grouped by morphotype (upversus downstream) irrespective of geographic origin, yet clustered separately within the haplotype network. Lake Rara snowtrout were morphologically but not genetically distinct, due to incomplete lineage sorting. Morphological and genetic variability in Schizothorax from Nepal represent a mosaic driven by isolation (= vicariance) and specialization (= adaptation), with taxonomy insufficiently reflecting diversity. Additional data are required to appropriately derive management and effective conservation plans.


INTRODUCTION
The Himalayas exhibit unparalleled ecological variation shaped by steep mountains and characterized by fast-flowing rivers with pristine habitats in isolated mountain chains (Pellissier et al. 2018). Ecological gradients along these tributaries are strong drivers of morphological and genetic diversification (He & Chen 2006, Bhatt et al. 2012, Liang et al. 2017. Separation of small communities in remote physiographic settings often promotes endemism in terrestrial species , and is even more pronounced in freshwater species, given more extensive riverine distances and the presence of dispersal barriers such as deep gorges, waterfalls, and interspersed rapids (Guo et al. 2016). In fact, morphological diversity among closely related Himalayan fishes has been associated with habitat type and its complexity (Li et al. 2009).
Despite this, Himalayan freshwater biodiversity remains largely understudied, and hence may be pro-foundly underestimated and not recognized taxonomically, thus impeding effective conservation across the region. Snowtrout (genus Schizothorax), a group of minnows (Cyprinidae) adapted to high-elevation Himalayan streams, epitomize this situation. Although species-rich, no comprehensive taxonomic study of this genus has been conducted, and patterns of diversification across the Qinghai-Tibetan Plateau (QTP) into eastern, southern, and western Asia remain enigmatic. Our study addresses one component of this knowledge gap by evaluating morphological variation of snowtrout in both riverine and lacustrine habitats of Nepal and contrasting these results in a phylogeographic context so as to explore potential drivers of diversification.
Traditionally, comparative morphology and anatomy are the foundations of taxonomy, and species are delimited by one (or more) non-overlapping morphological character(s) (Wiens 2007, Sullivan et al. 2014. Over the past several decades, taxonomic studies have increasingly incorporated molecular ap proaches, due primarily to the ease with which genetic data can be generated (Avise 2004). A distinct advantage of genetic data is their capacity to characterize evolutionary relationships across all taxonomic hierarchies, from distinct populations to species, families, and orders (Blaxter 2004). Although both morphological and molecular data convey evolutionary signals, their juxtaposition offers a unique perspective, albeit a slowly adopted one (Lefébure et al. 2006, Douglas et al. 2007.
Geometric morphometrics (GM) has gained popularity as a new analytical method in that it allows for morphological shape variation to be accurately quantified, such that comparisons with molecular data can be performed with elevated statistical power (Rohlf 1998, Zelditch et al. 2012, Regmi 2019. Body shape in fishes reflects adaptive convergence or divergence across various taxonomic groups in response to particular habitat conditions (Barlow et al. 1997, Jones 2008. Head and body shape in particular are critical for niche separation and reproductive isolation, especially at the population−species interface, and reflect adaptive divergence in feeding, mating, and locomotion (Russo et al. 2008, Ventura et al. 2017. Studies that assay morphological variation in Schizothorax spp. in the Himalaya are either rare and/or non-synthetic, and when executed, often seek to merely discriminate species, but are not concerned with discerning distinct ecological niches, such as upstream versus downstream habitats or lotic versus lacustrine systems. For example, anal fin length seemingly distinguished 2 snowtrout (S. richardsonii and S. progastus) in the Upper Ganges River of northwest India (Fig. 1) (Pandey & Nautiyal 1997). Similarly, a truss analysis involving 14 anatomical landmarks identified significant phenotypic variability in S. richardsonii distributed across central and western Himalayan rivers of India (Mir et al. 2013). A similar approach (Pradhan & Wagle 2013) was used to quantify differences within the same species across 3 tributaries of the Gandaki River, Central Himalayas.
Three Schizothorax species endemic to Lake Rara in northwestern Nepal (i.e. S. macrophthalmus, S. nepal -20 81.000°E 86.000°F ig. 1. Sampling sites for 528 Schizothorax specimens collected in Nepal and analyzed for morphological and genetic variation. Two riverine species were collected in 1984−1986 from 3 river drainages (Karnali, Gandaki, and Koshi), with green indicating upstream sites for S. richardsonii (RIC) and red indicating downstream sites for S. progastus (PRO). Three lacustrine species endemic to Lake Rara (black site, northwestern Nepal) were collected in 1996: S. macrophthalmus (MAC), S. nepalensis (NEP), and S. raraensis (RAR). Panel on right shows geographic location of Nepal, with China to the northeast and India to the southwest. Sample size for each species is provided in Table 1 ensis, S. raraensis) are recognized as a lacustrine species flock, with S. nepalensis and S. raraensis classified as Critically Endangered (Jha & Rayamajhi 2010a,b), yet lack a thorough quantitative morphological assessment. A re-analysis of published meristic data (Terashima 1984) could only discriminate the 3 species according to the number of outer gill rakers (B. Regmi unpubl.). Interestingly, these differences are congruent with intuited trophic divergence among Lake Rara-endemic species, as surmised from field observations of habitat and food preferences (Terashima 1984). Data are even more sparse for snow trout that are more broadly distributed, with descriptions largely generalized to indicate, for example, a 'blunt-nosed' phenotype in upstream versus a 'pointed-nosed' phenotype that occurs downstream (Petr & Swar 2002). The goal of our study was to conduct a morphological analysis of shape variation in snowtrout from Nepal (Fig. 2). To accomplish this, we acquired GM data across 5 putative species (referred to as operational taxonomic units, OTUs) using 18 anatomical landmarks that encapsulate head and body shape variation (Fig. 3). We then applied multivariate morphometric analyses as a means to potentially distinguish OTUs in pairwise comparisons. Samples represented 3 lacustrine OTUs (Lake Rara) as well as 2 riverine OTUs distributed in adjacent river systems (Fig. 1). To evaluate character congruence, shape variation was juxtaposed with molecular data, derived from mitochondrial (mt)DNA cytochrome b gene sequences, a standard approach to evaluate genetic variation within a geographic context (i.e. phylogeo graphy; Avise 2000). Combined, these an alyses quantified the extent of divergence among snowtrout in Lake Rara and contrasted these with delimitation of Schizothorax in adjacent river basins. In doing so, we addressed 2 fundamental questions: (1) Based on shape variation, how morphologically distinct are Schizothorax endemic to Lake Rara from riverine congenerics distributed in the Central Himalayas? (2) Is shape variation consistent with phylogeography?

Study species
Snowtrout (Cypriniformes, Cyprinidae, subfamily Schizothoracinae, genus Schizothorax) are minnows widely distributed across the Himalaya and trans-Himalaya. Contrary to the common name, snowtrout are not related to cold-water fishes generally identified as 'trout' (Salmoniformes). While both are superficially similar in shape, snowtrout lack an adipose fin and possess a subterminal (rather than a terminal) mouth. The common name reflects a distribution that  (Shrestha 1981, 1999, Terashima 1984, Edds 1989, He & Chen 2006. The genus is diverse, comprising ~127 species (Froese & Pauly 2019). Although Cyprinidae are recognized as the largest and most diverse family of fishes, few are considered species flocks (Dimmick & Edds 2002), a situation currently being reevaluated (Nagelkerke et al. 2015, Ayvazyan et al. 2019. The 3 lacustrine snowtrout restricted to Lake Rara in Nepal are a noteworthy exception (Terashima 1984). Two other snowtrout in this study, the more common riverine S. richardsonii and S. progastus, are distributed widely in Central Himalayan drainages (Tera shima 1984, Dimmick & Edds 2002, Edds 2007).

Specimens and tissue samples
Samples for the 5 Schizothorax species were obtained from specimens collected by D. R. Edds in 1984−1986 from 3 major drainages in Nepal (i.e. Karnali, Gandaki, and Koshi Rivers) (Edds 1989, Shrestha & Edds 2012, and in 1996 from Lake Rara (Dimmick & Edds 2002) (Fig. 1). Specimens were accessioned into the University of Kansas Natural History Museum and the Oklahoma State University Collection of Vertebrates. Voucher specimens (N = 528, Table S1) provided images for shape analyses, and ethanolpreserved tissue samples (N = 115, Table S2) were used for genetic analyses (Table 1).

Shape data
High-resolution digital images were captured from 528 museum voucher specimens (Fig. 2) so as to quantify variation across all 5 species (Table 1). Image capture was performed using a Ca non (T3 Rebel) digital camera with a 300 mm lens. A photographic box, de signed for optimal light and ex posure, housed the specimens as they were being photo graphed (LED Pro Box Plus 1419 Photography Lighting Studio, MK Digital Direct). The left lateral view was recorded by placing specimens into a V-shape photarium filled with water which allowed us to configure standardized body position. Dis torted specimens and juveniles were excluded. Determining homologous anatomical points (Type I landmarks) is an important step in GM analyses. Eighteen Type I landmarks ( Fig. 3) were digitized to obtain 2-dimensional Cartesian coordinates (TPSDIG; Rohlf 2004). Landmarks were selected to capture head morphology, i.e. snout elongation, head depth, and placement of the eyes and nostrils. Body shape was primarily defined by origin/insertion of fins (dorsal, pelvic, and pectoral), and upper/lower caudal fin bases, with the urostyle (i.e. posterior end of vertebral column) as the most posterior point (Fig. 3). A tps file (i.e. specific file format for storing morphometric data, Adams & Otárola-Castillo 2013) was generated using TPSUTIL (F. J. Rohlf, http:// life. bio. sunysb. edu/ ee/ rohlf/software.html), which derives the x-/ycoordinates for each image, extracts the list of images,  Table 1. Overview of 528 snowtrout (Schizothorax spp.) specimens collected from Nepal and analyzed for morphological and genetic variability. Samples were identified as 1 of 5 species based on phenotype and location: S. macrophthalmus, S. nepalensis, and S. raraensis are lacustrine forms collected from Lake Rara, while S. richardsonii and S. progastus are riverine forms collected across 3 drainages (Fig. 1). Vouchers: number of museum specimens examined for morphological variation using shape analysis; Tissues: number of DNA samples analyzed for sequence variation across the mitochondrial cytochrome b gene. Collection information is listed in Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ b030 p019 _ supp .pdf and appends the images. An R-based software package (GEOMORPH v3.05; Adams & Otárola-Castillo 2013) was employed to read and superimpose the coordinate data using generalized Procrustes analysis (GPA; Adams et al. 2013).

Potential biases due to digitization errors and biological sex
If shape variation is subtle then digitization errors can overwhelm the signals associated with species and/or geography. To gauge this potential, we digitized 9 individuals of each species and basin 10 separate times for comparative purposes (N = 90), then tested these data using Hotelling's t-square statistic. We first standardized t by the number of features (18 × 2 landmark co-ordinates) and replicates (number of specimens) so as to conform to an F-distribution.
The biological sex of the museum specimens was unknown, and could not be determined via internal anatomy or gonad development. To gauge the potential for sexual dimorphism, we classified individuals with breeding tubercles (Fig. 2) as 'male' (M) and those without as 'sex unknown' (SU). Anecdotal observations (D. Edds unpubl. data) indicate that female snowtrout do not develop breeding tubercules. We then implemented Hotelling's t-square test to evaluate the significance of shape variation between M and SU categorizations.

Molecular data
We employed sequences from the mtDNA cytochrome b gene to reconstruct phylogenetic relationships among OTUs, allowing for genetic variation to be evaluated within a geographic context (i.e. phylogeography; Avise 2000). The gene is suitable to gauge the relationship of recently evolved taxa given its high substitution rate at synonymous sites (Farias et al. 2001, He & Chen 2006, Khan et al. 2017, as is likely the case for Lake Rara endemics. Genomic DNA was extracted from 115 ethanolpreserved tissue samples using the Qiagen DNeasy kit, with minor modification to the manufacturer's protocol (Table 1). The entire cytochrome b gene was amplified using published primers and conditions (Unmack et al. 2009, Houston et al. 2012, visualized on agarose gels, and enzymatically purified. Sanger sequencing was conducted using BigDye (ver. 3.1) chemistry (Applied Biosystems [ABI]) and recorded on an ABI 3730 DNA Analyzer (W.M. Keck Center for Comparative and Functional Genomics, University of Illinois, Urbana/Champaign). Sequences were aligned manually using SEQUENCHER (ver. 5.0, Gene Codes). Five available sequences in GenBank were added (Yang et al. 2015). These sequences represented 2 Lake Rara species (S. macrophthalmus, N = 1, KP712253; and S. nepalensis, N = 2, KP712254 and NC031537) plus S. richardsonii from the Koshi River (N = 2; KP712220, KP712249).

Data analysis
We employed 4 fundamental approaches to analyze our data. We first explored the distribution of total morphological variance (i.e. across mutually orthogonal axes) using principal component analysis (PCA). Components represented morphological variation independent of predetermined hypotheses (e.g. taxonomic identification). We then tested the statistical significance of morphological variation in the context of taxonomy (species) or geography (basins) by employing a multivariate analysis of variance (MANOVA). Hierarchal clustering of species based on the Euclidean distance of mean group vectors (morphological distance matrix, MDM) was implemented as our third approach to estimate differences among species in morpho-space. Finally, we constructed a haplotype network based on the cyto chrome b sequence data, then assigned each specimen to a haplogroup as a means to compare morphological differences within the context of phylogeography (i.e. genetically distinct lineages reflecting taxonomy and/or geography). Although genetic data could only be obtained from archived tissue (N = 115), we assigned all specimens representing a particular OTU from a given location to the same haplogroup, so as to facilitate comparison with our GM analyses.
To visualize morphological variation in tangent space, we derived thin-plate splines from Procrustes analysis to generate deformation grids (DGs,) that visualize shape variation along x-/y-axes and incorporated multiple principal components (GEOMORPH v.3.1.0; Adams et al. 2019) (see Fig. 4). We examined the DGs to infer shape differences among species and basins, and to assess the manner by which each variable contributed to the principal components.
Implementing PCA across different OTU groupings allowed us a detailed visualization of how variables are loaded onto these components (see Fig. 5). The groupings were: (1) all specimens using taxonomic designation (= species) independent of geog-raphy (= basins); (2) only specimens representing the 3 Lake Rara species; (3) S. richardsonii from all 3 basins; and (4) S. progastus from all 3 basins. Because each grouping used a different set of specimens, we implemented GPA prior to visualizing each analysis.
Data were partitioned by species, basins, or both so as to test if body shape was attributable to taxonomy or geography. We employed a MANOVA on the GPA-aligned data using Goodall's F-test. Significance of the F-ratio was tested with 1000 permutations. Obtaining the morphological distance matrix based on OTUs (= putative species), and their geographic distribution, was necessary to determine relationships in standardized shape space. We used the 'mshape' function (Claude 2008) in GEOMORPH to calculate the mean vector of GPA-aligned coordinates for each OTU, and the 'dist' function in R (v.3.4.0) to estimate the Euclidian distance matrix. We used this MDM to hierarchically cluster 'species' and 'species-by-basin' in a dendrogram (see Fig. 6) based on their shape variation, using the single linkage function of 'hclust' in R.
Our haplotype network visualized genetic variation among OTUs within a geographical context. Cytochrome b gene sequences were aligned using MUSCLE V.3.8 (Edgar 2004), and sequences were collapsed into haplotypes in DNASP v.6 (Rozas et al. 2017). A haplotype network was constructed using median-joining (Bandelt et al. 1999), and the distribution of the 5 putative OTUs was visualized across haplogroups using POPART (Leigh & Bryant 2015).

RESULTS
Of the approximately 800 digital images taken of museum voucher specimens, a subsample of 528 was selected for GM analysis after distorted, illpreserved, and juvenile specimens were removed (Table 1, Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ b030 p019 _ supp .pdf). The 2 riverine OTUs were sampled from more locations and were consequently represented by larger sample sizes than the endemic Lake Rara species (Fig. 1).

Biases due to digitization errors and sex
Analysis of shape variation due to digitization errors showed greater variance among species or basins (i.e. between groups) than among individuals within groups (i.e. replicate digitization of the same images), indicating that variation due to digitization errors was negligible. In addition, Hotelling's t-test was not significant between 2 mean vectors of classes (M versus SU; p = 0.39, df = 36), indicating lack of sexual dimorphism among our samples.

Principal components of shape variation across basins and species
Our shape data comprised 2-dimensional Cartesian coordinates (x/y) representing 18 landmarks (Fig. 3) and yielded 36 principal components (PCs). Of these, only 32 had non-zero eigenvalues. However, the 4 zero-eigenvalue PCs were retained because their contribution to computational complexity was negligible, and more importantly, a landmark coordinate with computationally zero-eigenvalues for one axis may have non-zero eigenvalue on another.
Across all species, PC1 explained 30% of shape variation, and PC2 explained 14%, with the first 6 PCs loading 74% of the variation. The majority of variables loaded with equal proportion onto PC1, whereas those for PC2−PC6 associated significantly with a subset of traits. PC2 primarily reflected landmarks associated with head shape (i.e. tip of mouth, ventral end of operculum, anterior end of mouth opening). PC3 was dominated by traits associated with fin placement (i.e. origin and insertion of dorsal and pectoral fins). PC4 defined placement of the eye and the urostyle (= posterior end of vertebral column), whereas PC5 involved orbital edges, neurocranial end, and mouth opening. PC6 defined the insertion of the anal fin.
DGs visualize shape variation. An example is provided in Fig. 4 based on a pairwise comparison along PC1 and PC2 between a lacustrine species (Schizothorax nepalensis, Lake Rara) and a riverine species (S. progastus, Koshi River Basin). The top DG (PC1 = 40%) reflects shape variation along most landmarks (Fig. 4a), whereas the bottom DG (PC2 = 18%) shows more subtle variation associated with the upper end of the operculum, mouth opening, and placement of pectoral fins (Fig. 4b).
Plots contrasting shape variation along PC1 and PC2 among different sets of OTUs is shown in Fig. 5. A comparison of all species reflects broad overlap, with the 2 riverine species separating along PC1 (Fig. 5a). Contrasting only the 3 species endemic to Lake Rara indicates clear separation of S. raraensis, whereas S. macrophthalmus and S. nepalensis are intermingled (Fig. 5b). Within riverine forms, S. richardsonii from the Koshi River clearly separate from those in the Gandaki and Karnali Rivers (Fig. 5c). a b Fig. 4. Visualization of morphological variation between the lacustrine Schizothorax nepalensis from Lake Rara and the riverine S. progastus from the Koshi River, Nepal. Deformation grids (thin-plate spline) illustrate shape changes from the reference form (

Multivariate analysis of shape variation
To explain shape variation, we tested models that reflect hypotheses relating to taxonomy (species) and/or geography (basins). MANOVA revealed low but significant regression coefficients for both parameters. However, influence of 'species' on shape variation was stronger than that of 'basin' (Table 2).
Mean vectors of the Euclidian distance (average configuration; Claude 2008) were derived from pairwise comparisons of either groupings by 'species' alone, or partitioning riverine OTUs by 'basin' (Tables S3 & S4, respectively). Hierarchal clustering of these mean vectors is shown in a dendrogram (Fig. 6). The dendrogram of 'species' (Fig. 6a) clusters riverine S. richardsonii with lacustrine S. nepalensis, and riverine S. progastus with lacustrine S. macrophthalmus. Surprisingly, the third Lake Rara endemic, S. raraensis, clusters separately in shape space (Fig. 6a). A different pattern emerges when riverine S. richardsonii and S. progastus are also partitioned by 'basin,' indicating a geographic context of shape variation (Fig. 6b). This analysis tests for potential effects of geographic isolation among Himalayan drainages, and reflects a more nuanced perspective with regards to Koshi River samples.

Genetic variation
All cytochrome b sequences were 1140 base pairs in length, with nucleotide composition as: T = 29.1%, C = 27.2%, A = 27.0%, G = 16.7%. A total of 167 polymorphic sites were detected, with 20 representing singletons, whereas 147 were parsimonyinformative. Nucleo tide diversity was relatively high (π = 0.012). A total of 89 haplotypes (Table S2) segregated into 3 major haplogroups across the 3 basins (Fig. 7). All Lake Rara samples clustered with those from the Karnali River, i.e. the drainage within which the lake is located. Interestingly, riverine samples did not group in species-specific clusters, but rather associated instead with basinspecific haplogroups (Fig. 7), indicating the presence of a robust geographic structure. Importantly, one haplotype was shared among all lacustrine (Lake Rara) and riverine samples from the Karnali Basin.  (Fig. 3) aligned within each group using generalized Procrustes analyses. Models tested whether variation was partitioned by taxonomy (species only), geography (basin only), or both (species + basin). Geographic location of basins is shown in Fig. 1

DISCUSSION
Freshwater habitats in the Himalaya are globally unparalleled, as they follow rapidly elevating topographies (Pellissier et al. 2018), varying temperatures (Regmi & Lamichhane 2019), and variable rainfall (Shrestha et al. 2019), all within a short geographic range. Heterogeneous habitat conditions and geographic isolation provide the context for allo patric as well as sympatric speciation in remote high-elevation lakes and streams (Li et al. 2009). Our re search focused on 5 snowtrout species from the Central Himalaya, a region replete with geographically remote and isolated habitats. Specifically, our research evaluated the potential for morphological shape variation between riverine species in the Karnali, Gandaki, and Koshi River basins in contrast to an endemic species flock of 3 of lacustrine forms in Lake Rara (Karnali River basin). Given that Lake Rara is recognized as a RAMSAR site (i.e. a wetland designated of international importance under the Ramsar Convention), and the endemic species are considered endangered, quantifying their distinctiveness from riverine OTUs is relevant for appropriate management and conservation strategies.
Species flocks represent fascinating examples of adaptive radiation, and comprise sympatric, monophyletic groups of species restricted to a particular geographic area (Echelle & Kornfield 1984), with African cichlids being the most well-known and extensively-studied example (Ivory et al. 2016, Singh et al. 2017. Species flocks are recognized in other groups of fishes, including salmonids (Douglas et al. 1999(Douglas et al. , 2005, but are surprisingly rare in the Cyprinidae (minnows), the most speciose fish family (Kornfield & Carpenter 1984, Dimmick & Edds 2002). Yet, some notable exceptions have been recently recorded (de Graaf et al. 2008, Ayvazyan et al. 2019).

Drivers of morphological variation in fishes
Fishes exhibit high intra-specific morphological variation that often reflects ecological and/or evolutionary significance (Robinson & Wilson 1996, Clabaut et al. 2007. It is often associated with sexual dimorphism, and well-known examples are found in livebearing fishes (Bisazza & Marin 1995, Regmi et al. 2016) with pronounced differences in body shape between males and females (Hassell et al. 2012). Although sexual dimorphism has been described in Lake Rara Schizothorax (Terashima 1984), none of the seemingly divergent traits were related to head or body shape. We tested for patterns of sexual dimorphism but were unable to statistically discriminate sex based on body shape, and thus did not consider biological sex as a source of bias.
Another driver of morphological variation in fishes relates to trophic ecology, as documented in examples of adaptive radiation, including African cichlids, as well as endemic lacustrine species of salmonids: Coregonus (Turgeon et al. 1999, Douglas et al. 2005 and Arctic char Salvelinus alpinus (Brunner et al. 1998). Head morphology is often associated with feeding performance, and modifications to the shape of the head, mouth position, and other such features frequently drive trophic diver-  Table S2. Colors represent unsupervised (K = 5) clusters gence (Ornelas-García et al. 2018). Lake Rara snowtrout occupy different trophic niches, and hence are considered a species flock: S. macrophthalmus is planktivorous, S. nepalensis is herbivorous and S. raraensis is insectivorous (Tera shima 1984). This trophic diversification is reflected in head morphology: S. macrophthalmus and S. raraensis have an oblique mouth, whereas the mouth of S. nepalensis is subterminal, an adjustment for scraping sessile algae. However, these observations are based on qualitative assessments (Terashima 1984), and have not as of yet been quantified. Our statistical analysis of shape variation in snowtrout revealed that PC2 loadings were strongly associated with head shape (Fig. 4b). We uncovered a distinct morphological trajectory for S. raraensis, as compared to other snowtrout distributed across the Central Himalaya (Fig. 6a). In contrast, other Lake Rara endemics were similar to riverine species with regard to head shape, with S. nepalensis comparable to the blunt-nosed S. richardsonii, whereas S. macrophthalmus approximated S. progastus in that both demonstrated an elongated snout (Fig. 2). However, morphological similarity could be due to convergence, and does not necessarily reflect a phylogenetic signal.
The riverine S. richardsonii is distributed upstream in headwaters of many Himalayan rivers, and represents a detritophagic bottom feeder that harvests algae growing on gravel and boulders (Sharma 1984, Wagle et al. 2015. The other riverine species, S. progastus, is found in downstream reaches, and feeds on insects and diatoms suspended in the water column (Edds 1989, 1993, Shrestha 1999). These differences are clearly reflected in head shape, with similar diet preferences seemingly responsible for trophic morphology and head shape variation in Lake Rara species as well.
Other morphological characteristics that strongly associated with PC loadings were traits related to locomotion, such as the placement of dorsal and pectoral fins (Fig. 4). This suggests that morphological divergence of S. richardsonii versus S. progastus reflects adaptations to high-gradient streams in upper catchment areas versus gently flowing rivers in the foothills of the Himalaya. Although MANOVA of multivariate shape variation based on both taxonomy (species) and geography (basin) was significant, we found 'species' to be a statistically stronger predictor than 'basin' (Table 2), again arguing for vicariant (geographic isolation) effects within each species that supersede a taxonomic diagnosis based on morphology alone.

Concordance between morphological and molecular variation
To date, no study has combined morphological and molecular data to delimit freshwater fish diversity in isolated headwaters of the QTP, Himalaya, and trans-Himalaya. However, several have examined phylogeography of snowtrout in the region. He & Chen (2006) employed cytochrome b data to examine phylogenetic relationships of 30 species distributed in 9 rivers originating on the QTP. Genetic clades were found to be inconsistent with traditional morphological-based taxonomy, and instead reflected geographic associations within basins. This, in turn, suggests that incomplete lineage sorting may be responsible for the discrepancies between taxonomy and molecular data. However, the pattern of diversification presented by He & Chen (2006) was not biologically (or geomorphologically) interpretable, and thus the evolutionary history of Schizothorax remains rather opaque. A more localized study, based on sequence analysis of the mitochondrial D-loop region and flanking tRNAs (Dimmick & Edds 2002), indicated substantial genetic divergence of snowtrout among the Karnali, Gandaki, and Koshi river basins in Nepal.
Our study of morphological shape variation showed partial concordance with phylogeographic signals based on cytochrome b data. Strong geographic structuring of haplotypes (Fig. 7) supported previous findings that Schizothorax in isolated headwaters of the QTP and trans-Himalaya are genetically distinct according to drainage. Genetic divergence is most pronounced in the Koshi River basin, a pattern also reflected in morphological data (Fig. 6b).
In contrast, the 3 Lake Rara species cluster with the 2 Karnali River basin riverine OTUs (Fig. 7), a result inconsistent with our morphological results (Fig. 6b). However, incomplete lineage sorting (i.e. retention of ancestral polymorphism) is the most likely argument to explain this lack of concordance. This indicates that insufficient time has passed for distinct mtDNA lineages to emerge due to recent origin of Lake Rara, a well-established phenomenon in phylogeographic studies (Avise 2000, Douglas & Douglas 2010).

CONCLUSIONS
Our data suggest that multiple evolutionary mechanisms have driven diversification of Himalayan snowtrout, with the combination of geographic isola-tion and ecological specialization yielding a mosaic of morphological and genetic variability. This raises the potential that phenotypic similarities as seen herein may result from homoplasy, as driven by the strong selective pressures produced by harsh environments, such that morphologically similar forms have evolved repeatedly across different drainages (Douglas & Brunner 2002). If so, current Schizothorax taxonomy might drastically underestimate actual diversity, with obvious implications for species-centric conservation plans. To fully understand the evolutionary history of snowtrout in the Central Himalaya, a more in-depth genetic evaluation is warranted, using both morphological and genomic markers. One confounding issue is that Schizothorax species are polyploid (Yang et al. 2015), and this limits phylogenetic analyses to mtDNA data (given its haploid nature), or to nuclear genes, either cloned and isolated as single-copy markers, or genomic markers that can be bioinformatically phased.
Although endemic species in Lake Rara appear to be morphologically divergent, molecular data do not support their distinctiveness. However, given the geologically young age of Lake Rara, a likely scenario is rapid speciation resulting in distinct morphotypes that display shallow genetic divergence, and given this, the 3 endemic OTUs should be considered distinct species for conservation and management purposes.
This situation warrants further investigation, potentially using genomic approaches that can distinguish between hybridization, incomplete lineage sorting, or reticulated evolution as potential causes of the morphological and genetic incongruence found between lacustrine and riverine snowtrout in Nepal. This of course argues for polyploid genomic data to be analytically manipulated for subsequent analysis.
An important outcome of this study is the substantial differentiation of both Schizothorax species in the Koshi River basin from conspecifics in other drainages, which argues for their recognition as evolutionarily significant units (Avise 1989, Moritz 1994, Douglas & Brunner 2002, and should be so presented in future management plans to guide conservation efforts that will ensure the persistence of this intra-specific diversity uniquely adapted to Himalayan rivers.