Research Article

Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake

See allHide authors and affiliations

Science  18 Dec 2015:
Vol. 350, Issue 6267, pp. 1493-1498
DOI: 10.1126/science.aac9927

Cichlids diverge within a crater lake

It is not clear how populations diversify and new species form at the genomic level, especially when they coexist in the same location. Malinsky et al. investigated how two ecomorphs of cichlid fish in a small lake in Tanzania are diversifying relative to each other. Although there is gene flow between the two forms, major regions of genetic divergence, known as genomic islands, separate the populations. Within these islands, the authors found genes likely to be associated with mate choice, supporting the idea that genetic changes related to breeding preferences are the first to diverge during speciation.

Science, this issue p. 1493


The genomic causes and effects of divergent ecological selection during speciation are still poorly understood. Here we report the discovery and detailed characterization of early-stage adaptive divergence of two cichlid fish ecomorphs in a small (700 meters in diameter) isolated crater lake in Tanzania. The ecomorphs differ in depth preference, male breeding color, body shape, diet, and trophic morphology. With whole-genome sequences of 146 fish, we identified 98 clearly demarcated genomic “islands” of high differentiation and demonstrated the association of genotypes across these islands with divergent mate preferences. The islands contain candidate adaptive genes enriched for functions in sensory perception (including rhodopsin and other twilight-vision–associated genes), hormone signaling, and morphogenesis. Our study suggests mechanisms and genomic regions that may play a role in the closely related mega-radiation of Lake Malawi.

Understanding the causes and consequences of speciation, including at the genetic level, requires the investigation of taxa at different stages on the speciation continuum (1, 2). East African cichlids have repeatedly undergone rapid adaptive radiation (3). The Lake Malawi radiation has generated over 500 species in less than 5 million years, involving divergence in habitat, feeding apparatus, and nuptial color. Thus, these phenomena present an opportunity to observe hundreds of varied, recent, and in some cases ongoing, speciation events (4). However, the investigation of early stages of speciation in the large cichlid radiations of Malawi as well as Lakes Tanganyika and Victoria has been hampered by difficulties in identifying sister species relationships, in reconstructing past geographical situations, and in controlling for possible introgression from non-sister taxa (5).

During 2011, we conducted a survey (table S1) of fish fauna in six crater lakes in the Rungwe District of Tanzania (Fig. 1A and table S2). In all six lakes, we found endemic haplochromine cichlids of the genus Astatotilapia, closely related to A. calliptera (Fig. 1A), a species widely distributed in the rivers, streams, and shallow lake margins of the region. Thus, the Rungwe District Astatotilapia are close relatives of the of Lake Malawi endemic haplochromine cichlids (5).

Fig. 1 Cichlid radiation in the crater lakes of southern Tanzania.

(A) A phylogeny of the crater lake Astatotilapia based on reference-aligned RAD data (7906 SNPs across 5010 polymorphic RAD loci). It demonstrates reciprocal monophyly between the populations in each lake except for Itamba, and a close relationship to A. calliptera from rivers and from Lake Malawi. Within Lake Massoko, yellow symbols indicate the littoral morph, blue symbols indicate the benthic morph, and gray symbols denote small, phenotypically ambiguous, and thus unassigned individuals. Additional Astatotilapia individuals from other crater lakes are denoted by open circles. A. calliptera from rivers and Lake Malawi are denoted by black circles. Bootstrap values are displayed for nodes with >50% support. (B) Breeding males of the yellow littoral and blue benthic morphs of Lake Massoko. The symbols next to the photographs correspond to symbols used in (C) to (G). (C to F) Morphological divergence between the two morphs of Lake Massoko. Relative to the littoral morph, the benthic morph has relatively longer head and jaw (C), lower body mass (D), narrower “papilliform” pharyngeal teeth (E), and lighter lower pharyngeal jaws (F). The benthic fish have stable isotope ratios that tend to be more depleted in C13 than the littoral fish, indicative of a more offshore-planktonic diet (G).

In Lake Massoko (fig. S1), the benthic zone in deep waters (~20 to 25 m) is very dimly lit and populated by cichlids with phenotypes clearly different from those typical of shallow waters (≲5 m) close to the shore (littoral). Deep-water males are dark blue-black, whereas most males collected from the shallow waters are yellow-green, similar to riverine A. calliptera (Fig. 1B, movie S1, and table S3). We also collected small (<65 mm standard length) males that were not readily field-assigned to either ecomorph (6). The benthic and littoral morphs are reminiscent of the species pair of Pundamilia cichlids from Lake Victoria (7), but within a potentially simpler historical and geographical context. Lake Massoko is steep-sided, has a strong thermocline at ~15 m, and an anoxic boundary at ~25 m (8). The estimated time of lake formation is ~50,000 years ago (9).

