Research article
Hierarchical structure of ecological and non-ecological processes of differentiation shaped ongoing gastropod radiation in the Malawi Basin
Abstract
Ecological processes, non-ecological processes or a combination of both may cause reproductive isolation and speciation, but their specific roles and potentially complex interactions in evolutionary radiations remain poorly understood, which defines a central knowledge gap at the interface of microevolution and macroevolution. Here I examine genome scans in combination with phenotypic and environmental data to disentangle how ecological and non-ecological processes contributed to population differentiation and speciation in an ongoing radiation of Lanistes gastropods from the Malawi Basin. I found a remarkable hierarchical structure of differentiation mechanisms in space and time: neutral and mutation-order processes are older and occur mainly between regions, whereas more recent adaptive processes are the main driver of genetic differentiation and reproductive isolation within regions. The strongest differentiation occurs between habitats and between regions, i.e. when ecological and non-ecological processes act synergistically. The structured occurrence of these processes based on the specific geographical setting and ecological opportunities strongly influenced the potential for evolutionary radiation. The results highlight the importance of interactions between various mechanisms of differentiation in evolutionary radiations, and suggest that non-ecological processes are important in adaptive radiations, including those of cichlids. Insight into such interactions is critical to understanding large-scale patterns of organismal diversity.
1. Introduction
The completion of reproductive isolation between populations is a key step during speciation in sexually reproducing Eukaryota [1]. Such isolation is obtained once gene flow is fully interrupted, and it may be prezygotic and often reversible, e.g. related to geographical barriers or assortative mating [2], up to postzygotic, intrinsic and irreversible, e.g. because of Dobzhansky–Muller incompatibilities [1,3]. The processes involved determine to a large extent the strength and stability of the arising reproductive barrier(s) through time [3]. As such, insights into the nature of reproductive isolation and its genomic basis are prerequisites to understanding speciation [4,5].
Ecological speciation, i.e. the development of reproductive isolation between populations as a direct or indirect result of adaptation to different environments, has been examined extensively [3], but non-ecological ‘drivers’ of population differentiation and speciation may be common too [4]. Prolonged geographical isolation in ecologically similar environments may lead to the fixation of alternate, incompatible mutations in allopatric populations under uniform selection (mutation-order speciation), or genetic drift could facilitate differentiation and speciation, at least when founding populations are small [4,6–9]. Ecological and non-ecological processes can interact in intricate ways during speciation, both in space and time [10,11], and as to their net genomic effects, i.e. their relative contributions to changes in allele frequencies [4,12]. Isolating the contributions of neutral and selective processes to population differentiation is often difficult [1,4,12,13], certainly in retrospect as it requires separating reproductive barriers that initiated and eventually completed speciation from those that evolved afterwards [4,11]. Incipient species may be exceptionally informative in this respect [13,14], and those of ongoing radiations may additionally provide important insights into the roles of various speciation mechanisms in evolutionary radiation, which represents an important knowledge gap at the interface of microevolution and macroevolution [15,16].
Here I document mechanisms of population differentiation and speciation in an ongoing gastropod radiation from Lake Malawi. This study represents—to my knowledge—the first in depth examination of such processes in an invertebrate radiation that coincided with a renowned adaptive radiation of cichlid fishes [17–19]. Specifically, I performed integrative studies on genomic and morphological data in an ecological context to examine how ecological and non-ecological processes have affected population structure in the emergent species flock. I found a remarkable hierarchical structure of older, non-ecological processes on a large geographical scale and younger, ecological processes on a regional scale. Substantial differences exist in the relative contribution of these processes to the genome-wide differentiation in each speciation event. The work illustrates that various processes act simultaneously in evolutionary radiations and that their interaction in space and time is fundamental to the prospects of radiation. Overall, genetic drift and mutation-order processes have been essential in causing genome-wide differentiation in this gastropod radiation, although, the occurrence of ecological processes may give the impression that it represents a ‘full-blown’ adaptive radiation.
2. Study system
Lake Malawi in the East African Rift is one of the most species-rich freshwater systems in the world [18], and its climatological and environmental history is becoming increasingly well documented [20–22]. Here I focus on the ongoing Lanistes radiation in the Malawi Basin, which comprises five nominal species that exhibit differentiation at various levels of organization (electronic supplementary material, text and figure S1) [23–25]. Two of these species were previously assigned to lineages occurring elsewhere in Africa based on morphological resemblance, but they are genetically distinct [26], and I refer to them as L. sp. (ovum-like) and L. sp. (ellipticus-like). They occur in the lake on shallow sandy to muddy substrates in association with shoreline vegetation, but mainly in shallow fringing pools, swamps, satellite lakes and rivers [25]. The other three species, Lanistes solidus, Lanistes nyassanus and Lanistes nasutus, are restricted to soft substrates within Lake Malawi and have presumably evolved intra-lacustrine [23,27]. The habitats of L. solidus and nyassanus overlap in water depth from approximately 1–35 m, whereas L. nasutus is rare and lives at 40–90 m depth [23,25,27,28]. Behavioural, anatomical and morphological differences among these species have been interpreted as adaptations to (i) wave action and water energy (shell morphology), (ii) the presence of predators (shell thickness, variation in growth rates, burying behaviour, development of a nocturnal lifestyle), and (iii) substrate-specific food sources (differences in radulae) [23]. For certain traits, e.g. shell morphology, specimens with intermediate characteristics exist, which was attributed to the recent (several 10 000s up to 100 000 years), still ongoing differentiation in the Lanistes radiation [23]. This timeframe of diversification is substantiated by molecular data [26] and the rising water level after severe megadroughts approximately 75–135 ka ago [20,21], which has also triggered diversification in Malawian cichlids [21,29] and viviparid gastropods [30,31].
3. Material and methods
(a) Sampling
Specimens were collected in 2006, 2007 and 2008 from soft substrates in the Malawi Basin, as Lanistes does not occur on rocks [24,25] (additional sampling information is provided in the electronic supplementary material, text, tables S1 and S2). The samples contain all abovementioned species, except for L. nasutus, which is rare and presumably strongly declining in abundance [28]. Shell and tissue vouchers are stored at the Biodiversity and Systematics collection of the Justus Liebig University Giessen and the Natural History Museum London. Individuals were identified to the currently valid nominal Lanistes species of Lake Malawi following relevant literature descriptions [24,25]. Environmental conditions at collecting sites were characterized with five variables that are relevant in shaping ampullariid habitats [24], i.e. substrate composition, temperature and the concentrations of dissolved oxygen, phosphate and calcium. Habitat diversity was examined with non-metric multidimensional scaling (nmMDS) and model-based clustering using multivariate normal mixture models (see below).
(b) Survey of genome-wide polymorphism
Technical data on DNA extractions and amplified fragment length polymorphism (AFLP) protocols are provided in the electronic supplementary material, text. In total, 14 of the 60 primer combinations tested for selective PCR yielded clearly readable profiles (electronic supplementary material, table S3). AFLP bands were digitally captured with e-Seq v. 2.0, and automatically scored with SagaMX v. 3.2 (both from LI-COR). Scoring results for all loci were checked by eye without prior knowledge of population affiliation. Only unambiguous polymorphic bands ranging 100–350 bp were selected for subsequent analysis to mitigate the amount of potentially homoplasic bins. To estimate the quality of the generated data the entire AFLP protocol was repeated with a random subset of 10% of the specimens following ref. [32]. This test resulted in an overall reproducibility of greater than 95% (electronic supplementary material, table S4), indicating that the method is robust and the obtained dataset reliable.
(c) Genetic analyses
To examine the effects of selective processes, I compared the coefficients of population differentiation for all loci with those expected under neutral evolution in BayeScan v. 2.1 [33] (20 pilot runs, sampling = 100 000 generations, burn-in = 50 000 generations, thinning interval = 10). Loci that have significantly higher or lower differentiation than expected under neutral evolution were considered under selection, whereas the others were considered neutral for further analyses.
Pairwise population differentiation was quantified with Wright's (locus-specific) FST following ref. [34] for selectively neutral loci and all loci with AFLP-SURV v. 1.0 [35]. Null allele frequencies were calculated with a Bayesian approach using a non-uniform prior distribution [36], which tolerates moderate departures from Hardy–Weinberg equilibrium [37]. Nei's genetic diversity (Hj), which is equivalent to the expected heterozygosity across loci, and the global FST were also calculated with AFLP-SURV following ref. [34].
Patterns of molecular differentiation were identified from all markers with discriminant analysis of principal components (DAPC) using Adegenet v. 1.4 [38] in R v. 3.2.1 [39] and with Structure v. 2.3.4 [40]. DAPC was performed with the output of K-means clustering based on a Bayesian information criterion (BIC). As no prior information on the degree of gene flow among populations is available I used an admixture model with correlated allele frequencies in Structure, which allows for individuals having inherited a fraction of their genome from each of the K putative clusters. Analyses were conducted for K = 1–10, with 10 independent runs per K. Each Markov chain ran for 100 000 generations following a burn-in of 50 000 generations. The underlying genetic structure was examined with the ΔK statistic [41] in ‘Clumpak’ [42].
(d) Morphometric analyses
The specimens included in AFLP analyses were subjected to semilandmark analyses as described in ref. [43]. I used mainly adult specimens with well-preserved shells for molecular analyses, but in some cases juveniles or specimens with damaged shells were included to increase sampling sizes. Such specimens (n = 18; <10% of the dataset; electronic supplementary material, table S1) were omitted from morphometric analyses, because allometric differences and/or ambiguous reconstructions of shape would inflate the morphological variability of populations. Because Lanistes displays indeterminate growth (although growth is slower in adults) I further reduced the impact of allometric variation by focusing only on shape differences. Voucher images were compiled and digitized in tpsDig v. 2.16 [44] with 11 landmarks (LMs) and two curves of equally spaced semilandmarks (electronic supplementary material, figure S2). Individuals were aligned using a baseline (LMs 1 and 2) and subjected to generalized Procrustes superimposition in CoordGen6h [45]. Subsequently, semilandmarks were slid along their respective curve using the minimum Procrustes distance criterion in SemiLand6 [45]. These coordinates were imported in R, converted to a Euclidean distance matrix and subjected to nmMDS in two dimensions using Vegan v. 2.3-0 [46] and Mass v. 7.3-41 [47]. The goodness of fit of the resulting ordination was examined using Kruskal's and Clarke's rules of thumb [48,49]: stress values <5 indicate excellent ordination, those less than 15 good ordination, with limited chances of false inferences. I examined shape variation along the nmMDS axes in GeoMorph v. 2.1.6. [50], and the contribution of individual morphometric variables to the morphospace occupation with vegan's envfit function based on 1000 permutations. The morphospace occupation was subjected to model-based clustering using multivariate normal mixture models (n = 14) in Mclust v. 5.0.1 [51] without a priori classification [30]. These models make different assumptions as to the expected shape, volume and orientation of trait variance and covariance, and support is evaluated with a BIC. Statistically preferred clustering solutions may not be biologically realistic, and the reliability of the outcome depends on the representation of ‘groups’ in the dataset (discussed in ref. [30]).
Subsequently, I estimated the degree of phenotypic differentiation, PST, as an analogue to QST [52]. If PST serves as a reasonable approximation of QST, PSTs can be compared with FSTs calculated from the neutral loci to examine the roles of genetic drift and selection in the observed differentiation in phenotypic traits. Beyond genetic factors, PST can be influenced by environmental and non-additive genetic effects, and, hence, the accuracy with which QST can be approximated by PST requires evaluation [53].
I calculated PSTs for pairwise population comparisons directly from the morphospace occupation using the formula: 
(eqn. 3 in ref. [53]). The phenotypic variances between and within populations (
and 
, respectively) were calculated from the mean squares generated by an analysis of variance [54], which is unbiased by sampling effects and represents a conservative way of calculating PSTs. Two parameters influence the approximation: c, the proportion of total genetic variance owing to additive genetic effects across all populations, and heritability, h2. Values for these parameters are difficult to obtain for populations from the wild. Therefore, I evaluated the sensitivity of my analyses for a wide range of scenarios for c/h2 (i.e. from 0.0 to 2.0 in steps of 0.1). Specifically, I bootstrapped the original Procrustes superimposition coordinates and repeated ordinations in morphospace and PST calculations 1000× for each scenario of c/h2. The resulting 95% confidence interval (CI) on pairwise PST comparisons between populations was compared with the 95% CI on the global FST for all neutral loci [53]. Additionally, I tested at what values of c/h2 the corresponding PST distribution significantly exceeds the average locus-specific FSTs from all neutral loci with Wilcoxon rank-sum tests [55]. Shell morphology was considered to be one multidimensional trait for these tests (although it arguably fulfils several functions; see §2).
(e) Inference of mechanisms of differentiation
I examined processes of geographical differentiation in the Lanistes radiation by estimating the pairwise relatedness between individuals with a kinship coefficient, Fij, for dominant markers (based on all loci) [56]. These calculations require information on the inbreeding coefficient, FI, which is not available for Lanistes from the Malawi Basin. Hence, I repeated the calculations for scenarios of Hardy–Weinberg equilibrium up to substantial inbreeding (FI = 0.0 to 0.6, respectively, in steps of 0.2). This strategy is reasonable because deviations from the real, unknown value with up to 0.2 would bias kinship estimates less than 15% [56]. Kinship coefficients were regressed on the natural logarithm of spatial distances among localities in SPAGeDi, v. 1.4 [57] for the entire basin and for individual regions. Summary statistics are based on tests with 1000 permutations and jackknife procedures.
Genetic differentiation among populations (FST) may arise owing to selection, neutral processes or both. Given that PST reasonably approximates QST for the Lanistes dataset (see §4), I examined signals of adaptation by comparing pairwise population PSTs with locus-specific FSTs (for neutral loci), using the SPAGeDi results to account for geography. If the PST of a trait (nmMDS1) is significantly larger than the corresponding neutral FST, genetic drift alone is not sufficient to explain the observed differentiation, implying that diversifying selection occurred; if the observed PST is significantly lower than the neutral FST stabilizing selection is implied [55].
4. Results and discussion
(a) Environmental analysis indicates two habitat types
nmMDS of the sampling localities in environmental space resulted in very limited stress (=2.08), and permutation tests indicate that three variables (substrate composition, temperature and dissolved oxygen) contribute significantly to the ordination (figure 1). Model-based clustering without a priori group assignments reveals that localities cluster into two main groups or benthic habitats along nmMDS1, which relate strongly to differences in environmental stability and are labelled ‘stable’ and ‘fluctuating’ in what follows (more information on the examination of habitats is provided in the electronic supplementary material, text).
Figure 1. Non-metric multidimensional scaling (nmMDS) of the 15 sampling localities based on five environmental variables that characterize Lanistes habitats. Three of these variables, i.e. substrate composition, temperature and dissolved oxygen, contribute significantly to the environmental differentiation between localities. Model-based clustering without a priori group assignments indicates that sampling localities fall into two main groups or habitats (support of two-group model over models with 1, 3 or 4 groups: ΔBIC ≥ 1.69). Both habitats differ mainly in their degree of environmental stability, and are categorized as ‘stable’ and ‘fluctuating’. The mclust model identifiers in the inset are explained in ref. [51].
(b) Molecular diversity and differentiation
In total 201 AFLP loci were scored reliably for 182 Lanistes individuals from 15 populations (electronic supplementary material, table S1). BayeScan analysis indicates that 14 loci (6.97%) had α values indicative for diversifying (positive) selection (false-discovery rate (FDR) = 0.05) and these are subsequently referred to as outlier loci (electronic supplementary material, figure S3). Seven of these markers are particularly strong outliers (FDR < 0.001). As none of the examined loci show evidence for balancing selection, the other 187 loci are considered neutral. This proportion of outliers is very similar to what has been reported in studies of genome-wide differentiation associated with (local) adaptation and speciation (approx. 2–8%) [58–62].
Within-population genetic diversity ranges from 0.194 to 0.287 (for neutral markers; electronic supplementary material, table S1) and is significantly larger in the southern than in the northern region, and in stable than fluctuating habitats (Wilcoxon rank-sum test, W = 9, p = 0.029; W = 1, p < 0.001, respectively). The global FST of the dataset is 0.0975 (95% CI: 0.0313–0.1637) for neutral loci and 0.1232 (95% CI: 0.0732–0.1732) for all loci. These values are similar to those obtained in the study of speciation in Timemia and cichlids [12,63].
Examination of patterns of molecular diversity with DAPC and k-means clustering yields the best support for four molecular groups, each with an individual assignment probability of 1.00 (figure 2a). Other models have substantially higher BIC values (ΔBIC ≥ 1.26), except for a five-group model (ΔBIC = 0.05). Analysis in Structure indicates the highest ΔK for two groups (102.3 whereas ≤ 22.4 for other scenarios; electronic supplementary material, figure S4). Although different in the ‘preferred’ number of groups, DAPC and Structure provide almost identical attributions of specimens to clusters for scenarios of 2–5 groups. Both methods support a strong genetic separation between populations from the northern and southern regions of the Malawi Basin. The four-group DAPC model additionally recognizes two groups in each region. (The two-group Structure result can thus be reconstructed directly from the four-group DAPC model by regional lumping and, as such, I use DAPC for visualization.) Individuals from the two northern groups (I and II) are regularly found in sympatry, whereas those from the southern groups (III and IV) occur almost exclusively in allopatry (figure 2b). Groups I and II are very similar in their composition of (morphologically identified) nominal species, whereas groups III and IV consist almost exclusively of L. solidus and L. nyassanus and of L. sp. (ovum-like) and L. solidus, respectively (electronic supplementary material, figure S5). The differences between groups within regions are more limited than those between regions: L. sp. (ellipticus-like) and L. nyassanus occur only in the northern and southern region, respectively. In the five-group DAPC model the abovementioned group IV is further split into two parts: one part contains the specimens from Lake Chilingali (am0815), the other part the rest. Analyses of patterns of molecular diversity are important to provide biological meaning to the results, but I did not rely on specific grouping scenarios in subsequent analyses, beyond analysing patterns by region following the strong interregional molecular differentiation.
Figure 2. Genetic diversity and morphological disparity in the ongoing Lanistes radiation of the Malawi Basin. (a) Discriminant analysis of principal components (DAPC; n = 182) gives the best BIC support for four genetic groups, with a strong separation between the northern and southern regions of the basin. (b) Lake Malawi and surroundings with indication of sampling localities, colour-coded genetic diversity at these localities, bathymetry and shoreline substrates (after ref. [22]). In the north both genetic groups occur highly sympatric but in the south not. The dashed white outline at the centre of the rift indicates the approximate extent of the lake during Late Pleistocene low stands [20,21]. The morphospace occupation (n = 164) indicates an incomplete separation of molecular groups and of specimens from the northern and southern regions of the basin (c), but a strong separation as to habitats (d). (d) The nominal species occupy somewhat separate areas in morphospace although overlaps are substantial, especially between L. sp. (ellipticus-like) and L. sp. (ovum-like), which is inflated by allometric effects (young adults of L. sp. (ovum-like) resemble L. sp. (ellipticus-like)). Colours reflect habitats, symbols the nominal species; electronic supplementary material, figures S6 and S7 present additional representations of morphological variability.
(c) Morphological differentiation
Ordination of semilandmark data with nmMDS indicates a stress of 14.87, suggesting limited risk of drawing false inferences. The morphospace distribution (figure 2c,d and electronic supplementary material, figure S6) resembles a single ‘supercluster’ (ΔBIC ≥ 2.51 over model-based clustering solutions with two or more clusters), even though specimens occupy different parts of the morphospace depending on the nominal species to which they belong and the region and habitat they occupy. Substantial differences in morphospace occupation also exist between molecular groups, except for between groups I and II (electronic supplementary material, figure S7 details shape differences along both axes). If four clusters are enforced—given that the dataset comprises four valid, nominal species and following the results of the DAPC analysis—the resulting morphological groups coincide well with the nominal species (electronic supplementary material, figure S6) rather than with the molecular groups. Two-way multivariate analysis of variance (MANOVA) indicates significant differences in morphospace occupation by region (Wilks' λ = 0.81, F2,159 = 18.51, p < 0.001) and by habitat (Wilks’ λ = 0.31, F2,159 = 176.68, p < 0.001) but not by their interaction (Wilks' λ = 0.98, F2,159 = 1.84, p = 0.162). The morphospace differentiation between habitats is highly significantly larger than between regions (mean ± standard deviation: 0.0430 ± 0.0031 and 0.0419 ± 0.0035, respectively; Wilcoxon rank-sum test on 1000 bootstrap replicates: W = 404 500, p < 0.001), as can be observed by comparison of figure 2c,d.
(d) QST can be approximated by PST
For most of the relevant interval of c/h2 (i.e. 0.15–2.00) mean PSTs exceed the global FST from neutral loci, and the 95% CIs on PSTs do not overlap with the upper confidence limit on the global FST for c/h2 ≥ 0.45 (electronic supplementary material, figure S8a). Statistical tests indicate that PSTs significantly exceed the average locus-specific FSTs of neutral loci for values of c/h2 ≥ 0.5 (Wilcoxon rank-sum test for c/h2 = 0.5: W = 7872, p = 0.002) (electronic supplementary material, figure S8b). Therefore, the PSTs for this dataset are indicative of moderate-to-strong phenotypic divergence, and they can be used to approximate QST at c/h2 ≥ 0.5. To be conservative, I compared PSTs and average locus-specific FSTs for c/h2 = 0.5 rather than for the null assumption of c/h2 = 1.0 [53]. Population pairwise PSTs and FSTs are reported in electronic supplementary material, table S5.
(e) Patterns of population differentiation
The regression of kinship coefficients on spatial distance yields a significant negative correlation for the entire dataset (r2 = 0.216, p < 0.001; figure 3a), but this trend disappears when regions are analysed separately. A weak yet statistically significant regression is obtained in the north (r2 = 0.013, p = 0.023; figure 3b) but not in the south (r2 = 0.005, p = 0.237; figure 3c). These results are robust throughout the range of reasonable FI values (FI = 0.0–0.6).
Figure 3. Mechanisms of genetic differentiation in the Lanistes radiation. (a–c) Relationship of kinship coefficients (using all loci) on spatial distances for the complete Malawi Basin (a), the northern region (b), and the southern region (c). The black curve assumes Hardy–Weinberg equilibrium, but results are robust for progressively increasing inbreeding (FI = 0.2, 0.4 and 0.6; grey curves). The significant negative regression for the entire dataset does not persist within individual regions, suggesting the existence of a strong geographical barrier to gene flow coinciding with the centre of rifting. (d–g) Comparisons of phenotypic differentiation (PST) and neutral genetic differentiation (FST). (d) Comparisons coded according to regions; confidence intervals (CIs) are only shown for pairwise population comparisons (n = 70) for which the 95% CI on PSTs and on FSTs do not overlap with the line of FST−PST equality (dotted dash line), thus indicating adaptation or stabilizing selection. (e) The same comparisons, but coded as to habitat properties. Jointly (d,e) indicate that most adaptation occurs between regions and between habitats, but to a lesser extent also within regions and habitats. A positive and significant correlation between phenotypic differentiation and neutral genetic differentiation (blue regression lines), i.e. isolation by adaptation (IBA), exists for the entire dataset (e), but IBA occurs exclusively within regions (f) not between regions (g). The FSTs associated with comparisons between regions fall in narrow zones (blue vertical bars in g) that indicate a phase of genetic admixture in the past.
PSTs are significantly higher than the associated FSTs for 55 pairwise population comparisons (approx. 52%) and significantly lower in 15 other comparisons (approx. 14%). Most adaptation occurs between populations from different regions (figure 3d) and different habitats (figure 3e), but also within regions and habitats evidence for adaptation exists (figure 3d,e and electronic supplementary material, figure S9). A significant, positive correlation is observed between phenotypic divergence and neutral genetic differentiation for the entire dataset (r = 0.262, p = 0.007). Accounting for geographical effects (figure 3a–c), strong correlations are obtained for PST−FST comparisons within regions (south: n = 28; r = 0.596, p < 0.001; north: n = 21; r = 0.462, p = 0.035, all regional comparisons: n = 49; r = 0.488, p < 0.001; figure 3f), but not between regions (n = 56; r = −0.187; p = 0.168).
(f) Processes of population differentiation
Remarkable differences exist in patterns of differentiation within and between regions of the Malawi Basin. Below, I examine evidence for ecological and non-ecological mechanisms of differentiation separately before producing an integrative model of differentiation in the Lanistes radiation. Non-ecological mechanisms include neutral processes and divergence under uniform selection in separate regions. The role of neutral processes is illustrated by the significant negative correlation between gene flow (or kinship) and geographical distance for the entire dataset (figure 3a). If isolation by distance (IBD) were the main explanation for geographical differentiation, we would also expect significant negative correlations within regions and especially in the southern region given the overall larger distance between sampling localities than in the northern region. As correlations within regions are weak to non-existing, a single migration barrier explains the geographical differentiation in the Lanistes radiation from the Malawi Basin best, although weak IBD may overlay it. Geological and ecological data support this interpretation: Lanistes is restricted to shallow, soft-substrate habitats in the Malawi Basin, and, when lake levels are high, as during the last approximately 60 ka [20,21], such habitats occur only along the northernmost and southernmost shores of the basin. At the centre of rifting, border faults coincide with the shorelines of Lake Malawi, both east and west [64,65], and these steep, rocky shorelines (figure 2b) are uninhabitable for Lanistes [24,66]. In combination with the anoxic bottom waters [28] they form a strong migration barrier for molluscs. The molecular consequences of geographical separation since the water level rose are recorded in the pairwise population FSTs across the barrier (figure 3g): most of these FSTs fall in a narrow band at FST = ∼0.1, representing divergence since the last major event of genetic admixture, i.e. the last persistent low stand (figure 4). The occurrence of another narrow band at FST = ∼0.2 may provide evidence for an older period of admixture, however, all these comparisons involve satellite Lake Chilingali (am0815). The Lanistes specimens from this satellite lake form the most genetically differentiated population in the southern region (electronic supplementary material, figure S9), but they also have the lowest intrapopulational genetic diversity, which indicates that founder effects and genetic drift may have enhanced differentiation. Lake Chilingali is well separated from Lake Malawi and, despite its shallow depth (approx. 5 m), it also harbours a strongly differentiated Rhamphochromis population [67].
Figure 4. (a) Conceptual model of genetic and morphological differentiation in the ongoing Lanistes radiation of the Malawi Basin over time. Lake Malawi changed from a shallow, green to a deep, blue lake, most recently after megadroughts approximately 75–130 ka ago [20,21], which triggered diversification in Lanistes. Solid lines represent periods during which neutral processes predominate whereas dotted lines indicate phases of elevated selection and adaptation. Per unit of time larger genetic and morphological differentiation may occur under selection than under neutral processes. As the lake level rose, differentiation in Lanistes was initiated by geographical separation of populations in the northern and southern regions. Periods of more stable, high lake levels subsequently allowed adaptation and diversification, but these periods were interrupted occasionally by environmental perturbations. Perturbations may have caused a decrease in genetic and morphological diversity, e.g. by secondary contact, hybridization and reticulate evolution (arbitrarily represented in the grey lineage) and/or the extinction of genetic lineages (arbitrarily represented in the black lineage). The extant fauna (colour-dotted lines) displays a hierarchical structure of differentiation mechanisms: older, non-ecological mechanisms (geographical differentiation and parallel evolution) occur mainly between regions, and younger, ecological mechanisms (isolation by adaptation) within regions. (b) Hypothesis of the spatio-temporal constraints of water level fluctuations in the Malawi Basin (based on ref. [20]) on the observed molecular and morphological differentiation in Lanistes, indicating the last period of major genetic admixture, the separation between regions and adaptation within regions. SEC, shoreline elevation constraints.
Whereas molecular differentiation is strongest between regions (figure 2a), morphological differentiation occurs mainly between habitats (figure 2d), although highly significant differences in morphospace occupation occur also by region (figure 2c). Given that both habitats occur in each region, and that the common ancestor of the Lanistes radiation necessarily preferred one habitat type (adaptation within regions postdates molecular differentiation between regions (see below)), the molecular and morphological data jointly provide evidence for parallel evolution (electronic supplementary material, figure S10). Parallel evolution suggests that similar selective regimes existed in both regions and that mutation-order processes may have been influential in the differentiation between regions. Evidence for stabilizing selection is found in PST−FST comparisons between regions (red symbols in figure 3d) that correspond to comparisons between fluctuating or fluctuating and stable habitats (red and black symbols in figure 3e). This finding also suggests that at least some populations are near optima in the adaptive landscape. In recently separated lineages, as is the case for Lanistes, the probability that the same genes are involved in independent parallel adaptation between regions is high [68], which may increase the chance for rapid evolution of incompatibilities.
Ecological processes in the ongoing Lanistes radiation are evidenced by population pairwise PST−FST comparisons with PSTs that are significantly higher than the corresponding FSTs (figure 3d–g), and the differentiation in each region of quantitatively discernible ecotypes (figure 2c,d and electronic supplementary material, figure S10) (see refs [52,69]). The strongest evidence for adaptation occurs in pairwise comparisons across regions and across habitats (red symbols in figure 3d,e), but only within regions evidence for isolation by adaptation (IBA) exists, i.e. a positive correlation between adaptive phenotypic divergence and neutral genetic differentiation [58]. IBA is observed for all neutral AFLP loci, indicating it to be a strong driver of neutral genomic differentiation within regions. Between regions adaptation is occurring, but IBA does not contribute much to the associated neutral genetic differentiation. In summary, ecological processes are responsible for most of the neutral genetic differentiation within regions, whereas neutral genome-wide differentiation between regions is mainly caused by non-ecological processes.
(g) Spatio-temporal model of radiation
FSTs (of neutral loci) for comparisons between regions are significantly larger than those for comparisons within regions (Wilcoxon rank-sum test: W = 610, p < 0.001; compare figure 3f,g), which indicates that the geographical separation between regions predates IBA within each region. This finding suggests that only during prolonged high lake level phases the ‘stable’ and ‘fluctuating’ habitats persisted for a sufficient period for IBA to occur. The data reveal a strong structure of differentiation processes in time and space, tightly linked to environmental change in the Malawi Basin (figure 4): as the water level rose after the last persistent low stand geographical separation occurred between regions, and only later when lake levels were high and stable, adaptation to both habitat types occurred within each region. The hierarchical structure of various differentiation mechanisms, rather than their mere presence, has promoted evolutionary radiation by allowing multiple speciation events to unfold simultaneously, resulting in additive differentiation as mechanisms interacted synergistically in space and time. Interestingly, the ‘stage model’ of diversification in Lake Malawi cichlids [17], and the role of environmentally triggered habitat isolation and reconnection as a potential ‘species pump’ [67], suggest that similar interactions of differentiation mechanisms occurred in cichlids, however, currently these models focus on adaptive processes only. Cichlids would have primarily diversified during deep-lake conditions, whereas intermittent turbid, low lake level phases reduced rocky habitats, decreased visibility and opportunities for assortative mating, and promoted hybridization and the collapse of some clades [21,70–72]. Here I demonstrate that ecological and non-ecological processes of differentiation can interact in a structured way, especially in systems where geographical ranges and habitat diversity have been waxing and waning over time. The implications of my findings are probably important for cichlids and other adaptive radiations. As most of these radiations are more complex than that of Lanistes in the Malawi Basin, separating the mechanisms leading to speciation from those that followed afterwards will be more complicated, however. These complications are exemplified by cichlid radiations, for which it is often difficult to identify sister species relationships, to control for possible introgression from non-sister taxa, and to reconstruct the geographical context [73,74]. My results also suggest that as divergence progresses, interactions of various processes may cause the role of non-ecological processes to be underestimated as ecological processes are more effective in obliterating non-ecological signals, e.g. via selective sweeps or the progressive expansion of genomic islands of differentiation, than vice versa. Interactions of microevolutionary processes are critical to the potential for radiation, and documenting such interactions is therefore essential to understand the relative contribution of various processes to large-scale patterns of organismal diversity. Genomic approaches and experiments are required to test the genome-wide effects of the various processes involved and to unravel how selection acted in the differential adaptation in Lanistes.
Ethics
This study results from research cooperation KM/1/1.64 between Ghent University, Belgium and the Karonga Museum, Ministry of Tourism, Wildlife and Culture, Malawi. Fellowship 12N3915N of the FWO Vlaanderen received clearance from the ethical committee on animal experimentation of Ghent University.
Data accessibility
Datasets are deposited in Dryad at doi: http://dx.doi.org/10.5061/dryad.dn58p [75].
Competing interests
I have no competing interests.
Funding
Research was funded by postdoctoral fellowship 12N3915N of the FWO Vlaanderen and DFG grants SCHR 352/9-1, AL 1076/6-2, 7-1, 8-1.
Acknowledgements
I thank the Cultural and Museum Centre Karonga, Harrison Simfukwe and Friedemann Schrenk for support, Christian Albrecht, Tom Wilke and Jon Ablett for access to collections and Matthias Schrader and Silvia Nachtigall for laboratory work. Xavier Vekemans, Thomas Neubauer and an anonymous reviewer provided constructive advice. Editors: John Pandolfi and Sasha Dall.
References
- 1
- 2
Taylor EB, Boughman JW, Groenenboom M, Sniatynski M, Schluter D, Gow JL . 2006 Speciation in reverse: morphological and genetic evidence of the collapse of a three-spined stickleback (Gasterosteus aculeatus) species pair. Mol. Ecol. 15, 343–355. (doi:10.1111/j.1365-294X.2005.02794.x) Crossref, PubMed, Web of Science, Google Scholar - 3
- 4
Sobel JM, Chen GF, Watt LR, Schemske DW . 2010 The biology of speciation. Evolution 64, 295–315. (doi:10.1111/j.1558-5646.2009.00877.x) Crossref, PubMed, Web of Science, Google Scholar - 5
Seehausen O 2014 Genomics and the origin of species. Nat. Rev. Genet. 15, 176–192. (doi:10.1038/nrg3644) Crossref, PubMed, Web of Science, Google Scholar - 6
Gavrilets S, Hastings A . 1996 Founder effect speciation: a theoretical reassessment. Am. Nat. 147, 466–491. (doi:10.1086/285861) Crossref, Web of Science, Google Scholar - 7
Wessel A, Hoch H, Asche M, von Rintelen T, Stelbrink B, Heck V, Stone FD, Howarth FG . 2013 Founder effects initiated rapid species radiation in Hawaiian cave planthoppers. Proc. Natl Acad. Sci. USA 110, 9391–9396. (doi:10.1073/pnas.1301657110) Crossref, PubMed, Web of Science, Google Scholar - 8
Rosenblum EB, Römpler H, Schöneberg T, Hoekstra HE . 2010 Molecular and functional basis of phenotypic convergence in white lizards at White Sands. Proc. Natl Acad. Sci. USA 107, 2113–2117. (doi:10.1073/pnas.0911042107) Crossref, PubMed, Web of Science, Google Scholar - 9
Schluter D . 2001 Ecology and the origin of species. Trends Ecol. Evol. 16, 372–380. (doi:10.1016/S0169-5347(01)02198-X) Crossref, PubMed, Web of Science, Google Scholar - 10
Nosil P, Harmon LJ, Seehausen O . 2009 Ecological explanations for (incomplete) speciation. Trends Ecol. Evol. 24, 145–156. (doi:10.1016/j.tree.2008.10.011) Crossref, PubMed, Web of Science, Google Scholar - 11
Rundell RJ, Price TD . 2009 Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends Ecol. Evol. 24, 394–399. (doi:10.1016/j.tree.2009.02.007) Crossref, PubMed, Web of Science, Google Scholar - 12
Nosil P, Gompert Z, Farkas TE, Comeault AA, Feder JL, Buerkle CA, Parchman TL . 2012 Genomic consequences of multiple speciation processes in a stick insect. Proc. R. Soc. B 279, 5058–5065. (doi:10.1098/rspb.2012.0813) Link, Web of Science, Google Scholar - 13
Butlin R, Debelle A, Kerth C, Snook RR, Beukeboom LW , The_Marie_Curie_SPECIATION_Network. 2012 What do we need to know about speciation? Trends Ecol. Evol. 27, 27–39. (doi:10.1016/j.tree.2011.09.002) Crossref, PubMed, Web of Science, Google Scholar - 14
Malinsky M 2015 Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake. Science 350, 1493–1498. (doi:10.1126/science.aac9927) Crossref, PubMed, Web of Science, Google Scholar - 15
Schluter D . 2001 The ecology of adaptive radiation. Oxford, UK: Oxford University Press. Google Scholar - 16
Simpson GG . 1953 The major features of evolution. New York, NY: Columbia University Press. Crossref, Google Scholar - 17
Kocher TD . 2004 Adaptive evolution and explosive speciation: the cichlid fish model. Nat. Rev. Genet. 5, 288–298. (doi:10.1038/nrg1316) Crossref, PubMed, Web of Science, Google Scholar - 18
Salzburger W, Van Bocxlaer B, Cohen AS . 2014 Ecology and evolution of the African Great Lakes and their faunas. Annu. Rev. Ecol. Evol. Syst. 45, 519–545. (doi:10.1146/annurev-ecolsys-120213-091804) Crossref, Web of Science, Google Scholar - 19
Wagner CE, Harmon LJ, Seehausen O . 2012 Ecological opportunity and sexual selection together predict adaptive radiation. Nature 487, 366–369. (doi:10.1038/nature11144) Crossref, PubMed, Web of Science, Google Scholar - 20
Cohen AS 2007 Ecological consequences of early Late Pleistocene megadroughts in tropical Africa. Proc. Natl Acad. Sci. USA 104, 16 422–16 427. (doi:10.1073/pnas.0703873104) Crossref, Web of Science, Google Scholar - 21
Ivory SJ, Blome MW, King JW, McGlue MM, Cole JE, Cohen AS . 2016 Environmental change explains cichlid adaptive radiation at Lake Malawi over the past 1.2 million years. Proc. Natl Acad. Sci. USA 113, 11 895–11 900. (doi:10.1073/pnas.1611028113) Crossref, Web of Science, Google Scholar - 22
Lyons RP 2015 Continuous 1.3-million-year record of East African hydroclimate, and implications for patterns of evolution and biodiversity. Proc. Natl Acad. Sci. USA 112, 15 568–15 573. (doi:10.1073/pnas.1512864112) Crossref, Web of Science, Google Scholar - 23
Berthold T . 1990 Phylogenetic relationships, adaptations and biogeographic origin of the Ampullariidae (Mollusca, Gastropoda) endemic to Lake Malawi, Africa. Abh. Naturwiss. Vereins Hamburg NF. 31–32, 47–84. Google Scholar - 24
Brown DS . 1994 Freshwater snails of Africa and their medical importance. London, UK: Taylor & Francis. Crossref, Google Scholar - 25
Mandahl-Barth G . 1972 The freshwater Mollusca of Lake Malawi. Rev. Zool. Bot. Afr. 86, 257–289. Google Scholar - 26
Schultheiß R, Van Bocxlaer B, Wilke T, Albrecht C . 2009 Old fossils-young species: evolutionary history of an endemic gastropod assemblage in Lake Malawi. Proc. R. Soc. B 276, 2837–2846. (doi:10.1098/rspb.2009.0467) Link, Web of Science, Google Scholar - 27
Berthold T . 1990 Intralacustrine speciation and the evolution of shell sculpture in gastropods of ancient lakes—application of Günther's niche concept. Abh. Naturwiss. Vereins Hamburg NF. 31–32, 85–118. Google Scholar - 28
Van Bocxlaer B, Schultheiß R, Plisnier P-D, Albrecht C . 2012 Does the decline of gastropods in deep water herald ecosystem change in Lakes Malawi and Tanganyika? Freshwat. Biol. 57, 1733–1744. (doi:10.1111/j.1365-2427.2012.02828.x) Crossref, Web of Science, Google Scholar - 29
Genner MJ, Knight ME, Haesler MP, Turner GF . 2010 Establishment and expansion of Lake Malawi rock fish populations after a dramatic Late Pleistocene lake level rise. Mol. Ecol. 19, 170–182. (doi:10.1111/j.1365-294X.2009.04434.x) Crossref, PubMed, Web of Science, Google Scholar - 30
Van Bocxlaer B, Hunt G . 2013 Morphological stasis in an ongoing gastropod radiation from Lake Malawi. Proc. Natl Acad. Sci. USA 110, 13 892–13 897. (doi:10.1073/pnas.1308588110) Crossref, Web of Science, Google Scholar - 31
Schultheiß R, Wilke T, Jørgensen A, Albrecht C . 2011 The birth of an endemic species flock: demographic history of the Bellamya group (Gastropoda, Viviparidae) in Lake Malawi. Biol. J. Linn. Soc. 102, 130–143. (doi:10.1111/j.1095-8312.2010.01574.x) Crossref, Web of Science, Google Scholar - 32
Meudt HM, Clarke AC . 2007 Almost forgotten or latest practice? AFLP applications, analyses and advances. Trends Plant Sci. 12, 106–117. (doi:10.1016/j.tplants.2007.02.001) Crossref, PubMed, Web of Science, Google Scholar - 33
Foll M, Gaggiotti O . 2008 A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993. (doi:10.1534/genetics.108.092221) Crossref, PubMed, Web of Science, Google Scholar - 34
Lynch M, Milligan BG . 1994 Analysis of population genetic structure with RAPD markers. Mol. Ecol. 3, 91–99. (doi:10.1111/j.1365-294X.1994.tb00109.x) Crossref, PubMed, Web of Science, Google Scholar - 35
Vekemans X . 2002 AFLP-SURV: a program for genetic diversity analysis with AFLP (and RAPD) population data. Version 1.0. Bruxelles, Belgium: Université Libre de Bruxelles. Google Scholar - 36
Zhivotovsky LA . 1999 Estimating population structure in diploids with multilocus dominant DNA markers. Mol. Ecol. 8, 907–913. (doi:10.1046/j.1365-294x.1999.00620.x) Crossref, PubMed, Web of Science, Google Scholar - 37
Bonin A, Ehrich D, Manel S . 2007 Statistical analysis of amplified fragment length polymorphism data: a toolbox for molecular ecologists and evolutionists. Mol. Ecol. 16, 3737–3758. (doi:10.1111/j.1365-294X.2007.03435.x) Crossref, PubMed, Web of Science, Google Scholar - 38
Jombart T, Devillard S, Balloux F . 2010 Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11, e94. (doi:10.1186/1471-2156-11-94) Crossref, PubMed, Web of Science, Google Scholar - 39
Core Team R . 2015 R: a language and environment for statistical computing, version 3.2.1. Vienna, Austria: R Foundation for Statistical Computing. Google Scholar - 40
Pritchard JK, Stephens M, Donnelly P . 2000 Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Crossref, PubMed, Web of Science, Google Scholar - 41
Evanno G, Regnaut S, Goudet J . 2005 Detecting the number of clusters of individuals using the software structure: a simulation study. Mol. Ecol. 14, 2611–2620. (doi:10.1111/j.1365-294X.2005.02553.x) Crossref, PubMed, Web of Science, Google Scholar - 42
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I . 2015 Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. (doi:10.1111/1755-0998.12387) Crossref, PubMed, Web of Science, Google Scholar - 43
Van Bocxlaer B, Schultheiß R . 2010 Comparison of morphometric techniques for shapes with few homologous landmarks based on machine-learning approaches to biological discrimination. Paleobiology 36, 497–515. (doi:10.1666/08068.1) Crossref, Web of Science, Google Scholar - 44
Rohlf FJ . 2010 Tpsdig 2. Version 2.16. Stony Brook, NY: State University of New York at Stony Brook. Google Scholar - 45
Sheets HD . 2004 IMP: integrated morphometrics package. Buffalo, NY: Canisius College. Google Scholar - 46
Oksanen J 2015 vegan: Community Ecology Package. Version 2.3-0. See https://cran.r-project.org/package=vegan. Google Scholar - 47
Venables WN, Ripley BD . 2002 Modern applied statistics with S, 4th edn. New York, NY: Springer. Crossref, Google Scholar - 48
Clarke KR . 1993 Non-parametric multivariate analyses of changes in community structure. Aust. J. Ecol. 18, 117–143. (doi:10.1111/j.1442-9993.1993.tb00438.x) Crossref, Web of Science, Google Scholar - 49
Kruskal J . 1964 Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika 29, 1–27. (doi:10.1007/BF02289565) Crossref, Web of Science, Google Scholar - 50
Adams DC, Otarola-Castillo E . 2013 geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol. 4, 393–399. (doi:10.1111/2041-210X.12035) Crossref, Web of Science, Google Scholar - 51
Fraley C, Raftery AE, Murphy TB, Scrucca L . 2012 Mclust version 4 for R: normal mixture modeling for model-based clustering, classification, and density estimation In technical report no 597. Seattle, WA: Department of statistics, University of Washington. Google Scholar - 52
Leinonen T, Cano JM, Mäkinen H, Merilä J . 2006 Contrasting patterns of body shape and neutral genetic divergence in marine and lake populations of threespine sticklebacks. J. Evol. Biol. 19, 1803–1812. (doi:10.1111/j.1420-9101.2006.01182.x) Crossref, PubMed, Web of Science, Google Scholar - 53
Brommer JE . 2011 Whither PST? The approximation of QST by PST in evolutionary and conservation biology. J. Evol. Biol. 24, 1160–1168. (doi:10.1111/j.1420-9101.2011.02268.x) Crossref, PubMed, Web of Science, Google Scholar - 54
Lynch M . 1990 The rate of morphological evolution in mammals from the standpoint of the neutral expectation. Am. Nat. 136, 727–741. (doi:10.1086/285128) Crossref, Web of Science, Google Scholar - 55
Whitlock MC . 2008 Evolutionary inference from QST. Mol. Ecol. 17, 1885–1896. (doi:10.1111/j.1365-294X.2008.03712.x) Crossref, PubMed, Web of Science, Google Scholar - 56
Hardy OJ . 2003 Estimation of pairwise relatedness between individuals and characterization of isolation-by-distance processes using dominant genetic markers. Mol. Ecol. 12, 1577–1588. (doi:10.1046/j.1365-294X.2003.01835.x) Crossref, PubMed, Web of Science, Google Scholar - 57
Hardy OJ, Vekemans X . 2002 SPAGeDi: a versatile computer program to analyse spatial genetic structure at the individual or population levels. Mol. Ecol. Notes 2, 618–620. (doi:10.1046/j.1471-8286.2002.00305.x) Crossref, Google Scholar - 58
Nosil P, Egan SP, Funk DJ . 2008 Heterogeneous genomic differentiation between walking-stick ecotypes: ‘isolation by adaptation’ and multiple roles for divergent selection. Evolution 62, 316–336. (doi:10.1111/j.1558-5646.2007.00299.x) Crossref, PubMed, Web of Science, Google Scholar - 59
Wilding CS, Butlin RK, Grahame J . 2001 Differential gene exchange between parapatric morphs of Littorina saxatilis detected using AFLP markers. J. Evol. Biol. 14, 611–619. (doi:10.1046/j.1420-9101.2001.00304.x) Crossref, Web of Science, Google Scholar - 60
Fischer MC, Foll M, Excoffier L, Heckel G . 2011 Enhanced AFLP genome scans detect local adaptation in high-altitude populations of a small rodent (Microtus arvalis). Mol. Ecol. 20, 1450–1462. (doi:10.1111/j.1365-294X.2011.05015.x) Crossref, PubMed, Web of Science, Google Scholar - 61
Savolainen V 2006 Sympatric speciation in palms on an oceanic island. Nat. Rev. Genet. 441, 210–213. (doi:10.1038/nature04566) Google Scholar - 62
Bonin A, Taberlet P, Miaud C, Pompanon F . 2006 Explorative genome scan to detect candidate loci for adaptation along a gradient of altitude in the common frog (Rana temporaria). Mol. Biol. Evol. 23, 773–783. (doi:10.1093/molbev/msj087) Crossref, PubMed, Web of Science, Google Scholar - 63
Mallet J, Barton NH, Lamas M. G, Santisteban C. J, Eeley MH . 1990 Estimates of selection and gene flow from measures of cline width and linkage disequilibrium in Heliconius hybrid zones. Genetics 124, 921–936. Crossref, PubMed, Web of Science, Google Scholar - 64
Crossley R . 1984 Controls of sedimentation in the Malawi Rift Valley, Central Africa. Sediment. Geol. 40, 33–50. (doi:10.1016/0037-0738(84)90038-1) Crossref, Web of Science, Google Scholar - 65
Ring U, Betzler C . 1995 Geology of the Malawi Rift: kinematic and tectonosedimentary background to the Chiwondo Beds, northern Malawi. J. Hum. Evol. 28, 7–21. (doi:10.1006/jhev.1995.1003) Crossref, Web of Science, Google Scholar - 66
Fryer G . 1959 The trophic interrelationship and ecology of some littoral communities of Lake Nyasa, with especial reference to the fishes, and a discussion of the evolution of a group of rock-frequenting Cichlidae. Proc. Zool. Soc. Lond. 132, 153–281. (doi:10.1111/j.1469-7998.1959.tb05521.x) Crossref, Google Scholar - 67
Genner MJ, Nichols P, Carvalho GR, Robinson RL, Shaw PW, Smith A, Turner GF . 2007 Evolution of a cichlid fish in a Lake Malawi satellite lake. Proc. R. Soc. B 274, 2249–2257. (doi:10.1098/rspb.2007.0619) Link, Web of Science, Google Scholar - 68
Conte GL, Arnegard ME, Peichel CL, Schluter D . 2012 The probability of genetic parallelism and convergence in natural populations. Proc. R. Soc. B 279, 5039–5047. (doi:10.1098/rspb.2012.2146) Link, Web of Science, Google Scholar - 69
Wolf J.BW, Lindell J, Backström N . 2010 Speciation genetics: current status and evolving approaches. Phil. Trans. R. Soc. B 365, 1717–1733. (doi:10.1098/rstb.2010.0023) Link, Web of Science, Google Scholar - 70
Genner MJ, Turner GF . 2012 Ancient hybridization and phenotypic novelty within Lake Malawi's cichlid fish radiation. Mol. Biol. Evol. 29, 195–206. (doi:10.1093/molbev/msr183) Crossref, PubMed, Web of Science, Google Scholar - 71
Seehausen O . 2004 Hybridization and adaptive radiation. Trends Ecol. Evol. 19, 198–207. (doi:10.1016/j.tree.2004.01.003) Crossref, PubMed, Web of Science, Google Scholar - 72
Seehausen O 2008 Speciation through sensory drive in cichlid fish. Nature 455, 620–626. (doi:10.1038/nature07285) Crossref, PubMed, Web of Science, Google Scholar - 73
Joyce DA, Lunt DH, Genner MJ, Turner GF, Bills R, Seehausen O . 2011 Repeated colonization and hybridization in Lake Malawi cichlids. Curr. Biol. 21, R108–R109. (doi:10.1016/j.cub.2010.11.029) Crossref, PubMed, Web of Science, Google Scholar - 74
Martin CH, Cutler JS, Friel JP, Touokong CD, Coop G, Wainwright PC . 2015 Complex histories of repeated gene flow in Cameroon crater lake cichlids cast doubt on one of the clearest examples of sympatric speciation. Evolution 69, 1406–1422. (doi:10.1111/evo.12674) Crossref, PubMed, Web of Science, Google Scholar - 75
Van Bocxlaer B . 2017Data from: Hierarchical structure of ecological and non-ecological processes of differentiation shaped ongoing gastropod radiation in the Malawi Basin . Dryad Digital Repository. (http://dx.doi.org/10.5061/dryad.dn58p) Google Scholar
