Multiple species of grayling (Thymallus sp.) found in sympatry in a remote tributary of the Amur River

Large‐scale phenotypic and genetic studies of the salmonid genus Thymallus (grayling) in the Palaearctic suggest that most major phylogeographic lineages represent good biological species. Evaluating such a premise in areas involved in palaeo‐hydrological dynamics where multiple species are found in sympatry should serve to assess the level of reproductive isolation, the traditional sine qua non of species recognition. Molecular sequence (mtDNA) and microsatellite (nDNA, seven loci) analysis of grayling in the upper Bureya River watershed support the occurrence of three distinctive species of grayling living in sympatry in this large oligotrophic tributary of the Amur River. One of these lineages is primarily found throughout the Lena River basin and is recognized as Baikal‐Lena grayling Thymallus baicalolenensis; the second, the upper Amur grayling Thymallus grubii is found over large areas of the Amur catchment including the entire headwater region; and the third, the Bureya grayling Thymallus burejensis is endemic to the study area. A limited number of hybrids were identified, primarily between Baikal‐Lena grayling T. baicalolenensis and Bureya grayling T. burejensis with little to no signs of introgression among non‐hybrid individuals. Morphological distinctiveness among populations of these species living in sympatry was greater than between populations living in allopatry, suggesting character displacement. Divergence estimates among taxa range up to 6.2 MY, and allopatric origins for all three species’ is suggested. To our knowledge, this is the first data‐based confirmation of three species of grayling living in sympatry.

The Bureya River is the second largest left tributary of the Amur River (after the Zeya River) with a mean flow at its mouth of 943 m 3 /s, a length of 739 km (including its two headwater branches) and a catchment area of ca. 70,700 km 2 . It has a large (230 km long) impoundment built in 2003 in its middle reaches but is otherwise a relatively remote river with few settlements and a strict protected area (Bureinsky Zapovednik) in its headwaters. This nature reserve is over 350,000 ha in size, characterized by a subarctic climate with mountain tundra and taiga forest landscapes. At the time of our sampling (August 2007), two species were reported to occur in the Bureya River, T. grubii and T. burejensis. However, it was suspected via informal observations, that both T. tugarinae and Thymallus baicalolensis also occurred in the basin, the latter in the headwater reaches. We investigated this reach of the Bureya River with the intention of clarifying how many grayling species occur in sympatry in the upper Bureya River, and whether or not hybridization and introgression occurs among them.

| Sampling
A total of 81 graylings were sampled via angling just below the Nature Reserve in both the main stem of the Bureya River and a small right hand tributary (Umalta makit; Figure 1). In the field, adult specimens were identified based on their prominent dorsal fin and body colouration as described in Antonov (2004), Knizhin, Weiss, Antonov, and Froufe (2004) and Knizhin, Weiss, Bogdanov, Samarina, and Froufe (2006) as upper Amur, Bureya or Lena-Baikal grayling. Small fin-clips were preserved in 96% ethanol for genetic analysis, while a subset of full specimens were preserved in 4% formalin for subsequent morphological analysis in Irkutsk University.

| DNA extraction and sequencing
Whole genomic DNA was extracted using a standard high-salt protocol (Sambrook, Fritsch, & Maniatis, 1989). The complete mitochondrial (mtDNA) control region and portions of flanking tRNA gene regions were amplified in all individuals using the LRBT-25 and LRBT-1195 primers and PCR conditions described in Uiblein Jagsch Honsig-Erlenburg and Weiss (2001). Amplified DNA templates were purified and the forward half of the template sequenced using the forward PCR primer. The quality of the sequences was analysed with the software Chromas 2.6.5 (Technelysium Pty Ltd), and all new sequences have been deposited in GenBank database under the accession numbers MK746000-MK746080 (Table S1). The sequences were aligned with the MAFFT multiple sequence alignment algorithm (Katoh & Standley, 2013) and edited in the program Bioedit 7.0.5.3 (Hall, 1999), together with all available sequences on GenBank from the three species identified in the field as well as the yellow-spotted grayling T. flavomaculatus and a few representative sequences from T. tugarinae and the broad-ranged Thymallus arcticus (Table S1).

| Phylogenetic analyses
Phylogenetic analyses were performed using both maximum likelihood (ML) and Bayesian inference (BI) methods. An individual sequence from each of two other salmonid genera (Brachymystax and Hucho) were used as an outgroup. ML analyses were performed using RAxML-HPC2 Workflow on XSEDE version 8.2.10 (Stamatakis, 2016) with 1,000 bootstrap replicates and model GTR + G. jModelTest2 on XSEDE 2.1.6 (Darriba, Taboada, Doallo, & Posada, 2012) was used to calculate the best nucleotide substitution model, which revealed itself as the HKY + I model, and this was used in the BI analyses. BI analyses were performed in MrBayes on XSEDE (3.2.6) (Ronquist & Huelsenbeck, 2003) with two independent runs (10 7 generations with a sampling frequency of one tree for every 100 generations), each with four chains (three hot and one cold). The program Tracer version 1.7.1 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018) was then used to remove all the trees in a burn-in of 20,000. All the phylogenetic analyses have been performed on the portal CIPRES science gateway. Sequences were also joined in unrooted networks using a 95% connection limit criterion implemented in TCS version 1.21 (Clement, Posada, & Crandall, 2000). Sequence and nucleotide diversity measures were calculated in DnaSP version 6.12.01 (Rozas et al., 2017). Sequence divergences (uncorrected p distance) were assessed using MEGA version 10.0.5 (Kumar, Stecher, Li, Knyaz, & Tamura, 2018).