Ecomorph separation

To examine relationships between crater lake and riverine A. calliptera of southern Tanzania, we obtained restriction site–associated DNA (RAD) data from 30 fish from the Rungwe District and 11 outgroup Astatotilapia from the broader Lake Malawi catchment and a neighboring catchment (fig. S2 and table S4). A maximum-likelihood phylogeny constructed on the basis of these data (6) demonstrates monophyly of all specimens from Lake Massoko (Fig. 1A). Thus, the RAD phylogeny provides evidence that Massoko morphs might have evolved in primary sympatry, as proposed for crater lake cichlid radiations of Cameroon (10) and Nicaragua (11).

Morphological analyses of these two color morphs revealed significant differences in head and body shape, body mass, the shape of pharyngeal teeth, and pharyngeal jaw mass [Fig. 1, C to F, and table S5; analysis of covariance (ANCOVA) tests, all P < 0.001]. We also found significant differences in stable isotope ratios (Fig. 1G and table S5; ANCOVA test, P < 0.001), which are indicative of dietary differences. Together these results demonstrate ecomorph separation and adaptation to different ecological environments.

Whole-genome evidence

To study the genome-wide pattern of Massoko ecomorph divergence and to further clarify its geographical context, we obtained whole-genome sequence data at ~15× coverage for 6 individuals each of the yellow littoral and blue benthic ecomorphs and 16 additional A. calliptera from the wider Lake Malawi catchment and two adjacent catchments (fig. S2), supplemented by lower-coverage (~6×) data from 87 specimens from Lake Massoko (25 littoral, 32 benthic, and 30 small unassigned) and 30 individuals from Lake Itamba (Fig. 1A and table S6). Sequence data were aligned to the Metriaclima zebra reference assembly (12), from which divergence was 0.2 to 0.3%, and variants were called at 4,755,448 sites (1.2 to 1.6 million sites per individual).

A maximum-likelihood phylogeny built from whole-genome sequence data confirmed reciprocal monophyly of Astatotilapia within Lakes Massoko and Itamba and revealed the sister group of Massoko fish to be an A. calliptera population from the nearby Mbaka River (Fig. 2A). All specimens of the benthic ecomorph formed a monophyletic clade derived from the littoral ecomorph (Fig. 2A). Principal components analysis (PCA) showed strong population structure (Tracy-Widom statistics: P < 1 × 10−12), with benthic and littoral individuals separated by the first eigenvector and forming separate clusters (Fig. 2B). In contrast, within Lake Itamba, PCA did not reveal significant population structure (Tracy-Widom statistics: P = 0.11). Individuals from Massoko that were not field-assigned to either of the ecomorphs did not form a monophyletic clade in the phylogeny (Fig. 2A) or a distinct cluster in PCA (Fig. 2B).

Fig. 2 Whole-genome sequence data.

(A) A maximum-likelihood whole-genome phylogenetic tree. Black stars indicate nodes with 100% bootstrap support. The red star highlights the branch that separates all the Massoko benthic samples from the rest of the phylogeny (benthic ecomorph individuals are monophyletic in 50% of bootstrap samples). (B) PCA of genetic variation within Lake Massoko. (C to E) Genome-wide pattern of FST divergence in windows of 15 variants each. Darker color indicates greater density of data points. (C) Divergence between benthic and littoral ecomorphs within Massoko. (D) Divergence between combined Massoko and Itamba populations. (E) Absolute standard scores of Massoko-Itamba divergence (purple) overlaid on divergence between benthic and littoral ecomorphs (green).

Analysis of fine-scale genetic relationships with fineSTRUCTURE (13) supports the monophyly of the benthic ecomorph within the littoral, but also suggests that compared with the benthic population, the littoral population has greater coancestry with other A. calliptera; in particular with the Mbaka River sample (fig. S3). Therefore, we tested for evidence of secondary gene flow, as seen in cichlid populations from Cameroonian crater lakes (14). Under the null hypothesis of no differential gene flow into Massoko, A. calliptera from the Mbaka River should share derived alleles equally often with the littoral and with the benthic populations (15, 16). Instead, we found a small excess of shared derived alleles between A. calliptera from the Mbaka River and the littoral population, when compared with the benthic population (Patterson’s D = 1.1%; 4.86 SD from 0% or P < 5.8 × 10−7) (6). The proportion of admixture f with Mbaka was estimated at 0.9 ± 0.2%. This value is low, at a proportion that is approximately half of the Neanderthal introgression into non-African humans (15), and cross-coalescence rate analysis with the multiple sequentially Markovian coalescent (MSMC) (6, 17) indicates an average separation time of both Massoko ecomorphs from other A. calliptera samples (including the Mbaka River) approximately 10 times earlier than the split between the two ecomorphs (fig. S5). Thus, it is unlikely that a secondary invasion from the neighboring river systems (fig. S11B) contributed to the divergence of the ecomorphs.

We estimated individual ancestries for all Massoko and A. calliptera specimens with ADMIXTURE (6, 18) (fig. S4). Focusing on the Massoko samples, 11 of the 31 samples field-assigned as littoral were identified as admixed, with admixture fraction >25% from the benthic gene pool. No individuals identified as benthic were estimated to be admixed to the same extent; therefore, recent gene flow may be biased from deep to shallow waters. Ten of the 30 unassigned individuals were also identified as >25% admixed, whereas the remaining 20 unassigned samples appear to represent subadult individuals of both benthic and littoral ecomorphs (fig. S4A). When additional A. calliptera samples were included in ADMIXTURE analysis, a small amount of gene flow into Massoko was apparent with K = 2 ancestral populations (fig. S4C), consistent with the fineSTRUCTURE and Patterson’s D results described above. This analysis also suggests similar or even stronger gene flow out of Massoko and into the Mbaka River (fig. S4C).

Islands of speciation

Interestingly, there are no fixed differences between Massoko benthic and littoral ecomorphs. Genome-wide divergence FST is 0.038, and almost half (47.6%) of the variable sites have zero FST (table S7). Above the low background, a genome-wide FST profile shows clearly demarcated “islands” of high differentiation (Fig. 2, C and E). For single sites, the maximum FST is 13.6 SD above the mean, and 7543 sites have FST over 6 SD above the mean. In contrast, comparisons of the combined Massoko population and Itamba population revealed a pattern of consistently high FST across the genome (Fig. 2D) and large variance, making it impossible to detect statistical outliers (Fig. 2E). Similar results were obtained when varying the window size to include 15, 50, 100, or 500 variants (table S7 and figs. S6 and S7).

Comparing observed levels of divergence with neutral coalescent simulations (6) under a range of possible demographic models revealed that the top 1% of observed FST values (approximately FST ≥ 0.25) are always higher than the corresponding neutral FST values from simulations, consistent with divergent selection acting on approximately this top 1% of variant sites (figs. S8 and S9).

We identified genomic regions with observed benthic-littoral FST ≥ 0.25 (i.e., with FST above maximum levels seen in neutral simulations) (6). For this we used windows of 15 variants each, providing a balance between fine genomic resolution and reducing stochastic variation by averaging over variants. Small gaps arising from brief dips of FST below the threshold were eliminated by merging regions within 10 kb of one another. We found 344 such regions, with total length of 8.1 Mb (~1% of the genome). Next, to focus on the more significant outliers, we narrowed the list down to a set of 98 highly diverged regions (HDRs) for further characterization (table S8) by adding the requirement that at least one 10-kb window must have reached FST ≥ 0.3. The HDRs varied in length from 4.4 kb to 285 kb (median, 36.1 kb), with a total length of 5.5 Mb.

A key prediction of speciation with gene flow models is that loci participating in speciation should have both high relative divergence (FST) and high absolute sequence divergence (dXY) (19, 20). However, previous studies [examined in (20)] revealed low dXY and low nucleotide diversity (π) in regions of high FST. In contrast, we found that dXY in Massoko is significantly higher in HDRs relative to the rest of the genome (P < 2.2 × 10−16, two-tailed Mann-Whitney test; Fig. 3A). In the benthic ecomorph, π in HDRs is not significantly different from the rest of the genome (P = 0.34, two-tailed Mann-Whitney test; fig. S10A), and in the littoral ecomorph π in HDRs is elevated (P = 5.47 × 10−6, two-tailed Mann-Whitney test; fig. S10B). Individually, 55 HDRs have dXY above the 90th percentile of the genome-wide distribution. The convergence of FST and dXY measures in regions of normal π suggests that these 55 “islands of speciation” (table S9) may have been involved in reducing gene flow in sympatry and thereby directly causing speciation to progress. In contrast, loci involved in continuing local adaptation after the ecomorphs split sufficiently to constitute two largely separate gene pools (or during a period of allopatry) would be expected to have elevated FST but not dXY (20).

