Ancient hybridization fuels rapid cichlid fish adaptive radiations

Joana I. Meier; David A. Marques; Salome Mwaiko; Catherine E. Wagner; Laurent Excoffier; Ole Seehausen

doi:10.1038/ncomms14363

Ancient hybridization fuels rapid cichlid fish adaptive radiations

Meier, Joana I.; Marques, David A.; Mwaiko, Salome; Wagner, Catherine E.; Excoffier, Laurent; Seehausen, Ole 2017-02-10 00:00:00 ARTICLE Received 14 May 2016 | Accepted 20 Dec 2016 | Published 10 Feb 2017 DOI: 10.1038/ncomms14363 OPEN Ancient hybridization fuels rapid cichlid ﬁsh adaptive radiations 1,2,3 1,2,3 1,2 1,2,4 3,5 Joana I. Meier , David A. Marques , Salome Mwaiko , Catherine E. Wagner , Laurent Excofﬁer 1,2 & Ole Seehausen Understanding why some evolutionary lineages generate exceptionally high species diversity is an important goal in evolutionary biology. Haplochromine cichlid ﬁshes of Africa’s Lake Victoria region encompass 4700 diverse species that all evolved in the last 150,000 years. How this ‘Lake Victoria Region Superﬂock’ could evolve on such rapid timescales is an enduring question. Here, we demonstrate that hybridization between two divergent lineages facilitated this process by providing genetic variation that subsequently became recombined and sorted into many new species. Notably, the hybridization event generated exceptional allelic variation at an opsin gene known to be involved in adaptation and speciation. More generally, differentiation between new species is accentuated around variants that were ﬁxed differences between the parental lineages, and that now appear in many new combinations in the radiation species. We conclude that hybridization between divergent lineages, when coincident with ecological opportunity, may facilitate rapid and extensive adaptive radiation. 1 2 Aquatic Ecology and Evolution, Institute of Ecology and Evolution, University of Bern, 3012 Bern, Switzerland. Department of Fish Ecology and Evolution, Centre for Ecology, Evolution & Biogeochemistry, Eawag: Swiss Federal Institute of Aquatic Science and Technology, 6047 Kastanienbaum, Switzerland. 3 4 Computational and Molecular Population Genetics Lab, Institute of Ecology and Evolution, University of Bern, 3012 Bern, Switzerland. Biodiversity Institute & Department of Botany, University of Wyoming, Laramie Wyoming 82071, USA. Swiss Institute of Bioinformatics, 1015 Lausanne, Switzerland. Correspondence and requests for materials should be addressed to O.S. (email: [email protected]). NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 eterogeneity in diversiﬁcation rates among lineages is a of the major lakes of the region (Lakes Victoria, Edward, major factor shaping biodiversity. Yet, biological and Albert and Kivu). The largest of them is in Lake Victoria, Henvironmental factors underlying this variation are which has at least 500 endemic species that evolved in the 1 34–36 incompletely understood . Adaptive radiations are prime study past 15,000 years . Each radiation comprises enormous systems to learn about these factors as they are characterized by diversity in habitat occupation, trophic ecology, colouration a rapid origin of many species showing a diversity of ecological and behaviour. The high diversiﬁcation rate but also the high 34,37 adaptations. This process requires high levels of heritable nuclear genomic variation in the LVRS despite its young age variation in traits related to ecological and reproductive suggest that large amounts of standing genetic variation 2,3 isolation. However, adaptive radiations are often too rapid for must have been present at the onset of the radiation . the emergence of new relevant mutations between successive Previous work showing cytonuclear discordance in phylogenetic speciation events and are thus more likely to stem from reconstructions between LVRS and several riverine cichlid species standing variation. Hybridization between species can raised the possibility of ancient hybridization between divergent instantaneously boost genetic variation, which may facilitate species at the base of the radiation but could not demonstrate it . 3–9 speciation and adaptive radiation . Hybridization on secondary contact is not unlikely, as allopatric Hybridization among members of an adaptive radiation, has cichlid species, divergent by even millions of years, readily been suggested to potentially facilitate speciation events, an idea produce fertile offspring in the lab . known as the ‘syngameon hypothesis’ . Introgression of traits Using genomic data from riverine haplochromine cichlids involved in adaptation or reproductive isolation has been sampled from all major African drainage systems, and demonstrated among members of several adaptive radiations representative species from all lineages within the Lake Victoria (for example, traits related to host shift in Rhagoletis fruit ﬂies , region, we demonstrate here that the LVRS evolved from a hybrid 11,12 wing patterns in Heliconius butterﬂies or beak shape in swarm. All lake radiations show very similar proportions of Darwin’s ﬁnches ). In other radiations, the hybrid ancestry of mixed ancestry derived from two distantly related haplochromine some species has been inferred, but a direct link between lineages that had evolved in isolation from one another in introgressed traits and speciation awaits further testing (for different river systems for more than a million years before 14–16 17 example, cichlid ﬁshes of Lakes Tanganyika , Malawi , hybridizing in the Lake Victoria region. We ﬁnd evidence that 18,19 20 Victoria and Barombi Mbo ). this hybridization event facilitated subsequent adaptive radiation Another hypothesis for a perhaps more fundamental role of by providing genetic variation that has been recombined hybridization in adaptive radiation, distinct from the ‘syngameon and sorted into many new species. Variants that were ﬁxed hypothesis’, is the idea that hybridization between distinct between the parental lineages show accentuated differentiation lineages may seed the onset of an entire adaptive radiation . between young Lake Victoria species, but appear in many Such hybridization can be common when allopatric lineages new combinations in the different species. Notably, each of the come into secondary contact , and selection against hybrids may two major allele classes of an opsin gene involved in adaptation 39,40 be weak during colonization of new environments. In this and speciation among Lake Victoria cichlids is likely derived situation, the formation of a hybrid swarm, if coincident with from one of the two parental lineages. This indicates that a major ecological opportunity, may accelerate adaptive radiation by (i) part of the variation at this gene segregating in the LVRS stems providing functional genetic variation that can recombine into from hybridization between these lineages. Our results suggest novel trait combinations favoured by selection and mate choice, that hybridization between relatively distantly related species, and (ii) breaking genetic correlations that constrained the when coincident with ecological opportunity, may facilitate rapid evolvability of parental lineages . In addition, hybridization may adaptive radiation. Thus, hybridization, even in the distant facilitate speciation when multiple ﬁxed differences that confer past, may have important implications for understanding reproductive isolation between the two parental species decouple variation in extant species richness between lineages as well as and segregate in a hybrid swarm, such that selection against variation in recent rates of diversiﬁcation. incompatible gene combinations can generate more than two new 21–23 reproductively isolated species . This ‘hybrid swarm origin of adaptive radiation’ hypothesis has Results been more challenging to test. So far the only adaptive radiation Identifying the closest relatives of the LVRS. To identify the for which a hybrid origin has been robustly demonstrated is the closest relatives of the LVRS, we use comprehensive sampling of Hawaiian silverswords, which have radiated from an allopoly- haplochromine cichlids from the LVRS and from all major river ploid hybrid population between two North American tarweed systems harbouring haplochromines (Fig. 1, Supplementary species . Because gene and genome duplication are also Data 1). Maximum likelihood phylogenetic trees were recon- proposed to facilitate adaptive radiation , it is difﬁcult to structed from 3.15 Mb of concatenated restriction site associated distinguish between effects of hybridization per se and those of DNA (RAD) sequences (436,166 SNPs, 8.1% missing data, gene or genome duplication in this case. Evidence consistent with Supplementary Data 2), and from two mitochondrial markers a hybrid swarm origin of entire radiations has also been found in (1,897 bp). In the nuclear phylogeny, the entire LVRS forms Alpine whiteﬁsh , the ‘mbuna’ group of the Lake Malawi cichlid a well-supported clade that also includes cichlids from the 27 28,29 ﬁsh radiation , and allopolyploid Hawaiian endemic mints lower Nile in Egypt (Haplochromis spp. Egypt) and Haplochromis and possibly other polyploid plant radiations on Hawaii . sp. ‘Nyangara’ from the Rusizi River, the outﬂow of Lake However, it remains to be tested if, in these systems, Kivu and an inﬂow to Lake Tanganyika (orange labels in Fig. 1, hybridization occurred before or after the radiation had started, Supplementary Fig. 1, Supplementary Discussion). Astatotilapia and if hybridization-derived polymorphisms played a role in stappersi from the Kalambo River (an inﬂow to Lake Tanganyika, speciation and adaptive diversiﬁcation. Congo drainage, referred to as Haplochromis sp. ‘Chipwa’ in The Lake Victoria Region Superﬂock of cichlid ﬁsh (LVRS) Meyer et al. ), and an undescribed Astatotilapia species from is a group of 700 haplochromine cichlid species endemic to the middle Congo (A. sp. ‘Yaekama’, DRC), together hereafter the region around Lake Victoria and nearby western rift lakes in called the ‘Congolese lineage’ (red labels), form the superﬂock’s East Africa that started diversifying about 100–200 thousand closest relatives. The sister group to the LVRS plus the 31–34 years ago . It includes several adaptive radiations, one in each Congolese lineage taxa is a clade including lineages in 2 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE L. Albert ab L. Edward Mpanga R. / L. Saka L. Kivu Bootstrap values L. Edward L. Kivu >90 LVRS L. Kivu >80 L. Edward LVRS L. Victoria / Kagera L. L. Edward Haplochromis spp. Egypt Haplochromis sp. "Nyangara" Rusizi R. Astatotilapia stappersi Kalambo R. Congolese Astatotilapia sp. "Yaekama" Lineage Upper Nile Astatotilapia bloyeti Mkuzi R. / Ruaha Astatotilapia sparsidens L. Manyara Eastern Astatotilapia cf. paludinosa Malagarasi R. Astatotilapia sp. Wami R. "Haplochromis" gracilior L. Kivu Upper Nile Thoracochromis pharyngalis L. Edward Lineage Astatotilapia flavijosephi L. Kinneret Astatotilapia desfontainii Tunisia Rhamphochromis sp. Astatotilapia calliptera Lake Malawi Metriaclima zebra Haplochromines Protomelas spp. Congolese Haplochromis sp. Ruvu R. Astatotilapia cf. bloyeti Wami R. / kilossana Sonjo R. Astatotilapia sp. Rovuma Astatotilapia tweddlei L. Chilwa Astatotilapia burtoni Kalambo R. / L. Tanganyika Astatoreochromis straeleni Rusizi R. Astatoreochromis alluaudi L. Victoria / Saka Lobochilotes labiatus Lake Tanganyika Petrochromis sp. "kazumbe" Tropheini Ctenochromis horei Drainage systems Orthochromis sp. "red cheek red tail" L. Mweru Orthochromis polyacanthus L. Mweru East coast Lake Mweru Pseudocrenilabrus spp. L. Mweru Pseudocrenilabrus Haplochromis demeusii Congo Congo Serranochromis robustus L. Mweru Nile basin Lake Mweru Sargochromis spp. L. Mweru Sargochromis Haplochromis stigmatogenys Kasai/Luvua North Interior Cyclopharynx schwetzi L. Fwa Lake Fwa Endorheic drainage systems Thoracochromis callichromis L. Fwa Zambezi Ctenochromis pectoralis Pangani Figure 1 | Phylogenetic context of the Lake Victoria Region cichlid radiation. (a) Maximum likelihood phylogeny built from concatenated RADtag sequences of Lake Victoria Region Superﬂock (LVRS) cichlids and relatives including all known lineages of haplochromine cichlids (n ¼ 156). Radiations are indicated as grey triangles in the phylogenetic tree and multiple samples of a lineage are visually collapsed to a single terminal branch (full treein Supplementary Fig. 1). Members of the LVRS (including Haplochromis sp. ‘Nyangara’ from the Rusizi River and Haplochromis spp. Egypt, see Supplementary Discussion) are indicated with orange stars both in the tree and in the sampling map (b) and are labelled by lake (L) or river (R) they were sampled in. ‘Congolese lineage’ LVRS relatives are highlighted with red triangles, members of the ‘Upper Nile lineage’ with blue triangles, those from Eastern rivers with dark blue squares, and all other more distantly related lineages with black circles. (b) Sampling map. River drainage systems that we sampled are shown as coloured polygons. The radiation ancestor’s closest living relatives are shown in images: ‘Haplochromis’ gracilior/Thoracochromis pharyngalis from the Upper Nile lineage and Astatotilapia sp. ‘Yaekama’/A. stappersi from the Congolese lineage. The Lake Victoria cichlids shown in the grey triangle on the right represent some of the many and varied species that arose from the hybrid swarm (Photo credits: Ole Seehausen, Salome Mwaiko, Frans Witte, ‘Teleos’, Uli Schliewen, Adrian Indermaur, Oliver Selz; map adapted from http://www.worldwildlife.org/hydrosheds ). endorheic (that is, closed) drainage systems of the Eastern the LVRS, from which the radiations in the others lakes are Rift and in Indian Ocean drainage systems, east of the derived, and that connections between Lakes Edward and Kivu Lake Victoria region (‘Eastern taxa’, dark blue), as well as existed until recently. ‘Haplochromis’ gracilior from Lake Kivu and Thoracochromis pharyngalis from Lake Edward (light blue, Supplementary Discussion). Although the latter two species are sympatric Tests of hybridization. Phylogenetic reconstruction based on with LVRS members, they are clearly not part of the radiations concatenated sequence data assumes a single history across the themselves. We refer to H. gracilior and T. pharyngalis as genome, although recombination, incomplete lineage sorting, and the ‘Upper Nile lineage’ because all known members are conﬁned introgressive hybridisation can cause extensive genome-wide to the region that was historically the uppermost Nile variation in genealogy . To account for such genealogical 37,42 33,43 45 drainage . Consistent with previous publications ,we variation, we reconstructed a SNP based species tree . estimate the split of this lineage from the Congolese lineage The resulting tree supports the sister relationship of the LVRS to date to B1.6–5.8 million years ago, coincident with the with the Congolese taxa, but reveals incongruence between uplift of the ﬂanks of the Western branch of the East African this group and the Eastern clades and Upper Nile clade Rift truncating the paleo-river network and isolating the (Supplementary Fig. 2). To test for genetic admixture, we LVR (the new Upper Nile) from the Congo (Supplementary computed Patterson’s D statistics (ABBA-BABA test) .We Table 1; Supplementary Fig. 3). found that all LVRS radiations show strong signals of admixture Our genome-wide nuclear sampling provides unprecedented with the two Upper Nile species, H. gracilior and T. pharyngalis resolution of the relationships between radiations within (Fig. 2, Supplementary Table 2, Supplementary Fig. 4). Excess the LVRS (Fig. 1, Supplementary Fig. 1). The radiations in allele sharing between the LVRS and either species is equally Lake Victoria (including its major tributary, the Kagera River), strong (Supplementary Table 3, rows 2.1-2.4, Supplementary in Crater Lake Saka and associated Mpanga River, and in Lake Fig. 5). Importantly, neither H. gracilior nor T. pharyngalis Albert, each form monophyletic clades. In contrast, the show greater admixture with the LVRS species with which Lake Edward members of the LVRS are a paraphyletic group they live in sympatry than with allopatric members of the that includes taxa basal to the radiations of Lakes Victoria and LVRS from other lakes that have been geographically isolated for Saka, and others that form a clade together with species from many thousands of years (Fig. 2 and Supplementary Fig. 5). Lake Kivu. These ﬁndings are consistent with the hypothesis To test the directionality of gene ﬂow, we applied an extended that Lake Edward constitutes the oldest extant radiation within version of the partitioned D statistic (5 population test, NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 a b c L. Albert L. Edward L. Kivu 0.4 LVRS (P1) L. Saka g g t g t L. Victoria t g Kagera L. 0.2 t g g t P1 P2 P3 Outgroup A. stappersi Congo (P2) B A B A A. sp. "Yaekama" w w A. bloyeti l (b) ww 0.0 p p p p b s s s b s b b s b A. sparsidens (s) Eastern (P3) A. cf. paludinosa (p) −0.2 A. sp. "Wami" (w) T. pharyngalis (t) Upper Nile (P3) P1 P2 P3 Outgroup H. gracilior (g) Albert Edward Kivu Saka Victoria Kagera A B B A A. flavijosephi Outgroup P1 Figure 2 | Evidence for Congo-Nilotic hybridization in the ancestry of the LVRS. (a) Schematic genealogy with taxa used for D statistics (n ¼ 73 individuals, see Supplementary Data 2). Abbreviations used in other panels are given in parentheses and the color scheme is the same as in Figure 1. The inferred gene ﬂow edge is shown with an arrow (Note the directionality of gene ﬂow is inferred with the ﬁve-population test not shown in this ﬁgure). (b) D statistics to test for potential gene ﬂow between each Eastern and Upper Nile taxon (P3) separately (abbreviations given in a) and cichlids from each LVR lake radiation (P1) or the Congolese taxon A. stappersi (P2). Vertical bars correspond to three standard errors. Positive D values indicate gene ﬂow between P1 (LVR lake radiation) and P3 (Eastern or Upper Nile taxon), whereas negative D values indicate gene ﬂow between P2 (A. stappersi) and P3 (Eastern or Upper Nile taxon) as illustrated in (c). Exact values and more test results are given in Supplementary Table 2. ab 0.8 0.6 0.4 20 0.2 0.0 Conspe- Sister LVi-LVi LKi-LKi LVi-LKi cifics species Size of ancestry blocks (kb) Figure 3 | Congo-Nilotic ancestry blocks in LVRS genomes. (a) The size distribution of putative ancestry blocks shows mostly small ancestry blocks and slightly larger Congolese (red) than Upper Nile blocks (blue). The plot shows the counts of ancestry blocks in different size categories summed up for ﬁve LVRS radiation species across all scaffolds calculated with 3 kb windows. As most blocks do not span multiple windows of 3 kb and many blocks cannot be clearly allocated to Congolese or Upper Nile ancestry (grey) (Supplementary Fig. 8), the average ancestry block size is likely 3 kb or smaller, consistent with hybridization many thousands of years ago. (b) Correlation of ancestry blocks between whole-genome sequenced LVRS members is high overall but decreases with phylogenetic distance. The boxplots show correlation of f (ref. 78) in 10 kb windows between single individuals of conspeciﬁcs (Pundamilia individuals of the same species), sister species (Pu. pundamilia versus Pu. nyererei), more distantly related Lake Victoria (LVi) species (Paralabidochromis ﬂavus versus Pu. pundamilia and Pu. nyererei), Lake Kivu (LKi) species (Pa. paucidens versus Harpagochromis vittatus) and Lake Victoria against Lake Kivu species. This suggests that all radiation member species share the same hybridization event in their ancient history but vary in how long after that event they remained part of the same recombining population. It also suggests that some of the admixture variation still segregates within individual species (indicated by the deviation from an f correlation of 1 among conspeciﬁcs). Supplementary Table 3). We ﬁnd evidence for gene ﬂow from Importance of the ancient hybridization event. Many pheno- H. gracilior (g in Fig. 2b) and from T. pharyngalis (t in Fig. 2b) typic traits known to contribute to ecological adaptation and into each LVRS lake radiation, but no evidence for gene ﬂow in reproductive isolation have diverged in multiple speciation events the opposite direction (Supplementary Table 3, rows 1.1–1.4). in each of the LVRS radiations (for example, tooth shapes, male 48 39,40 Using F4-ratio tests , we estimated the Upper Nile ancestry nuptial colouration, opsin alleles) . Some of these traits are proportion to be B20% in all LVRS radiations (Supplementary also divergent between the Congolese Astatotilapia and Upper 37,42 Fig. 6 and Supplementary Table 4). Whole-genome sequencing Nile Thoracochromis species . Intrinsic incompatibilities of Lake Victoria and Lake Kivu radiation species, and have also been found in experimental hybrid crosses of A. stappersi, T. pharyngalis and H. gracilior, conﬁrms excess similarly divergent Astatotilapia species . In addition, some of allele sharing between LVRS members and the Upper Nile the allelic variation in the lake radiations has been shown to lineage, and reveals ancestry blocks of B3 kb, consistent with predate the origin of the LVRS (refs 49,50), perhaps most notably evidence that the admixture event occurred many thousands of that of the long-wavelength sensitive (LWS) opsin gene, one of generations ago in the common ancestor of all superﬂock the best-studied genes in Lake Victoria cichlids . This gene codes radiations in the different lakes (Fig. 3a, Supplementary Fig. 8). for the protein moiety of red-sensitive visual pigments in retinal The genome-wide signatures of admixture are correlated cones and is exceptionally diverse in Lake Victoria Region among different species from Lake Victoria and Lake cichlids . It plays a crucial role in adaptation on the Kivu (Fig. 3b), in line with a shared admixture event at the steep ambient light gradients associated with water depth and onset of the LVRS. turbidity gradients that are characteristic for the lakes in the Lake 4 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications >50 Counts (in thousands) D(BABA–ABBA)/(BABA+ABBA) f correlation d NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE Victoria region. More red-shifted LWS opsin variants are outside the radiation, LV outliers were not more often beneﬁcial in deep and murky water where the light spectrum is differentially ﬁxed in the control groups than non-outliers among relatively more red-shifted because of particulate matter category 4 SNPs (Supplementary Table 5, Supplementary 39,52 absorbing and scattering light of shorter wave lengths . This Discussion). The mosaic-like ancestry pattern is also corroborated opsin gene likely also plays a role in behavioural reproductive by our whole-genome sequence data (Supplementary Fig. 8). isolation because divergence in colour perception between species with different LWS opsin genotypes is often associated 39,49,54 with divergent male breeding colouration , which is an Discussion important mate choice cue . Here, we demonstrate that multiple large animal adaptive The LWS opsin gene is highly diverse in the LVRS with radiations arose from a hybrid swarm. Our results suggest two deeply divergent haplotype clades, often referred to as ancient admixture between two lineages of haplochromines, haplotype classes because the alleles within each clade are one from the Upper Congo and one from the Upper Nile functionally similar to each other . A substitution at amino acid drainage, at the origin of the entire Lake Victoria Region position 177 in the cichlid LWS gene shifts peak absorbance Superﬂock of cichlid ﬁsh. Contemporary Upper Nile lineage towards longer wavelength (red) in haplotype class I relative to introgression or earlier independent introgression events into class II by 7 nm (refs 49,51) and also other substitutions inﬂuence each lake radiation are unlikely given (i) the highly supported spectral sensitivity (Supplementary Discussion, Supplementary genomic monophyly of the entire Lake Victoria Region Data 3). The two haplotype classes are often associated Superﬂock (Fig. 1, Supplementary Fig. 1), (ii) similar Congo with different light environments, class I with shallow and and Upper Nile lineage ancestry proportions in all lake radiations clear water, and class II with deep and turbid water. For example, (Supplementary Fig. 6) and (iii) highly correlated Congolese and rocky shore algae scrapers of the genera Neochromis and Mbipia Upper Nile lineage ancestry tracts (Fig. 3). In addition, the fact and mud bottom detritivores of the genus Enterochromis, that neither H. gracilior (Lake Kivu) nor T. pharyngalis (Lake likely representing an early divergence event in the Lake Edward) show greater admixture into LVRS species that are Victoria radiation, occupy the opposite ends of the light sympatric with them than with allopatric members of the LVRS spectrum in Lake Victoria: clear and shallow versus murky and from other lakes (Fig. 2 and Supplementary Fig. 5) speaks against deep. All 438 LWS haplotypes previously sequenced from recent introgression. The radiations in Lakes Albert, Saka and 39,49 Neochromis and Mbipia algae scrapers are part of Victoria have been geographically separated from Lakes Edward the haplotype class I, whereas 11 of the 12 Enterochromis and Kivu, and hence from T. pharyngalis and H. gracilior, for 31,33,34 LWS haplotypes belong to class II (ref. 39). Within ecologically between 4,000 (Saka) and 100,000 (Victoria) years and variable lineages, young sister species occupying shallow versus show similar ancestry proportions. deep or clear versus turbid waters, often have different Highly similar amounts of allele sharing between the LVRS frequencies of these two LWS haplotype classes or species and both H. gracilior and T. pharyngalis (Supplementary recombinants between them, and LWS variation is consistently Table 3, rows 2.1–2.4, Supplementary Fig. 3) suggests that the associated with divergence in male nuptial coloration (Fig. 4, taxon that hybridized with the Congolese ancestor of the LVRS Supplementary Discussion). We sequenced the LWS opsin was a close relative of both Upper Nile species, or their ancestor. gene in Congolese and Upper Nile taxa, and found that the It is likely that the Congolese lineage colonized the Lake Victoria two haplotype classes in the LVRS are each shared exclusively region through capture of Malagarasi (Congo) tributaries during with just one of these parental lineages (Fig. 4, Supplementary a humid phase (for example, 145,000–120,000 years ago, Fig. 7, Supplementary Data 3, Supplementary Discussion). Supplementary Fig. 3) and encountered representatives of the All haplotypes in individuals from the Congolese lineage Upper Nile lineage that would by then have occupied the Lake take basal positions in the class I clade, whereas the Victoria region. The existence of large lakes in the Lake Victoria Upper Nile lineage haplotypes take basal positions in the class region at that time is likely to have provided ecological II clade. opportunity for diversiﬁcation, which could be exploited by To see if similar patterns of sorting of admixture variation a genetically diverse hybrid swarm. could be detected elsewhere in the genomes, we tested if The ﬁnding that the two LWS opsin haplotype classes of genomic SNPs that are highly differentiated among LVRS species the LVRS are each shared uniquely with one of the parental were commonly derived from admixture between the lineages (Fig. 4, Supplementary Fig. 7) suggests that the ancient two parental lineages. We inferred the origin of bi-allelic Congo-Nilotic admixture event was the source of the SNPs in six phenotypically diverse Lake Victoria species, the high functional variation among LVRS cichlids at the LWS opsin 39,51 insectivore Pundamilia pundamilia, the insectivore/zooplankti- gene . Recombination between the ancestral haplotypes within vore P. nyererei, the algae grazer Neochromis omnicaeruleus, the admixed populations further enhanced the functional the paedophage Lipochromis melanopterus, the insectivore variation of LWS haplotypes available for the LVRS radiation Paralabidochromis chilotes and the piscivore Harpagochromis (Supplementary Discussion, Supplementary Data 3). This gene is cf. serranus (Fig. 5). We found that at 30% of sites that particularly important for adaptation to the red-shifted end of are exceptionally strongly differentiated between Lake ambient aquatic light spectra that characterize the waters of 39,51 Victoria species (LV outliers), one of the alleles was indeed Lake Victoria . Bringing together the Congolese and Upper likely introduced into the ancestry of the radiation through Upper Nile opsin haplotypes thus appears to have facilitated adaptation Nile lineage introgression (categories 3 þ 4 in Fig. 5a, to an extreme range of light conditions and visual ecologies, with Supplementary Discussion). Sites at which the Congolese many transitions between them at both early as well as late stages and Upper Nile lineage taxa are ﬁxed for alternative alleles of the radiation. Divergence into major ecologically differentiated (category 4 in Fig. 5a) are enriched for LV outliers and show clades (Enterochromis vs Neochromis/Mbipia), but also major a pattern consistent with mosaic-like sorting of ancestral variants habitat shifts between young sister species are associated amongst Lake Victoria species (Fig. 5b). The enrichment for with recruitment of Congolese versus Upper Nile lineage LV outliers in category 4 SNPs does not seem to be because of derived haplotypes at the LWS locus (Supplementary Discussion). inherently increased ﬁxation probability of these loci, as in Our data suggest a much more general impact of the ancient pairwise comparisons of six control group cichlid species from admixture event on evolution in this set of young adaptive NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 Clear shallow water Shallow water algivores P. macrocephala 'blue' Class I Congolese lineage Deep water Turbid deep water P. macrocephala 'yellow' detritivores Class II Upper Nile lineage Figure 4 | High LWS opsin diversity likely because of the ancient hybridization event. The two major LWS opsin haplotype classes, I and II, in the LVRS (orange) are each shared exclusively with either the Congolese (A. stappersi and A. sp. ‘Yaekama’, red) or the Upper Nile lineage (H. gracilior and T. pharyngalis, blue), respectively (details in Supplementary Fig. 7; Supplementary Data 3). LWS haplotype class I is generally associated with cichlids living in shallow and clear water habitats, whereas class II is associated with deeper and more turbid habitats. Speciation by divergence in habitat type seems to have been accompanied by ﬁxation of alternative LWS haplotypes both at early and at late stages of the adaptive radiation. This is exempliﬁed by near ﬁxation of alternative haplotype classes between ecologically divergent genera such as shallow water rocky shore algae scrapers of Neochromis and Mbipia versus mud bottom detritivores of the genus Enterochromis, and by the young incipient species pair of Pundamilia macrocephala ‘blue’ (living very shallow) and ‘yellow’ (living deeper) which have predominantly alleles of haplotype class I and II, respectively (Photo credits: Ole Seehausen, Adrian Indermaur, ‘Teleos’, Oliver Selz, Uli Schliewen). radiations. Experimental crosses of cichlids have shown correlation of Congo-Nilotic ancestry block patterns decreases 38,56 57 that intrinsic incompatibilities and phenotypic novelty between species with increasing phylogenetic distance within the both increase with genetic distance between the crossed species. LVRS (Fig. 3b), it is likely that Congolese and Upper Nile lineage Both kinds of variation could be important in speciation. alleles have been segregating in the LVRS for a long time and that The divergence time of the Congolese and Upper Nile lineage they were progressively sorted among the radiation members taxa, interpreted in the context of this previous experimental during species diversiﬁcation. work, lets us predict that a hybrid swarm between these lineages The combined evidence from the LWS opsin gene and would contain both intrinsic incompatibilities and transgressive the genome-wide patterns of sorting of ancestral variation trait variation. In agreement with this prediction, we found is consistent with an important role of the ancestral hybridization that genomic sites with alternative alleles ﬁxed in the two for the evolution of the Lake Victoria Region Superﬂock. parental lineages are enriched for outlier loci likely involved in The resulting large genetic variation may help to explain how divergent adaptation and species differentiation in Lake Victoria multiple adaptive radiations in the Lake Victoria Region, together (Fig. 5). This suggests that speciation has commonly been forming 700 species, could arise in just 100–200 thousand 31,32 associated with sorting of alleles brought together in the radiation years , including 500 þ genetically and phenotypically well- ancestor by the admixture event. Sites showing negative epistatic differentiated species endemic to Lake Victoria which evolved in 32,35 interactions (Bateson Dobzhansky Muller (BDM) probably just 15,000 years . Hence, our data provide evidence incompatibilities) would be expected to be ﬁxed for alternative for the hypothesis that hybridization between divergent lineages, alleles in the parental lineages and thus fall in ancestry category 4 when coincident with novel ecological opportunity such as in Figure 5a. The strong enrichment of LV outliers in this colonisation of newly formed lakes, may facilitate rapid adaptive category is in line with differential sorting of BDM radiation through recombination and sorting of admixture- incompatibilities among species that emerged from the hybrid derived polymorphisms by natural and sexual selection. 21–23 swarm , which may have facilitated reproductive isolation That hybridization may fuel entire adaptive radiations has been among the members of the adaptive radiation. Given that the hypothesized for several systems, but robust evidence has so far 6 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NA NA NA NA NA U472 U38 LG23 LG22 LG15 LG14 LG11 LG10 LG9 LG7 LG4 LG2 NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE P<0.0001 P=0.028 P=0.0003 0.30 10,799 1,053 1,038 100 179 0.25 20% 0.20 0.15 10.1% 9.7% 0.10 7.9% 0.05 1.4% 0.00 1. 2. 3. 4. 5. Only one Both Both alleles Different Ancestral allele in alleles only if C + N alleles in LV maf C + N in C combined C and N 14-18% combined Upper Nile (N) b Congolese (C) Frequency of the Congolese allele 0.0 0.2 0.4 0.6 0.8 1.0 Pu. pundamilia N. omnicaeruleus Pu. nyererei L. melanopterus Pa. chilotes H. cf. serranus Figure 5 | Differential sorting of hybridization-derived variation in the LVRS. (a) Sites ﬁxed for alternative alleles in the Congolese (C) and Upper Nile (N) taxa are enriched for high global F outlier SNPs in Lake Victoria. Of the 12,890 biallelic SNPs among six sympatric Lake Victoria species (shown in ST (b)), 340 are outliers of high global F (LV outliers). We assigned SNPs to ﬁve different ancestry categories according to the presence or absence of the ST two alleles in the Congolese (C) and Upper Nile (N) lineage taxa. The grey bars show the proportion of LV outliers among all SNPs in each ancestry category. Total SNP counts in each category and P-values of two-sided Fisher’s exact tests are shown on top. Ancestry category (1) includes all SNPs for which only one of the two LV alleles was found in the Congolese and Upper Nile taxa together (novel LV allele or unsampled in parental lineages), (2) both LV alleles found in the Congolese taxa (polymorphic in LVRS even without Upper Nile hybridization), (3) only one allele found in Congolese but both alleles found in Upper Nile taxa (not available in LVRS without hybridization), and 4) Congolese and Upper Nile taxa each ﬁxed for alternative LV alleles (not available in LVRS without hybridization) potentially including Bateson–Dobzhansky–Muller incompatibilities. Category 5 includes sites with similar initial allele frequency in Lake Victoria (16%) than sites ﬁxed for alternative alleles in the parental lineages (category 4) to test if the enrichmentin category 4 could simply be because of high initial allele frequency. (b) Differential sorting of parental alleles between Lake Victoria cichlid species at LV outliers ﬁxed for alternative alleles in the Congolese and Upper Nile lineage taxa (mean global F among LV species ¼ 0.52). Each square represents ST a SNP coloured according to the allele frequency in that species ranging from red (ﬁxed for Congolese allele) to blue (ﬁxed for Upper Nile allele). All except two sites (2 þ 3 from the right) are located on different scaffolds of the Pundamilia nyererei reference genome. If known, chromosomal positions on the Oreochromis niloticus genome are shown below. At least 10 of the 22 chromosomes are involved in mosaic-like allele sorting between radiation species at loci that were ﬁxed for alternative alleles in the parental lineages of the ancestral hybrid swarm (Photo credits: Oliver Selz, Ole Seehausen, Adrian Indermaur, ‘Teleos’, Uli Schliewen). been conﬁned to the allopolyploid Hawaiian silverswords . performed restriction associated DNA (RAD) sequencing with haplochromine cichlid species from all major drainage systems that are either currently Here, we report strong evidence that hybrid ancestry connected to the Lake Victoria region, or have been in the past, and representatives fuelled diversiﬁcation in several large animal adaptive of all major radiations of the Lake Victoria Region Superﬂock. The closest relatives radiations. Thanks to the advent of high-throughput were inferred with maximum likelihood trees based on these data and on mito- sequencing, our power to test for ancient hybridization and to chondrial sequences of the same ﬁsh. D statistics, ﬁve population tests, and F4-ratio tests, were used to identify hybridization between lineages, determine the direction study its impact on subsequent evolution has increased of gene ﬂow and quantify ancestry proportions, respectively. Whole-genome enormously. Future studies will reveal if hybridization in the sequencing was then performed with ﬁsh of the lineages identiﬁed as ancestral to ancestry of major adaptive radiations is widespread, and whether the LVRS plus representatives of the superﬂock in order to corroborate the sig- its occurrence may explain some of the observed large variation in natures of hybridization and to estimate ancestry block sizes. the rates and volume of species diversiﬁcation among lineages. Taxonomic sampling. We sampled haplochromines from all major lakes and Methods several small lakes in the Lake Victoria region and from the Mpanga River that Experimental design. To identify the closest relatives of the Lake Victoria drains the Rwenzori Mountains in the drainage divide between Lakes Edward and Region Superﬂock (LVRS) and identify potential hybridization events, we Victoria (Supplementary Data 1, Fig. 1). Further, we collected or obtained NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 7 Proportion of LV outliers in each ‘ancestry category’ ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 the posterior distribution of the remaining 500 trees as consensus trees in from other collectors haplochromines belonging to all phylogenetic lineages Densitree . and from almost all African river drainages hosting haplochromines. We also included members of the haplochromine species ﬂocks of the other African Great Lakes Malawi, Tanganyika and Mweru (Fig. 1, Supplementary Data 1). Mitochondrial chronograms. Dated phylogenies were reconstructed based on Samples were collected under research and sample export permissions of the mitochondrial D-loop and ND2 (ref. 64) sequences using BEAST v. 2.3.0 (ref. 71) Tanzania Fisheries Research Institute and the Tanzanian Ministry of Agriculture, and four different sets of calibration nodes (Supplementary Methods). We Livestock and Fisheries; the Ugandan National Fisheries Resources Research caution that the mitochondrial tree only shows the phylogeny of the maternal Institute and Ministry of Agriculture, Animal Industry and Fisheries Uganda; line and that time estimates more recent than one million years are most likely and the Department of Fisheries Zambia. All samples were collected in compliance overestimates because of the increase of molecular rates towards the recent with applicable international and national guidelines for the use of animals, (Supplementary Methods). and ethical standards. To compare the splitting time between cichlid lineages with the major relevant geological events, we reconstructed paleogeographic maps at different time points based on previously published data and reviews (see Supplementary RAD sequencing. DNA was extracted from ﬁn clips or muscle tissue with Methods). a standard phenol-chloroform-isoamyl alcohol extraction method . Restriction- site Associated DNA sequencing (RADseq) was performed following a standard Patterson’s D statistics. To test for evidence of ancient admixture among protocol . Restriction digestion was done overnight (8–10 h) using the 46,72 restriction endonuclease HF-SbfI (NewEngland Biolabs) and 400–1,000 ng lineages, we computed Patterson’s D statistic (ABBA-BABA test), a method to detect admixture based on the frequencies of discordant SNP genealogies DNA per sample. P1 adaptors contained 5–8 bp long barcodes differing by at least two nucleotides from all other barcodes. The DNA was sheared with a Covaris in a four-taxon tree, with the software package ADMIXTOOLS v. 1.1 (ref. 48). Genotypes were discarded if they had less than 6 reads or a genotype quality Phred S220 Focused-Ultra sonicator and fragments of 300–600 bp length were extracted from an agarose gel. We performed 18 PCR cycles to amplify the RAD fragments scoreo20 (that is, error probability41%). Signiﬁcance of D statistics was assessed with a block jackknife procedure using a z score of three standard errors as (30 s 98 C, 18 (10 s 98 C, 30 s 65 C, 30 s 72 C), 5 min 72 C). All libraries a threshold . We used three individuals of Astatotilapia ﬂavijosephi from Lake were single-end sequenced on an Illumina HiSeq 2,500 sequencer. The reads were de-multiplexed and trimmed to 84 nucleotides (nt, after barcode removal) with Kinneret as the outgroup population and we tested for evidence of gene ﬂow from each Eastern and Upper Nile clade into the LVRS relative to allele sharing with the the process_radtags script from the Stacks pipeline , correcting single errors in the barcode and discarding reads with incomplete restriction sites. The length closest relative of the LVRS, the Congo drainage taxa A. stappersi from Zambia (Fig. 2) or A. sp. ‘Yaekama’ from the central Congo. For these analyses we used all of 84 nt results from removing the barcode (mostly 6 nt long) and trimming off the last 10 nt because of reduced quality at the read ends. The FastX toolkit individuals with at least 50% of the sites sequenced at a depth of at least 10 reads. To exclude the possibility that the D statistics are biased by the alignment of the (http://hannonlab.cshl.edu/fastx_toolkit) was used to remove all reads containing at least one base with a Phred quality score below 10 and reads with more than reads to the Lake Victoria species Pundamilia nyererei, we also aligned the reads to the Astatotilapia burtoni reference genome, which is an outgroup to all taxa used in 10% of bases with quality less than 30. The reads of each individual were then mapped to the Pundamilia nyererei the D statistics. To rule out that the choice of outgroup (A. ﬂavijosephi) biases the D statistics, we repeated the tests with A. burtoni or A. desfontainii as outgroups. reference genome using Bowtie2 (ref. 61) with the end-to-end alignment option. Single nucleotide polymorphisms (SNPs) and genotypes were called using All species combinations with number of individuals and SNPs included are given in Supplementary Table 2. GATK Uniﬁed Genotyper v. 3.1 (ref. 62). All sites were then ﬁltered with a custom- made Python script and vcftools v. 4.1 (ref. 63). SNPs within 5 nt from indels (insertions and deletions) were removed to avoid false SNPs because of Five population tests. To infer directionality of gene ﬂow between LVRS and the misalignment problems, and SNPs were required to have a quality value of Nilotic taxa, we used an extended version of the partitioned D statistic test at least 30. developed by Eaton & Ree (we call it a ‘ﬁve population test’). Similar to the D statistics, this test is based on discordant allele sharing patterns, but by considering ﬁve populations it allows one to infer the directionality of gene ﬂow. Mitochondrial sequencing. Two mitochondrial markers (NADH Dehydrogenase The ﬁve taxa with the topology ((P1,P2), (P3a,P3b)),O include the potential source Subunit 2 (ND2) using the primers ND2Met-F 5 -CAT ACC CCA AAC ATG of gene ﬂow (for example, P3a), the taxon receiving gene ﬂow (for example, P1), 0 0 0 TTG GT-3 and ND2Trp-R 5 -GTS GST TTT CAC TCC CGC TTA-3 and the a close relative of each of these two taxa (for example, P3b and P2) and an 65 0 mtDNA control region (D-loop) with the primers FISHL15926-F 5 -GAG CGC outgroup (O, see also Supplementary Table 3 for visualization). If genes had 0 0 CGG TCT TGT AA-3 and FISH12s-R 5 -TGC GGA GAC TTG CAT GTG TAA introgressed from for example, P3a into P1, a close relative (P3b) of the G-3 were ampliﬁed with PCR and Sanger sequenced for the same individuals introgression donor (P3a) would also show excess allele sharing with P1, but or downloaded from GenBank for the same species as those included in the to a lesser extent than P3a. This is because many derived alleles that introgressed RADseq dataset (Supplementary Data 1). The sequences were aligned in ClustalW from P3a into P1 will be shared by P3a and P3b because of their recent common implemented in BioEdit 7.2.5 (ref. 66) and manually curated for correct local ancestry. In contrast, we would not expect that a close relative (P2) of the receiver alignment. of gene ﬂow, P1, would show excess allele sharing with P3a and P3b. In the genomic data, this would be seen as an excess number of BABBA patterns, where P1 shares a derived allele (‘B’) with both P3a and P3b, whereas P2 has the outgroup Phylogenetic analyses. Phylogenies were reconstructed for the concatenated allele (‘A’), as compared with the number of ABBBA patterns (P2 shares a derived mitochondrial genes and the concatenated RAD sequences separately, including allele with both P3 taxa). On the other hand, if the direction of gene ﬂow was both variant and invariant sites using a maximum likelihood approach from P1 into P3a, P3b would not show excess allele sharing with P1 but instead, 67,68 (RAxML v. 7.7.7 and ExaML v. 1.0.4) . For the mitochondrial dataset, P2 would show excess allele sharing with P3a because of ancestrally shared a maximum likelihood tree was reconstructed using three partitions, one for Dloop, alleles between P1 and P2 that introgressed into P3a. Therefore, gene ﬂow from one for the ﬁrst and the second codon positions of ND2, and one for the third P1 into P3a would not affect the relative frequencies of BABBA and ABBBA codon position of ND2. For each dataset, we performed a RAxML analysis with patterns, but the pattern BBBAA (P1, P2 and P3a share a derived allele) would be 100 rapid bootstraps using the GTRGAMMA model of rate heterogeneity. For the more frequent than BBABA (P1 and P2 share a derived allele with P3b, see RAD-seq dataset, we used all concatenated sites with no more than 40 individuals Supplementary Table 3). Counting discordant allele sharing patterns thus allows us missing (25%) to reconstruct a maximum likelihood tree with RAxML (ref. 67) to infer the direction of gene ﬂow. and ExaML (ref. 68). Each of 100 bootstraps was performed by randomly sampling We computed the ﬁve population tests with a custom made script using the with replacement sites from the concatenated dataset to get a dataset of the original three individuals with least missing data for each LVRS group and the single size. The maximum likelihood tree was then inferred for each resampled individual with the most complete data for all other taxa. In contrast to the dataset with ExaML (ref. 68) using a GTRGAMMA model of rate heterogeneity, D statistics computed with ADMIXTOOLS, which are based on allele frequencies, as recommended in the RAxML-light manual . We calculated bootstrap support our ﬁve population tests are calculated from a single individual for each values based on these 100 topologies with RAxML (ref. 67). The nuclear and focal population, following Eaton and Ree . We tested each of the three mitochondrial trees were rooted with the reference genome of Oreochromis individuals of each LVRS group separately, and report the means for each niloticus and ladderized and plotted using the R-package Ape v. 3.1 (ref. 69). radiation. At heterozygous sites, one allele was chosen at random. For each We then used RAD-derived SNP data to infer the species tree for the LVRS combination of individuals tested, we counted the eight patterns needed to groups and their closest relatives with SNAPP (ref. 45). SNAPP bypasses gene compute the four D statistics (Supplementary Table 3) using all sites without trees and computes species trees directly from independently inherited markers missing data. We calculated z scores in units of s.d. from 100 bootstrapped datasets by integrating over all possible gene trees . We restricted this analysis to two 47 (sites resampled with replacement) as in Eaton and Ree . individuals per species (one for A. sp. ‘Yaekama’) and as SNAPP assumes no linkage among loci, we included only biallelic sites that were at least 500 kb apart from each other. The resulting data set contained 31 individuals and 1,817 sites. Estimation of ancestry proportions. The absolute value of D statistics does not We ran SNAPP for 1,000,000 iterations, sampling every 1,000th iteration using only depend on the admixture proportion but also on the demographic history, the default priors. We discarded the ﬁrst 50% of the trees as burn-in and visualized divergence between the hybrid parental lineages, and the genetic distance between 8 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE 73 inferred Upper Nile ancestry proportion for the Lake Victoria radiation). To the real source taxon of gene ﬂow and the sampled surrogate population . take the effect of this moderately high initial allele frequency into account in our Thus, the D statistic cannot be used directly to infer the magnitude of 48,73 subsequent analyses, we created a dataset of SNPs with comparable expected introgression. Instead, we applied the F4-ratio test to infer the Upper Nile ancestral allele frequency in Lake Victoria (category 5) by weighting the allele ancestry proportions in the different LVRS taxa using all genotypes with at least ten frequencies of the Congolese and Upper Nile lineage taxa at each site by 0.84 and reads and SNPs with a maximum of 10% missing data (72,443 SNPs) with 0.16, respectively, corresponding to the inferred ancestry proportions for the Lake ADMIXTOOLS v. 1.1 (ref. 48). The F4-ratio test is based on four populations with Victoria radiation. Of these, we extracted all sites with a weighted minor allele the genealogy (((A,B)C)O) and a ﬁfth potential hybrid population X that is tested frequency of 14–18% excluding sites that are divergently ﬁxed between the for ancestry proportions from B and C (Supplementary Fig. 6). We repeated the Congolese and Upper Nile taxa. As an example, if at a given bi-allelic site the test with different Eastern taxa used as population ‘A’, T. pharyngalis or H. gracilior frequency of one allele was 0.1 in the Congolese taxa and the frequency of the same as ‘B’, A. sp. ‘Yaekama’ as ‘C’ and A. ﬂavijosephi as ‘O’ as shown in Supplementary allele was 0.5 in the Upper Nile taxa, we calculated 0.1 0.84 þ 0.5 0.16 ¼ 0.164. Table 4. The Upper Nile ancestry proportions estimated with the combination of As this site would have a weighted minor allele frequency between 14% and 18%, it populations with the lowest s.d. were used for further analyses and to explain the would be included in the control set of SNPs. Note that such loci are not predicted F4-ratio test in Supplementary Fig. 6. to be involved in BDM incompatibilities because all their alleles were present within at least one of the two parental populations. The proportion of sites Long-wave sensitive (LWS) opsin gene. We sequenced the LWS opsin gene from that were outliers of high global F (Po0.05) among the Lake Victoria species ST exon 2 to 6 using the Carleton & Kocher primer combinations F2 (5 -TTT GAG (‘LV outliers’) was calculated for each dataset and compared between datasets 0 0 0 GGT CCC AAT TAC CA-3 ), R2 (5 -TCC ACA CAG CAA GGT AGC AC-3 ) and using a two-sided Fisher’s exact test. The robustness of the test to erroneously 0 0 0 F3 (5 -ACT GGC CTC ATG GAC TGA AG-3 ) and R4 (5 -TCC CAA AAT GGA inferred ancestry proportions was assessed by repeating the calculations GAA CAT GG-3 ) (http://cichlid.umd.edu/cichlidlabs/protocols/Basic/pcrpri- with different Upper Nile allele frequency weighting, ranging from 10 to mer.html) . The sequences were aligned together with all available Lake Victoria 30% (Supplementary Discussion). cichlid LWS sequences from GenBank using ClustalW implemented in BioEdit To check if LV outliers ﬁxed for alternative alleles in the parental lineages 7.2.5 (ref. 66) and then manually curated. For sequences used, see Supplementary inherently have a higher ﬁxation probability than non-outliers, we analysed ﬁxation Data 3. A maximum likelihood tree was reconstructed including both variant patterns in other closely related cichlid species that are not part of the radiation. and invariant sites using the GTRGAMMA model of rate heterogeneity with We used Astatotilapia bloyeti, A. sparsidens, A. ﬂavijosephi, A. calliptera, A. burtoni RAxML v. 7.7.7 (ref. 67). We used one partition with concatenated introns, one for and Astatoreochromis alluaudi as we have at least three sequenced individuals of the ﬁrst plus second codon position of concatenated exons, and one for the each. For each pairwise comparison between these species, we extracted sites with third codon positions. The tree was rooted with the Oreochromis niloticus at least three genotypes sequenced at a minimum depth of six reads in both species. sequence and nodal support values were drawn from 100 rapid bootstraps We then checked if sites ﬁxed for alternative alleles in the Congolese and performed with RAxML. The tree was ladderized and plotted using the R-package Upper Nile lineage taxa were also ﬁxed for alternative alleles between these other Ape v. 3.1 (ref. 69). species. Finally, we tested with Fisher’s exact tests if sites that were high global To visually assess the effects of recombination on the LWS opsin alleles in the F outliers among Lake Victoria species were also more often differentially ﬁxed ST Lake Victoria Region cichlids, we extracted all positions of the alignment that were between these other species than non-outlier sites (see also Supplementary polymorphic among the LVRS sequences and also differed between the Congolese Discussion). and Upper Nile taxa. If one of the alleles was only found in the Congolese taxa but To analyse the genomic distribution of LV outliers that are ﬁxed for not in the Upper Nile taxa, we coloured it red, whereas alleles found only in alternative alleles between the Congolese and Upper Nile taxa, we used liftOver the Upper Nile taxa were coloured blue. Many LVRS haplotypes could be explained (http://genome.ucsc.edu/) with a chain ﬁle from Brawand et al. to get positions by recombination between the two divergent allele classes occurring in the Congo on the Oreochromis niloticus genome which includes chromosomal information. and Upper Nile clades. The second LWS gene tree was reconstructed with all recombinant sequences falling between the two major allele classes in the phylogenetic tree removed from the alignment. Whole-genome sequencing data. To corroborate our ﬁndings from For information about the occurrence of LWS haplotypes across LVRS species, 51 RAD sequence data, and for the analysis of ancestry block sizes, we sequenced we compiled data from Tables 1 þ 2 in Terai et al. , Fig. 1 and Supplementary 49 39 whole genomes of nine ﬁsh from Lake Victoria, two from Lake Kivu, and of the Table 4 in Terai et al. , Supplementary Fig. 4b in Seehausen et al. , Figure S3 in 54 taxa that we found to be the closest extant representatives of the lineages directly Miyagi et al. and additional GenBank sequences by Carleton et al. (GenBank). ancestral to the LVRS, the Congolese Astatotilapia stappersi, and the Upper Nile Information sources and GenBank accession numbers are provided in Thoracochromis pharyngalis and ‘Haplochromis’ gracilior. Whole-genome Supplementary Data 3. sequencing data was generated using PCR-free library preparation and Illumina HiSeq 3000 paired-end sequencing for 11 individuals. Three additional individuals, Sorting of ancestral alleles. We created a RAD sequence dataset of our sequenced the same way on the same machine, were taken from McGee et al. Congolese and Upper Nile LVRS relatives and previously published sequences of (see Supplementary Data 1 for sample information). Local alignment against the six sympatrically occurring Lake Victoria species that are ecologically and Astatotilapia burtoni reference genome was performed with Bowtie 2 (ref. 61). phenotypically diverse (Pundamilia pundamilia (n ¼ 10, blue male nuptial For variant calling and genotyping we used Haplotype Caller (GATK v. 3.5) . coloration, benthic insectivore-omnivore), P. nyererei (n ¼ 11, red planktivore- Genotypes with fewer than ﬁve reads, multiallelic sites, and indels were removed omnivore), Harpagochromis cf. serranus (n ¼ 10, blue piscivore), Lipochromis using vcftools v. 4.1 (ref. 63). melanopterus (n ¼ 9, yellow paedophage), Paralabidochromis chilotes (n ¼ 9, blue Whole-genome D statistics were calculated with ADMIXTOOLS v. 1.1 (ref. 48) benthic insectivore) and Neochromis omnicaeruleus (n ¼ 12, blue algae scraper)). using all biallelic sites with a 1% minor allele frequency cutoff and a maximum The sites were ﬁltered for having at least 5 individuals covered in each species missing data proportion of 20% across all 14 genomes. The reference genome at a sequencing depth of at least six reads. We used Arlequin v. 3.514 (ref. 75) to (A. burtoni) was used as an outgroup. calculate global F values based on a hierarchical island model among the six Lake To study signatures of admixture along the genome, we used f statistics by ST d Victoria cichlid species. We calculated allele frequencies for the Congolese Martin et al. . In comparison to D statistics, f is more suited for small genomic (Astatotilapia stappersi from Zambia and A. sp. ‘Yaekama’ from DRC, n ¼ 5) and regions . ABBA and BABA pattern counts were calculated using allele frequencies the Upper Nile taxa (‘Haplochromis’ gracilior and Thoracochromis pharyngalis, by weighting each segregating site according to its ﬁt to the ABBA or BABA 46,72,78 n ¼ 9) to get estimates of the allele frequencies for the parental lineages of the LVRS pattern . As an example, if a site has derived allele frequencies of 0, 0.5, 1 and and to infer the likely ancestry of the alleles segregating in Lake Victoria. We 0 in P1, P2, P3 and the outgroup, respectively, it would count as half ABBA site. excluded sites monomorphic in the subset of Lake Victoria species or sequenced in We used A. stappersi as P1, different LVRS individuals as P2, H. gracilior and less than three individuals in either the Congolese or the Upper Nile lineage taxa. T. pharyngalis as P3 and A. burtoni as outgroup (Supplementary Fig. 8). Note that We extracted sites for which only one of the two alleles segregating in Lake Victoria P1 and P2 are switched as compared to D statistics, where LVRS are used as P1 and cichlids were found in the Congolese and Upper Nile lineage taxa together the Congolese lineage as P2. This difference is simply for consistency with the (category 1). Next, we extracted sites with both alleles segregating in the Congolese description of f in Martin et al. , and thus here ABBA represents sites with taxa representing ancestral standing genetic variation that would have been present a derived allele shared between LVRS and Upper Nile, whereas BABA represents in the Lake Victoria radiation also without ancient hybridization with the Upper sites with a derived allele shared between Congolese and Upper Nile lineage taxa. Nile lineage (category 2). f is calculated as the difference between ABBA and BABA patterns compared to We identiﬁed Lake Victoria polymorphic sites for which only one allele was the maximum possible difference where gene ﬂow between P2 and P3 would equal found in the Congolese taxa, and the second allele was found exclusively in the random mating . Positive f values indicate gene ﬂow between P2 (LVRS) and Upper Nile taxa, representing genetic variation that can be attributed to the ancient P3 (Upper Nile). We calculated f in non-overlapping sliding windows of 10 kb admixture event (category 3). We also identiﬁed the subset of sites ﬁxed for along the A. burtoni scaffolds using the Python script by Martin et al. . Windows alternative alleles in the Congolese and Upper Nile taxa (category 4). These are sites with less than ﬁve total ABBA and BABA patterns were excluded. As only positive that could be involved in epistatic BDM-like incompatibilities between Congolese f values are indicative of excess allele sharing between Upper Nile and LVRS, and and Upper Nile alleles. Sites ﬁxed for alternative alleles in the Congolese and are correctly standardized, values with negative D scores were set to 0 as in Martin Upper Nile lineage taxa are expected to have been present at the onset of the et al. . Pearson’s product-moment correlations of f values between different Lake Victoria radiation at a mean minor allele frequency of 16% (that is, the individuals were calculated with the stats R-package. NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 Putative Congolese or Upper Nile ancestry was assessed by comparing the 20. Schliewen, U. K. & Klee, B. Reticulate sympatric speciation in Cameroonian frequency of ABBA sites (LVRS shares the derived allele exclusively with crater lake cichlids. Front. Zool. 1, 5 (2004). Upper Nile taxa) with the frequency of BBAA sites (LVRS shares the derived 21. Seehausen, O. Conditions when hybridization might predispose populations for allele exclusively with the Congolese lineage representative A. stappersi)in adaptive radiation. J. Evol. Biol. 26, 279–281 (2013). non-overlapping sliding windows of 3 kb. Scaffolds with o100 kb data were 22. Schumer, M., Cui, R., Rosenthal, G. G. & Andolfatto, P. Reproductive isolation removed. As for the f statistic, ABBA and BBAA pattern counts were of hybrid populations driven by genetic incompatibilities. PLoS Genet. 11, calculated using allele frequencies by weighting each segregating site according e1005041 (2015). to its ﬁt to the ABBA or BBAA pattern. Only windows with a total ABBA and 23. Hermansen, J. S. et al. Hybrid speciation through sorting of parental BBAA pattern count exceeding one were used. Windows with a minimum ABBA incompatibilities in Italian sparrows. Mol. Ecol. 23, 5831–5842 (2014). proportion ((ABBA/(ABBA þ BBAA)) of 0.7 were deﬁned as candidate windows of 24. Barrier, M., Baldwin, B. G., Robichaux, R. H. & Purugganan, M. D. Interspeciﬁc Upper Nile lineage ancestry and are coloured in blue, whereas genomic regions hybrid ancestry of a plant adaptive radiation: alloploidy of the Hawaian with an ABBA proportion of 0.3 or less were deﬁned as putative Congolese lineage silversword alliance (Asteraceae) inferred from ﬂoral homeotic gene derived windows and are highlighted in red. Ancestry tracts were deﬁned as duplication. Mol. Biol. Evol. 16, 1105–1113 (1999). consecutive sliding windows of the same colour (red or blue) ignoring single 25. Lynch, M. & Conery, J. S. The evolutionary fate and consequences of duplicate sliding windows without data. Expected ancestry block sizes were calculated using a formula from Racimo et al. . Assuming B50,000–100,000 generations genes. Science 290, 1151–1155 (2000). since admixture, a recombination rate of 2.5 10 and an Upper Nile 26. Hudson, A. G., Vonlanthen, P. & Seehausen, O. Rapid parallel adaptive proportion of 20%, the expected mean length of the admixture tracts is radiations from a single hybridogenic ancestral population. Proc. R Soc. B 278, 8 1 (0.8 2.5 10 (50,000 to 100,000) 1)) ¼ 500 bp to 1 kb (ref. 79). 58–66 (2011). The breakup of ancestry blocks may have slowed as the original hybrid population 27. Joyce, D. A. et al. Repeated colonization and hybridization in Lake Malawi began to undergo genomic stabilization, speciated and formed geographically cichlids. Curr. Biol. 21, R108–R109 (2011). isolated radiations in separate lakes each of which underwent further stabilization 28. Lindqvist, C., Motley, T. J., Jeffrey, J. J. & Albert, V. A. Cladogenesis and independent of each other, likely associated with differential sorting of the ancestral reticulation in the Hawaiian endemic mints (Lamiaceae). Cladistics 19, 480–495 variation. (2003). 29. Lindqvist, C. & Albert, V. A. Origin of the Hawaiian endemic mints within North American Stachys (Lamiaceae). Am. J. Bot. 89, 1709–1724 Data availability. Mitochondrial and LWS opsin sequences are available on ð2002Þ: GenBank under the accession numbers KY366716-KY366843 for ND2, KY366844- 30. Baldwin, B. G. & Wagner, W. L. Hawaiian angiosperm radiations of North KY366970 for D-loop, and KY366971-KY366986 for LWS opsin sequences. American origin. Ann. Bot. 9, 849–879 (2010). RADseq and whole-genome sequencing reads generated in this study can be 31. Verheyen, E., Salzburger, W., Snoeks, J. & Meyer, A. Origin of the superﬂock downloaded from the NCBI Sequence Read Archive under Bioproject of cichlid ﬁshes from Lake Victoria, East Africa. Science 300, 325–329 PRJNA355227. The JAVA program for 5-population tests is publicly available on ð2003Þ: GitHub (https://github.com/joanam/scripts). 32. Seehausen, O. African cichlid ﬁsh: a model system in adaptive radiation research. Proc. R. Soc. B 273, 1987–1998 (2006). 33. Genner, M. J. et al. Age of cichlids: New dates for ancient lake ﬁsh radiations. References Mol. Biol. Evol. 24, 1269–1282 (2007). 1. Wagner, C. E., Harmon, L. J. & Seehausen, O. Ecological opportunity and 34. Bezault, E., Mwaiko, S. & Seehausen, O. Population genomic tests of models sexual selection together predict adaptive radiation. Nature 487, 366–369 of adaptive radiation in Lake Victoria Region cichlid ﬁsh. Evolution 65, (2012). 3381–3397 (2011). 2. Barrett, R. D. H. & Schluter, D. Adaptation from standing genetic variation. 35. Johnson, T. C., Kelts, K. & Odada, E. The holocene history of Lake Victoria. Trends Ecol. Evol. 23, 38–44 (2008). Ambio 29, 2–11 (2000). 3. Seehausen, O. Hybridization and adaptive radiation. Trends Ecol. Evol. 19, 36. Wagner, C. E. et al. Genome-wide RAD sequence data provide unprecedented 198–207 (2004). resolution of species boundaries and relationships in the Lake Victoria cichlid 4. Abbott, R. et al. Hybridization and speciation. J. Evol. Biol. 26, 229–246 adaptive radiation. Mol. Ecol. 22, 787–798 (2013). ð2013Þ: 37. Seehausen, O. et al. Nuclear markers reveal unexpected genetic variation and 5. Rieseberg, L. H. Hybrid origins of plant species. Ann. Rev. Ecol. Syst. 28, a Congolese-Nilotic origin of the Lake Victoria cichlid species ﬂock. Proc. 359–389 (1997). R. Soc. Lond. B Biol. Sci. 270, 129–137 (2003). 6. Mallet, J. Hybrid speciation. Nature 446, 279–283 (2007). 38. Stelkens, R. B., Young, K. A. & Seehausen, O. The accumulation of reproductive 7. Anderson, E. & Stebbins, G. L. Hybridization as an evolutionary stimulus. incompatibilities in African cichlid ﬁsh. Evolution 64, 617–632 (2010). Evolution 8, 378–388 (1954). 39. Seehausen, O. et al. Speciation through sensory drive in cichlid ﬁsh. Nature 8. Carlquist, S. Island Biology Vol. 581, 5279 (Columbia Univ Press, 1974). 455, 620–626 (2008). 9. Arnold, M. L. Natural Hybridization and Evolution (Oxford University Press, 40. Seehausen, O., van Alphen, J. & Witte, F. Can ancient colour polymorphisms 1997). explain why some cichlid lineages speciate rapidly under disruptive sexual 10. Feder, J. L. et al. Mayr, Dobzhansky, and Bush and the complexities of selection? Belg. J. Zool. 129, 43–60 (1999). sympatric speciation in Rhagoletis. Proc. Natl Acad. Sci. USA 102, 6573–6580 41. Meyer, B. S. et al. Back to Tanganyika: a case of recent trans-species-ﬂock (2005). dispersal in East African haplochromine cichlid ﬁshes. R. Soc. Open Sci. 2, 11. The Heliconius Genome Consortium. Butterﬂy genome reveals promiscuous 140498 (2015). exchange of mimicry adaptations among species. Nature 487, 94–98 42. Greenwood, P. H. Towards a phyletic classiﬁcation of the ‘genus’ Haplochromis ð2012Þ: (Pisces, Cichlidae) and related taxa. Part 1. Bull. Br. Mus. Nat. Hist. Zool. 35, 12. Pardo-Diaz, C. et al. Adaptive introgression across species boundaries in 265–322 (1979). Heliconius butterﬂies. PLOS Genet. 8, e1002752 (2012). 43. Hermann, C. M., Sefc, K. M. & Koblmu¨ller, S. Ancient origin and recent 13. Lamichhaney, S. et al. Evolution of Darwin’s ﬁnches and their beaks revealed by genome sequencing. Nature 518, 371–375 (2015). divergence of a haplochromine cichlid lineage from isolated water bodies in the East African Rift system. J. Fish Biol. 79, 1356–1369 (2011). 14. Salzburger, W., Baric, S. & Sturmbauer, C. Speciation via introgressive 44. Maddison, W. P. Gene trees in species trees. Syst. Biol. 46, 523–536 (1997). hybridization in East African cichlids? Mol. Ecol. 11, 619–625 (2002). 45. Bryant, D., Bouckaert, R., Felsenstein, J., Rosenberg, N. A. & RoyChoudhury, A. 15. Meyer, B. S., Matschiner, M. & Salzburger, W. Disentangling incomplete lineage sorting and introgression to reﬁne species-tree estimates for Lake Inferring species trees directly from biallelic genetic markers: bypassing Tanganyika cichlid ﬁshes. Syst. Biol. doi: 10.1093/sysbio/syw069 (2016). gene trees in a full coalescent analysis. Mol. Biol. Evol. 29, 1917–1932 16. Weiss, J. D., Cotterill, F. P. & Schliewen, U. K. Lake Tanganyika—A’Melting ð2012Þ: 46. Green, R. E. et al. A draft sequence of the Neandertal genome. Science 328, Pot’of Ancient and Young Cichlid Lineages (Teleostei: Cichlidae)? PLoS ONE 10, e0125043 (2015). 710–722 (2010). 17. Genner, M. J. & Turner, G. F. Ancient hybridization and phenotypic 47. Eaton, D. A. R. & Ree, R. H. Inferring phylogeny and introgression using novelty within Lake Malawi’s cichlid ﬁsh radiation. Mol. Biol. Evol. 29, 195–206 RADseq data: an example from ﬂowering plants (Pedicularis: Orobanchaceae). (2012). Syst. Biol. 62, 689–706 (2013). 18. Keller, I. et al. Population genomic signatures of divergent adaptation, gene 48. Patterson, N. et al. Ancient admixture in human history. Genetics 192, ﬂow and hybrid speciation in the rapid radiation of Lake Victoria cichlid ﬁshes. 1065–1093 (2012). Mol. Ecol. 22, 2848–2863 (2013). 49. Terai, Y. et al. Divergent selection on opsins drives incipient speciation in Lake 19. Meier, J. I. et al. Demographic modelling with whole-genome data reveals Victoria cichlids. PLoS Biol. 4, 2244–2251 (2006). 50. Brawand, D. et al. The genomic substrate for adaptive radiation in African parallel origin of similar Pundamilia cichlid species after hybridization. Mol. Ecol, 26, 123–141 (2017). cichlid ﬁsh. Nature 513, 375–381 (2014). 10 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE 51. Terai, Y., Mayer, W. E., Klein, J., Tichy, H. & Okada, N. The effect of selection 76. Kozarewa, I. et al. Ampliﬁcation-free Illumina sequencing-library preparation on a long wavelength-sensitive (LWS) opsin gene of Lake Victoria cichlid ﬁshes. facilitates improved mapping and assembly of (G þ C)-biased genomes. Nat. Methods 6, 291–295 (2009). Proc. Natl Acad. Sci. USA 99, 15501–15506 (2002). 77. McGee, M. D., Neches, R. Y. & Seehausen, O. Evaluating genomic divergence 52. Okullo, W. et al. Parameterization of the inherent optical properties of Murchison Bay, Lake Victoria. Appl. Opt. 46, 8553–8561 (2007). and parallelism in replicate ecomorphs from young and old cichlid adaptive 53. Maan, M. E., Hofker, K. D., van Alphen, J. J. M. & Seehausen, O. Sensory drive radiations. Mol. Ecol. 25, 260–268 (2016). 78. Martin, S. H., Davey, J. W. & Jiggins, C. D. Evaluating the use of ABBA–BABA in cichlid speciation. Am. Nat. 167, 947–954 (2006). statistics to locate introgressed loci. Mol. Biol. Evol. 32, 244–257 (2015). 54. Miyagi, R. et al. Correlation between nuptial colors and visual sensitivities 79. Racimo, F., Sankararaman, S., Nielsen, R. & Huerta-Sanchez, E. Evidence tuned by opsins leads to species richness in sympatric Lake Victoria cichlid for archaic adaptive introgression in humans. Nat. Rev. Genet. 16, 359–371 ﬁshes. Mol. Biol. Evol. 29, 3281–3296 (2012). (2015). 55. Selz, O. M., Pierotti, M. E. R., Maan, M. E., Schmid, C. & Seehausen, O. Female 80. Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from preference for male color is necessary and sufﬁcient for assortative mating in spaceborne elevation data. EOS 89, 93–94 (2008). 2 cichlid sister species. Behav. Ecol. 25, 612–626 (2014). 56. Stelkens, R. B., Schmid, C. & Seehausen, O. Hybrid breakdown in cichlid ﬁsh. PLoS ONE 10, e0127207 (2015). Acknowledgements 57. Stelkens, R. B., Schmid, C., Selz, O. & Seehausen, O. Phenotypic novelty in We thank all collaborators who provided tissue samples over the years: Sigal Balshine- experimental hybrids is predicted by the genetic distance between species of Earn (McMaster), Roger Bills (SAIAB), Julia Day (UCL), Yves Fermon (Paris), John Friel cichlid ﬁsh. BMC Evol. Biol. 9, 283–295 (2009). (CUMV), Martin Genner (U Bristol), Sylvain Piry (Paris), Lukas Ru¨ber (Bern Natural 58. Sambrook, J. & Russell, D. W. Molecular Cloning: A Laboratory Manual 3rd History Museum), Walter Salzburger (U Basel), Uli Schliewen (Zoologische Staats- edn (Cold Spring Harbor Laboratory Press, 2001). sammlung Munich), Erwin Schraml (Augsburg), Jos Snoeks (Africa Museum Tervuren), 59. Baird, N. A. et al. Rapid SNP discovery and genetic mapping using sequenced Melanie Stiassny (AMNH), George Turner (U Bangor), Sylvester Wandera (NAFIRI; RAD markers. PLoS ONE 3, e3376 (2008). Uganda) and Marco Welss (Kressberg). We also thank Matt McGee for help with 60. Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: DNA extraction, Keith Harshman of the Lausanne Genomic Technologies Facility and an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140 (2013). Cord Drogemu¨ller, Tosso Leeb, Muriel Fragnie`re and Miche`le Ackermann of the 61. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. NGS platform of the University of Bern for Illumina sequencing support, Aria Minder Methods 9, 357–U354 (2012). and Stefan Zoller of the Genetic Diversity Center (GDC) at ETH Zu¨rich and Irene Keller 62. McKenna, A. et al. The Genome Analysis Toolkit: A MapReduce framework for of the Interfaculty Bioinformatics Unit at the University of Bern for lab and bioinfor- analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 matics support and discussion, and Andy Cohen for comments on the paleogeographic (2010). maps. This research was supported by the Swiss National Science Foundation grant 63. Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, PDFMP3 134657 to O.S. and L.E. 2156–2158 (2011). 64. Kocher, T. D., Conroy, J. A., McKaye, K. R., Stauffer, J. R. & Lockwood, S. F. Author contributions Evolution of NADH dehydrogenase subunit 2 in east African cichlid ﬁsh. Mol. O.S., J.I.M. and C.E.W. designed the study; O.S. gathered and identiﬁed the cichlid Phylogenet. Evol. 4, 420–432 (1995). samples; J.I.M., S.M. and D.A.M. performed the lab work (DNA extraction, mtDNA and 65. Kocher, T. D. et al. Dynamics of mitochondrial DNA evolution in animals: LWS sequencing and RAD library preparation); J.I.M. conducted the analyses, with ampliﬁcation and sequencing with conserved primers. Proc. Natl Acad. Sci. assistance from O.S., C.E.W., D.A.M. and L.E.; J.I.M. prepared the manuscript together USA 86, 6196–6200 (1989). with O.S. and C.E.W. and L.E. and D.A.M. contributed to writing. 66. Hall, T. A. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acid Symp. 41, 95–98 (1999). Additional information 67. Stamatakis, A. RAxML-VI-HPC: maximum likelihood-based phylogenetic Supplementary Information accompanies this paper at http://www.nature.com/ analyses with thousands of taxa and mixed models. Bioinformatics 22, naturecommunications 2688–2690 (2006). 68. Stamatakis, A. et al. RAxML-Light: a tool for computing terabyte phylogenies. Competing ﬁnancial interests: The authors declare no competing ﬁnancial interests. Bioinformatics 28, 2064–2066 (2012). 69. Popescu, A. A., Huber, K. T. & Paradis, E. ape 3.0: new tools for distance-based Reprints and permission information is available online at http://npg.nature.com/ phylogenetics and evolutionary analysis in R. Bioinformatics 28, 1536–1537 reprintsandpermissions/ (2012). How to cite this article: Meier, J. I. et al. Ancient hybridization fuels rapid cichlid ﬁsh 70. Bouckaert, R. R. DensiTree: making sense of sets of phylogenetic trees. adaptive radiations. Nat. Commun. 8, 14363 doi: 10.1038/ncomms14363 (2017). Bioinformatics 26, 1372–1373 (2010). 71. Bouckaert, R. et al. BEAST 2: a software platform for Bayesian evolutionary Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in analysis. PLoS Comput. Biol. 10, e1003537 (2014). published maps and institutional afﬁliations. 72. Durand, E. Y., Patterson, N., Reich, D. & Slatkin, M. Testing for ancient admixture between closely related populations. Mol. Biol. Evol. 28, 2239–2252 (2011). This work is licensed under a Creative Commons Attribution 4.0 73. Reich, D., Thangaraj, K., Patterson, N., Price, A. L. & Singh, L. Reconstructing International License. The images or other third party material in this Indian population history. Nature 461, 489–494 (2009). article are included in the article’s Creative Commons license, unless indicated otherwise 74. Carleton, K. L. & Kocher, T. D. Cone opsin genes of African cichlid ﬁshes: in the credit line; if the material is not included under the Creative Commons license, tuning spectral sensitivity by differential gene expression. Mol. Biol. Evol. 18, users will need to obtain permission from the license holder to reproduce the material. 1540–1550 (2001). To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ 75. Excofﬁer, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567 (2010). r The Author(s) 2017 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 11 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nature Communications Springer Journals http://www.deepdyve.com/lp/springer-journals/ancient-hybridization-fuels-rapid-cichlid-fish-adaptive-radiations-AAzvql8Bcs

Loading next page...

References (96)

R. Bouckaert, Joseph Heled, D. Kühnert, T. Vaughan, Chieh-Hsi Wu, Dong Xie, M. Suchard, A. Rambaut, A. Drummond (2014)
BEAST 2: A Software Platform for Bayesian Evolutionary Analysis
PLoS Computational Biology, 10
(2002)
The effect of selection on a long wavelength-sensitive (LWS) opsin gene of Lake Victoria cichlid fishes
Proceedings of the National Academy of Sciences of the United States of America, 99
M. McGee, R. Neches, O. Seehausen (2016)
Evaluating genomic divergence and parallelism in replicate ecomorphs from young and old cichlid adaptive radiations
Molecular Ecology, 25
U. Schliewen, Barbara Klee (2004)
Reticulate sympatric speciation in Cameroonian crater lake cichlids
Frontiers in Zoology, 1
T. Hall (1999)
BIOEDIT: A USER-FRIENDLY BIOLOGICAL SEQUENCE ALIGNMENT EDITOR AND ANALYSIS PROGRAM FOR WINDOWS 95/98/ NT
C. Wagner, I. Keller, S. Wittwer, O. Selz, S. Mwaiko, L. Greuter, A. Sivasundar, O. Seehausen (2013)
Genome‐wide RAD sequence data provide unprecedented resolution of species boundaries and relationships in the Lake Victoria cichlid adaptive radiation
Molecular Ecology, 22
N. Patterson, Priya Moorjani, Yontao Luo, Swapan Mallick, N. Rohland, Yiping Zhan, Teri Genschoreck, Teresa Webster, D. Reich (2012)
Ancient Admixture in Human History
Genetics, 192
E. Durand, N. Patterson, D. Reich, M. Slatkin (2011)
Testing for ancient admixture between closely related populations.
Molecular biology and evolution, 28 8
C. Wagner, L. Harmon, O. Seehausen (2012)
Ecological opportunity and sexual selection together predict adaptive radiation
Nature, 487
R. Barrett, D. Schluter (2008)
Adaptation from standing genetic variation.
Trends in ecology & evolution, 23 1
J. Mallet (2007)
Hybrid speciation
Nature, 446
O. Seehausen, E. Koetsier, M. Schneider, L. Chapman, C. Chapman, M. Knight, G. Turner, J. Alphen, R. Bills (2003)
Nuclear markers reveal unexpected genetic variation and a Congolese-Nilotic origin of the Lake Victoria cichlid species flock
Proceedings of the Royal Society of London. Series B: Biological Sciences, 270
Kurt Haberyan, R. Hecky (1987)
The late Pleistocene and Holocene stratigraphy and paleolimnology of Lakes Kivu and Tanganyika
Palaeogeography, Palaeoclimatology, Palaeoecology, 61
M. Genner, B. Ngatunga, S. Mzighani, ALAN Smith, G. Turner (2015)
Geographical ancestry of Lake Malawi's cichlid fish diversity
Biology Letters, 11
Y. Terai, O. Seehausen, Takeshi Sasaki, Kazuhiko Takahashi, Shinji Mizoiri, T. Sugawara, Tetsu Sato, Masakatsu Watanabe, N. Konijnendijk, H. Mrosso, H. Tachida, H. Imai, Y. Shichida, N. Okada (2006)
Divergent Selection on Opsins Drives Incipient Speciation in Lake Victoria Cichlids
PLoS Biology, 4
C. Lindqvist, T. Motley, J. Jeffrey, V. Albert (2003)
Cladogenesis and reticulation in the Hawaiian endemic mints (Lamiaceae)
Cladistics, 19
Simon Martin, J. Davey, C. Jiggins (2014)
Evaluating the Use of ABBA–BABA Statistics to Locate Introgressed Loci
Molecular Biology and Evolution, 32
T. Johnson, K. Kelts, E. Odada (2000)
The Holocene History of Lake Victoria
, 29
S. Koblmüller, U. Schliewen, Nina Duftner, K. Sefc, Cyprian Katongo, C. Sturmbauer (2008)
Age and spread of the haplochromine cichlid fishes in Africa.
Molecular phylogenetics and evolution, 49 1
B. Lehner, K. Verdin, A. Jarvis (2008)
New Global Hydrography Derived From Spaceborne Elevation Data
Eos, Transactions American Geophysical Union, 89
Molly Schumer, Rongfeng Cui, G. Rosenthal, P. Andolfatto (2014)
Reproductive Isolation of Hybrid Populations Driven by Genetic Incompatibilities
PLoS Genetics, 11
Britta Meyer, A. Indermaur, Xenia Ehrensperger, B. Egger, G. Banyankimbona, J. Snoeks, W. Salzburger (2015)
Back to Tanganyika: a case of recent trans-species-flock dispersal in East African haplochromine cichlid fishes
Royal Society Open Science, 2
Thomas Kocher, Janet Conroy, K. McKaye, Jay Stauffer, Samuel Lockwood (1995)
Evolution of NADH dehydrogenase subunit 2 in east African cichlid fish.
Molecular phylogenetics and evolution, 4 4
Britta Meyer, Michael Matschiner, W. Salzburger (2016)
Disentangling Incomplete Lineage Sorting and Introgression to Refine Species‐Tree Estimates for Lake Tanganyika Cichlid Fishes
Systematic Biology, 66
M. Genner, O. Seehausen, D. Lunt, D. Joyce, P. Shaw, G. Carvalho, G. Turner (2007)
Age of cichlids: new dates for ancient lake fish radiations.
Molecular biology and evolution, 24 5
J. Feder, Xianfa Xie, J. Rull, S. Vélez, A. Forbes, B. Leung, Hattie Dambroski, Kenneth Filchak, M. Aluja (2005)
Mayr, Dobzhansky, and Bush and the complexities of sympatric speciation in Rhagoletis
Proceedings of the National Academy of Sciences of the United States of America, 102
C. Lindqvist, V. Albert (2002)
Origin of the Hawaiian endemic mints within North American Stachys (Lamiaceae).
American journal of botany, 89 10
Joana Meier, V. Sousa, D. Marques, O. Selz, C. Wagner, L. Excoffier, O. Seehausen (2017)
Demographic modelling with whole‐genome data reveals parallel origin of similar Pundamilia cichlid species after hybridization
Molecular Ecology, 26
J. Catchen, Paul Hohenlohe, S. Bassham, A. Amores, W. Cresko (2013)
Stacks: an analysis tool set for population genomics
Molecular Ecology, 22
Bin Ma, Ming Li, Louxin Zhang (2000)
From Gene Trees to Species Trees
SIAM J. Comput., 30
R. Stelkens, Corinne Schmid, O. Selz, O. Seehausen (2009)
Phenotypic novelty in experimental hybrids is predicted by the genetic distance between species of cichlid fish
BMC Evolutionary Biology, 9
S Carlquist (1974)
5279
Island Biology, Vol. 581
M. Arnold (1997)
Natural Hybridization and Evolution
Friedman M, Keck Bp, M. Friedman, T. Near, Benjamin Keck, A. Dornburg, R. Eytan, C. Martin, C. Hulsey, P. Wainwright (2013)
Molecular and fossil evidence place the origin of cichlid fishes long after Gondwanan rifting
Proceedings of the Royal Society B: Biological Sciences, 280
Ben Langmead, S. Salzberg (2012)
Fast gapped-read alignment with Bowtie 2
Nature Methods, 9
L. Excoffier, H. Lischer (2010)
Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows
Molecular Ecology Resources, 10
O. Seehausen, J. Alphen, F. Witte (1999)
Can ancient colour polymorphisms explain why some cichlid lineages speciate rapidly under disruptive sexual selection
Belgian Journal of Zoology, 129
I. Kozarewa, Z. Ning, M. Quail, M. Sanders, M. Berriman, D. Turner (2009)
Amplification-free Illumina sequencing-library preparation facilitates improved mapping and assembly of GC-biased genomes
Nature methods, 6
J. Hermansen, F. Haas, Cassandra Trier, R. Bailey, A. Nederbragt, A. Marzal, Glenn-Peter Sætre (2014)
Hybrid speciation through sorting of parental incompatibilities in Italian sparrows
Molecular Ecology, 23
A. McKenna, M. Hanna, E. Banks, A. Sivachenko, K. Cibulskis, A. Kernytsky, K. Garimella, D. Altshuler, S. Gabriel, M. Daly, M. DePristo (2010)
The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.
Genome research, 20 9
Deren Eaton, R. Ree (2013)
Inferring Phylogeny and Introgression using RADseq Data: An Example from Flowering Plants (Pedicularis: Orobanchaceae)
Systematic Biology, 62
D. Reich, K. Thangaraj, N. Patterson, A. Price, L. Singh (2009)
Reconstructing Indian Population History
Nature, 461
D. Damme, M. Pickford (2003)
The late Cenozoic Thiaridae (Mollusca, Gastropoda, Cerithioidea) of the Albertine Rift Valley (Uganda-Congo) and their bearing on the origin and evolution of the Tanganyikan thalassoid malacofauna
Hydrobiologia, 498
R. Stelkens, K. Young, O. Seehausen (2010)
THE ACCUMULATION OF REPRODUCTIVE INCOMPATIBILITIES IN AFRICAN CICHLID FISH
Evolution, 64
O. Selz, M. Pierotti, M. Maan, Corinne Schmid, O. Seehausen (2014)
Female preference for male color is necessary and sufficient for assortative mating in 2 cichlid sister species
Behavioral Ecology, 25
M. Barrier, B. Baldwin, R. Robichaux, M. Purugganan (1999)
Interspecific hybrid ancestry of a plant adaptive radiation: allopolyploidy of the Hawaiian silversword alliance (Asteraceae) inferred from floral homeotic gene duplications.
Molecular biology and evolution, 16 8
I. Keller, C. Wagner, L. Greuter, S. Mwaiko, O. Selz, A. Sivasundar, S. Wittwer, O. Seehausen (2012)
Population genomic signatures of divergent adaptation, gene flow and hybrid speciation in the rapid radiation of Lake Victoria cichlid fishes
Molecular Ecology, 22
É. Bezault, S. Mwaiko, O. Seehausen (2011)
POPULATION GENOMIC TESTS OF MODELS OF ADAPTIVE RADIATION IN LAKE VICTORIA REGION CICHLID FISH
Evolution, 65
C. Hermann, K. Sefc, S. Koblmüller (2011)
Ancient origin and recent divergence of a haplochromine cichlid lineage from isolated water bodies in the East African Rift system.
Journal of fish biology, 79 5
E. Verheyen, W. Salzburger, J. Snoeks, A. Meyer (2003)
Origin of the Superflock of Cichlid Fishes from Lake Victoria, East Africa
Science, 300
O. Seehausen (2004)
Hybridization and adaptive radiation.
Trends in ecology & evolution, 19 4
M. Lynch, J. Conery (2000)
The evolutionary fate and consequences of duplicate genes.
Science, 290 5494
F. Racimo, S. Sankararaman, R. Nielsen, Emilia Huerta-Sánchez (2015)
Evidence for archaic adaptive introgression in humans
Nature Reviews Genetics, 16
K. Carleton, T. Kocher (2001)
Cone opsin genes of african cichlid fishes: tuning spectral sensitivity by differential gene expression.
Molecular biology and evolution, 18 8
J. Sambrook, E. Fritsch, T. Maniatis (2001)
Molecular Cloning: A Laboratory Manual
B. Baldwin, W. Wagner (2010)
Hawaiian angiosperm radiations of North American origin.
Annals of botany, 105 6
W. Okullo, T. Ssenyonga, B. Hamre, O. Frette, K. Sørensen, J. Stamnes, A. Steigen, K. Stamnes (2007)
Parameterization of the inherent optical properties of Murchison Bay, Lake Victoria.
Applied optics, 46 36
Alexandros Stamatakis, A. Aberer, Christian Goll, Stephen Smith, S. Berger, F. Izquierdo-Carrasco (2012)
RAxML-Light: a tool for computing terabyte phylogenies
Bioinformatics, 28
W. Salzburger, S. Barić, C. Sturmbauer (2002)
Speciation via introgressive hybridization in East African cichlids?
Molecular Ecology, 11
(2015)
Evolution of Darwin’s finches and their beaks revealed by genome sequencing
Nature, 518
J. Weiss, F. Cotterill, U. Schliewen (2015)
Lake Tanganyika—A 'Melting Pot' of Ancient and Young Cichlid Lineages (Teleostei: Cichlidae)?
PLoS ONE, 10
O. Seehausen (2013)
Conditions when hybridization might predispose populations for adaptive radiation
Journal of Evolutionary Biology, 26
W. Salzburger, B. Bocxlaer, A. Cohen (2014)
Ecology and Evolution of the African Great Lakes and Their Faunas
Annual Review of Ecology, Evolution, and Systematics, 45
S. Yokoyama, F. Radlwimmer (1998)
The "five-sites" rule and the evolution of red and green color vision in mammals.
Molecular biology and evolution, 15 5
D. Joyce, D. Lunt, M. Genner, G. Turner, R. Bills, O. Seehausen (2011)
Repeated colonization and hybridization in Lake Malawi cichlids
Current Biology, 21
D. Brawand, C. Wagner, Yang Li, M. Malinsky, I. Keller, Shaohua Fan, O. Simakov, A. Ng, Z. Lim, É. Bezault, Jason Turner-Maier, Jeremy Johnson, Rosa Alcazar, H. Noh, P. Russell, Bronwen Aken, Jessica Alföldi, C. Amemiya, N. Azzouzi, J. Baroiller, F. Barloy-Hubler, A. Berlin, R. Bloomquist, K. Carleton, Matthew Conte, H. D'cotta, O. Eshel, L. Gaffney, F. Galibert, H. Gante, S. Gnerre, L. Greuter, R. Guyon, Natalie Haddad, W. Haerty, Rayna Harris, H. Hofmann, Thibaut Hourlier, G. Hulata, D. Jaffe, M. Lara, Alison Lee, Iain Maccallum, S. Mwaiko, M. Nikaido, Hidenori Nishihara, C. ozouF-cosTaz, D. Penman, Dariusz Przybylski, M. Rakotomanga, S. Renn, Filipe Ribeiro, M. Ron, W. Salzburger, L. Sánchez-Pulido, M. Santos, S. Searle, Ted Sharpe, Ross Swofford, F. Tan, Louise Williams, Sarah Young, Shuangye Yin, N. Okada, T. Kocher, E. Miska, E. Lander, B. Venkatesh, R. Fernald, A. Meyer, C. Ponting, J. Streelman, K. Lindblad-Toh, O. Seehausen, F. Palma (2014)
The genomic substrate for adaptive radiation in African cichlid fish
Nature, 513
G. Fryer (1975)
The Inland Waters of Tropical Africa. An Introduction to Tropical Limnology.
C. Scholz, T. Johnson, A. Cohen, J. King, J. Peck, J. Overpeck, M. Talbot, E. Brown, L. Kalindekafe, Philip Amoako, R. Lyons, T. Shanahan, I. Castañeda, C. Heil, S. Forman, L. McHargue, Kristina Beuning, J. Gómez, J. Pierson (2007)
East African megadroughts between 135 and 75 thousand years ago and bearing on early-modern human origins
Proceedings of the National Academy of Sciences, 104
S. Ho, G. Larson (2006)
Molecular clocks: when times are a-changin'.
Trends in genetics : TIG, 22 2
R. Green, J. Krause, Adrian Briggs, T. Maricic, U. Stenzel, Martin Kircher, N. Patterson, Heng Li, W. Zhai, Markus Fritz, Nancy Hansen, E. Durand, Anna-Sapfo Malaspinas, J. Jensen, Tomàs Marquès-Bonet, C. Alkan, Kay Prüfer, M. Meyer, H. Burbano, J. Good, R. Schultz, Ayinuer Aximu-Petri, Anne Butthof, Barbara Höber, Barbara Höffner, M. Siegemund, A. Weihmann, C. Nusbaum, E. Lander, C. Russ, Nathaniel Novod, J. Affourtit, M. Egholm, C. Verna, P. Rudan, D. Brajković, Ž. Kućan, I. Gusic, V. Doronichev, L. Golovanova, C. Lalueza-Fox, M. Rasilla, J. Fortea, A. Rosas, R. Schmitz, Philip Johnson, E. Eichler, D. Falush, E. Birney, J. Mullikin, M. Slatkin, R. Nielsen, J. Kelso, M. Lachmann, D. Reich, S. Pääbo (2010)
A Draft Sequence of the Neandertal Genome
Science, 328
A. Stamatakis (2006)
RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models
Bioinformatics, 22 21
Ryutaro Miyagi, Y. Terai, Mitsuto Aibara, T. Sugawara, H. Imai, H. Tachida, S. Mzighani, Takashi Okitsu, A. Wada, N. Okada (2012)
Correlation between nuptial colors and visual sensitivities tuned by opsins leads to species richness in sympatric Lake Victoria cichlid fishes.
Molecular biology and evolution, 29 11
M. Genner, G. Turner (2012)
Ancient hybridization and phenotypic novelty within Lake Malawi's cichlid fish radiation.
Molecular biology and evolution, 29 1
R. Abbott, D. Albach, S. Ansell, J. Arntzen, S. Baird, N. Bierne, J. Boughman, A. Brelsford, C. Buerkle, R. Buggs, R. Butlin, U. Dieckmann, F. Eroukhmanoff, A. Grill, S. Cahan, J. Hermansen, G. Hewitt, A. Hudson, C. Jiggins, Julia Jones, B. Keller, T. Marczewski, J. Mallet, P. Martínez-Rodríguez, M. Möst, S. Mullen, R. Nichols, A. Nolte, C. Parisod, K. Pfennig, A. Rice, M. Ritchie, Burkhardt Seifert, C. Smadja, R. Stelkens, J. Szymura, R. Väinölä, J. Wolf, D. Zinner (2013)
Hybridization and speciation
Journal of Evolutionary Biology, 26
R. Stelkens, Corinne Schmid, O. Seehausen (2015)
Hybrid Breakdown in Cichlid Fish
PLoS ONE, 10
L. Rieseberg (1997)
Hybrid Origins of Plant Species
Annual Review of Ecology, Evolution, and Systematics, 28
Margaret Blome, A. Cohen, C. Tryon, A. Brooks, J. Russell (2012)
The environmental context for the origins of modern human diversity: a synthesis of regional variability in African climate 150,000-30,000 years ago.
Journal of human evolution, 62 5
M. Maan, Kees Hofker, J. Alphen, O. Seehausen (2006)
Sensory Drive in Cichlid Speciation
The American Naturalist, 167
R. Bouckaert (2010)
DensiTree: making sense of sets of phylogenetic trees
Bioinformatics, 26 10
A. Cohen, J. Stone, Kristina Beuning, L. Park, P. Reinthal, D. Dettman, C. Scholz, T. Johnson, J. King, M. Talbot, E. Brown, S. Ivory (2007)
Ecological consequences of early Late Pleistocene megadroughts in tropical Africa
Proceedings of the National Academy of Sciences, 104
S. Ho, R. Lanfear, L. Bromham, M. Phillips, Julien Soubrier, A. Rodrigo, A. Cooper (2011)
Time‐dependent rates of molecular evolution
Molecular Ecology, 20
P. Danley, M. Husemann, Baoqing Ding, Lyndsay DiPietro, E. Beverly, D. Peppe (2012)
The Impact of the Geologic History and Paleoclimate on the Diversification of East African Cichlids
International Journal of Evolutionary Biology, 2012
T. Kocher, W. Thomas, A. Meyer, S. Edwards, S. Pääbo, F. Villablanca, A. Wilson (1989)
Dynamics of mitochondrial DNA evolution in animals: amplification and sequencing with conserved primers.
Proceedings of the National Academy of Sciences of the United States of America, 86 16
A. Hudson, P. Vonlanthen, O. Seehausen (2011)
Rapid parallel adaptive radiations from a single hybridogenic ancestral population
Proceedings of the Royal Society B: Biological Sciences, 278
P. Greenwood (1979)
Towards a phyletic classification of the ‘genus’ Haplochromis (Pisces, Cichlidae) and related taxa. Part 1
Bulletin of the British Museum of Natural History, 35
Carolina Pardo-Díaz, C. Salazar, C. Salazar, S. Baxter, C. Mérot, W. Figueiredo-Ready, M. Joron, W. McMillan, C. Jiggins, C. Jiggins (2012)
Adaptive Introgression across Species Boundaries in Heliconius Butterflies
PLoS Genetics, 8
Andrei-Alin Popescu, K. Huber, E. Paradis (2012)
ape 3.0: New tools for distance-based phylogenetics and evolutionary analysis in R
Bioinformatics, 28 11
John McCormack, James Maley, S. Hird, E. Derryberry, Gary Graves, R. Brumfield (2009)
Next-generation sequencing reveals phylogeographic structure and a species tree for recent bird divergences.
Molecular phylogenetics and evolution, 62 1
Benjamin Vernot, J. Akey (2014)
Resurrecting Surviving Neandertal Lineages from Modern Human Genomes
Science, 343
(2010)
Sequence analysis Advance Access publication June 7, 2011 The variant call format and VCFtools
E. Anderson, G. Stebbins (1954)
HYBRIDIZATION AS AN EVOLUTIONARY STIMULUS
Evolution, 8
O. Seehausen, Y. Terai, I. Magalhaes, K. Carleton, H. Mrosso, Ryutaro Miyagi, I. Sluijs, M. Schneider, M. Maan, H. Tachida, H. Imai, N. Okada (2008)
Speciation through sensory drive in cichlid fish
Nature, 455
C. Sturmbauer, S. Barić, W. Salzburger, L. Rüber, E. Verheyen (2001)
Lake level fluctuations synchronize genetic divergences of cichlid fishes in African lakes.
Molecular biology and evolution, 18 2
N. Baird, Paul Etter, T. Atwood, M. Currey, A. Shiver, Zachary Lewis, E. Selker, W. Cresko, Eric Johnson (2008)
Rapid SNP Discovery and Genetic Mapping Using Sequenced RAD Markers
PLoS ONE, 3
K. Dasmahapatra, James Walters, A. Briscoe, J. Davey, A. Whibley, Nicola Nadeau, A. Zimin, D. Hughes, Laura Ferguson, Simon Martin, C. Salazar, James Lewis, S. Adler, Seung-Joon Ahn, D. Baker, S. Baxter, Nicola Chamberlain, R. Chauhan, B. Counterman, T. Dalmay, L. Gilbert, K. Gordon, D. Heckel, H. Hines, K. Hoff, P. Holland, E. Jacquin‐Joly, F. Jiggins, Robert Jones, Durrell Kapan, P. Kersey, G. Lamas, D. Lawson, Daniel Mapleson, L. Maroja, Arnaud Martin, Simon Moxon, W. Palmer, R. Papa, A. Papanicolaou, Y. Pauchet, D. Ray, N. Rosser, S. Salzberg, Megan Supple, A. Surridge, Ayşe Tenger-Trolander, H. Vogel, P. Wilkinson, Derek Wilson, J. Yorke, Furong Yuan, Alexi Balmuth, Cathlene Eland, K. Gharbi, Marian Thomson, R. Gibbs, Yi Han, J. Jayaseelan, C. Kovar, T. Mathew, D. Muzny, Fiona Ongeri, Ling-Ling Pu, J. Qu, Rebecca Thornton, K. Worley, Yuanqing Wu, M. Linares, M. Blaxter, R. Constant, M. Joron, M. Kronforst, S. Mullen, R. Reed, S. Scherer, S. Richards, J. Mallet, W. McMillan, C. Jiggins (2012)
Butterfly genome reveals promiscuous exchange of mimicry adaptations among species
Nature, 487
O. Seehausen (2006)
African cichlid fish: a model system in adaptive radiation research
Proceedings of the Royal Society B: Biological Sciences, 273

Publisher: Springer Journals
Copyright: Copyright © 2017 by The Author(s)
Subject: Science, Humanities and Social Sciences, multidisciplinary; Science, Humanities and Social Sciences, multidisciplinary; Science, multidisciplinary
eISSN: 2041-1723
DOI: 10.1038/ncomms14363
Publisher site: See Article on Publisher Site

Abstract

ARTICLE Received 14 May 2016 | Accepted 20 Dec 2016 | Published 10 Feb 2017 DOI: 10.1038/ncomms14363 OPEN Ancient hybridization fuels rapid cichlid ﬁsh adaptive radiations 1,2,3 1,2,3 1,2 1,2,4 3,5 Joana I. Meier , David A. Marques , Salome Mwaiko , Catherine E. Wagner , Laurent Excofﬁer 1,2 & Ole Seehausen Understanding why some evolutionary lineages generate exceptionally high species diversity is an important goal in evolutionary biology. Haplochromine cichlid ﬁshes of Africa’s Lake Victoria region encompass 4700 diverse species that all evolved in the last 150,000 years. How this ‘Lake Victoria Region Superﬂock’ could evolve on such rapid timescales is an enduring question. Here, we demonstrate that hybridization between two divergent lineages facilitated this process by providing genetic variation that subsequently became recombined and sorted into many new species. Notably, the hybridization event generated exceptional allelic variation at an opsin gene known to be involved in adaptation and speciation. More generally, differentiation between new species is accentuated around variants that were ﬁxed differences between the parental lineages, and that now appear in many new combinations in the radiation species. We conclude that hybridization between divergent lineages, when coincident with ecological opportunity, may facilitate rapid and extensive adaptive radiation. 1 2 Aquatic Ecology and Evolution, Institute of Ecology and Evolution, University of Bern, 3012 Bern, Switzerland. Department of Fish Ecology and Evolution, Centre for Ecology, Evolution & Biogeochemistry, Eawag: Swiss Federal Institute of Aquatic Science and Technology, 6047 Kastanienbaum, Switzerland. 3 4 Computational and Molecular Population Genetics Lab, Institute of Ecology and Evolution, University of Bern, 3012 Bern, Switzerland. Biodiversity Institute & Department of Botany, University of Wyoming, Laramie Wyoming 82071, USA. Swiss Institute of Bioinformatics, 1015 Lausanne, Switzerland. Correspondence and requests for materials should be addressed to O.S. (email: [email protected]). NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 eterogeneity in diversiﬁcation rates among lineages is a of the major lakes of the region (Lakes Victoria, Edward, major factor shaping biodiversity. Yet, biological and Albert and Kivu). The largest of them is in Lake Victoria, Henvironmental factors underlying this variation are which has at least 500 endemic species that evolved in the 1 34–36 incompletely understood . Adaptive radiations are prime study past 15,000 years . Each radiation comprises enormous systems to learn about these factors as they are characterized by diversity in habitat occupation, trophic ecology, colouration a rapid origin of many species showing a diversity of ecological and behaviour. The high diversiﬁcation rate but also the high 34,37 adaptations. This process requires high levels of heritable nuclear genomic variation in the LVRS despite its young age variation in traits related to ecological and reproductive suggest that large amounts of standing genetic variation 2,3 isolation. However, adaptive radiations are often too rapid for must have been present at the onset of the radiation . the emergence of new relevant mutations between successive Previous work showing cytonuclear discordance in phylogenetic speciation events and are thus more likely to stem from reconstructions between LVRS and several riverine cichlid species standing variation. Hybridization between species can raised the possibility of ancient hybridization between divergent instantaneously boost genetic variation, which may facilitate species at the base of the radiation but could not demonstrate it . 3–9 speciation and adaptive radiation . Hybridization on secondary contact is not unlikely, as allopatric Hybridization among members of an adaptive radiation, has cichlid species, divergent by even millions of years, readily been suggested to potentially facilitate speciation events, an idea produce fertile offspring in the lab . known as the ‘syngameon hypothesis’ . Introgression of traits Using genomic data from riverine haplochromine cichlids involved in adaptation or reproductive isolation has been sampled from all major African drainage systems, and demonstrated among members of several adaptive radiations representative species from all lineages within the Lake Victoria (for example, traits related to host shift in Rhagoletis fruit ﬂies , region, we demonstrate here that the LVRS evolved from a hybrid 11,12 wing patterns in Heliconius butterﬂies or beak shape in swarm. All lake radiations show very similar proportions of Darwin’s ﬁnches ). In other radiations, the hybrid ancestry of mixed ancestry derived from two distantly related haplochromine some species has been inferred, but a direct link between lineages that had evolved in isolation from one another in introgressed traits and speciation awaits further testing (for different river systems for more than a million years before 14–16 17 example, cichlid ﬁshes of Lakes Tanganyika , Malawi , hybridizing in the Lake Victoria region. We ﬁnd evidence that 18,19 20 Victoria and Barombi Mbo ). this hybridization event facilitated subsequent adaptive radiation Another hypothesis for a perhaps more fundamental role of by providing genetic variation that has been recombined hybridization in adaptive radiation, distinct from the ‘syngameon and sorted into many new species. Variants that were ﬁxed hypothesis’, is the idea that hybridization between distinct between the parental lineages show accentuated differentiation lineages may seed the onset of an entire adaptive radiation . between young Lake Victoria species, but appear in many Such hybridization can be common when allopatric lineages new combinations in the different species. Notably, each of the come into secondary contact , and selection against hybrids may two major allele classes of an opsin gene involved in adaptation 39,40 be weak during colonization of new environments. In this and speciation among Lake Victoria cichlids is likely derived situation, the formation of a hybrid swarm, if coincident with from one of the two parental lineages. This indicates that a major ecological opportunity, may accelerate adaptive radiation by (i) part of the variation at this gene segregating in the LVRS stems providing functional genetic variation that can recombine into from hybridization between these lineages. Our results suggest novel trait combinations favoured by selection and mate choice, that hybridization between relatively distantly related species, and (ii) breaking genetic correlations that constrained the when coincident with ecological opportunity, may facilitate rapid evolvability of parental lineages . In addition, hybridization may adaptive radiation. Thus, hybridization, even in the distant facilitate speciation when multiple ﬁxed differences that confer past, may have important implications for understanding reproductive isolation between the two parental species decouple variation in extant species richness between lineages as well as and segregate in a hybrid swarm, such that selection against variation in recent rates of diversiﬁcation. incompatible gene combinations can generate more than two new 21–23 reproductively isolated species . This ‘hybrid swarm origin of adaptive radiation’ hypothesis has Results been more challenging to test. So far the only adaptive radiation Identifying the closest relatives of the LVRS. To identify the for which a hybrid origin has been robustly demonstrated is the closest relatives of the LVRS, we use comprehensive sampling of Hawaiian silverswords, which have radiated from an allopoly- haplochromine cichlids from the LVRS and from all major river ploid hybrid population between two North American tarweed systems harbouring haplochromines (Fig. 1, Supplementary species . Because gene and genome duplication are also Data 1). Maximum likelihood phylogenetic trees were recon- proposed to facilitate adaptive radiation , it is difﬁcult to structed from 3.15 Mb of concatenated restriction site associated distinguish between effects of hybridization per se and those of DNA (RAD) sequences (436,166 SNPs, 8.1% missing data, gene or genome duplication in this case. Evidence consistent with Supplementary Data 2), and from two mitochondrial markers a hybrid swarm origin of entire radiations has also been found in (1,897 bp). In the nuclear phylogeny, the entire LVRS forms Alpine whiteﬁsh , the ‘mbuna’ group of the Lake Malawi cichlid a well-supported clade that also includes cichlids from the 27 28,29 ﬁsh radiation , and allopolyploid Hawaiian endemic mints lower Nile in Egypt (Haplochromis spp. Egypt) and Haplochromis and possibly other polyploid plant radiations on Hawaii . sp. ‘Nyangara’ from the Rusizi River, the outﬂow of Lake However, it remains to be tested if, in these systems, Kivu and an inﬂow to Lake Tanganyika (orange labels in Fig. 1, hybridization occurred before or after the radiation had started, Supplementary Fig. 1, Supplementary Discussion). Astatotilapia and if hybridization-derived polymorphisms played a role in stappersi from the Kalambo River (an inﬂow to Lake Tanganyika, speciation and adaptive diversiﬁcation. Congo drainage, referred to as Haplochromis sp. ‘Chipwa’ in The Lake Victoria Region Superﬂock of cichlid ﬁsh (LVRS) Meyer et al. ), and an undescribed Astatotilapia species from is a group of 700 haplochromine cichlid species endemic to the middle Congo (A. sp. ‘Yaekama’, DRC), together hereafter the region around Lake Victoria and nearby western rift lakes in called the ‘Congolese lineage’ (red labels), form the superﬂock’s East Africa that started diversifying about 100–200 thousand closest relatives. The sister group to the LVRS plus the 31–34 years ago . It includes several adaptive radiations, one in each Congolese lineage taxa is a clade including lineages in 2 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE L. Albert ab L. Edward Mpanga R. / L. Saka L. Kivu Bootstrap values L. Edward L. Kivu >90 LVRS L. Kivu >80 L. Edward LVRS L. Victoria / Kagera L. L. Edward Haplochromis spp. Egypt Haplochromis sp. "Nyangara" Rusizi R. Astatotilapia stappersi Kalambo R. Congolese Astatotilapia sp. "Yaekama" Lineage Upper Nile Astatotilapia bloyeti Mkuzi R. / Ruaha Astatotilapia sparsidens L. Manyara Eastern Astatotilapia cf. paludinosa Malagarasi R. Astatotilapia sp. Wami R. "Haplochromis" gracilior L. Kivu Upper Nile Thoracochromis pharyngalis L. Edward Lineage Astatotilapia flavijosephi L. Kinneret Astatotilapia desfontainii Tunisia Rhamphochromis sp. Astatotilapia calliptera Lake Malawi Metriaclima zebra Haplochromines Protomelas spp. Congolese Haplochromis sp. Ruvu R. Astatotilapia cf. bloyeti Wami R. / kilossana Sonjo R. Astatotilapia sp. Rovuma Astatotilapia tweddlei L. Chilwa Astatotilapia burtoni Kalambo R. / L. Tanganyika Astatoreochromis straeleni Rusizi R. Astatoreochromis alluaudi L. Victoria / Saka Lobochilotes labiatus Lake Tanganyika Petrochromis sp. "kazumbe" Tropheini Ctenochromis horei Drainage systems Orthochromis sp. "red cheek red tail" L. Mweru Orthochromis polyacanthus L. Mweru East coast Lake Mweru Pseudocrenilabrus spp. L. Mweru Pseudocrenilabrus Haplochromis demeusii Congo Congo Serranochromis robustus L. Mweru Nile basin Lake Mweru Sargochromis spp. L. Mweru Sargochromis Haplochromis stigmatogenys Kasai/Luvua North Interior Cyclopharynx schwetzi L. Fwa Lake Fwa Endorheic drainage systems Thoracochromis callichromis L. Fwa Zambezi Ctenochromis pectoralis Pangani Figure 1 | Phylogenetic context of the Lake Victoria Region cichlid radiation. (a) Maximum likelihood phylogeny built from concatenated RADtag sequences of Lake Victoria Region Superﬂock (LVRS) cichlids and relatives including all known lineages of haplochromine cichlids (n ¼ 156). Radiations are indicated as grey triangles in the phylogenetic tree and multiple samples of a lineage are visually collapsed to a single terminal branch (full treein Supplementary Fig. 1). Members of the LVRS (including Haplochromis sp. ‘Nyangara’ from the Rusizi River and Haplochromis spp. Egypt, see Supplementary Discussion) are indicated with orange stars both in the tree and in the sampling map (b) and are labelled by lake (L) or river (R) they were sampled in. ‘Congolese lineage’ LVRS relatives are highlighted with red triangles, members of the ‘Upper Nile lineage’ with blue triangles, those from Eastern rivers with dark blue squares, and all other more distantly related lineages with black circles. (b) Sampling map. River drainage systems that we sampled are shown as coloured polygons. The radiation ancestor’s closest living relatives are shown in images: ‘Haplochromis’ gracilior/Thoracochromis pharyngalis from the Upper Nile lineage and Astatotilapia sp. ‘Yaekama’/A. stappersi from the Congolese lineage. The Lake Victoria cichlids shown in the grey triangle on the right represent some of the many and varied species that arose from the hybrid swarm (Photo credits: Ole Seehausen, Salome Mwaiko, Frans Witte, ‘Teleos’, Uli Schliewen, Adrian Indermaur, Oliver Selz; map adapted from http://www.worldwildlife.org/hydrosheds ). endorheic (that is, closed) drainage systems of the Eastern the LVRS, from which the radiations in the others lakes are Rift and in Indian Ocean drainage systems, east of the derived, and that connections between Lakes Edward and Kivu Lake Victoria region (‘Eastern taxa’, dark blue), as well as existed until recently. ‘Haplochromis’ gracilior from Lake Kivu and Thoracochromis pharyngalis from Lake Edward (light blue, Supplementary Discussion). Although the latter two species are sympatric Tests of hybridization. Phylogenetic reconstruction based on with LVRS members, they are clearly not part of the radiations concatenated sequence data assumes a single history across the themselves. We refer to H. gracilior and T. pharyngalis as genome, although recombination, incomplete lineage sorting, and the ‘Upper Nile lineage’ because all known members are conﬁned introgressive hybridisation can cause extensive genome-wide to the region that was historically the uppermost Nile variation in genealogy . To account for such genealogical 37,42 33,43 45 drainage . Consistent with previous publications ,we variation, we reconstructed a SNP based species tree . estimate the split of this lineage from the Congolese lineage The resulting tree supports the sister relationship of the LVRS to date to B1.6–5.8 million years ago, coincident with the with the Congolese taxa, but reveals incongruence between uplift of the ﬂanks of the Western branch of the East African this group and the Eastern clades and Upper Nile clade Rift truncating the paleo-river network and isolating the (Supplementary Fig. 2). To test for genetic admixture, we LVR (the new Upper Nile) from the Congo (Supplementary computed Patterson’s D statistics (ABBA-BABA test) .We Table 1; Supplementary Fig. 3). found that all LVRS radiations show strong signals of admixture Our genome-wide nuclear sampling provides unprecedented with the two Upper Nile species, H. gracilior and T. pharyngalis resolution of the relationships between radiations within (Fig. 2, Supplementary Table 2, Supplementary Fig. 4). Excess the LVRS (Fig. 1, Supplementary Fig. 1). The radiations in allele sharing between the LVRS and either species is equally Lake Victoria (including its major tributary, the Kagera River), strong (Supplementary Table 3, rows 2.1-2.4, Supplementary in Crater Lake Saka and associated Mpanga River, and in Lake Fig. 5). Importantly, neither H. gracilior nor T. pharyngalis Albert, each form monophyletic clades. In contrast, the show greater admixture with the LVRS species with which Lake Edward members of the LVRS are a paraphyletic group they live in sympatry than with allopatric members of the that includes taxa basal to the radiations of Lakes Victoria and LVRS from other lakes that have been geographically isolated for Saka, and others that form a clade together with species from many thousands of years (Fig. 2 and Supplementary Fig. 5). Lake Kivu. These ﬁndings are consistent with the hypothesis To test the directionality of gene ﬂow, we applied an extended that Lake Edward constitutes the oldest extant radiation within version of the partitioned D statistic (5 population test, NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 a b c L. Albert L. Edward L. Kivu 0.4 LVRS (P1) L. Saka g g t g t L. Victoria t g Kagera L. 0.2 t g g t P1 P2 P3 Outgroup A. stappersi Congo (P2) B A B A A. sp. "Yaekama" w w A. bloyeti l (b) ww 0.0 p p p p b s s s b s b b s b A. sparsidens (s) Eastern (P3) A. cf. paludinosa (p) −0.2 A. sp. "Wami" (w) T. pharyngalis (t) Upper Nile (P3) P1 P2 P3 Outgroup H. gracilior (g) Albert Edward Kivu Saka Victoria Kagera A B B A A. flavijosephi Outgroup P1 Figure 2 | Evidence for Congo-Nilotic hybridization in the ancestry of the LVRS. (a) Schematic genealogy with taxa used for D statistics (n ¼ 73 individuals, see Supplementary Data 2). Abbreviations used in other panels are given in parentheses and the color scheme is the same as in Figure 1. The inferred gene ﬂow edge is shown with an arrow (Note the directionality of gene ﬂow is inferred with the ﬁve-population test not shown in this ﬁgure). (b) D statistics to test for potential gene ﬂow between each Eastern and Upper Nile taxon (P3) separately (abbreviations given in a) and cichlids from each LVR lake radiation (P1) or the Congolese taxon A. stappersi (P2). Vertical bars correspond to three standard errors. Positive D values indicate gene ﬂow between P1 (LVR lake radiation) and P3 (Eastern or Upper Nile taxon), whereas negative D values indicate gene ﬂow between P2 (A. stappersi) and P3 (Eastern or Upper Nile taxon) as illustrated in (c). Exact values and more test results are given in Supplementary Table 2. ab 0.8 0.6 0.4 20 0.2 0.0 Conspe- Sister LVi-LVi LKi-LKi LVi-LKi cifics species Size of ancestry blocks (kb) Figure 3 | Congo-Nilotic ancestry blocks in LVRS genomes. (a) The size distribution of putative ancestry blocks shows mostly small ancestry blocks and slightly larger Congolese (red) than Upper Nile blocks (blue). The plot shows the counts of ancestry blocks in different size categories summed up for ﬁve LVRS radiation species across all scaffolds calculated with 3 kb windows. As most blocks do not span multiple windows of 3 kb and many blocks cannot be clearly allocated to Congolese or Upper Nile ancestry (grey) (Supplementary Fig. 8), the average ancestry block size is likely 3 kb or smaller, consistent with hybridization many thousands of years ago. (b) Correlation of ancestry blocks between whole-genome sequenced LVRS members is high overall but decreases with phylogenetic distance. The boxplots show correlation of f (ref. 78) in 10 kb windows between single individuals of conspeciﬁcs (Pundamilia individuals of the same species), sister species (Pu. pundamilia versus Pu. nyererei), more distantly related Lake Victoria (LVi) species (Paralabidochromis ﬂavus versus Pu. pundamilia and Pu. nyererei), Lake Kivu (LKi) species (Pa. paucidens versus Harpagochromis vittatus) and Lake Victoria against Lake Kivu species. This suggests that all radiation member species share the same hybridization event in their ancient history but vary in how long after that event they remained part of the same recombining population. It also suggests that some of the admixture variation still segregates within individual species (indicated by the deviation from an f correlation of 1 among conspeciﬁcs). Supplementary Table 3). We ﬁnd evidence for gene ﬂow from Importance of the ancient hybridization event. Many pheno- H. gracilior (g in Fig. 2b) and from T. pharyngalis (t in Fig. 2b) typic traits known to contribute to ecological adaptation and into each LVRS lake radiation, but no evidence for gene ﬂow in reproductive isolation have diverged in multiple speciation events the opposite direction (Supplementary Table 3, rows 1.1–1.4). in each of the LVRS radiations (for example, tooth shapes, male 48 39,40 Using F4-ratio tests , we estimated the Upper Nile ancestry nuptial colouration, opsin alleles) . Some of these traits are proportion to be B20% in all LVRS radiations (Supplementary also divergent between the Congolese Astatotilapia and Upper 37,42 Fig. 6 and Supplementary Table 4). Whole-genome sequencing Nile Thoracochromis species . Intrinsic incompatibilities of Lake Victoria and Lake Kivu radiation species, and have also been found in experimental hybrid crosses of A. stappersi, T. pharyngalis and H. gracilior, conﬁrms excess similarly divergent Astatotilapia species . In addition, some of allele sharing between LVRS members and the Upper Nile the allelic variation in the lake radiations has been shown to lineage, and reveals ancestry blocks of B3 kb, consistent with predate the origin of the LVRS (refs 49,50), perhaps most notably evidence that the admixture event occurred many thousands of that of the long-wavelength sensitive (LWS) opsin gene, one of generations ago in the common ancestor of all superﬂock the best-studied genes in Lake Victoria cichlids . This gene codes radiations in the different lakes (Fig. 3a, Supplementary Fig. 8). for the protein moiety of red-sensitive visual pigments in retinal The genome-wide signatures of admixture are correlated cones and is exceptionally diverse in Lake Victoria Region among different species from Lake Victoria and Lake cichlids . It plays a crucial role in adaptation on the Kivu (Fig. 3b), in line with a shared admixture event at the steep ambient light gradients associated with water depth and onset of the LVRS. turbidity gradients that are characteristic for the lakes in the Lake 4 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications >50 Counts (in thousands) D(BABA–ABBA)/(BABA+ABBA) f correlation d NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE Victoria region. More red-shifted LWS opsin variants are outside the radiation, LV outliers were not more often beneﬁcial in deep and murky water where the light spectrum is differentially ﬁxed in the control groups than non-outliers among relatively more red-shifted because of particulate matter category 4 SNPs (Supplementary Table 5, Supplementary 39,52 absorbing and scattering light of shorter wave lengths . This Discussion). The mosaic-like ancestry pattern is also corroborated opsin gene likely also plays a role in behavioural reproductive by our whole-genome sequence data (Supplementary Fig. 8). isolation because divergence in colour perception between species with different LWS opsin genotypes is often associated 39,49,54 with divergent male breeding colouration , which is an Discussion important mate choice cue . Here, we demonstrate that multiple large animal adaptive The LWS opsin gene is highly diverse in the LVRS with radiations arose from a hybrid swarm. Our results suggest two deeply divergent haplotype clades, often referred to as ancient admixture between two lineages of haplochromines, haplotype classes because the alleles within each clade are one from the Upper Congo and one from the Upper Nile functionally similar to each other . A substitution at amino acid drainage, at the origin of the entire Lake Victoria Region position 177 in the cichlid LWS gene shifts peak absorbance Superﬂock of cichlid ﬁsh. Contemporary Upper Nile lineage towards longer wavelength (red) in haplotype class I relative to introgression or earlier independent introgression events into class II by 7 nm (refs 49,51) and also other substitutions inﬂuence each lake radiation are unlikely given (i) the highly supported spectral sensitivity (Supplementary Discussion, Supplementary genomic monophyly of the entire Lake Victoria Region Data 3). The two haplotype classes are often associated Superﬂock (Fig. 1, Supplementary Fig. 1), (ii) similar Congo with different light environments, class I with shallow and and Upper Nile lineage ancestry proportions in all lake radiations clear water, and class II with deep and turbid water. For example, (Supplementary Fig. 6) and (iii) highly correlated Congolese and rocky shore algae scrapers of the genera Neochromis and Mbipia Upper Nile lineage ancestry tracts (Fig. 3). In addition, the fact and mud bottom detritivores of the genus Enterochromis, that neither H. gracilior (Lake Kivu) nor T. pharyngalis (Lake likely representing an early divergence event in the Lake Edward) show greater admixture into LVRS species that are Victoria radiation, occupy the opposite ends of the light sympatric with them than with allopatric members of the LVRS spectrum in Lake Victoria: clear and shallow versus murky and from other lakes (Fig. 2 and Supplementary Fig. 5) speaks against deep. All 438 LWS haplotypes previously sequenced from recent introgression. The radiations in Lakes Albert, Saka and 39,49 Neochromis and Mbipia algae scrapers are part of Victoria have been geographically separated from Lakes Edward the haplotype class I, whereas 11 of the 12 Enterochromis and Kivu, and hence from T. pharyngalis and H. gracilior, for 31,33,34 LWS haplotypes belong to class II (ref. 39). Within ecologically between 4,000 (Saka) and 100,000 (Victoria) years and variable lineages, young sister species occupying shallow versus show similar ancestry proportions. deep or clear versus turbid waters, often have different Highly similar amounts of allele sharing between the LVRS frequencies of these two LWS haplotype classes or species and both H. gracilior and T. pharyngalis (Supplementary recombinants between them, and LWS variation is consistently Table 3, rows 2.1–2.4, Supplementary Fig. 3) suggests that the associated with divergence in male nuptial coloration (Fig. 4, taxon that hybridized with the Congolese ancestor of the LVRS Supplementary Discussion). We sequenced the LWS opsin was a close relative of both Upper Nile species, or their ancestor. gene in Congolese and Upper Nile taxa, and found that the It is likely that the Congolese lineage colonized the Lake Victoria two haplotype classes in the LVRS are each shared exclusively region through capture of Malagarasi (Congo) tributaries during with just one of these parental lineages (Fig. 4, Supplementary a humid phase (for example, 145,000–120,000 years ago, Fig. 7, Supplementary Data 3, Supplementary Discussion). Supplementary Fig. 3) and encountered representatives of the All haplotypes in individuals from the Congolese lineage Upper Nile lineage that would by then have occupied the Lake take basal positions in the class I clade, whereas the Victoria region. The existence of large lakes in the Lake Victoria Upper Nile lineage haplotypes take basal positions in the class region at that time is likely to have provided ecological II clade. opportunity for diversiﬁcation, which could be exploited by To see if similar patterns of sorting of admixture variation a genetically diverse hybrid swarm. could be detected elsewhere in the genomes, we tested if The ﬁnding that the two LWS opsin haplotype classes of genomic SNPs that are highly differentiated among LVRS species the LVRS are each shared uniquely with one of the parental were commonly derived from admixture between the lineages (Fig. 4, Supplementary Fig. 7) suggests that the ancient two parental lineages. We inferred the origin of bi-allelic Congo-Nilotic admixture event was the source of the SNPs in six phenotypically diverse Lake Victoria species, the high functional variation among LVRS cichlids at the LWS opsin 39,51 insectivore Pundamilia pundamilia, the insectivore/zooplankti- gene . Recombination between the ancestral haplotypes within vore P. nyererei, the algae grazer Neochromis omnicaeruleus, the admixed populations further enhanced the functional the paedophage Lipochromis melanopterus, the insectivore variation of LWS haplotypes available for the LVRS radiation Paralabidochromis chilotes and the piscivore Harpagochromis (Supplementary Discussion, Supplementary Data 3). This gene is cf. serranus (Fig. 5). We found that at 30% of sites that particularly important for adaptation to the red-shifted end of are exceptionally strongly differentiated between Lake ambient aquatic light spectra that characterize the waters of 39,51 Victoria species (LV outliers), one of the alleles was indeed Lake Victoria . Bringing together the Congolese and Upper likely introduced into the ancestry of the radiation through Upper Nile opsin haplotypes thus appears to have facilitated adaptation Nile lineage introgression (categories 3 þ 4 in Fig. 5a, to an extreme range of light conditions and visual ecologies, with Supplementary Discussion). Sites at which the Congolese many transitions between them at both early as well as late stages and Upper Nile lineage taxa are ﬁxed for alternative alleles of the radiation. Divergence into major ecologically differentiated (category 4 in Fig. 5a) are enriched for LV outliers and show clades (Enterochromis vs Neochromis/Mbipia), but also major a pattern consistent with mosaic-like sorting of ancestral variants habitat shifts between young sister species are associated amongst Lake Victoria species (Fig. 5b). The enrichment for with recruitment of Congolese versus Upper Nile lineage LV outliers in category 4 SNPs does not seem to be because of derived haplotypes at the LWS locus (Supplementary Discussion). inherently increased ﬁxation probability of these loci, as in Our data suggest a much more general impact of the ancient pairwise comparisons of six control group cichlid species from admixture event on evolution in this set of young adaptive NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 Clear shallow water Shallow water algivores P. macrocephala 'blue' Class I Congolese lineage Deep water Turbid deep water P. macrocephala 'yellow' detritivores Class II Upper Nile lineage Figure 4 | High LWS opsin diversity likely because of the ancient hybridization event. The two major LWS opsin haplotype classes, I and II, in the LVRS (orange) are each shared exclusively with either the Congolese (A. stappersi and A. sp. ‘Yaekama’, red) or the Upper Nile lineage (H. gracilior and T. pharyngalis, blue), respectively (details in Supplementary Fig. 7; Supplementary Data 3). LWS haplotype class I is generally associated with cichlids living in shallow and clear water habitats, whereas class II is associated with deeper and more turbid habitats. Speciation by divergence in habitat type seems to have been accompanied by ﬁxation of alternative LWS haplotypes both at early and at late stages of the adaptive radiation. This is exempliﬁed by near ﬁxation of alternative haplotype classes between ecologically divergent genera such as shallow water rocky shore algae scrapers of Neochromis and Mbipia versus mud bottom detritivores of the genus Enterochromis, and by the young incipient species pair of Pundamilia macrocephala ‘blue’ (living very shallow) and ‘yellow’ (living deeper) which have predominantly alleles of haplotype class I and II, respectively (Photo credits: Ole Seehausen, Adrian Indermaur, ‘Teleos’, Oliver Selz, Uli Schliewen). radiations. Experimental crosses of cichlids have shown correlation of Congo-Nilotic ancestry block patterns decreases 38,56 57 that intrinsic incompatibilities and phenotypic novelty between species with increasing phylogenetic distance within the both increase with genetic distance between the crossed species. LVRS (Fig. 3b), it is likely that Congolese and Upper Nile lineage Both kinds of variation could be important in speciation. alleles have been segregating in the LVRS for a long time and that The divergence time of the Congolese and Upper Nile lineage they were progressively sorted among the radiation members taxa, interpreted in the context of this previous experimental during species diversiﬁcation. work, lets us predict that a hybrid swarm between these lineages The combined evidence from the LWS opsin gene and would contain both intrinsic incompatibilities and transgressive the genome-wide patterns of sorting of ancestral variation trait variation. In agreement with this prediction, we found is consistent with an important role of the ancestral hybridization that genomic sites with alternative alleles ﬁxed in the two for the evolution of the Lake Victoria Region Superﬂock. parental lineages are enriched for outlier loci likely involved in The resulting large genetic variation may help to explain how divergent adaptation and species differentiation in Lake Victoria multiple adaptive radiations in the Lake Victoria Region, together (Fig. 5). This suggests that speciation has commonly been forming 700 species, could arise in just 100–200 thousand 31,32 associated with sorting of alleles brought together in the radiation years , including 500 þ genetically and phenotypically well- ancestor by the admixture event. Sites showing negative epistatic differentiated species endemic to Lake Victoria which evolved in 32,35 interactions (Bateson Dobzhansky Muller (BDM) probably just 15,000 years . Hence, our data provide evidence incompatibilities) would be expected to be ﬁxed for alternative for the hypothesis that hybridization between divergent lineages, alleles in the parental lineages and thus fall in ancestry category 4 when coincident with novel ecological opportunity such as in Figure 5a. The strong enrichment of LV outliers in this colonisation of newly formed lakes, may facilitate rapid adaptive category is in line with differential sorting of BDM radiation through recombination and sorting of admixture- incompatibilities among species that emerged from the hybrid derived polymorphisms by natural and sexual selection. 21–23 swarm , which may have facilitated reproductive isolation That hybridization may fuel entire adaptive radiations has been among the members of the adaptive radiation. Given that the hypothesized for several systems, but robust evidence has so far 6 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NA NA NA NA NA U472 U38 LG23 LG22 LG15 LG14 LG11 LG10 LG9 LG7 LG4 LG2 NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE P<0.0001 P=0.028 P=0.0003 0.30 10,799 1,053 1,038 100 179 0.25 20% 0.20 0.15 10.1% 9.7% 0.10 7.9% 0.05 1.4% 0.00 1. 2. 3. 4. 5. Only one Both Both alleles Different Ancestral allele in alleles only if C + N alleles in LV maf C + N in C combined C and N 14-18% combined Upper Nile (N) b Congolese (C) Frequency of the Congolese allele 0.0 0.2 0.4 0.6 0.8 1.0 Pu. pundamilia N. omnicaeruleus Pu. nyererei L. melanopterus Pa. chilotes H. cf. serranus Figure 5 | Differential sorting of hybridization-derived variation in the LVRS. (a) Sites ﬁxed for alternative alleles in the Congolese (C) and Upper Nile (N) taxa are enriched for high global F outlier SNPs in Lake Victoria. Of the 12,890 biallelic SNPs among six sympatric Lake Victoria species (shown in ST (b)), 340 are outliers of high global F (LV outliers). We assigned SNPs to ﬁve different ancestry categories according to the presence or absence of the ST two alleles in the Congolese (C) and Upper Nile (N) lineage taxa. The grey bars show the proportion of LV outliers among all SNPs in each ancestry category. Total SNP counts in each category and P-values of two-sided Fisher’s exact tests are shown on top. Ancestry category (1) includes all SNPs for which only one of the two LV alleles was found in the Congolese and Upper Nile taxa together (novel LV allele or unsampled in parental lineages), (2) both LV alleles found in the Congolese taxa (polymorphic in LVRS even without Upper Nile hybridization), (3) only one allele found in Congolese but both alleles found in Upper Nile taxa (not available in LVRS without hybridization), and 4) Congolese and Upper Nile taxa each ﬁxed for alternative LV alleles (not available in LVRS without hybridization) potentially including Bateson–Dobzhansky–Muller incompatibilities. Category 5 includes sites with similar initial allele frequency in Lake Victoria (16%) than sites ﬁxed for alternative alleles in the parental lineages (category 4) to test if the enrichmentin category 4 could simply be because of high initial allele frequency. (b) Differential sorting of parental alleles between Lake Victoria cichlid species at LV outliers ﬁxed for alternative alleles in the Congolese and Upper Nile lineage taxa (mean global F among LV species ¼ 0.52). Each square represents ST a SNP coloured according to the allele frequency in that species ranging from red (ﬁxed for Congolese allele) to blue (ﬁxed for Upper Nile allele). All except two sites (2 þ 3 from the right) are located on different scaffolds of the Pundamilia nyererei reference genome. If known, chromosomal positions on the Oreochromis niloticus genome are shown below. At least 10 of the 22 chromosomes are involved in mosaic-like allele sorting between radiation species at loci that were ﬁxed for alternative alleles in the parental lineages of the ancestral hybrid swarm (Photo credits: Oliver Selz, Ole Seehausen, Adrian Indermaur, ‘Teleos’, Uli Schliewen). been conﬁned to the allopolyploid Hawaiian silverswords . performed restriction associated DNA (RAD) sequencing with haplochromine cichlid species from all major drainage systems that are either currently Here, we report strong evidence that hybrid ancestry connected to the Lake Victoria region, or have been in the past, and representatives fuelled diversiﬁcation in several large animal adaptive of all major radiations of the Lake Victoria Region Superﬂock. The closest relatives radiations. Thanks to the advent of high-throughput were inferred with maximum likelihood trees based on these data and on mito- sequencing, our power to test for ancient hybridization and to chondrial sequences of the same ﬁsh. D statistics, ﬁve population tests, and F4-ratio tests, were used to identify hybridization between lineages, determine the direction study its impact on subsequent evolution has increased of gene ﬂow and quantify ancestry proportions, respectively. Whole-genome enormously. Future studies will reveal if hybridization in the sequencing was then performed with ﬁsh of the lineages identiﬁed as ancestral to ancestry of major adaptive radiations is widespread, and whether the LVRS plus representatives of the superﬂock in order to corroborate the sig- its occurrence may explain some of the observed large variation in natures of hybridization and to estimate ancestry block sizes. the rates and volume of species diversiﬁcation among lineages. Taxonomic sampling. We sampled haplochromines from all major lakes and Methods several small lakes in the Lake Victoria region and from the Mpanga River that Experimental design. To identify the closest relatives of the Lake Victoria drains the Rwenzori Mountains in the drainage divide between Lakes Edward and Region Superﬂock (LVRS) and identify potential hybridization events, we Victoria (Supplementary Data 1, Fig. 1). Further, we collected or obtained NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 7 Proportion of LV outliers in each ‘ancestry category’ ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 the posterior distribution of the remaining 500 trees as consensus trees in from other collectors haplochromines belonging to all phylogenetic lineages Densitree . and from almost all African river drainages hosting haplochromines. We also included members of the haplochromine species ﬂocks of the other African Great Lakes Malawi, Tanganyika and Mweru (Fig. 1, Supplementary Data 1). Mitochondrial chronograms. Dated phylogenies were reconstructed based on Samples were collected under research and sample export permissions of the mitochondrial D-loop and ND2 (ref. 64) sequences using BEAST v. 2.3.0 (ref. 71) Tanzania Fisheries Research Institute and the Tanzanian Ministry of Agriculture, and four different sets of calibration nodes (Supplementary Methods). We Livestock and Fisheries; the Ugandan National Fisheries Resources Research caution that the mitochondrial tree only shows the phylogeny of the maternal Institute and Ministry of Agriculture, Animal Industry and Fisheries Uganda; line and that time estimates more recent than one million years are most likely and the Department of Fisheries Zambia. All samples were collected in compliance overestimates because of the increase of molecular rates towards the recent with applicable international and national guidelines for the use of animals, (Supplementary Methods). and ethical standards. To compare the splitting time between cichlid lineages with the major relevant geological events, we reconstructed paleogeographic maps at different time points based on previously published data and reviews (see Supplementary RAD sequencing. DNA was extracted from ﬁn clips or muscle tissue with Methods). a standard phenol-chloroform-isoamyl alcohol extraction method . Restriction- site Associated DNA sequencing (RADseq) was performed following a standard Patterson’s D statistics. To test for evidence of ancient admixture among protocol . Restriction digestion was done overnight (8–10 h) using the 46,72 restriction endonuclease HF-SbfI (NewEngland Biolabs) and 400–1,000 ng lineages, we computed Patterson’s D statistic (ABBA-BABA test), a method to detect admixture based on the frequencies of discordant SNP genealogies DNA per sample. P1 adaptors contained 5–8 bp long barcodes differing by at least two nucleotides from all other barcodes. The DNA was sheared with a Covaris in a four-taxon tree, with the software package ADMIXTOOLS v. 1.1 (ref. 48). Genotypes were discarded if they had less than 6 reads or a genotype quality Phred S220 Focused-Ultra sonicator and fragments of 300–600 bp length were extracted from an agarose gel. We performed 18 PCR cycles to amplify the RAD fragments scoreo20 (that is, error probability41%). Signiﬁcance of D statistics was assessed with a block jackknife procedure using a z score of three standard errors as (30 s 98 C, 18 (10 s 98 C, 30 s 65 C, 30 s 72 C), 5 min 72 C). All libraries a threshold . We used three individuals of Astatotilapia ﬂavijosephi from Lake were single-end sequenced on an Illumina HiSeq 2,500 sequencer. The reads were de-multiplexed and trimmed to 84 nucleotides (nt, after barcode removal) with Kinneret as the outgroup population and we tested for evidence of gene ﬂow from each Eastern and Upper Nile clade into the LVRS relative to allele sharing with the the process_radtags script from the Stacks pipeline , correcting single errors in the barcode and discarding reads with incomplete restriction sites. The length closest relative of the LVRS, the Congo drainage taxa A. stappersi from Zambia (Fig. 2) or A. sp. ‘Yaekama’ from the central Congo. For these analyses we used all of 84 nt results from removing the barcode (mostly 6 nt long) and trimming off the last 10 nt because of reduced quality at the read ends. The FastX toolkit individuals with at least 50% of the sites sequenced at a depth of at least 10 reads. To exclude the possibility that the D statistics are biased by the alignment of the (http://hannonlab.cshl.edu/fastx_toolkit) was used to remove all reads containing at least one base with a Phred quality score below 10 and reads with more than reads to the Lake Victoria species Pundamilia nyererei, we also aligned the reads to the Astatotilapia burtoni reference genome, which is an outgroup to all taxa used in 10% of bases with quality less than 30. The reads of each individual were then mapped to the Pundamilia nyererei the D statistics. To rule out that the choice of outgroup (A. ﬂavijosephi) biases the D statistics, we repeated the tests with A. burtoni or A. desfontainii as outgroups. reference genome using Bowtie2 (ref. 61) with the end-to-end alignment option. Single nucleotide polymorphisms (SNPs) and genotypes were called using All species combinations with number of individuals and SNPs included are given in Supplementary Table 2. GATK Uniﬁed Genotyper v. 3.1 (ref. 62). All sites were then ﬁltered with a custom- made Python script and vcftools v. 4.1 (ref. 63). SNPs within 5 nt from indels (insertions and deletions) were removed to avoid false SNPs because of Five population tests. To infer directionality of gene ﬂow between LVRS and the misalignment problems, and SNPs were required to have a quality value of Nilotic taxa, we used an extended version of the partitioned D statistic test at least 30. developed by Eaton & Ree (we call it a ‘ﬁve population test’). Similar to the D statistics, this test is based on discordant allele sharing patterns, but by considering ﬁve populations it allows one to infer the directionality of gene ﬂow. Mitochondrial sequencing. Two mitochondrial markers (NADH Dehydrogenase The ﬁve taxa with the topology ((P1,P2), (P3a,P3b)),O include the potential source Subunit 2 (ND2) using the primers ND2Met-F 5 -CAT ACC CCA AAC ATG of gene ﬂow (for example, P3a), the taxon receiving gene ﬂow (for example, P1), 0 0 0 TTG GT-3 and ND2Trp-R 5 -GTS GST TTT CAC TCC CGC TTA-3 and the a close relative of each of these two taxa (for example, P3b and P2) and an 65 0 mtDNA control region (D-loop) with the primers FISHL15926-F 5 -GAG CGC outgroup (O, see also Supplementary Table 3 for visualization). If genes had 0 0 CGG TCT TGT AA-3 and FISH12s-R 5 -TGC GGA GAC TTG CAT GTG TAA introgressed from for example, P3a into P1, a close relative (P3b) of the G-3 were ampliﬁed with PCR and Sanger sequenced for the same individuals introgression donor (P3a) would also show excess allele sharing with P1, but or downloaded from GenBank for the same species as those included in the to a lesser extent than P3a. This is because many derived alleles that introgressed RADseq dataset (Supplementary Data 1). The sequences were aligned in ClustalW from P3a into P1 will be shared by P3a and P3b because of their recent common implemented in BioEdit 7.2.5 (ref. 66) and manually curated for correct local ancestry. In contrast, we would not expect that a close relative (P2) of the receiver alignment. of gene ﬂow, P1, would show excess allele sharing with P3a and P3b. In the genomic data, this would be seen as an excess number of BABBA patterns, where P1 shares a derived allele (‘B’) with both P3a and P3b, whereas P2 has the outgroup Phylogenetic analyses. Phylogenies were reconstructed for the concatenated allele (‘A’), as compared with the number of ABBBA patterns (P2 shares a derived mitochondrial genes and the concatenated RAD sequences separately, including allele with both P3 taxa). On the other hand, if the direction of gene ﬂow was both variant and invariant sites using a maximum likelihood approach from P1 into P3a, P3b would not show excess allele sharing with P1 but instead, 67,68 (RAxML v. 7.7.7 and ExaML v. 1.0.4) . For the mitochondrial dataset, P2 would show excess allele sharing with P3a because of ancestrally shared a maximum likelihood tree was reconstructed using three partitions, one for Dloop, alleles between P1 and P2 that introgressed into P3a. Therefore, gene ﬂow from one for the ﬁrst and the second codon positions of ND2, and one for the third P1 into P3a would not affect the relative frequencies of BABBA and ABBBA codon position of ND2. For each dataset, we performed a RAxML analysis with patterns, but the pattern BBBAA (P1, P2 and P3a share a derived allele) would be 100 rapid bootstraps using the GTRGAMMA model of rate heterogeneity. For the more frequent than BBABA (P1 and P2 share a derived allele with P3b, see RAD-seq dataset, we used all concatenated sites with no more than 40 individuals Supplementary Table 3). Counting discordant allele sharing patterns thus allows us missing (25%) to reconstruct a maximum likelihood tree with RAxML (ref. 67) to infer the direction of gene ﬂow. and ExaML (ref. 68). Each of 100 bootstraps was performed by randomly sampling We computed the ﬁve population tests with a custom made script using the with replacement sites from the concatenated dataset to get a dataset of the original three individuals with least missing data for each LVRS group and the single size. The maximum likelihood tree was then inferred for each resampled individual with the most complete data for all other taxa. In contrast to the dataset with ExaML (ref. 68) using a GTRGAMMA model of rate heterogeneity, D statistics computed with ADMIXTOOLS, which are based on allele frequencies, as recommended in the RAxML-light manual . We calculated bootstrap support our ﬁve population tests are calculated from a single individual for each values based on these 100 topologies with RAxML (ref. 67). The nuclear and focal population, following Eaton and Ree . We tested each of the three mitochondrial trees were rooted with the reference genome of Oreochromis individuals of each LVRS group separately, and report the means for each niloticus and ladderized and plotted using the R-package Ape v. 3.1 (ref. 69). radiation. At heterozygous sites, one allele was chosen at random. For each We then used RAD-derived SNP data to infer the species tree for the LVRS combination of individuals tested, we counted the eight patterns needed to groups and their closest relatives with SNAPP (ref. 45). SNAPP bypasses gene compute the four D statistics (Supplementary Table 3) using all sites without trees and computes species trees directly from independently inherited markers missing data. We calculated z scores in units of s.d. from 100 bootstrapped datasets by integrating over all possible gene trees . We restricted this analysis to two 47 (sites resampled with replacement) as in Eaton and Ree . individuals per species (one for A. sp. ‘Yaekama’) and as SNAPP assumes no linkage among loci, we included only biallelic sites that were at least 500 kb apart from each other. The resulting data set contained 31 individuals and 1,817 sites. Estimation of ancestry proportions. The absolute value of D statistics does not We ran SNAPP for 1,000,000 iterations, sampling every 1,000th iteration using only depend on the admixture proportion but also on the demographic history, the default priors. We discarded the ﬁrst 50% of the trees as burn-in and visualized divergence between the hybrid parental lineages, and the genetic distance between 8 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE 73 inferred Upper Nile ancestry proportion for the Lake Victoria radiation). To the real source taxon of gene ﬂow and the sampled surrogate population . take the effect of this moderately high initial allele frequency into account in our Thus, the D statistic cannot be used directly to infer the magnitude of 48,73 subsequent analyses, we created a dataset of SNPs with comparable expected introgression. Instead, we applied the F4-ratio test to infer the Upper Nile ancestral allele frequency in Lake Victoria (category 5) by weighting the allele ancestry proportions in the different LVRS taxa using all genotypes with at least ten frequencies of the Congolese and Upper Nile lineage taxa at each site by 0.84 and reads and SNPs with a maximum of 10% missing data (72,443 SNPs) with 0.16, respectively, corresponding to the inferred ancestry proportions for the Lake ADMIXTOOLS v. 1.1 (ref. 48). The F4-ratio test is based on four populations with Victoria radiation. Of these, we extracted all sites with a weighted minor allele the genealogy (((A,B)C)O) and a ﬁfth potential hybrid population X that is tested frequency of 14–18% excluding sites that are divergently ﬁxed between the for ancestry proportions from B and C (Supplementary Fig. 6). We repeated the Congolese and Upper Nile taxa. As an example, if at a given bi-allelic site the test with different Eastern taxa used as population ‘A’, T. pharyngalis or H. gracilior frequency of one allele was 0.1 in the Congolese taxa and the frequency of the same as ‘B’, A. sp. ‘Yaekama’ as ‘C’ and A. ﬂavijosephi as ‘O’ as shown in Supplementary allele was 0.5 in the Upper Nile taxa, we calculated 0.1 0.84 þ 0.5 0.16 ¼ 0.164. Table 4. The Upper Nile ancestry proportions estimated with the combination of As this site would have a weighted minor allele frequency between 14% and 18%, it populations with the lowest s.d. were used for further analyses and to explain the would be included in the control set of SNPs. Note that such loci are not predicted F4-ratio test in Supplementary Fig. 6. to be involved in BDM incompatibilities because all their alleles were present within at least one of the two parental populations. The proportion of sites Long-wave sensitive (LWS) opsin gene. We sequenced the LWS opsin gene from that were outliers of high global F (Po0.05) among the Lake Victoria species ST exon 2 to 6 using the Carleton & Kocher primer combinations F2 (5 -TTT GAG (‘LV outliers’) was calculated for each dataset and compared between datasets 0 0 0 GGT CCC AAT TAC CA-3 ), R2 (5 -TCC ACA CAG CAA GGT AGC AC-3 ) and using a two-sided Fisher’s exact test. The robustness of the test to erroneously 0 0 0 F3 (5 -ACT GGC CTC ATG GAC TGA AG-3 ) and R4 (5 -TCC CAA AAT GGA inferred ancestry proportions was assessed by repeating the calculations GAA CAT GG-3 ) (http://cichlid.umd.edu/cichlidlabs/protocols/Basic/pcrpri- with different Upper Nile allele frequency weighting, ranging from 10 to mer.html) . The sequences were aligned together with all available Lake Victoria 30% (Supplementary Discussion). cichlid LWS sequences from GenBank using ClustalW implemented in BioEdit To check if LV outliers ﬁxed for alternative alleles in the parental lineages 7.2.5 (ref. 66) and then manually curated. For sequences used, see Supplementary inherently have a higher ﬁxation probability than non-outliers, we analysed ﬁxation Data 3. A maximum likelihood tree was reconstructed including both variant patterns in other closely related cichlid species that are not part of the radiation. and invariant sites using the GTRGAMMA model of rate heterogeneity with We used Astatotilapia bloyeti, A. sparsidens, A. ﬂavijosephi, A. calliptera, A. burtoni RAxML v. 7.7.7 (ref. 67). We used one partition with concatenated introns, one for and Astatoreochromis alluaudi as we have at least three sequenced individuals of the ﬁrst plus second codon position of concatenated exons, and one for the each. For each pairwise comparison between these species, we extracted sites with third codon positions. The tree was rooted with the Oreochromis niloticus at least three genotypes sequenced at a minimum depth of six reads in both species. sequence and nodal support values were drawn from 100 rapid bootstraps We then checked if sites ﬁxed for alternative alleles in the Congolese and performed with RAxML. The tree was ladderized and plotted using the R-package Upper Nile lineage taxa were also ﬁxed for alternative alleles between these other Ape v. 3.1 (ref. 69). species. Finally, we tested with Fisher’s exact tests if sites that were high global To visually assess the effects of recombination on the LWS opsin alleles in the F outliers among Lake Victoria species were also more often differentially ﬁxed ST Lake Victoria Region cichlids, we extracted all positions of the alignment that were between these other species than non-outlier sites (see also Supplementary polymorphic among the LVRS sequences and also differed between the Congolese Discussion). and Upper Nile taxa. If one of the alleles was only found in the Congolese taxa but To analyse the genomic distribution of LV outliers that are ﬁxed for not in the Upper Nile taxa, we coloured it red, whereas alleles found only in alternative alleles between the Congolese and Upper Nile taxa, we used liftOver the Upper Nile taxa were coloured blue. Many LVRS haplotypes could be explained (http://genome.ucsc.edu/) with a chain ﬁle from Brawand et al. to get positions by recombination between the two divergent allele classes occurring in the Congo on the Oreochromis niloticus genome which includes chromosomal information. and Upper Nile clades. The second LWS gene tree was reconstructed with all recombinant sequences falling between the two major allele classes in the phylogenetic tree removed from the alignment. Whole-genome sequencing data. To corroborate our ﬁndings from For information about the occurrence of LWS haplotypes across LVRS species, 51 RAD sequence data, and for the analysis of ancestry block sizes, we sequenced we compiled data from Tables 1 þ 2 in Terai et al. , Fig. 1 and Supplementary 49 39 whole genomes of nine ﬁsh from Lake Victoria, two from Lake Kivu, and of the Table 4 in Terai et al. , Supplementary Fig. 4b in Seehausen et al. , Figure S3 in 54 taxa that we found to be the closest extant representatives of the lineages directly Miyagi et al. and additional GenBank sequences by Carleton et al. (GenBank). ancestral to the LVRS, the Congolese Astatotilapia stappersi, and the Upper Nile Information sources and GenBank accession numbers are provided in Thoracochromis pharyngalis and ‘Haplochromis’ gracilior. Whole-genome Supplementary Data 3. sequencing data was generated using PCR-free library preparation and Illumina HiSeq 3000 paired-end sequencing for 11 individuals. Three additional individuals, Sorting of ancestral alleles. We created a RAD sequence dataset of our sequenced the same way on the same machine, were taken from McGee et al. Congolese and Upper Nile LVRS relatives and previously published sequences of (see Supplementary Data 1 for sample information). Local alignment against the six sympatrically occurring Lake Victoria species that are ecologically and Astatotilapia burtoni reference genome was performed with Bowtie 2 (ref. 61). phenotypically diverse (Pundamilia pundamilia (n ¼ 10, blue male nuptial For variant calling and genotyping we used Haplotype Caller (GATK v. 3.5) . coloration, benthic insectivore-omnivore), P. nyererei (n ¼ 11, red planktivore- Genotypes with fewer than ﬁve reads, multiallelic sites, and indels were removed omnivore), Harpagochromis cf. serranus (n ¼ 10, blue piscivore), Lipochromis using vcftools v. 4.1 (ref. 63). melanopterus (n ¼ 9, yellow paedophage), Paralabidochromis chilotes (n ¼ 9, blue Whole-genome D statistics were calculated with ADMIXTOOLS v. 1.1 (ref. 48) benthic insectivore) and Neochromis omnicaeruleus (n ¼ 12, blue algae scraper)). using all biallelic sites with a 1% minor allele frequency cutoff and a maximum The sites were ﬁltered for having at least 5 individuals covered in each species missing data proportion of 20% across all 14 genomes. The reference genome at a sequencing depth of at least six reads. We used Arlequin v. 3.514 (ref. 75) to (A. burtoni) was used as an outgroup. calculate global F values based on a hierarchical island model among the six Lake To study signatures of admixture along the genome, we used f statistics by ST d Victoria cichlid species. We calculated allele frequencies for the Congolese Martin et al. . In comparison to D statistics, f is more suited for small genomic (Astatotilapia stappersi from Zambia and A. sp. ‘Yaekama’ from DRC, n ¼ 5) and regions . ABBA and BABA pattern counts were calculated using allele frequencies the Upper Nile taxa (‘Haplochromis’ gracilior and Thoracochromis pharyngalis, by weighting each segregating site according to its ﬁt to the ABBA or BABA 46,72,78 n ¼ 9) to get estimates of the allele frequencies for the parental lineages of the LVRS pattern . As an example, if a site has derived allele frequencies of 0, 0.5, 1 and and to infer the likely ancestry of the alleles segregating in Lake Victoria. We 0 in P1, P2, P3 and the outgroup, respectively, it would count as half ABBA site. excluded sites monomorphic in the subset of Lake Victoria species or sequenced in We used A. stappersi as P1, different LVRS individuals as P2, H. gracilior and less than three individuals in either the Congolese or the Upper Nile lineage taxa. T. pharyngalis as P3 and A. burtoni as outgroup (Supplementary Fig. 8). Note that We extracted sites for which only one of the two alleles segregating in Lake Victoria P1 and P2 are switched as compared to D statistics, where LVRS are used as P1 and cichlids were found in the Congolese and Upper Nile lineage taxa together the Congolese lineage as P2. This difference is simply for consistency with the (category 1). Next, we extracted sites with both alleles segregating in the Congolese description of f in Martin et al. , and thus here ABBA represents sites with taxa representing ancestral standing genetic variation that would have been present a derived allele shared between LVRS and Upper Nile, whereas BABA represents in the Lake Victoria radiation also without ancient hybridization with the Upper sites with a derived allele shared between Congolese and Upper Nile lineage taxa. Nile lineage (category 2). f is calculated as the difference between ABBA and BABA patterns compared to We identiﬁed Lake Victoria polymorphic sites for which only one allele was the maximum possible difference where gene ﬂow between P2 and P3 would equal found in the Congolese taxa, and the second allele was found exclusively in the random mating . Positive f values indicate gene ﬂow between P2 (LVRS) and Upper Nile taxa, representing genetic variation that can be attributed to the ancient P3 (Upper Nile). We calculated f in non-overlapping sliding windows of 10 kb admixture event (category 3). We also identiﬁed the subset of sites ﬁxed for along the A. burtoni scaffolds using the Python script by Martin et al. . Windows alternative alleles in the Congolese and Upper Nile taxa (category 4). These are sites with less than ﬁve total ABBA and BABA patterns were excluded. As only positive that could be involved in epistatic BDM-like incompatibilities between Congolese f values are indicative of excess allele sharing between Upper Nile and LVRS, and and Upper Nile alleles. Sites ﬁxed for alternative alleles in the Congolese and are correctly standardized, values with negative D scores were set to 0 as in Martin Upper Nile lineage taxa are expected to have been present at the onset of the et al. . Pearson’s product-moment correlations of f values between different Lake Victoria radiation at a mean minor allele frequency of 16% (that is, the individuals were calculated with the stats R-package. NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 Putative Congolese or Upper Nile ancestry was assessed by comparing the 20. Schliewen, U. K. & Klee, B. Reticulate sympatric speciation in Cameroonian frequency of ABBA sites (LVRS shares the derived allele exclusively with crater lake cichlids. Front. Zool. 1, 5 (2004). Upper Nile taxa) with the frequency of BBAA sites (LVRS shares the derived 21. Seehausen, O. Conditions when hybridization might predispose populations for allele exclusively with the Congolese lineage representative A. stappersi)in adaptive radiation. J. Evol. Biol. 26, 279–281 (2013). non-overlapping sliding windows of 3 kb. Scaffolds with o100 kb data were 22. Schumer, M., Cui, R., Rosenthal, G. G. & Andolfatto, P. Reproductive isolation removed. As for the f statistic, ABBA and BBAA pattern counts were of hybrid populations driven by genetic incompatibilities. PLoS Genet. 11, calculated using allele frequencies by weighting each segregating site according e1005041 (2015). to its ﬁt to the ABBA or BBAA pattern. Only windows with a total ABBA and 23. Hermansen, J. S. et al. Hybrid speciation through sorting of parental BBAA pattern count exceeding one were used. Windows with a minimum ABBA incompatibilities in Italian sparrows. Mol. Ecol. 23, 5831–5842 (2014). proportion ((ABBA/(ABBA þ BBAA)) of 0.7 were deﬁned as candidate windows of 24. Barrier, M., Baldwin, B. G., Robichaux, R. H. & Purugganan, M. D. Interspeciﬁc Upper Nile lineage ancestry and are coloured in blue, whereas genomic regions hybrid ancestry of a plant adaptive radiation: alloploidy of the Hawaian with an ABBA proportion of 0.3 or less were deﬁned as putative Congolese lineage silversword alliance (Asteraceae) inferred from ﬂoral homeotic gene derived windows and are highlighted in red. Ancestry tracts were deﬁned as duplication. Mol. Biol. Evol. 16, 1105–1113 (1999). consecutive sliding windows of the same colour (red or blue) ignoring single 25. Lynch, M. & Conery, J. S. The evolutionary fate and consequences of duplicate sliding windows without data. Expected ancestry block sizes were calculated using a formula from Racimo et al. . Assuming B50,000–100,000 generations genes. Science 290, 1151–1155 (2000). since admixture, a recombination rate of 2.5 10 and an Upper Nile 26. Hudson, A. G., Vonlanthen, P. & Seehausen, O. Rapid parallel adaptive proportion of 20%, the expected mean length of the admixture tracts is radiations from a single hybridogenic ancestral population. Proc. R Soc. B 278, 8 1 (0.8 2.5 10 (50,000 to 100,000) 1)) ¼ 500 bp to 1 kb (ref. 79). 58–66 (2011). The breakup of ancestry blocks may have slowed as the original hybrid population 27. Joyce, D. A. et al. Repeated colonization and hybridization in Lake Malawi began to undergo genomic stabilization, speciated and formed geographically cichlids. Curr. Biol. 21, R108–R109 (2011). isolated radiations in separate lakes each of which underwent further stabilization 28. Lindqvist, C., Motley, T. J., Jeffrey, J. J. & Albert, V. A. Cladogenesis and independent of each other, likely associated with differential sorting of the ancestral reticulation in the Hawaiian endemic mints (Lamiaceae). Cladistics 19, 480–495 variation. (2003). 29. Lindqvist, C. & Albert, V. A. Origin of the Hawaiian endemic mints within North American Stachys (Lamiaceae). Am. J. Bot. 89, 1709–1724 Data availability. Mitochondrial and LWS opsin sequences are available on ð2002Þ: GenBank under the accession numbers KY366716-KY366843 for ND2, KY366844- 30. Baldwin, B. G. & Wagner, W. L. Hawaiian angiosperm radiations of North KY366970 for D-loop, and KY366971-KY366986 for LWS opsin sequences. American origin. Ann. Bot. 9, 849–879 (2010). RADseq and whole-genome sequencing reads generated in this study can be 31. Verheyen, E., Salzburger, W., Snoeks, J. & Meyer, A. Origin of the superﬂock downloaded from the NCBI Sequence Read Archive under Bioproject of cichlid ﬁshes from Lake Victoria, East Africa. Science 300, 325–329 PRJNA355227. The JAVA program for 5-population tests is publicly available on ð2003Þ: GitHub (https://github.com/joanam/scripts). 32. Seehausen, O. African cichlid ﬁsh: a model system in adaptive radiation research. Proc. R. Soc. B 273, 1987–1998 (2006). 33. Genner, M. J. et al. Age of cichlids: New dates for ancient lake ﬁsh radiations. References Mol. Biol. Evol. 24, 1269–1282 (2007). 1. Wagner, C. E., Harmon, L. J. & Seehausen, O. Ecological opportunity and 34. Bezault, E., Mwaiko, S. & Seehausen, O. Population genomic tests of models sexual selection together predict adaptive radiation. Nature 487, 366–369 of adaptive radiation in Lake Victoria Region cichlid ﬁsh. Evolution 65, (2012). 3381–3397 (2011). 2. Barrett, R. D. H. & Schluter, D. Adaptation from standing genetic variation. 35. Johnson, T. C., Kelts, K. & Odada, E. The holocene history of Lake Victoria. Trends Ecol. Evol. 23, 38–44 (2008). Ambio 29, 2–11 (2000). 3. Seehausen, O. Hybridization and adaptive radiation. Trends Ecol. Evol. 19, 36. Wagner, C. E. et al. Genome-wide RAD sequence data provide unprecedented 198–207 (2004). resolution of species boundaries and relationships in the Lake Victoria cichlid 4. Abbott, R. et al. Hybridization and speciation. J. Evol. Biol. 26, 229–246 adaptive radiation. Mol. Ecol. 22, 787–798 (2013). ð2013Þ: 37. Seehausen, O. et al. Nuclear markers reveal unexpected genetic variation and 5. Rieseberg, L. H. Hybrid origins of plant species. Ann. Rev. Ecol. Syst. 28, a Congolese-Nilotic origin of the Lake Victoria cichlid species ﬂock. Proc. 359–389 (1997). R. Soc. Lond. B Biol. Sci. 270, 129–137 (2003). 6. Mallet, J. Hybrid speciation. Nature 446, 279–283 (2007). 38. Stelkens, R. B., Young, K. A. & Seehausen, O. The accumulation of reproductive 7. Anderson, E. & Stebbins, G. L. Hybridization as an evolutionary stimulus. incompatibilities in African cichlid ﬁsh. Evolution 64, 617–632 (2010). Evolution 8, 378–388 (1954). 39. Seehausen, O. et al. Speciation through sensory drive in cichlid ﬁsh. Nature 8. Carlquist, S. Island Biology Vol. 581, 5279 (Columbia Univ Press, 1974). 455, 620–626 (2008). 9. Arnold, M. L. Natural Hybridization and Evolution (Oxford University Press, 40. Seehausen, O., van Alphen, J. & Witte, F. Can ancient colour polymorphisms 1997). explain why some cichlid lineages speciate rapidly under disruptive sexual 10. Feder, J. L. et al. Mayr, Dobzhansky, and Bush and the complexities of selection? Belg. J. Zool. 129, 43–60 (1999). sympatric speciation in Rhagoletis. Proc. Natl Acad. Sci. USA 102, 6573–6580 41. Meyer, B. S. et al. Back to Tanganyika: a case of recent trans-species-ﬂock (2005). dispersal in East African haplochromine cichlid ﬁshes. R. Soc. Open Sci. 2, 11. The Heliconius Genome Consortium. Butterﬂy genome reveals promiscuous 140498 (2015). exchange of mimicry adaptations among species. Nature 487, 94–98 42. Greenwood, P. H. Towards a phyletic classiﬁcation of the ‘genus’ Haplochromis ð2012Þ: (Pisces, Cichlidae) and related taxa. Part 1. Bull. Br. Mus. Nat. Hist. Zool. 35, 12. Pardo-Diaz, C. et al. Adaptive introgression across species boundaries in 265–322 (1979). Heliconius butterﬂies. PLOS Genet. 8, e1002752 (2012). 43. Hermann, C. M., Sefc, K. M. & Koblmu¨ller, S. Ancient origin and recent 13. Lamichhaney, S. et al. Evolution of Darwin’s ﬁnches and their beaks revealed by genome sequencing. Nature 518, 371–375 (2015). divergence of a haplochromine cichlid lineage from isolated water bodies in the East African Rift system. J. Fish Biol. 79, 1356–1369 (2011). 14. Salzburger, W., Baric, S. & Sturmbauer, C. Speciation via introgressive 44. Maddison, W. P. Gene trees in species trees. Syst. Biol. 46, 523–536 (1997). hybridization in East African cichlids? Mol. Ecol. 11, 619–625 (2002). 45. Bryant, D., Bouckaert, R., Felsenstein, J., Rosenberg, N. A. & RoyChoudhury, A. 15. Meyer, B. S., Matschiner, M. & Salzburger, W. Disentangling incomplete lineage sorting and introgression to reﬁne species-tree estimates for Lake Inferring species trees directly from biallelic genetic markers: bypassing Tanganyika cichlid ﬁshes. Syst. Biol. doi: 10.1093/sysbio/syw069 (2016). gene trees in a full coalescent analysis. Mol. Biol. Evol. 29, 1917–1932 16. Weiss, J. D., Cotterill, F. P. & Schliewen, U. K. Lake Tanganyika—A’Melting ð2012Þ: 46. Green, R. E. et al. A draft sequence of the Neandertal genome. Science 328, Pot’of Ancient and Young Cichlid Lineages (Teleostei: Cichlidae)? PLoS ONE 10, e0125043 (2015). 710–722 (2010). 17. Genner, M. J. & Turner, G. F. Ancient hybridization and phenotypic 47. Eaton, D. A. R. & Ree, R. H. Inferring phylogeny and introgression using novelty within Lake Malawi’s cichlid ﬁsh radiation. Mol. Biol. Evol. 29, 195–206 RADseq data: an example from ﬂowering plants (Pedicularis: Orobanchaceae). (2012). Syst. Biol. 62, 689–706 (2013). 18. Keller, I. et al. Population genomic signatures of divergent adaptation, gene 48. Patterson, N. et al. Ancient admixture in human history. Genetics 192, ﬂow and hybrid speciation in the rapid radiation of Lake Victoria cichlid ﬁshes. 1065–1093 (2012). Mol. Ecol. 22, 2848–2863 (2013). 49. Terai, Y. et al. Divergent selection on opsins drives incipient speciation in Lake 19. Meier, J. I. et al. Demographic modelling with whole-genome data reveals Victoria cichlids. PLoS Biol. 4, 2244–2251 (2006). 50. Brawand, D. et al. The genomic substrate for adaptive radiation in African parallel origin of similar Pundamilia cichlid species after hybridization. Mol. Ecol, 26, 123–141 (2017). cichlid ﬁsh. Nature 513, 375–381 (2014). 10 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14363 ARTICLE 51. Terai, Y., Mayer, W. E., Klein, J., Tichy, H. & Okada, N. The effect of selection 76. Kozarewa, I. et al. Ampliﬁcation-free Illumina sequencing-library preparation on a long wavelength-sensitive (LWS) opsin gene of Lake Victoria cichlid ﬁshes. facilitates improved mapping and assembly of (G þ C)-biased genomes. Nat. Methods 6, 291–295 (2009). Proc. Natl Acad. Sci. USA 99, 15501–15506 (2002). 77. McGee, M. D., Neches, R. Y. & Seehausen, O. Evaluating genomic divergence 52. Okullo, W. et al. Parameterization of the inherent optical properties of Murchison Bay, Lake Victoria. Appl. Opt. 46, 8553–8561 (2007). and parallelism in replicate ecomorphs from young and old cichlid adaptive 53. Maan, M. E., Hofker, K. D., van Alphen, J. J. M. & Seehausen, O. Sensory drive radiations. Mol. Ecol. 25, 260–268 (2016). 78. Martin, S. H., Davey, J. W. & Jiggins, C. D. Evaluating the use of ABBA–BABA in cichlid speciation. Am. Nat. 167, 947–954 (2006). statistics to locate introgressed loci. Mol. Biol. Evol. 32, 244–257 (2015). 54. Miyagi, R. et al. Correlation between nuptial colors and visual sensitivities 79. Racimo, F., Sankararaman, S., Nielsen, R. & Huerta-Sanchez, E. Evidence tuned by opsins leads to species richness in sympatric Lake Victoria cichlid for archaic adaptive introgression in humans. Nat. Rev. Genet. 16, 359–371 ﬁshes. Mol. Biol. Evol. 29, 3281–3296 (2012). (2015). 55. Selz, O. M., Pierotti, M. E. R., Maan, M. E., Schmid, C. & Seehausen, O. Female 80. Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from preference for male color is necessary and sufﬁcient for assortative mating in spaceborne elevation data. EOS 89, 93–94 (2008). 2 cichlid sister species. Behav. Ecol. 25, 612–626 (2014). 56. Stelkens, R. B., Schmid, C. & Seehausen, O. Hybrid breakdown in cichlid ﬁsh. PLoS ONE 10, e0127207 (2015). Acknowledgements 57. Stelkens, R. B., Schmid, C., Selz, O. & Seehausen, O. Phenotypic novelty in We thank all collaborators who provided tissue samples over the years: Sigal Balshine- experimental hybrids is predicted by the genetic distance between species of Earn (McMaster), Roger Bills (SAIAB), Julia Day (UCL), Yves Fermon (Paris), John Friel cichlid ﬁsh. BMC Evol. Biol. 9, 283–295 (2009). (CUMV), Martin Genner (U Bristol), Sylvain Piry (Paris), Lukas Ru¨ber (Bern Natural 58. Sambrook, J. & Russell, D. W. Molecular Cloning: A Laboratory Manual 3rd History Museum), Walter Salzburger (U Basel), Uli Schliewen (Zoologische Staats- edn (Cold Spring Harbor Laboratory Press, 2001). sammlung Munich), Erwin Schraml (Augsburg), Jos Snoeks (Africa Museum Tervuren), 59. Baird, N. A. et al. Rapid SNP discovery and genetic mapping using sequenced Melanie Stiassny (AMNH), George Turner (U Bangor), Sylvester Wandera (NAFIRI; RAD markers. PLoS ONE 3, e3376 (2008). Uganda) and Marco Welss (Kressberg). We also thank Matt McGee for help with 60. Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: DNA extraction, Keith Harshman of the Lausanne Genomic Technologies Facility and an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140 (2013). Cord Drogemu¨ller, Tosso Leeb, Muriel Fragnie`re and Miche`le Ackermann of the 61. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. NGS platform of the University of Bern for Illumina sequencing support, Aria Minder Methods 9, 357–U354 (2012). and Stefan Zoller of the Genetic Diversity Center (GDC) at ETH Zu¨rich and Irene Keller 62. McKenna, A. et al. The Genome Analysis Toolkit: A MapReduce framework for of the Interfaculty Bioinformatics Unit at the University of Bern for lab and bioinfor- analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 matics support and discussion, and Andy Cohen for comments on the paleogeographic (2010). maps. This research was supported by the Swiss National Science Foundation grant 63. Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, PDFMP3 134657 to O.S. and L.E. 2156–2158 (2011). 64. Kocher, T. D., Conroy, J. A., McKaye, K. R., Stauffer, J. R. & Lockwood, S. F. Author contributions Evolution of NADH dehydrogenase subunit 2 in east African cichlid ﬁsh. Mol. O.S., J.I.M. and C.E.W. designed the study; O.S. gathered and identiﬁed the cichlid Phylogenet. Evol. 4, 420–432 (1995). samples; J.I.M., S.M. and D.A.M. performed the lab work (DNA extraction, mtDNA and 65. Kocher, T. D. et al. Dynamics of mitochondrial DNA evolution in animals: LWS sequencing and RAD library preparation); J.I.M. conducted the analyses, with ampliﬁcation and sequencing with conserved primers. Proc. Natl Acad. Sci. assistance from O.S., C.E.W., D.A.M. and L.E.; J.I.M. prepared the manuscript together USA 86, 6196–6200 (1989). with O.S. and C.E.W. and L.E. and D.A.M. contributed to writing. 66. Hall, T. A. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acid Symp. 41, 95–98 (1999). Additional information 67. Stamatakis, A. RAxML-VI-HPC: maximum likelihood-based phylogenetic Supplementary Information accompanies this paper at http://www.nature.com/ analyses with thousands of taxa and mixed models. Bioinformatics 22, naturecommunications 2688–2690 (2006). 68. Stamatakis, A. et al. RAxML-Light: a tool for computing terabyte phylogenies. Competing ﬁnancial interests: The authors declare no competing ﬁnancial interests. Bioinformatics 28, 2064–2066 (2012). 69. Popescu, A. A., Huber, K. T. & Paradis, E. ape 3.0: new tools for distance-based Reprints and permission information is available online at http://npg.nature.com/ phylogenetics and evolutionary analysis in R. Bioinformatics 28, 1536–1537 reprintsandpermissions/ (2012). How to cite this article: Meier, J. I. et al. Ancient hybridization fuels rapid cichlid ﬁsh 70. Bouckaert, R. R. DensiTree: making sense of sets of phylogenetic trees. adaptive radiations. Nat. Commun. 8, 14363 doi: 10.1038/ncomms14363 (2017). Bioinformatics 26, 1372–1373 (2010). 71. Bouckaert, R. et al. BEAST 2: a software platform for Bayesian evolutionary Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in analysis. PLoS Comput. Biol. 10, e1003537 (2014). published maps and institutional afﬁliations. 72. Durand, E. Y., Patterson, N., Reich, D. & Slatkin, M. Testing for ancient admixture between closely related populations. Mol. Biol. Evol. 28, 2239–2252 (2011). This work is licensed under a Creative Commons Attribution 4.0 73. Reich, D., Thangaraj, K., Patterson, N., Price, A. L. & Singh, L. Reconstructing International License. The images or other third party material in this Indian population history. Nature 461, 489–494 (2009). article are included in the article’s Creative Commons license, unless indicated otherwise 74. Carleton, K. L. & Kocher, T. D. Cone opsin genes of African cichlid ﬁshes: in the credit line; if the material is not included under the Creative Commons license, tuning spectral sensitivity by differential gene expression. Mol. Biol. Evol. 18, users will need to obtain permission from the license holder to reproduce the material. 1540–1550 (2001). To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ 75. Excofﬁer, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567 (2010). r The Author(s) 2017 NATURE COMMUNICATIONS | 8:14363 | DOI: 10.1038/ncomms14363 | www.nature.com/naturecommunications 11

Journal

Nature Communications – Springer Journals

Published: Feb 10, 2017

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Ancient hybridization fuels rapid cichlid fish adaptive radiations

Ancient hybridization fuels rapid cichlid fish adaptive radiations

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Ancient hybridization fuels rapid cichlid fish adaptive radiations

Ancient hybridization fuels rapid cichlid fish adaptive radiations

References (96)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies