- Research article
The evolutionary history of sharp- and blunt-snouted lenok (Brachymystax lenok (Pallas, 1773)) and its implications for the paleo-hydrological history of Siberia
BMC Evolutionary Biologyvolume 8, Article number: 40 (2008)
Broad-scale phylogeographic studies of freshwater organisms provide not only an invaluable framework for understanding the evolutionary history of species, but also a genetic imprint of the paleo-hydrological dynamics stemming from climatic change. Few such studies have been carried out in Siberia, a vast region over which the extent of Pleistocene glaciation is still disputed. Brachymystax lenok is a salmonid fish distributed throughout Siberia, exhibiting two forms hypothesized to have undergone extensive range expansion, genetic exchange, and multiple speciation. A comprehensive phylogeographic investigation should clarify these hypotheses as well as provide insights on Siberia's paleo-hydrological stability.
Molecular-sequence (mtDNA) based phylogenetic and morphological analysis of Brachymystax throughout Siberia support that sharp- and blunt-snouted lenok are independent evolutionary lineages, with the majority of their variation distributed among major river basins. Their evolutionary independence was further supported through the analysis of 11 microsatellite loci in three areas of sympatry, which revealed little to no evidence of introgression. Phylogeographic structure reflects climatic limitations, especially for blunt-snouted lenok above 56° N during one or more glacial maxima. Presumed glacial refugia as well as interbasin exchange were not congruent for the two lineages, perhaps reflecting differing dispersal abilities and response to climatic change. Inferred demographic expansions were dated earlier than the Last Glacial Maximum (LGM). Evidence for repeated trans-basin exchange was especially clear between the Amur and Lena catchments. Divergence of sharp-snouted lenok in the Selenga-Baikal catchment may correspond to the isolation of Lake Baikal in the mid-Pleistocene, while older isolation events are apparent for blunt-snouted lenok in the extreme east and sharp-snouted lenok in the extreme west of their respective distributions.
Sharp- and blunt-snouted lenok have apparently undergone a long, independent, and demographically dynamic evolutionary history in Siberia, supporting their recognition as two good biological species. Considering the timing and extent of expansions and trans-basin dispersal, it is doubtful that these historical dynamics could have been generated without major rearrangements in the paleo-hydrological network, stemming from the formation and melting of large-scale glacial complexes much older than the LGM.
Our knowledge on the evolutionary history of north temperate fishes has been fundamentally altered due to the advent and application of broad-scale phylogeography [1–4]. Phylogeographic investigations of freshwater fishes in Europe are numerous and inferences drawn on the history of intraspecific lineages often relate to how river courses and their accompanying catchment basins dynamically change through several glacial epochs [e.g., [5, 6]]. For cold tolerant fishes such inferences can be complex. Genetic lineages can be distributed mosaically among basins, reflecting repeated population expansions and contractions across the shifting colonization corridors that have resulted from river capture events, the formation and dynamics of pro-glacial lakes and fluctuating levels and salinities of seas [7–9]. Despite relatively sound knowledge of European glaciation and attempts to find common patterns, phylogeographic scenarios are often species specific.
There are few similar studies in Siberia and far less certainty concerning the extent of glaciation and paleohydrological stability . One of the first broad-scale phylogeographic studies in Siberia reported that genetic lineages of grayling (genus Thymallus), corresponded to major Siberian river systems (e.g. Amur, Lena, Enisei) . The study also supported that grayling had been extirpated from Lake Baikal during the early to mid-Pleistocene as the result of some climate-induced environmental perturbation. Subsequently, grayling were able to recolonize Lake Baikal when its waters over spilled forming a new outlet into the Enisei basin, 110,000 to 450,000 years ago . The authors speculated that this event might relate to highly controversial hypotheses concerning the paleo-climate in Siberia. Most geologists consider Siberian glaciation to have been rather limited based on the modeling of sparse precipitation during the Pleistocene (minimum model) . However, field evidence supports extensive glaciation along the polar continental shelves and coastal Pacific lowlands (maximum model) . Such ice sheets would have blocked north flowing rivers and created a series of pro-glacial lakes. Evidence for such blockage has been presented for the Ob and Enisei systems [14, 15].
Furthermore, interior mountain regions (e.g. Trans-Baikalian) were glaciated perhaps above 1000–1200 m. However, many potential refugia for cold tolerant organisms must have existed in central and east Siberia, north of interior mountain systems, as supported by phylogeographic patterns found in grayling from the Lena basin . Siberian glacial scenarios, however, are much in dispute, especially for the last glacial maximum (LGM) . Recent studies reflect an appreciation for the region's paleohydrological dynamics and its effects on organismal history [6, 18–21]. Nonetheless, no study has of yet covered the majority of Siberia where four of the world's ten largest rivers occur (Ob, Lena, Enisei, and Amur).
The Asian endemic salmonid fish Brachymystax lenok occurs in all major Siberian river systems (Figure 1) and thus can serve as a phylogeographic model for assessing paleohydrological events. Lenoks occur in two morphological forms, differing in the length and shape of their snouts as well as a number of external morphological and osteological characters. The so-called sharp- and blunt-snouted lenoks are viewed as either one species complex B. lenok represented by two infraspecies , or as two nominal species. Morphological variation within each form exists among major basins [22, 23] and in zones of sympatry F1 hybrids are found [[24–26], Additional File 1]. Interestingly, longitudinal clines in morphological characters led to the hypothesis of countercurrent dispersal, with sharp-snouted lenok expanding from the west and blunt-snouted from the east . Such dispersal, combined with character displacement in contact zones is thought to have resulted in the formation of clines [22, 23].
Analysis of allozyme variation supported limited gene flow between the two forms [26, 27]. Several enzyme systems can distinguish sharp- and blunt-snouted lenoks within basins, particularly in sympatry, but no alleles are consistently diagnostic throughout Siberia. Moreover, the range-wide pattern of allelic distribution for di-allelic loci is complex with no consistency as to which alleles are fixed or found in high frequencies within a form among all river basins. These results were interpreted as evidence of extensive hybridization and gene exchange in the past, with the present structure of the genus formed via "multiple hybrid speciation" , though such interpretation received criticism .
Considering the interesting evolutionary scenario of sharp- and blunt-snouted lenoks as well as their widespread distribution in Siberia, we proposed a combined genetic and morphological approach to test the existing hypotheses of counter-current dispersal, hybridization and gene flow between the two forms. Concomitantly, we aimed to use lenok as a phylogeographic model for further understanding the paleohydrological dynamics of the Pleistocene in Siberia.
The final alignment included 494 bp of the mtDNA control region (CR) (N = 151), 987 bp of the NADH-1 gene (ND1) (N = 142), and 1481 bp (N = 114) with both genes combined. There was no significant deviation from base frequency homogeneity across taxa for either gene. In lenok the transition/transversion ratio was 2.56 for the CR, which revealed one 2-bp and four 1-bp indels, with outgroup taxa included. For the coding ND1 there were no indels nor amino acid changes. Neither transitions nor transversions were saturated in either gene segment, including third codon positions. There were 33 variable sites for the CR, 30 of which were parsimony informative (excl. indels) and 208 variable sites (109 parsimony informative) for the ND1. Pairwise sequence divergence within B. lenok ranged from 0 to 3.0% (CR) to 0 to 6.2% (ND1). In all analyses two monophyletic groups are identified corresponding to blunt- and sharp-snouted lenoks (Figure 2). Net divergence between groups ranged from 1.7% (CR) to 4.7% (ND1), while within lineage divergences ranged from 0.2% to 0.7%. Net divergences between outgroup taxa and lenok varied with the two genes: 4.3% (CR) or 7.7% (ND1) for H. hucho, 5.4%(CR) or 7.4% (ND1) for H. taimen, and 7.6% (CR) or 13.3% (ND1) for P. perryi.
The phenogram based on external morphological and osteological characters parallels the genetic results revealing two clusters corresponding to sharp- and blunt-snouted lenoks (Figure 3). Similarly, a plot of the first two factors of a Principal Component Analysis (PCA) reveals two clusters clearly representing the two forms (Figure 4). Only the position of sharp-snouted lenok from the Ob basin, intermediate along the first factor (PC1) between blunt-snouted lenok and all other sharp-snouted individuals prevents 100% diagnosis of all individuals to a form. Size effects were assumed to be minimal as there was no correlation between PC1 and size (Additional File 2) and mean total length of individuals analyzed in both forms was nearly identical (sharp-snouted 41.2 cm; blunt-snouted 41.6 cm). Application of a Discriminant Function (DF) produced from a 25% random sample of the data set resulted in 99.9% and 99.6% correct identification for the blunt- and sharp-snouted lenok, respectively. DF values for the misclassified individuals (two Ob basin sharps and 5 Primor'e region blunts) were intermediate to the ranges of values defining each form.
The CR network revealed 19 haplotypes, 10 in blunt- and 9 in sharp-snouted lenok (Figure 5; Additional File 3). Maximum pairwise divergence was considerably less within blunt-snouted (six steps) compared to sharp-snouted (twelve steps) lenok (Figure 5). Both forms revealed at least one haplotype shared between the Lena and Amur drainages supporting dispersal of some kind between these basins. Both forms also revealed private haplotypes at the eastern (Uda) as well as western (Ob) edge of their distributions, supporting isolation of these drainages.
The addition of 15 GenBank sequences (from China and Korea) to the CR network revealed a significant finding. While several haplotypes from mainland China clustered with both our blunt- and sharp-snouted clades, the remaining haplotypes from Korea and the Paleo-Yellow River formed a distinct cluster, a minimum of seven steps divergent from our sharp-snouted clade (Figure 5). Despite the lack of morphological data it is likely that the first two groups of haplotypes represent blunt-snouted and sharp-snouted lenoks, respectively.
Parsimony analysis of the ND1 gene revealed nearly twice the number of haplotypes (N = 35) as the CR, suggesting a considerably higher substitution rate for this coding gene seqment. Two separate networks were revealed, in addition to one highly divergent sharp-snouted haplotype from the Ob basin (Figures 6A and 6B; Additional File 4). Both networks revealed at least one star-like cluster of haplotypes, reflecting demographic expansions of populations primarily from the Amur (blunt-snouted) and the Amur and Lena (sharp-snouted) basins. Each network also revealed one or more highly divergent groups, spanning 12 steps in blunt- and 13 steps in the sharp-snouted lenok.
The blunt-snouted lenok ND1 network (18 haplotypes) contained no shared haplotypes between the Lena and Amur catchments in contrast to analysis with the more slowly evolving CR. Some Amur haplotypes were shared with the Tugur River and/or one of the island populations, whereby the Shantar Islands also revealed a unique haplotype reflecting a degree of long-term isolation. Considerably longer periods of isolation were suggested for both the extreme eastern (Ob) and western (Sea of Japan) regions of the blunt-snouted lenok's distribution. The Ob basin revealed five unique haplotypes that were a minimum of two steps from either the Lena's or Amur's most common haplotype. In the Sea of Japan basin, two very distinct haplotypes (ND15 & ND16) were found, both a minimum of eight steps from any other in the network.
The sharp-snouted ND1 network (17 haplotypes) revealed one presumably ancestral haplotype (ND17) shared among populations throughout the Lena, as well as locations in the Indigirka, Kolyma and Amur basins. This haplotype is central to a star-like cluster of 11 haplotypes, each one or two base pairs divergent from ND17. Unique haplotypes were found in the Amur, Tugur, Enisei, and Uda basins. Three haplotypes from the Selenga-Baikal basin (ND31, 32, 33) form a highly divergent group, suggesting a refuge or isolation, not reflected in data sets from other salmonid fishes . The Ob basin is fixed for one, highly divergent haplotype (ND34) beyond the 95% parsimony limit defining the network, which, together with the results of the blunt-snouted lenok, demonstrate long-term isolation of this basin.
Among basin comparisons (Amova & pairwise differences)
For both lenok forms, the among group variance (ΦCT), representing differences between major ocean basins (Arctic and Pacific) was minimal (or negative) and statistically non-significant, while the largest and highly signficant component (ΦSC) represented the among river basin variance (Additional File 5). A overview this among basin variability can be provided with a table of pairwise differences, which demonstrates relatively large average differences among most pairs of basins, except those that share some haplotypes (Additional File 6), such as the Amur and Lena (both forms) or Lena and Enisey (sharp only). Within basin variation is mostly much lower, except for those basins harboring multiple divergent haplotypes such as the Primor'ye and Ob basins for blunt-snouted lenok, and the Enisei or Tugur/Uda region for sharp-snouted lenok.
The CR pairwise mismatch distribution was uni-modal for blunt-, and bi-modal for sharp-snouted lenok (Additional File 7). Removal of regionally restricted haplotypes (Ob and Uda basins) in sharp-snouted lenok, reflecting population subdivision allowed analysis of uni-modal distributions for both forms, and these distributions both differed significantly from those expected for stable population sizes (K-S tests; P < 0.0001). Additionally, using the least-squares approach, there was no significant difference between observed and simulated data under an expansion model for either form (blunts: SSd, P = 0.108; Harpending's raggedness index, P = 0.116; sharps: SSd, P = 0.1680; Harpending's raggedness index, P = 0.1000).
The ND1 mismatch distribution was multi-modal for both forms (Additional File 8). Similar to the CR analysis, removal of geographically distinct haplotype groups allowed analysis of uni-modal distributions, representing the Amur and Uda basins for blunt-snouted, and the Amur and Lena basins for sharp-snouted lenok. Again, observed mismatch distributions were significantly different from that expected under stable population models (K-S tests; P < 0.0001), but the data did not fit the expansion model either.
Using estimates for the expansion parameter tau, along with divergence rate estimates ranging from 0.5 to 3% per million years for the CR and 1.5 to 6% for the ND1, we estimated the mean age of expansion for both genes in both forms (Figure 7A–D). Mean estimates of expansion times, regardless of the gene or presumed substitution rate clearly relate to periods in the mid- to late Pleistocene (50,000 to 400,000 year ago) but in all cases earlier than the LGM (18,000 years ago).
Microsatellite loci revealed significant deviations from HWE within populations across loci (primarily positive FIS values), most likely due to the sampling of spawning and post-spawning aggregates of individuals. A test of deviation from HWE across all loci and populations was also significant. The population tree reveals two groups (sharp- and blunt-snouted forms) supported by a moderate D AS bootstrap value (73%) and a lower value (56%) for Nei's genetic distance (Figure 8).
The Bayesian population structure analysis for the entire data set resulted in the highest delta log-likelihood between K = 1 and K = 2, representing the two groups of individuals corresponding to sharp- and blunt-snouted lenoks. Strong structure was further revealed within forms (data not shown), especially for blunt-snouted lenok corresponding to different sampling dates, providing the most plausible cause of deviations from HWE. That is, a Wahlund effect is assumed caused by the pooling of distinct population units into single tests for HWE.
We further tested for introgression between forms by running the Bayesian structure analysis on each sympatric sampling site individually, with K set at 2. Using this analysis there was little to no signal for a mixed ancestory for individual genotypes (Additional File 9).
Correlation between morphological and genetic distances
Genetic (mtDNA) and morphological distances among populations were highly correlated in sharp-snouted lenok (Mantel's r = 0.75, P < 0.001) but weakly correlated in blunt-snouted lenok (Mantel's r = 0.36, P = 0.014). Morphological analysis corresponds with genetic data as populations are almost exclusively grouped according to major river basins, except for divergence within the Selenga basin, a phenomenon common to the genetic data set (Figure 3).
Phylogenetic relationships and taxonomic inferences
Extensive mtDNA screening and comprehensive morphological analysis of Brachymystax throughout its Siberian range support that the two forms represent two independent evolutionary lineages, and analysis of 11 microsatellite loci in three sympatric populations (N = 130) revealed no signs of introgression. Thus, the hypothesis of "multiple hybrid speciation"  is not supported with our selection of genetic markers.
The complete concordance between morphology and mtDNA haplotype as well as no evidence of admixture in sympatric populations analyzed with microsatellite loci strongly suggests that introgression between the two forms is low or non-existent. Nonetheless, hybrids exist [[24–26], Additional File 1], some evidence of introgression has been reported  and in captivity, F1 hybrids and their progeny (backcrosses and F2) are viable [Alekseev S., unpublished data]. Thus, one or more mechanisms must be operating to prevent higher levels of introgression, such as selection against hybrids in nature. However, a pre-zygotic barrier, such as spatial segregation of spawning may be the primary isolating mechanism. Despite widespread co-occurrence in the Amur  and Lena basins  blunt-snouted lenok are reported to prefer cooler water during summer months and use distinct spawning grounds. Blunt-snouted lenok spawned both in mountain tributaries (at 3–4°C) and in the main channel (at 5–7°C) of a second-order Amur tributary whereas sharp-snouted lenok spawned in the main channel only, typically in its more downstream reaches .
The lack of correlation between genetic and morphological distances in blunt-snouted lenok may be due to the fact that they exist in smaller and more fragmented populations, which would result in higher genetic drift and more stochastic evolutionary change. Likewise, strong genetic drift among basins, as well as fragmentation within basins could explain the range-wide incongruities in allozyme allele distributions between the two forms [26, 27, 30], and thus more complex models of hybridization and introgression need not be invoked.
The data from China and Korea are difficult to incorporate into our conclusions, although these populations may prove pivotal for understanding the evolutionary history of the genus. The divergent mtDNA lineage in the Yellow River basin for example, is only weakly supported as a sister group to the sharp-snouted clade (data not shown). Lenok from this region are known as B. lenok tsinlingensis , but  and  place them together with Siberian blunt-snouted lenoks based on external morphology.  assigned the name B. lenok to the sharp-snouted form, and B. savinovi to the blunt-snouted, but the latter name was shown to be invalid . Alternatively the name B. tumensis  first used for Tumen River (Korea) lenok is currently used for the blunt-snouted lenok . If the Yellow River and Korean lineage can be definitively defined as blunt-snouted lenok, then the strong correspondence between genetic lineage and form presented in our work would be eroded.
The taxonomy of sharp- and blunt-snouted lenok remains unresolved. The genetic proximity of the Tumen River haplotype to the sharp-snouted lineage (Figure 5) casts serious doubt on the applicability of the name B. tumensis for blunt-snouted lenok. Thus, taxonomic harmonization of Brachymystax can only be achieved when both genetic and morphological (especially craniological) characters of Chinese and Korean lenoks are compared to material from the Siberian range.
Comparative phylogeography of blunt- and sharp-snouted lenok
Both lenok forms revealed strong phylogeographic structure with some common patterns such as the existence of highly divergent mtDNA lineages below 56° N and no such variation above 56° N. This suggests that extensive regions of contemporary lenok habitat in Siberia may have been uninhabitable during one or more Pleistocene glacial maxima. In contrast, high within or among basin diversity in more southern latitudes supports the existence of long-term refugia. Interestingly, there has apparently been a complex history within several of these refugia, analogous to the concept of "refugia-within-refugia" promoted for the high lineage diversity of numerous taxa within the Iberian Peninsula in Europe . This "within-refugia" pattern can be seen for the Ob and Primor'e basin blunt-snouted lenoks as well as the Enisey sharp-snouted lenoks. Only the Ob sharp-snouted lenoks run counter to this pattern, as only a single haplotype has been thus far found.
Furthermore, as the most divergent groups within each form are not geographically congruent, different regions apparently served as refugia for each form during unfavorable paleohydrological conditions. For example, relatively divergent haplotypes in sharp-snouted lenok are found in the Selenga/Baikal portion of the Enisei basin where blunt-snouted lenok are wholly absent. Additionally, in contrast to the broadly distributed sharp-snouted lenok, blunt-snouted lenok above 56° N are only found in the Lena basin, a very large system where at least its upstream regions are known to have served as a refuge for an endemic lineage of grayling Thymallus . These patterns imply that the two lenok forms, perhaps related to their spatial or temperature preferences differ in both dispersal ability and response to climatic change.
Whereas the majority of genetic variance in both forms was distributed among major basins more specific mechanisms generating phylogeographic structure were inferred. Both long distance dispersal (Uda) and contiguous range expansion (Shantar Islands) underscore the importance of the Amur region as a center of radiation for blunt-snouted lenok. Additionally, the morphological phenogram groups blunt-snouted lenoks from both islands together with those from the Amur. That blunt-snouted lenoks from the Tugur carry the most widespread Amur haplotype, suggests that this river was a corridor between the Amur and the Shantar islands. This supports earlier conclusions of  and  and is concordant with geological data indicating that Shantar and Sakhalin islands were connected to the continent during glacial maxima. As both genetic and morphological divergences are greater between the Uda and Amur, than the Amur and islands suggests that the Uda may have been colonized during an earlier marine regression than the islands, which were presumably colonized during the LGM. This perspective was supported by the mismatch analysis, revealing expansion from the Amur basin. The age estimate of this event clearly reaches back into the mid-Pleistocene (Figure 7), based on moderate substitution rates of 1%/MY for the CR, or moderately faster rates for the ND-1. The earliest split within blunt-snouted lenok, however, is clearly between populations of the Primor'e region and all others. Coupled with their morphological distinctiveness, Primor'e region blunts have clearly experienced a long isolation from the Amur basin.
The demographic history of sharp-snouted lenoks is more difficult to describe. Mismatch analysis supports exponential growth for the Lena, Amur and Enisei data combined, and long distance dispersal was documented from the Lena to the Amur drainage, as well as continuous range expansion within the Amur and into the Enisei basin. The allopatric fragmentation in sharp-snouted lenoks involving the Selenga-Baikal lineage may correspond to the prior isolation of the Baikal basin and subsequent connection to the Enisei basin as discussed in relation to Thymallus , or about 500,000 years ago according to geological data . Yet longer periods of isolation are clear for sharp-snouted lenok in the extreme east (Uda) and west (Ob).
In both forms, the CR network reveals haplotypes shared between Lena and Amur basins, as previously described  and thus supporting paleohydrological exchange between these major basins. Interestingly, this basin sharing is seen with the ND1 gene in sharp-snouted but not blunt-snouted lenok perhaps reflecting a longer period of allopatry between the two basins for blunts, allowing sufficient time for mutations to occur.
While data from China and Korea adds some confusion to phylogeographic summary of Siberia, a scenario suggesting that the Chinese portion of the lenok's range was colonized by expansion from Siberia southwards  is unlikely, or at best oversimplified. Chinese lenok populations contain diverse lineages, perhaps reflecting both long-term relictism as well as admixture from southward expansion.
The primary demographic inference was that sudden expansions in both forms relate to paleo-climatic conditions much older than the LGM, and thus are better reflected in the more slowly evolving CR. Moreover, signals of sudden expansions are limited to the Amur, Lena, or Enisei, and were not evident in the extreme eastern or western areas of occurrence. These results appear incongruent with countercurrent dispersal , as no signal of expansion from Western or Central Siberia can be seen for sharp-snouted lenoks, but rather from the aggregated distributions across major basins. Nonetheless, the occurrence of highly divergent lineages in the east for blunts and the west for sharps supports long-term residence, but there is a lack of genetic resolution for the deep past, so no inference concerning origins of either lineage can be made.
Although major basins account for a substantial percentage of both morphological and genetic variation, underscoring isolation, the inferred trans-basin expansions support repeated availability of paleohydrological dispersal (active or passive) corridors. This implies a dynamic paleo-landscape with refugia perhaps in the form of periglacial lakes and shifting drainage patterns stemming from headwater captures and river flow reversals. More specific details concerning Siberia's paleohistory remain highly controversial, and can not be addressed based on the moleulcar phylogeography of Brachymystax alone. Glacial reconstructions range from nearly complete coverage of northern Siberia by a marine-based ice-sheet [13, 42, 43] to individual ice caps centered on arctic archipelagos that advanced onto adjacent shelves . It remains to be seen whether or not phylogeographic models, and especially comparative phylogeography could be used to support or refute one of these contrasting scenarios of paleo-Siberian landscapes. While prior studies [11, 16, 20, 21] have invoked highy dynamic paleo-climatic events to explain the distribution of salmonid mtDNA lineages in various Siberian basins, we offer no further speculations on these scenarios. It is clear, however, that the level of resolution in understanding the effect of paleohydrology on current genetic diversity and distribution of freshwater fishes in Siberia lags behind that of Europe and North America.
To date, there have been no Siberian-wide phylogeographic studies of aquatic organisms, and as well, an underestimation of the effects of paleo-hydrological dynamics on their current distributions. Our data clearly identify a "northern" gradient in lineage diversity particularly evident across 56° N, multiple genetic imprints of inter-basin exchange, long-distance dispersal and sudden expansions across broad expanses of Siberia in two evolutionarily independent lineages of the salmonid fish Brachymystax lenok. Genetic signals for these events clearly date to climatic periods prior to the LGM, and thus were better reflected by the more slowly evolving mtDNA control region rather than the faster mutating NADH-1 gene segment.
Several previous hypotheses concerning the evolution of lenok are not supported by our data. Namely, multiple hybrid speciation , countercurrent dispersal , and a Siberian source for lenok in the Chinese portion of their current range . While no inference can be drawn on the potential origins of the two forms of lenok, both exhibit highly divergent lineages located at opposite end's of the composite distribution range (to the extreme west for the sharp-snouted form, and the extreme east for the blunt-snouted form). A simple taxonomic division of two genetic lineages of lenok, corresponding to two phenotypically recognizable forms is compromised by the existence of a third mtDNA lineage in China and Korea where phenotypic data is lacking. Our data support species recognition of the two forms, but taxonomic harmonization rests on the ability to integrate both genetic and phenotypic data from the China and Korea.
Blunt-snouted (n = 663) and sharp-snouted (n = 1028) lenoks from 78 locations were sampled in 1975–2005. In 22 locations both forms were collected, whereas 23 sites yielded exclusively blunt-and 33 exclusively sharp-snouted lenoks (Table 1, Figure 1). Several of the most closely related taxa were added for comparative purposes to the genetic analyses including one Hucho hucho, one Hucho taimen, and one Sakhalin taimen Parahucho perryi; and for the morphological analysis 101 H. taimen from various basins.
Amplification and sequencing
Whole genomic DNA was extracted using a standard high-salt protocol. Two mtDNA fragments, the control region (CR) and NADH-1 subunit (ND1) were amplified using the polymerase chain reaction (PCR). The complete CR (including segments of flanking tRNA) was amplified in 97 B. lenok and three outgroup individuals using the primers LRBT-25 and LRBT-1195 (8). The remaining CR sequences were taken from previously published research (GenBank accession no. AY230451–AY230472). The ND1 primers BlNDF and BlNDR  were used to amplify 142 B. lenok as well as outgroup individuals. PCR conditions (25 μl reactions) were identical to those described in , and sequences of the left-domain of the CR (486 bp) and complete ND1 (987 bp) were produced on an ABI-3100 genotyper. When necessary, primers designed to amplify shorter segments were used to amplify degraded DNA (see Additional File 10). New mtDNA sequences have been deposited under accession numbers [GenBank: EU395714–EU395769].
Sequence alignment and phylogenetic analysis
All sequences were easily aligned by eye including an alignment incorporating one sequence from South Korea (GenBank accession no. AF125519), and 14 from China . Quantitative assessment was limited to the sequences produced in our laboratory.
Sequences were imported into PAUP*4.0b10  for phylogenetic analyses, pairwise sequence divergence (uncorrected p distances) and the number of transitions and transversions. Saturation in ND1 was assessed by plotting the number of transitions and transversions against uncorrected p distances for each codon position. A chi-square (χ2) test was used to evaluated base composition homogeneity in the ND1 gene for each codon position.
To evaluate relationships among closely related haplotypes, unrooted networks were constructed with a 95% parsimony criterion using TCS ver 1.13 . Between-group variation was calculated using net nucleotide divergence (Da) in MEGA version 2.1 . Haplotype or clade divergence was also calculated using Da distances between groups using the Kimura two-parameter model. Uncorrected p distances were used for divergence estimates between in- and outgroup taxa.
Maximum parsimony (MP), Maximum-Likelihood (ML) and Neighbor-Joining (NJ) were used for phylogenetic reconstruction. Modeltest 3.0  was used to choose models of nucleotide evolution. To estimate the most likely topology for ML and MP methodologies, heuristic searches (10 replicates) started with stepwise addition trees, with each replicate beginning with a random order of sequences. Bootstrap analysis was used to estimate node support with 10000 (NJ and MP) or 1000 replicates (ML). Full heuristic search algorithms were applied for the MP and the "fast" stepwise addition method for the ML analysis.
Genetic variation among and within major basins was evaluated with an analysis of molecular variance (AMOVA) using Arlequin 3.1 . The AMOVA structure was defined by major ocean basins (Arctic/Pacific), then by large river basins or regions (Ob, Enisei, Lena including Indigirka and Kolyma, Amur, Islands, Primor'e, and Uda and Tugur) (see Table 1). A more detailed overview of within and among basin differentiation is provided with a table of average (and corrected) pairwise differences also done with Arlequin.
The demographic signature of mtDNA haplotype variation was evaluated with the pairwise mismatch distribution . The goodness-of-fit of the observed data to a simulated model of expansion was tested with the sum of squared deviations and Harpending's raggedness index . The age of expansion was estimated with the formula τ = 2 μt, where τ is drawn from the mismatch distribution, μ equals the aggregate substitution rate across all nucleotides per generation (5 years for lenok) and t is the expansion time in generations, graphically displayed in years. The aggregate substituion rate was based on a plausible range of divergence rates for mtDNA in salmonid fishes [see Discussions in [8, 11, 20, 53] and references therin]. Thus, estimates of the age of expansion were plotted across a range of substitution rates, whereby it is assumed that the substitution rate of the coding ND1 gene is higher than the control region, which is common for salmonid fishes [45, 54], and supported by our haplotype diversities (i.e., twice the number of haplotypes for the ND-1 gene compared to the CR). A Kolmogorov-Smirnov two-sample test was used to test the distribution of observed values against those expected under the null hypothesis of a stable population using SPSS ver. 12.0.
To assess introgression between sharp- and blunt-snouted lenoks, we applied bi-parentally inherited microsatellites to three sets of sympatric populations from the Amur basin: the Anui (n = 56), Khor (n = 48), and Sukpai (n = 26) rivers.
Four tri-nucleotide and three tetra-nucleotide microsatellite loci , and four unpublished di-nucleotide microsatellite loci were analyzed (see Supplementary Materials). All forward primers were fluorescently labeled and PCR and genotyping were performed as described in , (see Additional File 10). Exact probability tests for deviations from Hardy-Weinberg equilibrium (HWE) across populations (within loci) and loci (within populations), exact tests for deviations from genotypic linkage equilibrium (LE) across populations, and tests for genic differentiation among populations were performed with Genepop 3.2a . Corrections for multiple significance tests were performed using a sequential Bonferroni-type correction .
To estimate the proportion of each individual's genome originating in each parental species and the patterns of intraspecific variability among sharp- and blunt snouted lenoks we used the Structure . This software implements a Bayesian model-based clustering algorithm that reveals population structure in a data set by resolving clusters of individuals that minimize Hardy-Weinberg and linkage disequilibrium. The parameter settings included the assumption of admixture and a correlated allele frequencies model. In exploratory runs we did not provide the software with prior information regarding the origin of individuals, but in final runs, and when analyzing the three sympatric populations separately, individuals from the two forms were assumed to represent pure blunt- or sharp-snouted lenoks, and used as a proxy for determining the degree of admixture of the individuals within each sympatric population. Structure was run for 100,000 steps, of which the first 10,000 where discarded as burn-in and we conducted five independent replicates of the MCMC for each value of k. We also performed analyses using the independent allele frequencies model to test for robustness of our conclusions to the violation of prior assumptions because of recent suggestions that the choice of the model might strongly influence the outcome of the clustering algorithm [59, 60].
External morphology was analyzed using fresh fish or heads fixed with salt. Measurements were taken with dividers to the nearest mm. We used a modified morphometric scheme  for salmonids [see ], with reductions. Osteological characters were assessed according to  (for details see ). To adequately reflect shape variation, body measurements were expressed as percent of fork length (FL), head measurements as percent of head length (c), skull measurements as percent of skull base length (Lcr) and bone widths (depths) as percent of bone lengths. Allometric growth is known to be weak in lenok above 20 cm FL. Therefore, to reduce potential allometric effects only fish larger than 25 cm FL were used. A total of 38 indexes and 8 counts were analyzed (see Additional File 10).
To explore the relationship in morphology between forms and among populations within forms both principal component analysis (PCA) and cluster analysis were done using the NTSYS-pc package, v2.0 . For PCA, indexes and counts of individual fish were used. Eigenvectors were calculated from the variance-covariance matrix and the eigenvector loadings were scaled so that the length of the vectors equaled the square root of their eigenvalues. For cluster analysis, index and count means were standardized and used to calculate taxonomic distances between populations. These distances served as raw data for construction of a UPGMA phenogram. To further control for potential allometric effects, the first factor in the PCA was regressed with body size.
To quantify the morphological differentiation between forms we used a canonical discriminant analysis (CDA), using the first three PCA factors as input variables, and an equal probability of assignment of each individual to a form (blunt or sharp). A discriminant function (DF) was derived from a 25% random sample of the data set, and this DF was then applied to assign individuals to a form.
Comparison of morphological and genetic data
We assessed the correlation between genetic and morphological data across the range of the two lenok forms using Mantel tests, done in the NTSYS-pc package. To produce a pairwise genetic (mtDNA) distance matrix, between group variation (corrected for within-group variation) was calculated using the net nucleotide distances (DA) (Kimura two-parameter model). For the morphological distance matrix, the same taxonomic distances generated for the clustering were used. As not all populations had matching genetic and morphological matrices, only 22 populations of sharp-snouted, and 17 populations of blunt-snouted lenoks were used. The significance of Mantel's r was evaluated with 9999 permutations.
Nesbo CL, Fossheim T, Vollestad A, Jakobsen KS: Genetic divergence and phylogeographic relationships among European perch (Perca fluviatilis) populations reflect glacial refugia and postglacial colonization. Mol Ecol. 1999, 8: 1387-1404.
Bernatchez L: The evolutionary history of the brown trout (Salmo trutta L.) inferred from phylogeographic, nested clade, and mistmatch analysis of mitochondrial DNA variation. Evolution Int J Org Evolution. 2001, 55: 351-379.
Volkaert FAM, Hänfling B, Hellman B, Carvalho GR: Timing of the population dynamics of bullhead Cottus gobio (Teleostei: Cottidae) during the Pleistocene. J Evol Biol. 2002, 15: 930-944.
Van Houdt JK, Hellemans B, Volckaert FAM: Phylogenetic relationships among Palearctic and Nearctic burbot (Lota lota): Pleistocene extinctions and recolonization. Mol Phylogenet Evol. 2003, 29: 599-612.
Englebrecht CC, Freyhof J, Nolte A, Rassman K, Schliewen U, Tautz D: Phylogeography of the bullhead Cottus gobio (Pisces: Teleostei: Cottidae) suggests a pre-Pleistocene origin of the major central European populations. Mol Ecol. 2000, 9: 709-722.
Brunner PC, Douglas MR, Osinov A, Wilson CC, Bernatchez L: Holarctic phylogeography of Arctic charr (Salvelinus alpinus L.) inferred from mitochondrial DNA sequences. Evolution. 2001, 55 (3): 573-586.
Kontula T, Väinöla R: Postglacial colonization of Northern Europe by distinct phylogeographic lineages of the bullhead, Cottus gobio. Mol Ecol. 2001, 10: 1983-2002.
Weiss S, Persat H, Eppe R, Schlötterer C, Ublein F: Complex patterns of colonization and refugia revealed for European grayling Thymallus thymallus, based on complete sequencing of the mtDNA control region. Mol Ecol. 2002, 11: 1393-1407.
Šlechtova V, Bohlen J, Freyhof J, Persat H, Delmastro G: The Alps as barrier to dispersal in cold-adapted freshwater fishes? Phylogeographic history and taxonomic status of the bullhead in the Adriatic freshwater drainage. Mol Phylogenet Evol. 2004, 33: 225-239.
Dawson AG: Ice Age earth: Late Quaternary Geology and Climate. 1992, Rutledge
Koskinen MT, Knizhin I, Primmer CR, Schlötterer C, Weiss S: Mitochondrial and nuclear DNA phylogeography of Thymallus spp. (grayling) provides evidence of ice-age mediated environmental perturbations in the world's oldest body of freshwater, Lake Baikal. Mol Ecol. 2002, 11: 2599-2611.
Rozenbaum GE, Shpolyanskii NA: Late Cenozoic history of the cryolithozone of the Arctic and the tendencies of its further development. 2000, Moscow: Nauchnyi Mir publishing house, 104-[In Russian].
Grosswald MG: New Approach to the Ice Age Paleohydrology of Northern Eurasia. Paleohydrology and Environmental Change. Edited by: Benito G, Baker VR, Gregory KJ. 1998, Wiley and Sons, Chichester, England, 199-214.
Yamskikh AF: Late Pleistocene and Holocene Siberian River Valley Geomorphogenesis as a Result of Palaeogeographical cyclic changes. Paleohydrology and Environmental Change. Edited by: Benito G, Baker VR, Gregory KJ. 1998, Wiley and Sons, Chichester, England, 111-124.
Mangerud J, Astakhov V, Jakobsson M, Svendsen JI: Huge ice-age lakes in Russia. J Quat Sci. 2001, 16 (8): 773-777.
Weiss S, Knizhin I, Kirillov A, Froufe E: Phenotypic and genetic differentiation of two major phylogeographic lineages of arctic grayling Thymallus arcticus in the Lena River, and surrounding Arctic draining basins. Biol J Lin Sci. 2006, 88: 511-525.
Brigham-Grette J: New perspectives on Beringian Quaternary paleogeography, stratigraphy, and glacial history. Quat Sci Review. 2001, 20: 15-24.
Politov DV, Bickham JW, Patton JC: Molecular Phylogeography of Palearctic and Nearctic Ciscoes. Annales Zoologici Fennici. 2004, 41: 13-23.
Radchenko OA: Variability of Mitochondrial DNA in Populations of Lake Chars from Genus Salvelinus in the Far East and Siberia. J Ichth. 2003, 43 (7): 539-547.
Froufe E, Alekseyev S, Knizhin I, Alexandrino P, Weiss S: Comparative phylogeography of salmonid fishes (Salmonidae) reveals late to post-Pleistocene exchange between three now-disjunct rivers basins in Siberia. Diver and Distrib. 2003, 9: 269-282.
Froufe E, Knizhin I, Koskinen M, Primmer CR, Weiss S: Identification of reproductively isolation lineages of Amur grayling (Thymallus grubii Dybowski 1869): concordance between phenotypic and genotypic variation. Mol Ecol. 2003, 12: 2345-2355.
Mina MV: Microevolution of fishes: evolutionary aspects of phenetic diversity. 1991, Oxonian Press Pvt. Ltd., New Delhi, Calcutta, 215-
Alekseyev SS, Mina MV, Kondrashov AS: Parallel clines as the result of countercurrent dispersion and character displacement with the special reference to the genus Brachymystax (Salmoniformes, Salmonidae). Zoologicheskii Zhurnal. 1986, 65 (2): 227-234. [In Russian].
Alekseyev SS: Morpho-ecological characteristics of lenoks (Brachymystax, Salmonidae, Salmoniformes) from the Amur River basin and from the Uda River. Zoologicheskii Zhurnal. 1983, 62 (7): 1057-1068. [In Russian].
Alekseyev SS, Kirillov AF, Samusenok VP: Distribution and morphology of the sharp-snouted and the blunt-snouted lenoks of the genus Brachymystax (Salmonidae) of East Siberia. J Ichth. 2003, 43 (3): 311-333.
Osinov AG, Il'in II, Alekseyev S: Forms of lenok, Brachymystax (Salmoniformes, Salmonidae) delineated by genetic analysis. J Ichth. 1990, 30 (5): 138-153.
Osinov AG: Countercurrent Dispersal, Secondary Contacts, and Speciation in Lenoks of the Genus Brachymystax (Salmonidae, Salmoniformes). Genetika. 1993, 29 (4): 654-669.
Shed'ko SV, Ginatulina LK, Parpura IZ, Ermolenko AV: Evolutionary and taxonomic relationships among Far-Eastern salmonid fishes inferred from mitochondrial DNA divergence. J Fish Biol. 1996, 49: 815-829.
Besednov LN, Kucherov AI: On the Systematic Position of the Genus Brachymystax from the Iman River. Zoologicheskie problemy Sibiri (Zoological Problems of Siberia). Edited by: Novosibirsk Nauka. 1972, 220-221. [In Russian].
Alekseyev SS, Osinov AG: Blunt-snouted lenoks (genus Brachymystax: Salmoniformes, Salmonidae) from the Ob' basin: new data on morphology and allozyme variation. J Ichth. 2006, 46 (7): 500-516.
Li SZ: Discussion on the geographical distribution of the salmonid fishes in Chin. Annual Bulletin of the Freshwater Fish Protection Association. Freshwater Fisheries. 1985, 11: 89-93. [In Chinese].
Kifa MI: Morphology of two forms of lenok (genus Brachymystax, fam. Salmonidae) from the Amur basin and their systematic position. Zoogeografiya i sistematika ryb (Zoogeography and systematics of fishes). Edited by: Skarlato. 1976, Zoological Institute of the Academy of Sciences of the USSR Publishing House, Leningrad [In Russian], 142-156.
Shed'ko SV: List of cyclostomes and fishes of coastal fresh water of Primor'e. Vladimir Ya. Levanidov's biennial Memorial meetings. 2001, Vladivostok, Dal'nauka, 1: 229-249. [In Russian].
Vasil'eva ED, Mina MV: Comparative Analysis of Morphological Characters of Lenoks from Different Regions of the Range of the Genus Brachymystax (Salmoniformes, Salmonidae). Zool Zhur. 1980, 59 (1): 79-90. [In Russian].
Mori T: On the fresh water fishes from the Tumen River, Korea, with descriptions of new species. J Chosen Nat Hist Soc. 1930, 11: 39-49. pl. 3.
Shedko SV, Shedko MB: New data on freshwater ichthyofauna of the south of the Russian Far East. Vladimir Ya. Levanidov's biennial memorial meetings. 2003, Vladivostok, Dal'nauka, 319-336. [In Russian]., 2
Gómez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of Southern European Refugia. Edited by: Weiss S, Ferrand N. 2007, Springer, Dordrecht, Netherlands, 155-188.
Alekseyev SS, Dudnik YuI: The lenok Brachymystax lenok from the rivers of Sakhalin Island and its phenetic relationship to lenok from waters of the Far Eastern part of the continent. J Ichth. 1989, 29 (5): 145-147.
Alekseyev SS, Gruzdeva MA, Skopets MB: Fish fauna of the Shantar islands. J Ichth. 2004, 44 (1): 36-51.
Mats VD, Ufimtsev GF, Mandelbaum MM, Alakshin AM, Pospeev AV, Shimarev MN, Khlystov OM: Cenozoic era of the Baikal rift depression. Structure and geological history. 2001, Novosibirsk:SO RAN (Geo) Publishers, [In Russian].
Xia Y, Sheng Y, Chen Y: DNA sequence variation in the mitochondrial control region of lenok (Brachymystax lenok) populations in China. Bio Sci. 2006, 14 (1): 48-54.
Lambeck K: Constraints on the Late Weichselian ice sheet over the Barents Sea from observations of raised shorelines. Quat Sci Rev. 1995, 14: 1-16.
Peltier WR: Mantle viscosity and ice-age ice sheet topography. Science. 1996, 273: 1359-1364.
Velichko AA, Kononov YM, Faustova MA: The last glaciation of Earth: Size and volume of ice sheets. Quat Inter. 1997, 41–42: 43-52.
Froufe E, Alekseyev S, Knizhin I, Weiss S: Comparative mtDNA sequence (control region, ATPase 6 and NADH-1) divergence in Hucho taimen (Pallas, 1773) across four Siberian river basins. J Fish Biol. 2005, 67: 1040-1053.
Swofford DL: PAUP* ver 4.0.b10. Phylogenetic analysis using parsimony and other methods. 2002, Sinauer Associates, Sunderland, Sunderland, MA
Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9: 1657-1660.
Kumar S, Tamura K, Jakobsen IB, Nei M: MEGA2: molecular evolutionary genetics analysis software. Bioinf. 2001, 17: 1244-1245.
Posada D, Crandall KA: Modeltest: Testing the model of DNA substitution. Bioinf. 1998, 14: 817-818.
Schneider S, Roessli D, Excoffier L: Arlequin ver. 3.1: A software for population genetic data analysis. 2000, Genetics and Biometry Laboratory, University of Geneva. Geneva, Switzerland
Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9: 552-569.
Harpending H: Signature of ancient population growth in a low resolution mitochondrial DNA mismatch distribution. Human Biol. 1994, 66: 591-600.
Smith GR: Introgression in fishes: significance for paleontology, cladistics, and evolutionary rates. Sys Bio. 1992, 57: 41-57.
Froufe E, Knizhin I, Weiss S: Phylogenetic analysis of the genus Thymallus (grayling) based on mtDNA control region and ATPase 6 genes, with inferences on control region constraints and broad-scale Eurasian phylogeography. Mol Phylogenet Evol. 2005, 34: 106-117.
Froufe E, Sefc KM, Alexandrino P, Weiss S: Isolation and characterization of Brachymystax lenok microsatellite loci and cross-species amplification in Hucho spp. and Parahucho perryi. Mol Ecol Not. 2004, 4 (2): 150-152.
Raymond M, Rousset F: Genepop (version 1.2): population genetics software for exact tests and ecumenicism. J Heredity. 1995, 86: 248-249.
Rice WR: Analyzing tables for statistical tests. Evolution. 1989, 43: 223-225.
Pritchard JK, Stephens M, Donnelly P: Inference of population structure from multilocus genotype data. Genetics. 2000, 155: 945-959.
Serre D, Päabo S: Evidence for gradients of human genetic diversity within and among continents. Gen Res. 2004, 14: 1679-1685.
Rosenberg NA, Mahajan S, RAmachandran S, Chengfeng Z, Pritchard JK, Feldman MW: Clines, clusters and the effect of study design on the inference of human population structure. PLoS Genet. 2005, 1: e70-
Pravdin IF: Rukovodstvo Po Izucheniyu Rib. 1966, Pischevaya Promyshlennost Press, Moscow, [in Russian].
Alekseyev SS, Samusenok VP, Matveev AN, Pichugin My: Diversification, sympatric speciation, and trophic polymorphism of Arctic charr (Salvelinus alpinus complex) in Transbaikalia. Environ Biol Fish. 2002, 64 (1–3): 97-114.
Savvaitova KA, Maksimov VA, Medvedeva ED: Davatchan Salvelinus alpinus erythrinus (Georgi). Vopr Ikhtiol. 1977, 203-219.
Rohlf FJ: NTSYS-pc: Numerical taxonomy and multivariate analysis system, v. 2.0. Exeter Software. 1998, Setauket, New York
Financial support for this study was provided by the Portuguese Ministry of Science and Technology, "Fundação para a Ciência e Tecnologia" grant POCTI/BSE/33364/99 to P. Alexandrino and S. Weiss, as well as a PhD grant to E. Froufe (SFRH/BD/11377/2002), grants from the Russian Foundations for Basic Research (project 06-04-48352), and from the Program of Fundamental Studies of the Presidium of the Russian Academy of Sciences "Biodiversity and the Dynamics of Gene Pools". We are extremely grateful for the assistance of numerous colleagues in providing samples and/or field assistance (see Additional File 12). We thank M.V. Mina for valuable comments on the manuscript.
EF, SA, PA, and SW participated in the design and coordination of the project. EF and SW carried out the molecular genetic studies whereas SA carried out the morphological work and collected most of the samples. SW wrote the manuscript and all authors read, revised and approved the final version.