In the rivers: Multiple adaptive radiations of cyprinid fishes (Labeobarbus) in Ethiopian Highlands

Multiple repeated patterns of adaptive radiation were revealed in cyprinid fish inhabiting the compact geographic region of the Ethiopian Highlands. We found four independently evolved radiations in the evolutionary hexaploid (2n = 150) Labeobarbus lineage based on matrilineal relationships of >800 individuals. Each radiation displayed similar patterns of mouth phenotype diversification, and included ecomorphs/species of the generalized, lipped, scraping (one or two), and large-mouthed (one to three) types. All radiations were detected in geographically isolated rivers, and originated from different ancestral populations. This is the first documented case in which numerous parallel radiations of fishes occurred–via different ways–in a riverine environment. Some radiations are very recent and monophyletic, while others are older and include ecomorphs that originated in separate mini flocks and later combined into one. The diversification bursts among Ethiopian Labeobarbus were detected in the mid-upper reaches of rivers (1050-1550 m above sea level), which likely offer ecological opportunities that include diverse habitats yet poor fauna (i.e. lower competition and relaxed selection). This promising example of parallel evolution of adaptive radiation warrants further investigation.

The origin of biodiversity is one of the most intriguing questions in evolutionary biology and ecology. Along with the neutral divergence of geographically isolated lineages, adaptive radiation-the emergence of ecological and phenotypic diversity in rapidly diversifying lineages 1 -is considered to be the main mode that biodiversity is generated. Adaptive radiations usually result in biodiversity bursts in geographically restricted areas 2,3 . They often accompany the colonization of novel environments such as islands or lakes, or follow a massive species extinction when new ecological opportunities or niches become available [4][5][6] .

for reviews.
Although examples of adaptive radiation in riverine fish are not particularly frequent, they have substantially increased during the last decades. Such evolutionary phenomena have been revealed by molecular data in the African mormyrids, mochokid catfishes and cyprinids [16][17][18][19][20] , as well as in the African and South American cichlids [21][22][23][24][25][26][27] . For example, cichlids of the genus Crenicichla Heckel 1840 from Paraná and Uruguay Rivers demonstrate bright morpho-ecological diversification and trophic resource partitioning. In particular, five co-occurring ecomorphs diverged in trophic ecology, and are of recent origin based on genetic data 27 . The cyprinids of the Labeobarbus Rüppell 1835 from the Genale River, East Africa, constitute a sympatric assemblage of six ecomorphs, five of which are divergent in trophic ecology. All riverine forms of the Labeobarbus from the Genale River have intra-basin origins 20 . Apart these cases, there is evidence of several other putative riverine adaptive radiations in the African and Asian cyprinids [28][29][30][31][32] , though they have yet to be confirmed with genetic data.