Fig. 3 Islands of speciation between benthic and littoral ecomorphs.

(A) Elevated dXY in HDRs. (B) Clustering of putative speciation islands on five linkage groups. (C) Nucleotide diversity (π) within HDRs (red points) and outside HDRs (blue with shading corresponding to density). Each point corresponds to a 10-kb window (therefore, there may be multiple points per HDR). Overall, 95% of observations lie between the two curves (y = x ± 4.1 × 10−4). Putative sweeps in the benthic ecomorph are in the top left corner, and putative sweeps in the littoral are in the bottom right corner. (D) Patterns of FST, dXY, and πdiff in a speciation cluster on scaffold 15.

Another key prediction of speciation with gene flow models is that loci causing speciation should be located in relatively few linked clusters within the genome (2, 20, 21). Instead of a large number of scattered islands, the theory predicts a smaller number of clusters that grow in size because of the “divergence hitchhiking” process. We tested this prediction using a recently generated linkage map (22) and found that at least 27 out of the 55 putative speciation islands are co-localized on five linkage groups (LGs), with 26 of them clustered within their respective LGs (Fig. 3B and table S9). These potential speciation clusters extended for approximately 25 cM on LG5, 40 cM on LG7, 30 cM on LG12, 5 cM on LG20, and 45 cM on LG 23. In total, these regions account for under 7% of the genome, suggesting that divergence hitchhiking may play a role in shaping the observed pattern of genomic differentiation.

Although genomic islands within these clusters are often separated only by a few hundred kilobases, FST divergence between HDRs generally drops to background levels (Fig. 3D), with one exception on scaffold 88 ,where a broader “continent” of divergence has formed (fig. S11).

Further support for sympatric divergence

We next tested whether the HDRs correlated with the signal of gene flow into Lake Massoko, as identified using the sample from the nearby Mbaka River. Compared with the rest of the genome, the HDRs do not have elevated values of Patterson’s D (P = 0.22, two-tailed Mann-Whitney test; fig. S12C), nor elevated f statistics, which were recently proposed as an means by which one could identify introgressed loci (6, 23) (P = 0.08, two-tailed Mann-Whitney test; fig. S12D). These results suggest that introgression from the Mbaka River did not play a major role in generating the HDRs between the benthic and littoral ecomorphs within Lake Massoko (fig. S12), and strengthen the evidence that the ecomorph divergence has taken place within the lake.

Divergent SNPs associated with mate choice

Many recently diverged taxa, particularly those not geographically isolated, show stronger premating isolation than postmating isolation (1, 24, 25). We carried out laboratory experiments to test for reproductive isolation resulting from direct mate choice between the Massoko ecomorphs (6) (movie S2). Fifty Massoko females were given a choice from 16 males representing the variety of male phenotypes. In parallel, we designed a single-nucleotide polymorphism (SNP) assay with 117 polymorphic sites representing 44 HDRs identified from the first 12 genomes sequenced (table S10).

We genotyped a reference sample of 18 benthic and 16 littoral males, demonstrating that the SNP assay can reliably separate the ecomorphs along the first principal component (PC1) in PCA (Fig. 4A, top). We then genotyped all females and males participating in the mate choice experiments (Fig. 4A, bottom) and calculated an average of the PC1 distances between each female and the males she mated with during the experiment, as assayed by microsatellite paternity analysis (6). Compared with expectation under random mating (6), females had a moderate, but significant (P = 4.3 × 10−5, paired t test), preference for mating with males more genetically similar to themselves (i.e., closer to them along PC1) (Fig. 4B), demonstrating direct association between HDR variants and mate choice. Assortative mating by genotype was strong among females with positive (littoral) PC1 scores (P = 5.9 × 10−9, paired t test), whereas no assortative mating was detected among females with negative (benthic) PC1 scores (Fig. 4B).

Fig. 4 Mate choice trials.

(A) PCA based on 117 genotyped SNPs. (Top) The first axis of variation (PC1) in PCA reliably separates benthic and littoral males in a reference sample. (Bottom) PC1 positions of females (n = 50) and males (n = 16) participating in mate- hoice trials. (B) Results: Each point compares the average of absolute PC1 distances between a female and the males she mated with (observed PC1 distance) and all males she could have mated with (expected PC1 distance). Points are colored according to the PC1 score of the female. Females below and to the right of the dashed diagonal line on average mate with males more like themselves in terms of PC1 score than would be true if they mated at random.