| Microsatellite analyses
Allelic variation at seven microsatellite loci was analysed for population genetic analyses. The loci were all chosen from F I G U R E 1 Map of the ranges of all three species investigated in this study, and close-up of the Bureya River catchment, where all three species were found in sympatry. Green: Thymallus grubii; red: Thymallus burejensis; black: Thymallus baicalolenensis. The red rectangle represents the sampling area of the present study. The dots in the enlarged box represent records of occurrence reported in Knizhin (2011, 2014), Antonov and Mikheev (2016) and Ma and Jiang (2011) [Colour figure can be viewed at wileyonlinelibrary.com] published reports: BFRO009(Thy18) (Sušnik, Snoj, & Dovč, 1999a), (BFRO018(Thy54), BFRO015(Thy72) (Sušnik, Snoj, & Dovč, 1999b), BFRO004(Thy2/2) (Snoj, Sušnik, Pohar, & Dovč, 1999), Cocl23 (Rogers, Marchand, & Bernatche, 2004) and Tar104, Tar108 from (Diggs & Ardren, 2008). Loci were amplified in duplex or triplex PCRs using the same conditions as in Weiss et al. (2013). Allelic richness and diversity, observed (H O ) and expected (H E ) heterozygosity, and inbreeding coefficients (F IS ) were estimated using Genetix 4.03 (Belkhir, 2004), FSTAT v 2.9.3.2 (Goudet, 2001) and GENEPOP online version 4.2 (Raymond & Rousset, 1995). To estimate genetic differentiation among species with the microsatellite markers, we used Rst (Slatkin, 1995), as it incorporates allele size differences, which we expect to be significant when analysing different species. To evaluate the genetic relationship among all individuals as well as potential hybridization and introgression between the presumed species present, we applied the model-based approach implemented in STRUCTURE v.2.3.3 (Pritchard, Stephens, & Donnelly, 2000). The posterior probabilities of K were estimated assuming uniform prior values of K between one and four (the number of presumed species + 1) with 15 replicates of 500,000 iterations following a burn-n period of 100,000 iterations. A non-admixture model and an independent allele frequency in each population were assumed.
An alternative non-model-based clustering approach, discriminant analysis of principle components (DAPC) (Jombart, Devillard, & Balloux, 2010) was carried out using the adegenet (Jombart, 2008) and ade4 packages (Dray & Dufour, 2007) written for R. A K-means clustering algorithm of principle components was first carried together with the reiteriative find. clusters function. We ran the K-means clustering 100 times with a maximum of 100,000 iterations and setting the maximum number of clusters to 15. The Bayesian information criterion (BIC) was used to choose the optimal model (i.e. number of clusters, K). Within the DAPC, the optima.score function was used to choose the optimal number of components to be retained for discriminant analysis (Jombart et al., 2010).

| Hybrid analysis
Analysis of the hybrids was carried out with the software HybridLab v. 1.0 (Nielsen, Bach, & Kotlicki, 2006) and New Hybrids v. 1.1 (Anderson, 2002). HybridLab was used to generate hypothetical hybrids for each pairwise combination of the three taxa, and New Hybrids was used to test for membership of genotyped individuals to several simulated hybrid classes. For each potential combination of species, a simulated population of 30 hybrids of F1, F2 and first backcross (BX1) was generated.

| Morphological analysis
A total of 52 individuals were preserved in 4% formalin and transported to Irkutsk State University for morphological analysis. A total of 30 morphometric and 12 meristic measurements as described in Knizhin and Weiss (2009) (see also variable list in Table S2) were taken. These morphometric measurements followed Pravdin (1966) with the following minor modifications: we added the measurements head depth at the centre of the eye, interorbital width, upper jaw depth, body thickness and length of the posterior and anterior edges of the dorsal fin, and we left out total length, length of the middle part of the head and the longest ray of the dorsal fin. Additional meristic characters included branched rays of the pectoral and pelvic fins and separate counts of the branched and unbranched (as well as total) rays. We note that total number of vertebrae were counted excluding the urostyle, and the total number of gill rakers including rudimentary gill rakers. All measurements were done on adult fish and were standardized by fork length following, for example, Knizhin and Weiss (2009). Based on a-priori field diagnosis, the sample included 31 Bureya grayling T. burejensis, 17 upper Amur grayling T. grubii and four Baikal-Lena grayling T. baicalolenensis based on the key of Knizhin (2009), largely relating to body and dorsal fin coloration. A stepwise discriminant analysis was carried out in SPSS version 22.0 using default criteria for variable entry and removal (alpha probabilities). Prior probabilities were computed from group sizes, and a within-group covariance matrix was generated.

| Phylogenetic relationships
The final CR alignment was 582 bp long and contained 200 sequences. There were 145 variable (polymorphic) sites and 121 parsimony informative sites. As the topologies of ML and BI phylogenetic analysis were congruent, only the results of phylogenetic BI analyses are shown ( Figure 2a). Individual sequences of the four potential species occurring in the Amur drainage, along with the widespread Arctic grayling T. arcticus, were collapsed into triangles for visual efficiency. Haplotypes from three of the four grayling species (all but lower Amur grayling, T. tugarinae) potentially occurring in the Bureya River were recovered from samples in the study area. With individual haplotypes of each species collapsed into triangles, the first split in the tree following the outgroup is between the lower Amur grayling T. tugarinae and the remaining four grayling taxa. The next split is between upper Amur grayling T. grubii and the remaining taxa. Finally, the Bureya grayling T. burejensis forms a sister group to the clade containing both T. arcticus and T. baicalolenensis.
Pairwise p-distances ranged from a low of 1.6% between T. baicalolenensis and T. grubii to a high of 6.2% between T. grubii and T. baicalolenensis (Table 1). The three unrooted networks contrast in overall diversity. A single haplotype is shown for T. burejensis and is thus far the only one reported (Figure 2b). Only two of the 10 displayed haplotypes from T. baicalolenensis stem from this study site (Figure 2c). A total of 29 haplotypes is shown for T. grubii with six haplotypes coming from this study and others coming from distant locations in both Russia and Mongolia (Figure 2d). Connected to the central haplotype of T. grubii with a minimum of seven mutational steps is a cluster of haplotypes (shown in yellow) representing T. flavomaculatus (Figure 2d).

| Microsatellites
The Delta K approach to the STRUCTURE analysis gave an unambiguous signal of K = 3, and all individuals except for five showing potentially admixed group membership, corresponded to their presumed species status based on either phenotype or mtDNA clade assignment or both (Figure 3).
Removing the five presumed admixed genotypes (based on a cut-off Q-value of 0.95) from the data set, the expected heterozygosity ranged from 0.709 in T. burejensis (N = 52)  and 0.720 in T. grubii (N = 34) to 0.643 in the eight samples of T. baicalolenensis (Table S3). F is values were all positive and for T. burejensis highly significant (p < .001; Table S3). The number of private alleles for each species was relatively high ranging from 19 (T. baicalolenensis) to 41 and 45 for T. grubii and T. burejensis, respectively (Table S3). Pairwise Rst's between species were high ranging from 0.724 to 0.846 (Table S4). Retaining the putative admixed genotypes in the data set did not result in any significant change in the overall levels of allelic variation and species differentiation (data not shown). For the five presumed admixed genotypes, four are presumably hybrids between T. burejensis and T. baicalolenensis and one between T. burejensis and T. grubii. For three of these putative hybrids, both mtDNA and microsatellite genotypes were available, and in all cases, the female was T. burejensis (the male was T. baicalolenensis in two of these individuals and T. grubii in one). For the remaining two individuals, only microsatellite data was available, and the highest contribution was from T. baicalolenensis across all replicates, with some instability concerning the per cent contribution from T. burejensis. All individuals, with the exception of the five putative hybrids, were assigned to parental species with a probability of > 99%. The four putative hybrids between T. burejensis and T. baicalolenensis were all assigned to the F2 class with a mean probability of 84%, while the remaining putative hybrid between T. burejensis and T. grubii was assigned to the F2 class with a probability of 66%.

| Morphological analysis
Nine of the 30 morphological variables and three of the 12 meristic variables were retained in the discriminant analysis (Table S5). These included least body depth, total dorsal rays, length of anal fin base, the number of pores along the lateral line, the post-dorsal distance, the length of the posterior edge of the dorsal fin, the maximum body depth, the length of the caudal peduncle, the number of pyloric caecae, the ventro-anal distance and the number of ventral fin rays. The canonical correlations for both of the two functions were very high (0.998 and 0.979), and both were highly significant (p < .001), with the first function distinguishing T. burejensis from the other two species being strongly influenced by the minimum body depth, the length of the anal fin and the number of dorsal and ventral fin rays ( Figure S1). The second function distinguishes T. grubii and T. baicalolenensis from one another and was most strongly influenced by ante-dorsal distance, the maximum body depth and the number of ventral fin rays ( Figure S1 and Table S5). All individuals were assigned to their a-priori diagnosed species with a probability >99.9%.

| Additional analyses
Having identified the species of grayling living in the upper Bureya River, we augmented our phenotypic analysis in order gain perspective on the inter-and intraspecific variation in these taxa. Examining the morphological distinction of species that live both in sympatry and allopatry further allows for evaluation of potential competitive pressure or character displacement (Brown & Wilson, 1956;Pfennig & Pfennig, 2009). We note that our analysis must be limited to identifying the pattern of character displacement and not the ecological processes potentially underlying these patterns (Schluter, 2000). The additional data contained the same measurements as reported above for 237 individuals of Amur River drainage grayling (four species) as well as 352 individuals of T. baicalolenensis. The measurements stem from reports in Knizhin et al. (2004); ; Knizhin, Weiss, Bogdanov, and Kopun (2008); ; and Weiss et al. (2006) with some unpublished additions (Table S6). Figure S3 shows the locations of all sample sites for this additional analysis. With these data, two separate discriminant analyses were carried out. The first, targeting interspecific differences, selected a large number (31 of 52) of characters (Table S7) to derive four functions in discriminating all five species analysed. Correct assignment ranged from 93% to 100% for each species. We especially note the relatively high rate of successful classification with some overlap, particularly between T. baicalolenensis and T. grubii ( Figure S4). Correct classification of T. baicalolenensis, despite its occurrence in three drainage basins, was 98.6% with all five misclassified individuals coming from small lacustrine habitats. The second analysis aimed at evaluating the intraspecific variation of T. baicalolensis, applying an a-priori assignment to the three different drainage basins (Amur, Enisei and Lena). Here, the discriminant analysis showed clear differentiation among the drainages ( Figure S5) based on two functions and 15 (of 52) selected characters (Table S8). Here again, we find the same individuals mentioned above from the Gramna-Goudzhekit-Tiya misassigned, in this case to the Lena basin, although they are found in a Baikal tributary, meaning the Enisey catchment. The overall correct assignment is 100% for individuals from the Amur basin, over 99% for 337 individuals from the Lena basin and over 90% for the 40 individuals from the Enisey basin.

| DISCUSSION
We verified the sympatric occurrence of three species of the genus Thymallus using mtDNA, bi-parentally inherited molecular markers and morphological analysis. This is also a rare example of three congeneric salmonids living in close sympatry in a riverine system. Most reports of multiple species of Eurasian salmonids in sympatry are from lacustrine environments (Douglas, Brunner, & Douglas, 2005;Duguid, Ferguson, & Prodohl, 2006;Klemetsen, 2010;Vonlanthen et al., 2012), although the Zeta River in Montenegro (Mrdak, Nikolić, Tošić, & Simonović, 2012) and the Neretva River in Bosnia Herzegovina (Razpet, Sušnik, Jug, & Snoj, 2007) harbour multiple species of Salmo. For our data set, as no pair of the three species reported reflect a phylogenetic sister relationship (Froufe, Knizhin, & Weiss, 2005; and Figure 2 herein), it is clear that these species did not evolve in sympatry but rather reflect contact of lineages that have primarily evolved in allopatry. Historically, grayling throughout the interior Lena basin was referred to as T. arcticus pallasi, but this has shown to be based on a misconception discussed in Knizhin, Kirillov, et al. (2006) and Weiss et al. (2006). Currently, the small-bodied grayling found throughout the Lena basin, excluding the Delta region, is recognized as T. baicalolenensis stemming from a subspecies description in Matveev et al. (2005) and formal recognition in Bogutskaya Naseka Shedko Vasil'eva and Chereshnev (2008). This species, commonly known as Baikal-Lena grayling, corresponds to the lineage reported in (Koskinen Piironen and Primmer (2001) as T. arcticus (Lena basin) and subsequently described in Weiss et al. (2006); Knizhin, Kirillov, et al. (2006); ; and Knizhin et al. (2008) as the upper Lena grayling. While Knizhin et al. (2008) recognize that the Baikal-Lena grayling also occurs in upper tributaries of Lake Baikal, the species occurrence in some upper tributaries of the Amur basin, as reported here and in sympatry with one or more Amur basin species has not been as well investigated (but see Antonov & Mikheev, 2016). The palaeo-hydrological scenario leading to the distribution of the Baikal-Lena grayling in the Lena and Baikal (i.e. Enisey) basins is clear, as Baikal's major outflow during the Pliocene or early to mid-Pleistocene was the Lena River via the palaeo-Manzurka outflow (Arzhannikov et al., 2018;Florensov, 1978;Ivanov et al., 2016;Mats, 1993). Subsequent tectonic activity as well as shifting water levels led to a series of outflows and/or mega-floods connecting Lake Baikal to the Enisey basin. These outflows are thought to have first occurred through the Irkut River and then through the Primorkiy Ridge, eventually forming the now contemporary lake outlet and tributary to the Enisey, the Angara River (Arzhannikov et al., 2018;Ivanov et al., 2016;Mats, 1993; see also Koskinen, Piironen, & Primmer, 2001, related to colonization of Thymallus). Thus, the Baikal-Lena grayling is a relict in the Baikal drainage, confined to headwater tributaries and perhaps displaced in Lake Baikal itself by the Baikal grayling Thymallus baicalensis. The presence of the Baikal-Lena grayling in the Amur basin is probably based on less dramatic circumstances as various headwater tributaries of the Amur and Lena basins south-east of the Sea of Okhotsk come into close contact. Here, several river capture events undoubtedly took place throughout the Pleistocene epoch if not into the Holocene (Antonov, 2012;Antonov & Mikheev, 2016;Groswald, 2009). In the Bureya River, presence of the upper Amur grayling has been long reported and is not unusual, as this species is widespread throughout the vast upper Amur catchment including all tributaries (Antonov & Mikheev, 2016;Knizhin et al., 2004). Most interesting is the phylogenetic position of the Bureya grayling, which to date is reported as an endemic to the upper Bureya River (Antonov, 2004;Antonov & Mikheev, 2016). The most plausible explanation for the phylogenetic placement of this unique taxon is that it is most closely related to a common ancestor of T. arcticus and T. baicalolenensis but evolved in isolation in the upper Bureya catchment or a nearby isolated catchment that has since been fused with the Bureya drainage (see also Antonov & Knizhin, 2011). When calibrated estimates of divergence rates for the mtDNA CR are applied (Froufe et al., 2005;Koskinen et al., 2002), the Bureya grayling diverged from its most closely related congeneric 3-4 million years ago. Thus, it is likely that the isolation and evolution of the Bureya grayling occurred during the late Pliocene. Independent evidence supports that at least portions of the Bureya basin have been isolated and served as a glacial refuge. In addition to several endemic plant species, for example Ribes burejensis, Erigeron burejensis and Taraxacum lineare, recently an endemic chironomid Abiskomyia korbokhon (Makarchenko, 2017) and a vole of the genus Alexandromys was reported (Sheremetyeva, Kartavtseva, Vasiljeva, & Frisman, 2017). The basin experienced two Pleistocene glaciations during which many periglacial lakes were formed as well as river network rearrangements (Antonov & Knizhin, 2011;Groswald, 2009;Ivashinnikov, 2010;Komatsu et al., 2016;Savrasov, 1949). Currently, the Bureya River drainage experiences both perennial and seasonal permafrost up to 500 m in depth (RSWSU, 1966), with the river considered to be the coldest Amur River tributary (current records are not available, but historical data indicate a mean summer temperature in the lower reaches of Bureya of 11 degrees C). Small lakes located in several depressions in the region are inhabited by fish (Rhynchocypris percnurus, Percottus glenii, Misgurnus mohity, Carassius gibelio), whose populations are likely to be significantly isolated from the main areas in the lower reaches of the Bureya and the Amur Valley. Blunt-snouted lenok Brachmystax tumensis was found in Korbokhon Lake, with ecological and morphological differences from those populations inhabiting the Bureya River (Antonov, 2017a). Additionally, a unique morphological form of the lake minnow R. percnurus was reported from a single lake in the Bureya basin (Antonov, 2017b).
We speculate that the evolution of the three species of Thymallus living in the upper Bureya system evolved in allopatry. For most individuals analysed in this study, there was no evidence (either from morphology or microsatellitebased STRUCTURE analysis) of introgression, although the hybrid analysis allocated a few individuals to at least an F2 generation. This either reflects real introgression at a very low level, or poor resolution due to the limited number of samples and loci analysed. Nonetheless, at least one hybrid was identified due to phenotypic colouration, as the species are very easily diagnosed in the field ( Figure S2). We also note that there was a general trend in our samples for a lower observed versus expected heterozygosity in all three species, and this difference was highly significant in T. burejensis (Fis = 0.190, p < .001, Table S3). This again could be a reflection of the application of a relatively low number of markers to three relatively divergent species, or alternatively population structure in the sample.
The augmented phenotypic analysis focusing on both inter-and intraspecific variation provided some interesting results concerning evolution in sympatry and allopatry, as well as concordance with palaeo-hydrological views. First, four of the five misclassified T. baicalolensis stem from the Gramna-Goudzhekit-Tiya tributary system on the north shore of Lake Baikal (Knizhin et al., 2008).
Little data are available from this river, but its headwaters rise in a diffuse glacial valley with connection to the Lena River headwaters, most likely promoting Holocene exchange between the river drainages. The intraspecific analysis of T. baicalolenensis likewise assigning these individuals to the Lena River drainage supported this inference. Concerning other species, nine of 113 individuals a-priori identified as T. grubii were assigned to T. baicalolenensis, with all stemming from a single population (Zehlikhen) in the Ingoda catchment of the upper Amur River draiange, where recent reports of T. baicalolenensis have been made (Antonov & Mikheev, 2016). No population genetic data are available from this region, and thus, this may be an additional hybrid zone between T. grubii and T. baicalolenensis. The lower Amur grayling, which is found in sympatry with T. flavomaculatus (Froufe et al., 2003), was classified with 100% success (N = 74). In summary, this shows that interspecific phenotypic differentiation in these grayling is much stronger ( Figure S1) when species are in sympatry, compared to when they are in allopatry ( Figure S5), the expected pattern of character displacement. Despite the high correct assignment of T. baicalolenensis to individual catchments, the characters resulting in this intraspecific differentiation were different than those supporting interspecific variation. For T. baicalolenensis, the upper jaw length, eye diameter and anal fin length were the most heavily weighted characters (Table S8) in the discriminant analysis, characters that did not weigh heavily in the interspecific comparisons.
This pattern may result from the fact that Lena basin populations of T. baicalolenensis are primarily found in the absence of a congeneric competitor, whereas the Enisey and Amur populations share habitat with congenerics. With the Enisey populations of T. baicalolenensis (i.e. populations found in the tributaries of Lake Baikal), this sharing has not yet been confirmed to be explicitly sympatric, as no specimens have been recorded directly in the lake, where T. baicalensis is commonly found.
Similar to the morphometric analysis limited to the Bureya River drainage, there was no single diagnostic character for differentiating any of these species, nor the river drainage populations of T. baicalolensis. The raw data, or distribution ranges of all characters of all species analysed in this study, can be made available for taxonomic research upon request to the authors.
The sympatric occurrence of multiple species of Thymallus in Asian rivers is relatively unique among salmonid fishes, yet the mechanisms maintaining this diversity are not yet clear and might differ from the typical trophic niche differentiation often reported among salmonid taxa in lacustrine environments. Certainly, the lack of specific diagnostic morphometric characters for these species may underscore the lack of trophic niche partitioning. We suggest that a combination of mate choice, habitat niche differentiation and perhaps genetic incompatibility is the candidate mechanisms underlying this unique riverine diversity.

ACKNOWLEDGEMENTS
We dedicate this manuscript to our dear friend and colleague Dr. Igor Borisovich Knizhin, who accompanied us on this field expedition to the Bureya River and took all the morphometric measurements presented in this study. Igor was Head of the Department of Invertebrate Zoology and Hydrobiology at Irkutsk State University in Irkutsk, Russia. He was renowned for his work on grayling systematics and taxonomy and was responsible for many changes in our overview of grayling diversity throughout Eurasia, including the first descriptions of two new species (T. tugarinae and Thymallus svetovidovi) and one subspecies (T. grubii flavomaculatus) that has been subsequently treated as a species. He was author or co-author of over 100 scientific publications and was a talented teacher who was popular among students. His dedication, humour, integrity and enthusiasm for furthering our scientific understanding of freshwater fishes were infectious and helped promote the careers of several of the co-authors on this manuscript. Igor passed away in January of 2016 after a brief illness and is dearly missed by all. He left me (the first author) with the plea to continue with research on grayling and complete our unfinished goals. We thank Duarte V. Gonçalves for assisting with the map construction, Sergey Alekseev and Nina Bogutskaya for help understanding the morphometric scheme of Pravdin (1966), Theodore Kopun for production of some of the raw genetic data and Gernot Englmaier for reading earlier versions of the text. We thank one anonymous reviewer and N. Bogutskaya for their constructive comments on earlier versions of this manuscript. This work was supported by the Foundation for Science and Technology (SFRH/BD/139069/2018) to G. S. and by Strategic Funding UID/Multi/04423/2019 and by the Institute of Biology at the University of Graz.