Results
phylogeny of Labeobarbus in the ethiopian Highlands. Both Bayesian inference (BI) and maximum likelihood (ML) analyses supported the division of Ethiopian Labeobarbus into two clades representing the eastern (L. gananensis/L. jubae complex) and western (L. intermedius supercomplex) parts of the Ethiopian Highlands, with L. ethiopicus (Zolezzi 1939) as an outgroup (Fig. 3). The western clade also includes L. altianalis (Boulenger 1900) from the Lake Victoria basin in Kenya, which is a sister lineage to all Labeobarbus in this clade. Ethiopian Labeobarbus in the western clade are further subdivided into A) northern, and B) southern lineages. Lineage A includes populations from the Blue Nile basin with Lake Tana, the Atbara-Tekeze basin (the Nile), the mainly northern localities from the Ethiopian Rift Valley (the Awash system, Lakes Langano, Awassa and several representatives from Lakes Abaya and Chamo). Surprisingly, this lineage also includes one population from the White Nile basin (Sore River). Two riverine (the Didessa and Sore) and one lacustrine (Lake Tana) radiations were nested within northern lineage. The southern lineage is comprised of populations from the White Nile basin, with the exception of the Sore population, and populations from the Omo-Turkana and the southern part of the Ethiopian Rift Valley (Lakes Awassa, Abaya, Chamo, and Chew Bahir basin) (Fig. 3). This group includes Labeobarbus' diversification detected in the Gojeb River, a tributary of the Gibe River (Omo-Turkana system). These northern and southern lineages correspond to lineages A and B suggested by Beshera & Harris 72 , but include more basins and localities due to our wider geographic coverage. The remaining L. ethiopicus and L. beso (Rüppell 1835) that are endemic to Ethiopia but placed in an outgroup position to all other Ethiopian Labeobarbus.
Although the intragroup divergence was often weakly resolved or not resolved, the monophyly was supported by the Bayes factor test for the Gojeb (log BF = 7.82) and Sore (log BF = 8.10) assemblages, confirming their species flock status. However, there was no such support for the Didessa (log BF = 51.33 in favor of unconstrained topology) and Genale (log BF = 394.6) assemblages.
Cytb haplotype network/phylogeography. Among the 769 individuals of Ethiopian Labeobarbus sampled, 185 haplotypes were detected. The haplotype network is complex (Fig. 4), composed of five main haplogroups with a central haplotype represented by an individual from the Ethiopian Rift Valley (Fig. 4). Haplogroup (1) includes the most diversified labeobarbs of the L. gananensis/L. jubae complex from the Genale River and other tributaries of the ancient Juba-Wabe-Shebelle (JWS) drainage, with 24 mutational steps to the nearest haplotype from the Chamo Lake in the southern Ethiopian Rift Valley. Haplogroup (2) is also diverse and includes the lacustrine radiation of Labeobarbus from Lake Tana, along with populations from the Blue Nile (including the Didessa Labeobarbus, represented by two divergent compact haplogroups), the Nile (Tekeze basin), and populations mainly from the northern part of the Ethiopian Rift Valley. Haplogroup (3) comprises the riverine radiation of the Gojeb barbs within the Gibe River, as well as populations from the southern part of the Ethiopian Rift Valley including Lake Turkana. The portion of the Didessa radiation that is mainly composed of forms with the large-mouthed phenotype is very close to the Omo-Turkana haplogroup. Haplogroups (4) and (5) are much more compact and both represent the tributaries of the White Nile basin. While haplogroup (4) includes most samples from the Baro River and its tributaries, haplogroup (5) includes only the population from the Sore River, a radiation of Labeobarbus reported in the current study. These haplogroups were separated by 23 mutation steps and joined via a central haplotype.
For ease of visualization, Fig. 5 highlights only the riverine adaptive radiations. Based on both the cytb phylogenetic trees and haplotype network, all four radiations appear to have evolved independently. Each originated within its riverine basin, except for the Didessa radiation (see below): (I). The Genale radiation occurred within the Juba and Wabe-Shebelle basin (Indian Ocean drainage); it displays the deepest divergence and highest haplotype diversity (Table 2), which is consistent with previous results from a recent study 20 . (II). The Gojeb radiation occurred in the Gojeb-Gibe (Omo-Turkana basin). Some Gojeb individuals share haplotypes with or are genetically similar to a population from the southern part of the Ethiopian Rift Valley. (III). The Sore radiation was realized within the smallest haplogroup specific for the Sore River. Although the Sore River belongs to the Baro riverine system (White Nile), its haplotypes are distant to any others within the same basin (17 mutations to the closest haplotype). Here the Sore barbs can be considered as its own evolutionary lineage within White Nile basin. It is likely that this lineage is old (it has early divergence in the phylogenetic tree - Fig. 3), but the radiation is recent, as it has the lowest genetic and morphological diversity (unpubl. data) compared to the other three riverine radiations. (IV). The origin of the Didessa radiation is more complicated than that of the other three radiations, which occurred within their respective basins. The scraper form L. beso, which is phenotypically analogous to the scraping forms in other radiations, is distant to any of the other Ethiopian barbs and branched significantly earlier (in outgroup position in Fig. 3). Most of the haplotypes of the three large-mouthed forms belong to the Omo-Turkana haplogroup, although they were separated from the radiation detected in the Omo-Turkana by 6 mutational steps. Only few individuals shared haplotypes with the co-occurring generalized form from the Didessa, most likely due to introgression, considering the large genetic distance between these two (Figs. 3-4, S3). The group of the Didessa large-mouthed forms within Omo-Turkana haplogroup does not share haplotypes with any other forms, and was also distant to other sympatric radiations like Gojeb and Sore. This most likely represents a remnant of the past flocks like the scraper L. beso. Therefore, the Didessa radiation is not of sympatric origin (except for the generalised-lipped phenotypes; see below for discussion). The generalised, scraper and large-mouthed forms originated from different ancestors. Interestingly, despite this complex history, the forms in the Didessa are very similar to the www.nature.com/scientificreports www.nature.com/scientificreports/ forms from other basins in which sympatric origin was detected. This raises the question as to whether the riverine environment strongly and predictably shapes the adaptive landscape of the Labeobarbus? Genetic diversity within adaptive radiations and intra-radiation relationships. The highest haplotype diversity was detected in the Genale (h = 0.91) and Didessa (h = 0.80) radiations, which are non-monophyletic in origin (Levin et al. 20 and this study), compared to the Sore (h = 0.65) and Gojeb (h = 0.52), which had much lower nucleotide diversity and average number of nucleotide substitutions ( Table 2).
The genetic differentiation between sympatric forms was revealed within the Genale and Didessa riverine radiations in F ST values (Table S4). No subdivision or incomplete lineage sorting between sympatric phenotypes was detected within the Gojeb and Sore riverine radiations. However, some phenotypes in these rivers had haplotypes that were divergent from sympatrically occurring phenotypes. In particular, the large-mouthed phenotype differed from the generalized, lipped and scraper1 in the Gojeb River, the while in the Sore River the scraper phenotype differed from the generalized form (Table S4).

Discussion
The most important result of our study is the evidence for repeatedly evolved similar patterns of ecological diversification (namely in mouth phenotypes) in a compact geographic region. We found that cyprinid fish Labeobarbus exhibit parallel adaptive radiation in the four rivers in the Ethiopian Highlands that belong to independent riverine drainages. Although this phenomenon is commonly observed in lake-dwelling fishes (e.g. 1,38,42,43,46,76,77 ), it is extremely rare in river-dwelling fish assemblages. Thus far, only one case of repeated patterns composed of two riverine adaptive radiations has been documented within the South-American cichlids of the genus Crenicichla from the Paraná and Uruguay Rivers 25,27 . The four riverine radiations of Ethiopian Labeobarbus originated from different ancestral populations. Moreover, our mitochondrial data suggests that the repeated patterns of phenotypic diversification also have different origins, which we discuss below. Two of the four assemblages-the Gojeb (Omo-Turkana basin) and Sore (White Nile)-are monophyletic in origin without complete haplotype sorting between sympatric forms. This finding can most likely be explained by their recent origin. Nevertheless, some sympatric phenotypes in these radiations differed significantly in haplotype frequency, which may indicate reproductive isolation (Figures S1-S3; Table S4). These two can be considered as species flocks (sensu 78 ) and their sympatric origin is a plausible hypothesis. Both the Gojeb and Sore species flocks have significantly lower haplotype and nucleotide diversity compared to the Genale and Didessa assemblages. The Sore radiation differs from the Gojeb in that it has fewer forms (4 vs. 5). Moreover, the low F ST values suggest a younger origin or slower diversification rate of the Sore radiation compared to the Gojeb, which is supported by our morphological data (unpublished) that also revealed shallower phenotypic divergence among the sympatric Sore forms. Notably, the Sore lineage showed early branching in the phylogenetic tree, which would suggest that this relatively old lineage has only recently started to diversify.
The origin of the other two radiations (Genale and Didessa) is more complex. Our data indicate that their assemblages are the result of a combination of sympatric and allopatric speciation with secondary contact due to geological events during ancient periods. A recent study highlighted that the diversification of the Labeobarbus from the Genale River is a result of three mini species flocks that originated in different channels of an ancient riverine net with a dynamic geological history 20 . Two of the specialized Genale forms, large-mouthed and one of the scrapers, originated in two geographically distant parts of the riverine basin and later colonized/migrated to the Genale River where they combined with local forms (generalised, lipped, and short form) 20 . Therefore, the Genale adaptive radiation is of mixed origin fueled by both sympatric and allopatric speciation (see Levin et al. 20    www.nature.com/scientificreports www.nature.com/scientificreports/ for details). For such a diversification that results from the combination of products from different mini flocks of closely related lineages, we have suggested a term multiflock 20 .
Similar to the Genale, diversity of the Didessa Labeobarbus assemblage is a result of the combination of sympatric and allopatric speciation (generalised and lipped forms) during different geological time. The scraper form of the Didessa River, L. beso, is a morphologically highly specialized algae scraper 67 with more ancient origin than any Labeobarbus in the Ethiopian Highlands based on its phylogenetic position 66,68 . It is possible that L. beso is a relic remnant of a former (extinct) species flock that existed before the recent lineages of the Labeobarbus colonized the Ethiopian Highlands. This is plausible in light of the fact that the Ethiopian Highlands is a tectonically active territory with a recent Plio-Pleistocene volcanism [79][80][81] . An ancient species flock hypothesis was recently proposed to explain the origin of diversity within another lineage of algae scrapers of the genus Capoeta Valenciennes 1842 in the Armenian Highland 82 .
The Labeobarbus earliest fossil records were found in the Ethiopian Rift Valley and dated back to the late-Miocene 54 . The Ethiopian Highlands are a volcanic massif of flood and shield volcano basalts 0.5-3.0 km thick that form spectacular trap topography (1500-4500 m) flanking the Main Ethiopian Rift 81 . The geological history of the Ethiopian Highlands was very dynamic and rich in volcanic episodes; four main periods have been detected from Oligocene to Pleistocene time 81 . The volcanic activity has been severe enough to deleteriously affect the biota and cause major disruptions in ecosystems. This could likely explain i) the rare representatives of older branches among Ethiopian Labeobarbus (L. beso and L. ethiopicus), and ii) the young, Pleistocene origin of the majority of Labeobarbus species and populations 53 .
We uncovered three large-mouthed phenotypes in the Didessa River with no intra-divergence in mtDNA; these are genetically distant from sympatric generalised and lipped phenotypes, as well as any individuals from the Blue Nile-Lake Tana system. In a haplotype network, the Didessa large-mouthed forms had an intermediate position between populations from the Gojeb (Omo-Turkana) and the central haplotype (Ethiopian Rift Valley). We hypothesize that the large-mouthed phenotypes have diverged within an ancient flock that is undetected now. Accordingly, co-existence of sympatric ecomorphs in the Didessa would be the result of both sympatric speciation and ancient secondary contacts. Of course, diversification of the three closely-related large-mouthed (piscivorous) phenotypes in the riverine environment is an intriguing phenomenon for cyprinids and should be studied further. We cannot exclude the scenario that these phenotypes (and possibly others as well) diversified upon secondary contact and as a result of genetic admixture. This scenario was successfully tested to explain the origin of functional novelties in some cichlids 83 , and also the contemporary ecological speciation of threespined stickleback in Lake Constance 84 . The ancient hybridization between divergent lineages of cichlids was also reasonably hypothesized to fuel mega-radiations in African Great Lakes [85][86][87] .
To summarize, the very similar patterns of adaptive radiations in riverine drainages of the Ethiopian Highlands were achieved via different processes: • Sympatric ecological speciation, when all members of the adaptive radiation are of monophyletic origin, i.e. represent species flock and sympatric speciation (Gojeb and Sore). • Secondary contact of closely-related but reproductively isolated specialized phenotypes, which originated in different parts of the same riverine basin. This is essentially established for the Genale River multiflock 20 . • Secondary contact of closely-related but genetically isolated pools of phenotypes, which could initiate ecological speciation upon hybridisation (this scenario should be tested further using more comprehensive genomic approaches). Indeed, a similar scenario was recently suggested for Labeobarbus from the Lower Congo 69 .
Parallel origin of the very similar mouth polymorphism of Labeobarbus spp. suggests that these phenotypes are determined by several factors. For example, such propensity to produce different mouth phenotypes is apparently explained by a complex hybrid genome. The Labeobarbus is an evolutionary hexaploid fish of allopolyploid origin 47,48,50 . Its maternal lineage was from the tetraploid Arabibarbus Borkenhagen 2014 (=Tor) lineage 50 , distributed in the Middle East. Mainly representatives of the Arabibarbus have a generalised mouth with moderately developed lips, while few ones possess hypertrophied rubber lips 88,89 . The paternal lineage of the Labeobarbus is diploid Cyprinion, distributed in South Asia and the Middle East. Most Cyprinion species are specialized algae scrapers with a well-developed horny sheath on the lower jaw 89 . Hereby, a mouth polymorphism of Labeobarbus is apparently based on pre-existing genetic templates, a legacy of ancestors, and realized numerously as replicated pattern in plethora of isolated lineages distributed throughout Africa. Meanwhile, the novel phenotypes such as scraper2 in the Ethiopian Highlands, papillated mouth phenotype in West Africa, and numerous diversified large-mouthed phenotypes (e.g. 59 ) are likely the results of the new genomic combinations during further evolution.  www.nature.com/scientificreports www.nature.com/scientificreports/ Such a genomic prerequisite is highly important but not an exclusive explanation for the phenotypic diversification of Ethiopian Labeobarbus. Adaptive radiation is often linked to colonization of a new environment that is somehow isolated from further colonists like islands or lakes, which provide new ecological opportunities 1,90 . The significant elevation of the Ethiopian Highlands formed due to several episodes of volcanisms with one of the largest explosive volcanic events in Earth's history 81 , which probably had deleterious effects on local biota. The last major episode of volcanism in the Ethiopian Highlands was in the Pliocene-Quaternary period, however it was very recent in certain regions (33 kya Blue Nile basaltic blockade formed Tis-Isat waterfall 81 ). Indeed, the oldest fossils of Labeobarbus (late Miocene) were found in the Northern Ethiopian Rift Valley 54 but the main diversification among the Ethiopian Labeobarbus is of more recent origin (Pleistocene) in Africa 53 . Based on published phylogenetic trees 53 (and our study), the ancestor lineage that colonized the Ethiopian Highlands was from the Kenyan water bodies, like L. altianalis widely distributed in East African Rift from the northern part of Lake Tanganyika northward to Lakes Victoria and Kyoga.
The Labeobarbus is distributed along the rivers of the Ethiopian Highlands that vary gradually in several ways. For example, the elevation gradient from our sampled locations ranged from 175 m to ca. 2000 m above sea level. Interestingly, the diversification bursts were only detected in the mid-upper segments of the rivers varying between 1050 (Gojeb) and 1550 m (Sore) above sea level. This altitude effect can be attributed to a combination of two factors. First, the fauna in the mid-upper reaches is commonly poorer compared to that in the lower reaches (e.g. 15,20 ), where more diversified fauna-including highly specialized competitors-can be found. Hence, natural selection in the fauna-poor upper reaches might be relaxed due to lowered competition. Relaxed selection can help to increase morphological and ecological variability, leading to subsequent diversification. The second factor that likely contributes to the increased diversification in the middle reaches is the increased availability of ecological niches (i.e. more diverse ecotopes/habitats in this section offer more ecological opportunities). We suggest that continued diversification bursts among Labeobarbus are possible in both Ethiopian Highlands and elsewhere in its range. DNA extraction, PCR amplification and sequencing. DNA was extracted from a small piece of the fin or muscle using a standard salt method 91 or a BioSprint 15 kit for tissue and blood (Qiagen). A 1037 bp fragment of the mtDNA cytb gene was amplified by PCR using the following primers: GluDg: 5′-TGACTTGAARAACCAYCGTTG-3′ 92 and H16460: 5′-CGAYCTTCGGATTAACAAGACCG-3′ 93 Table S1 for details).