Stronger mating discrimination by ancestral populations as compared to derived ones has been previously found in Drosophila and sticklebacks, possibly because low population density after a founder event favors less choosy individuals (26). However, it is also possible that the benthic ecomorph mates assortatively only in the deep-water environment; given that our experiments used wide-spectrum lighting characteristic of shallow water. Overall, the moderate assortative mating suggests a role for sexual selection in ecomorph divergence but does not indicate that it is a primary force causing population-wide divergence.

Signals of adaptation

A reduced level of genetic polymorphism in one subpopulation may be indicative of a recent selective sweep. Overall, the magnitude of difference in nucleotide diversity (π) between benthic and littoral ecomorphs (πdiff) is significantly higher in the HDRs than in the rest of the genome (P < 2.2 × 10−16, two-tailed Mann-Whitney test; fig. S13A) (6). Individually 46 HDRs have πdiff above the 95th percentile of the genome-wide distribution and are likely to have been under recent positive selection in one of the two ecomorphs. There is a significant overlap between HDRs with high dXY (putative speciation islands) and HDRs with high πdiff (putative recent selective sweeps): 35 of 55 high-dXY islands also have high πdiff (fig. S13B; P = 3 × 10−5, hypergeometric test). On the other hand, the 11 putative sweeps that did not lead to elevated dXY are indicative of adaptation not directly involved in reproductive isolation. Reduced nucleotide diversity in high πdiff regions, indicative of selective sweeps, was significantly more prevalent in the benthic ecomorph (36 of 46; P < 1.6 × 10−4, two-tailed binomial test; Fig. 3C, top left, and fig. S13C) (6), consistent with the benthic ecomorph being derived and undergoing more extensive adaptation. Nevertheless, there are also a few strong outliers, suggesting selective sweeps in the littoral ecomorph (Fig. 3C, bottom right).

Functions of adaptation

To explore the function of candidate adaptive genes, we performed gene ontology (GO) enrichment analysis (6) on three sets: (i) genes in candidate islands of speciation ±50 kb (enriched terms in table S11); (ii) genes in all HDRs ±10 kb (table S12); (iii) genes in all HDRs ±50 kb (table S13). Combining results of all three analyses in a network (Fig. 5A) connecting GO terms with high overlap (i.e., they share many genes) revealed clear clusters of enriched terms related to (i) morphogenesis (e.g., cartilage and pharyngeal system development or fin morphogenesis), consistent with morphological differentiation; (ii) sensory systems (e.g., photoreceptor cell differentiation), consistent with previous studies showing the role of cichlid vision in adaptation and speciation (7, 27); and (iii) (steroid) hormone signaling.

Fig. 5 Characterizing function of genes in HDRs.

(A) Enrichment map for significantly enriched GO terms. The level of overlap between GO enriched terms is indicated by the thickness of the edge between them. The size of the node indicates the best P value for the term, and the color of the node indicates the gene group for which the term was found significant (i.e., has P < 0.05 in candidate speciation islands ±50 kb, blue; in all HDRs ±10 kb or ±50 kb, red; or in both groups, green). Broad functional groupings (morphogenesis, sensory systems, etc.) were derived using automatic clustering followed by manual editing. (B) The absorption spectrum of the H5 allele of rho, more prevalent in the benthic ecomorph, is shifted toward blue wavelengths. (C) The joint roles of rho, rdh5, and rp1l1b in photoreceptor rod cells. rdh5 produces the chromophore 11-cis-retinal that binds rho, whereas rp1l1b, located at the axoneme of the outer segment and connecting cilia, also contributes to photosensitivity.

We examined in more detail the functions of candidate genes involved in photoreceptor function (table S14) and two highly diverged alleles of the rhodopsin (rho) gene in Lake Massoko (alleles H4 and H5, separated by four amino acid changes; FST = 0.39; fig. S14). Blue-shifted rhodopsin absorption spectra are known to play a role in deep-water adaptation (27). Therefore, we expressed rhodopsins from H4 and H5 alleles and reconstructed them with 11-cis-retinal, measured their absorption spectra (6), and demonstrated that the H5 allele, associated with the deep-water benthic ecomorph, has a blue-shifted absorption spectrum (Fig. 5B). The retina-specific retinol dehydrogenase rdh5 (table S14) produces 11-cis-retinal, the visual pigment binding partner of rhodopsin (28) and thus is likely to have a direct role in dark adaptation. Finally, a mouse ortholog of rp1l1b affects photosensitivity and morphogenesis of the outer segment (OS) of rod photoreceptor cells, locating to the axoneme of the OS and of the connecting cilia (29) (Fig. 5C). Together, these results suggest that divergent selection on rho, rdh5, and rp1l1b may facilitate the adaptation of scotopic (twilight) vision to the darker conditions experienced by the benthic ecomorph.

Comparisons to other systems

Overall, our results suggest a pair of incipient species undergoing divergence with gene flow within crater lake Massoko. Their overall level of divergence (FST = 0.038) is low compared with background FST observed in other recent studies of speciation with gene flow in Anopheles mosquitoes (S and M forms; FST = 0.21) (20), Ficedula flycatchers (FST = 0.36) (30), and Heliconius butterflies (FST = 0.18) (31), highlighting that we are looking at an early stage of divergence. The MSMC analysis suggests that median effective divergence occurred within the past 500 to 1000 years (~200 to 350 generations), after the separation of lake fish from the Mbaka River population around 10,000 years ago (fig. S5). However, divergence may have started considerably earlier than these times and be masked by subsequent gene flow.

Among populations at similar levels of divergence to that of the Lake Massoko ecomorphs are Timema stick insects (FST = 0.015 for adjacent and FST = 0.03 for geographically isolated population pairs), in which thousands of regions of moderately elevated divergence were found all across the genome (32), and German carrion and Swedish hooded crows (FST = 0.017), which have strongly diverged with fixed differences, but at fewer than five loci (33). In Massoko, we observe an intermediate pattern between these two extremes, with a few dozen moderately elevated islands clustering within the genome, indicating close linkage and no fixed differences. A genome-wide pattern with multiple loci of moderate divergence suggests a genomic architecture similar to the ecological divergence of a sympatric threespine stickleback pair in Paxton Lake, Canada (34), and the sympatric divergence of dune-specialist sunflowers, Helianthus (35).

The ecomorphs of Lake Massoko show clear differences in traits normally associated with adaptive radiation in cichlid fishes, including body shape, pharyngeal jaw morphology, diet, microhabitat preference, retinal pigment sensitivity, male color, and mate preference (3, 4, 11, 12, 27). Therefore, our study suggests processes and specific genomic regions to investigate to determine whether they are involved in speciation events within the great cichlid radiations of Lakes Malawi, Victoria, and Tanganyika.

Supplementary Materials

Materials and Methods

Figs. S1 to S16

Tables S1 to S16

References (3667)

Movies S1 and S2

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: This work was funded by Royal Society–Leverhulme Trust Africa Awards AA100023 and AA130107 (M.J.G., B.P.N. and G.F.T.), Wellcome Trust Ph.D. studentship grant 097677/Z/11/Z (M.M.), Wellcome Trust grant WT098051 (S.S. and R.D.), Wellcome Trust and Cancer Research UK core support and a Wellcome Trust Senior Investigator Award (E.A.M.), Leverhulme Trust Research Fellowship RF-2014-686 (M.J.G.), a University of Bristol Research Committee award (M.G.), a Bangor University Anniversary Ph.D. studentship (A.M.T.), and a Fisheries Society of the British Isles award (G.F.T.). Raw sequencing reads Sequence Read Archive: RAD sequencing (BioProject: PRJNA286304; accession nos. SAMN03768857 to SAMN03768912) and whole-genome sequencing (BioProject PRJEB1254: sample accessions listed in table S16). The RAD-based phylogeny and alignments have been deposited in TreeBase (TB2:S18241). Whole-genome variant calls in the Variant Call Format (VCF), phylogenetic trees, and primer sequences for Sequenom genotyping are available from the Dryad Digital Repository ( R.D. declares his interests as a founder and non-executive director of Congenica Ltd., that he owns stock in Illumina from previous consulting, and is a scientific advisory board member of Dovetail Inc. We thank R. Schley for generating pharyngeal jaw data; S. Mzighani, J. Kihedu, and the staff of the Tanzanian Fisheries Research Institute for logistical support; A. Smith, H. Sungani, A. Shechonge, P. Parsons, J. Swanstrom, G. Cooke, and J. Bridle for contributions to sampling and aquarium maintenance; the Sanger Institute sequencing core for DNA sequencing; and H. Imai (Kyoto University) for the use of the spectrometer in his laboratory.
View Abstract

Navigate This Article