Material and
Sequence alignment and phylogenetic reconstructions. All reliable sequences from additional populations and species of Labeobarbus from Ethiopia and Kenya (n = 138) were retrieved from GenBank for comprehensive analyses (Table S2). Arabibarbus grypus (KF876026), L. beso (AF180862), L. nelspruitensis (Gilchrist & Thompson 1911) (AF180866) as well as L. ethiopicus (AF180828) were selected as outgroups. All sequences were aligned and edited using Clustal X 94 as implemented in MEGA v. 7.0 95 . The final dataset included 833 cytb sequences from all main drainages of Ethiopia, including Lake Tana.
The sequences were collapsed into common haplotypes using ALTER software 96 . DAMBE software 97 was used to analyze substitution saturation by calculating the entropy-based index of substitution saturation (Iss) and its critical value (Iss.c) 98 . The whole dataset was tested first, and then again with only the third codon positions. In both cases, the Iss values were significantly lower than the Iss.c (p < 0.0001), indicating no or little substitution saturation, hence the data were suitable for phylogenetic inference.
Bayesian phylogenetic inference (BI) was performed using MrBayes v.3.2.6 99 . Two simultaneous runs with four Markov chains each were run for 1 × 10 7 generations, sampled every 500 generations. The first 25% of runs were discarded as burn-in. Convergence of the runs was assessed by examining the average standard deviation of split frequencies and the potential scale reduction factor. In addition, stationarity was confirmed by examining posterior probability, log likelihood, and all model parameters by the effective sample sizes (ESSs) in the program Scientific RepoRtS | (2020) 10:7192 | https://doi.org/10.1038/s41598-020-64350-4 www.nature.com/scientificreports www.nature.com/scientificreports/ Tracer v1.6 100 . The best-fit model of molecular evolution by codon position for BI was estimated via the Bayesian information criterion using PartitionFinder v. 2.1.1 [101][102][103] . The BI model used the following: 1st codon position K80 + I + G, 2nd codon position HKY, 3rd codon position GTR + G. The maximum likelihood (ML) search was performed using IQ-TREE 1.6.12 104 . The best partition scheme for ML was selected in ModelTest 105 as implemented in IQ-TREE. Node support values were inferred using 1 000 bootstrap replicates. The ML model used the following: 1st codon position TIM3e + R2, 2nd codon position HKY + F, 3rd codon position GTR + F + G4.
The phylogenetic trees from the ML and BI analyses were visualized and edited using FigTree v.1.4.4 106 . In addition, a haplotype network was built with PopART 1.7 software 107 using the median joining algorithm 108 .
Bayes factor (BF) comparisons of constrained and unconstrained tree topologies were used to test for monophyly of the riverine radiations. MrBayes was used to calculate the harmonic mean estimator of marginal likelihood for trees with a hard constraint of monophyly of the examined group. Bayes factors were calculated as the difference of harmonic mean estimators of the two models (constrained vs. unconstrained) in log units. According to Kass & Raftery 109 a log difference of 3-5 is strong evidence, and a difference of > 5 is very strong evidence in favor of the better model. BI of constrained topologies were run with the same settings as described above for unconstrained phylogenetic analysis.
Genetic diversity and structure. We calculated the number of haplotypes (H), haplotype diversity (h), nucleotide diversity (π), and the average number of nucleotide substitutions (K) for each supposed adaptive radiation and geographical population (by basin) using DnaSP v.5.10 110 . Genetic differentiation among sampling locations was tested in Arlequin v. 3.5.2.2 111 using the fixation index F ST 112 based on pairwise differences or haplotype frequencies. In addition, an exact test of sample differentiation (global and pairwise), which tests the non-random distribution of haplotypes into population samples under the hypothesis of panmixia, was calculated in Arlequin v. 3.5.2.2.