Rugged adaptive landscapes shape a complex, sympatric radiation

Strong disruptive ecological selection can initiate speciation, even in the absence of physical isolation of diverging populations. Species evolving under disruptive ecological selection are expected to be ecologically distinct but, at least initially, genetically weakly differentiated. Strong selection and the associated fitness advantages of narrowly adapted individuals, coupled with assortative mating, are predicted to overcome the homogenizing effects of gene flow. Theoretical plausibility is, however, contrasted by limited evidence for the existence of rugged adaptive landscapes in nature. We found evidence for multiple, disruptive ecological selection regimes that have promoted divergence in the sympatric, incipient radiation of ‘sharpfin’ sailfin silverside fishes in ancient Lake Matano (Sulawesi, Indonesia). Various modes of ecological specialization have led to adaptive morphological differences between the species, and differently adapted morphs display significant but incomplete reproductive isolation. Individual fitness and variation in morphological key characters show that disruptive selection shapes a rugged adaptive landscape in this small but complex incipient lake fish radiation.

The fitness consequences of phenotypic variation in a given environment can be conceptualized with adaptive landscapes [1,11,12]. Phenotypes that are well adapted with respect to key evolutionary traits tend to occupy peaks in the adaptive landscape, whereas poorly adapted individuals occur in lower regions of the landscape [1,2,[11][12][13]. Adaptive landscape topologies may range from uniform (lack of selection) or unimodal (stabilizing selection) to bi-or multimodal (disruptive selection) distributions [2,12,13]. In the simplest case, such a landscape is based on fitness estimation over variation in one single phenotypic trait, such as the beak in birds [14,15]. However, it is likely that selection frequently affects more than a single structure, and adaptive landscapes can also be used to explore more complex, multivariate axes of phenotypic trait variation [1,11,13,16].
Adaptive radiations provide examples for the rapid evolution of conspicuous adaptations to alternative food resources and habitats [1]. Disruptive ecological selection may affect variation in multiple key traits, and gain additional momentum with the evolution of increasing species diversity [11]. The adaptive landscape of an adaptive radiation can therefore be expected to be complex and rugged, sensu Nosil [2, p. 58] (i.e. with multiple distinct peaks and deep valleys) [1,11 -13]. Empirical evidence for the theoretical prediction that rugged adaptive landscapes facilitate the evolution of species flocks remains still limited, certainly for adaptive radiations in which taxa evolved in sympatry. A notable exception is a study [17] of the small Cyprinodon pupfishes radiation on San Salvador islands, for which a direct link between phenotype and fitness is established, and which supports the role of rugged adaptive landscapes in adaptive radiation. However, this example is based on hybrids that were bred in cages or the laboratories, and from species that rarely interbreed in nature [18]. Therefore, evidence that rugged adaptive landscapes promote adaptive radiation in nature remains largely lacking.
Here we study a radiation of sailfin silversides that constitutes a promising system to examine the relation between the shape of adaptive landscapes and adaptive radiation. More precisely, we test whether ecologically based disruptive selection can drive the evolution of rugged adaptive landscapes under sympatric conditions.
The small, but morphologically complex fish radiation of 'sharpfins' (Telmatherinidae: Telmatherina) from ancient graben-Lake Matano (Central Sulawesi, Indonesia) is the most diverse lineage of sailfin silversides in the Malili Lakes System [19]. The species flock is endemic to Lake Matano [20] and consists of several ecological specialists with conspicuous morphological and behavioural adaptations to various ecological niches [21][22][23]. However, specific morphotypes do not form well-sorted clades in phylogenetic analyses based on mitochondrial markers [24]. High levels of interspecific gene flow [25] and direct evidence for hybridization [24,26] have been observed, but morphotypes nevertheless mate largely assortatively [25,27]. As such, reproductive isolation among Lake Matano's Telmatherina can be considered to be substantial but incomplete. The onset of sharpfin diversifications commenced approximately 0.9 Ma [28].
Using a random sample of sharpfins, we perform integrated analyses of an estimator for individual fitness, an array of morphological traits that influence ecological specialization, the habitat context and stomach contents. We estimate the adaptive landscape underlying the sharpfin radiation and present signatures of ecologically based disruptive selection in the light of population structure, habitat use and trophic ecology.

Material and methods (a) Sampling
We pre-defined a total of 160 sampling locations, distributed equally (approx. every 500 m) all around Lake Matano's shoreline (electronic supplementary material, figure S1a). Sampling locations were fished in a random order between May and August 2010, by placing multimesh gillnets at depths of 0.75, 1.5, 3, 4.5 and 6 m. Gillnets were placed consistently for 20 min; fishes were euthanized in MS222, and the gutted body and liver were weighed with a precision of 0.001 g. Subsequently, a muscle tissue sample was stored in RNAlater (Qiagen), whereas specimens were transferred to 80% ethanol. Only mature individuals (determined by fin shape and colour in males, and gonad state in females) were included in the analyses (n ¼ 1060).

(b) Ecomorphological traits
We analysed variation in overall body shape and head shape with geometric morphometrics from digital X-ray images [29]. Homologous landmarks (15 and 13 for the body and head, respectively; electronic supplementary material, figure S1b) were digitized in TPSDIG2 [30]. GLS procrustes superimposition coordinates [29] were calculated for each dataset in PAST [31] to reduce effects of size and position. The first right gill arches were dissected and measured with IMAGEJ [32] (electronic supplementary material, figure S1b). Principal component analyses (PCA) were carried out independently for each geometric morphometric dataset and the relative lengths of the measured gill traits (residuals of a linear regression of measured distances by the length of upper arch).

(c) Phenotypic structure
Principal components (PCs) of the shape analyses and the gill trait measurements were scored to their true rank, which is defined here as the number of PCs, explaining together 99% of the total variance in each dataset. The datasets were then scored by multiplying each PC's residual by the square root of the rank divided by the variance of the dataset [33]. We combined the PCs of the shape analyses and the gill trait measurements in a k-mean clustering approach, and in a discriminant principal component analysis (DAPC; R package adegenet [34,35]), to examine how the various morphotypes in the radiation are structured in morphospace. Overfitting during the DAPC was avoided by retaining less than one-third of the original number of PCs (18 of 56) [34].
Identifying the actual number of clusters with a statistical procedure (e.g. BIC support values for various scenarios) is difficult [34]. Here, the BIC values indicate that 5-11 clusters can appropriately summarize the data. To estimate the actual number of clusters within the obtained range we applied the NbClust function [36,37] to the dataset. This function uses 22 criteria to determining the number of clusters and indicated that five clusters (morphotypes) summarize the data most adequately (electronic supplementary material, table S1).

(d) Stomach content analysis
Individual stomach contents were analysed following Pfaender et al. [23,38]. Food items of the gastrointestinal tract between oesophagus and pylorus were determined to the most accurate taxonomic level feasible and their relative volumetric proportions (in %) were estimated. The 933 individuals containing food items in their stomach were included in the stomach content analysis. Strong negative correlations between some types of content were observed, especially between molluscs, terrestrial insects, copepods and fish eggs. Therefore, a PCA was applied to the correlation matrix of the proportional stomach contents to reduce multicollinearity. We applied statistical tests (ANOVAs and Tukey's post hoc tests when variances were equal; Welch tests with Dunnet T3 post hoc tests when they were unequal) on individual PCs to test for pairwise differentiation between morphotypes in feeding ecology.

(e) Habitat differentiation
We examined the habitats where fish were caught to identify habitat parameters that may explain differences in the abundance of sharpfin morphotypes. Four habitat parameters were included in the canonical correspondence analysis (CCA; R package vegan [39]): dominant substrate at the general sampling location (rock, gravel, sand); substrate directly at the gillnet (soft to hard in %); bathymetric slope where the gillnet was positioned (58 steps); and sedimentation, which is related to surrounding canopy coverage (low, medium, high). Significance of the habitat-morphotype associations and individual canonical axes were calculated using the anova.cca function [39] with 11 000 permutations per step. Each placed gillnet was treated as a single catch event.
(f ) Amplified fragment length polymorphism procedure DNA was extracted from muscle tissue (NucleoSpin 96; Machery and Nagel) of the 1060 individuals, following the manufacturer's protocol. Quality of DNA was checked on 1% agarose gels; DNA concentrations were measured with a NanoDrop 200 (Thermo Scientific) and adjusted to 25 ng ml 21 . We applied a slightly modified AFLP method described by Vos [40] (see the electronic supplementary material, notes and [25,24] for details). In total, 1218 polymorphic AFLP markers were retained.

(g) Genetic differentiation and population structure
To test how the identified phenotypic clusters relate to genetic population structure, a DAPC and reassignment test (see detailed description above) were applied to the retained AFLP markers. Analysis of molecular variance (AMOVA [41]), implemented in ARLEQUIN v. 3.5.1.2 [42], was used to investigate the structure of the genetic variation among phenotypic clusters.

(h) Loci under selection
Allele distribution models are a powerful method to identify alleles under selection from AFLP datasets [43,44]. They use multiple univariate logistic regressions to test the relation between the frequency of AFLP loci and phenotypic variables. We applied ADMs in MATSAM [43,44] on all retained loci and the PC residuals from the ecomorphological trait analyses, centroid size, sex and sample location to test for correlation. p-values resulting from multiple tests have been adjusted with the sequential goodness-of-fit correction (SGoF [45]). The correction controls for the family-wise error rate, which is suitable given that a large number of tests are performed [45].

(i) Geographical differentiation
Pairwise least-cost paths between the sampling locations were calculated from a resistance map with CIRCUITSCAPE [46]. As sharpfins inhabit the littoral zones of Lake Matano [19,24], we clipped out areas deeper than 50 m to avoid misleading interpretations of the shortest possible way across the lake for these fishes. We examined correlations between geographical distance (pairwise least-cost paths) and genetic (Euclidean) distance with Mantel tests (R package ade4 [47]) within morphotypes.

( j) Fitness proxies
Estimating fitness is challenging, especially in large natural populations. For fishes, various proxies have been applied to quantify lifetime reproduction and survivorship [48] (see also [8,48,49]). Here, we use the relative liver weight as proxy for individual fitness, because the liver is an important fat-storage organ in fishes [50] and is frequently used as an estimator for individual condition in fishes (e.g. [51,52]). Moreover, individual condition is tightly correlated with reproductive capacity in other fish species (e.g. [53][54][55][56]). Hence, it is reasonable to assume that relative liver mass is indicative for the individual energy status in sailfin silversides, and that an enhanced energy status enables an individual to produce larger numbers of viable offspring, thus leading to greater individual fitness. Relative liver weight was calculated as the residuals of a generalized linear model constructed with data on sex, body size and the ratio between log-transformed liver weight and gutted body weight. We tested the reliability of our fitness proxy by investigating the correlation between relative liver mass and reproductive capacity (the ratio between log-transformed gonad weight and gutted body weight) with a subset of 127 female sharpfins. Because female sharpfins spawn promiscuous numerous times per day (J.P. & F.H. 2012 -2014, personal observation), without distinct spawning season female gonad mass appears a reliable estimator of reproductive capacity [27].
The use of relative liver mass as a general proxy for individual fitness does not imply that all potential fitness components are adequately covered by it. Absolute size and phenotypic proportions such as body depth vary strongly among the morphotypes (see results below and [23]). To test for potential differences in how our proxy represents fitness for various morphotypes, we applied an ANOVA to the residuals of log liver weight versus log body weight; the test result did not show significant differences (F 4,1032 ¼ 1.902, p ¼ 0.108).

(k) Adaptive landscapes
To test for a signal of disruptive selection along major axes of morphological variation, we applied non-parametric regressions without an a priori assumption about the shape of the fitness surface. Therefore, we used the generalized additive models (GAM) routine of R package mgcv [34,57], with the fitness proxy as dependent variable and the linear discriminant axes (LDs) obtained from the linear discriminant analysis (implemented in the DAPC) of the ecomorphological PCs as predictor variables. The GAM models were calculated without splines, with cubic splines, and with thin plate splines. The final model was chosen based on the lowest Akaike information criterion values (electronic supplementary material, table S3f ). Splines were smoothed via generalized cross-validation and an approximation of the significance of smoothing terms was calculated following Wood [58]. Trait axes with significant smoothing according to the GAM analyses were visualized as pairwise combinations of two axes using the Tps function and its incorporated plot function (R package fields [59]). The smoothing parameter for the plot function minimizes the generalized cross-validation score, and incorporates the estimated degrees of freedom from the GAM analyses. To test whether sharpfin population means plot in regions of higher fitness than expected by chance, we compared the observed fitness with random simulations of mean fitness. Using the runif function in R [34], 1000 population means and 10 000 individuals were drawn from the adaptive landscape. Individual fitness values were obtained from the landscape, using the raster and extract functions (R package raster) [60]. A pairwise comparison of each population and randomly generated means was carried out, with an ANOVA and Tukey's post hoc test with individuals falling in a range of a standard error of 0.75 around random means. To adjust p-values given the large number of post hoc tests SGoF [45] was applied. To test whether a morphotype's mean fitness differed significantly from the population of random fitness means, a bimodal test was carried out in R [34]. Tests for which observed fitness mean was significantly greater than the pool of random fitness means were treated as positives and others as negatives. To identify whether morphotype means occupy an average peaks in the adaptive landscape, we compared (using methods described above) fitness values from the regions in the landscape that are occupied by two morphotypes with the values obtained in intermediate regions.

(a) Phenotypic diversity and ecology
Cluster analyses of the observed morphological variation in ecomorphological candidate traits (electronic supplementary material, figure S1b) [23,61 -63] support five morphotype groups (figure 1; electronic supplementary material, table S1; total assignment probability: 83.28%), which largely match the nominal sharpfin species [19,64]. The morphotype not confirm the existing hypothesis that two, fish-feeding sharpfins exist (T. abendanoni and T. sp. 'elongated') [19]. However, this result may be an artefact, because the fishfeeding morphs are rare in nature (and hence rare in our sampling), and because both morphs are very similar in comparison to the greater differences between various feeding types [19]. Therefore, further investigation is required in this respect. Because our samples contained mainly the 'elongated' fish-feeding morph we refer to this fifth morphotype as T. sp. 'elongated'.
As expected for an emerging adaptive radiation [1,2,65], the five morphotypes differ significantly in trophic ecology and habitat use ( figure 1; electronic [19], and hence occupies predominantly sandy habitats that serve as spawning grounds for these roundfins. The slightly smaller T. opudi feeds on insects and zooplankton; it prefers a hard substrate of inshore habitats that provide shelter, preferably covered by terrestrial plants. Predominantly fish-feeding sharpfins also prey on shrimps at sand and gravel habitats. Telmatherina sp. 'thicklip' forage nearly exclusively on shrimps on gravel fields, whereas T. wahjui inhabit the shallows of muddy areas (figure 1; electronic supplementary material, figures S2 and S3), where they feed on insects, molluscs and zooplankton.
Morphological differentiation among sharpfin morphotypes (figures 2a and 3) is consistent with adaptations predicted by their respective ecological niches [23,38,62,63,66]. These include short and widely spaced gill rakers in fish and shrimp predators, and long and closely spaced gill rakers in sharpfins feeding on fish eggs, zooplankton or insects (figures 2a and 3). Deep bodies, known to increase manoeuverability in fishes [63], are characteristic for the shrimp-and fish-egg-feeders, whereas slender shapes occur predominantly in fish predators and zooplankton-feeders (figures 1, 2a and 3). Pointed snouts are typically found in 'pick-feeding' fish species, as fish-egg-, zooplankton-and mollusc-feeders. The jaws of shrimp-feeders are, by contrast, compact and the lips are pronounced. The rather gracile and elongated jaws of fish-feeders indicate a modification for a fast movement, suggesting a ram-feeding mode of foraging (figures 1, 2a and 3) [38,62].

(b) Genotypic differentiation in sympatry
Sympatric speciation predicts weak genetic differentiation among diverging populations, restricted only to small proportions of the genome [67,68]. The morphotypes described above are generally well supported by our predominantly nuclear AFLP data as the probability for assignment to the correct morphotype is approximately 71%. However, this percentage also indicates that mismatches occur frequently, which relates to substantial overlap between the species in morphology and genetics ( table S2c). Removal of these putatively selective loci decreases the success of assignment to the correct morphotype markedly (52.29%). This finding suggests that sharpfin genetic population structure is maintained by comparatively small proportions of nuclear loci. We did not detect any spatial effects in population structure (electronic supplementary material, table S2d) and all morphotypes except for T. sp. 'thicklip' were found evenly distributed around Lake Matano's shoreline (electronic supplementary material, figure S1a).

(c) Adaptive landscapes
Disruptive selection is predicted to shape rugged adaptive landscapes during speciation in the face of gene flow [1,2,69]. As expected, relative liver weight is correlated with female gonad mass, and thus reproductive capacity in sharpfins (n ¼ 127C; t 9.42; R 2 41.5%; p , 0.001). The finding indicates that our earlier assumption that relative liver weight is a good proxy for individual fitness in sharpfins much as in other fishes [53][54][55][56] is justified. Upon comparing relative liver weight with variation in the size-independent morphological dataset in the GAM (R 2 5.85%; effective degrees of freedom: 11.57; generalized cross-validation: 0.191; electronic supplementary material, table S3), three LDs (LD1, LD2 and LD4) revealed a significant correlation. Resulting landscapes (figure 3; electronic supplementary material, figure S4) show multiple local maxima. Individual species generally plot near local fitness peaks in the landscape constructed with GAMs from the two main axes of morphological variation (LD1 versus LD2) and the fitness proxy (figure 3a). Some individuals do not appear to be particularly well adapted, but the mean fitness of three of the five species significantly exceeds random expectations (figure 3c). As mentioned, the shrimp-feeding T. sp. 'thicklip' is morphologically very different from the other species, and its fitness peak is separated from the region of enhanced fitness occupied by the four other species by a deep valley (figure 3a; electronic supplementary material, S5a). In the landscape LD2 versus LD4 constructed from the complete dataset, the macroinvertebratefeeder T. opudi occupies a significant fitness peak (electronic supplementary material, figure   'thicklip'), the number of individuals with fitness proxies exceeding neutral expectations is higher than that of sharpfins with below-average fitness; these cases are, however, not supported in rigorous binomial tests.
Differences in gill raker length dominate the morphological variance, explained by LD1, whereas changes in body and head shape are summarized in the following axes (LD2 and LD4; figure 3a,b; electronic supplementary material, figure  S4a,b). Therefore, disruptive selection is significant in gill traits, where species means tend to peak in the adaptive landscape with either short and wide-spaced or long and narrow-spaced gill rakers ( figure 3a,b). Signal for disruptive selection in head shape is especially prominent along axes of variation in snout proportions ( figure 3a,b). Finally, different kinds of body shapes are associated to high fitness proxy values: deep-bodied fishes with a short caudal peduncle, individuals with conspicuously oval body shapes and conspicuously slender sailfins with either an elongated or a short caudal peduncle (figure 3a,b).

Discussion
Although ecological speciation is considered a major driver of diversification [2,65], our knowledge on how ecological factors have contributed to lineage divergence and speciation in complex radiations remains limited [2,70]. Here, we document that ecological opportunities (such as different food sources in sharpfin silversides) promote radiation by creating complex, rugged adaptive landscapes to which populations adapt. Adaptation to ecological opportunities is driven by selection, which can be disruptive in rugged adaptive landscapes.
Our sampling of Lake Matano's sharpfins reveals the presence of five, mainly sympatric morphotypes (figure 2; electronic supplementary material, figure S1). Variation in morphological key traits, including the shape of the head, body and gill rakers, is directly related to different feeding strategies, behaviour and habitat use (figures 1 and 2; electronic supplementary material, figures S2 and S3). All morphotypes relate to ecologically and genetically largely differentiated species, from deep-bodied shrimp-feeders, inhabiting gravel fields, to smaller near-shore insectivores and slender fish-feeders dwelling with shoaling juveniles (figure 1). Genome-wide surveys of polymorphism demonstrate that the variation in ecomorphological traits is correlated with variation in a restricted number of genetic markers (electronic supplementary material, table S2c). These findings clearly meet predictions for an ongoing radiation, evolving sympatrically under ecological selection pressure [68,70]. The incipient nature of this adaptive radiation, however, also results in overlaps in morphology and genetic variability between four of the five morphotypes, and intermediate, usually less fit individuals remain relatively common (figures 2 and 3; electronic supplementary material, figure S5).
A central aspect of this study is the significant signal of disruptive selection along major axes of morphological variation (electronic supplementary material, table S3; figure 3), resulting in complex adaptive landscapes with multiple peaks (figure 3; electronic supplementary material, figure  S4). As predicted [2,12], the five morphotype means occupy several regions of high fitness ( figure 3). Interestingly, valleys between these fitness peaks are, with exception of T. sp. 'thicklip', not very deep, but frequently reflect regions of significant lower fitness (figure 3; electronic supplementary material, figure S5). This suggests a general lower fitness of individuals intermediate to the population means (i.e. hybrids) compared with the population means. Incomplete speciation, possibly in combination with adaptation that arises from combined divergence in several traits, or traits that are not quantified here, provides a possible explanation for a lack of significant valleys between some of the phenotypic less distinct sharpfins. Nevertheless, sharpfins clearly fulfil predictions for 'rugged' adaptive landscapes in an emerging adaptive radiation, driven by disruptive ecological selection.
Selective regimes may vary over time, as reported, for example, from postglacial stickleback systems [71]. We also cannot exclude that the adaptive landscapes reconstructed here, and the selective regimes that shape divergence of the sharpfin flock, may change through time. However, there are no apparent hints for this assumption. The environment of Lake Matano is considered long-term stable [20]. Furthermore, correlations among traits and sources of nutrition remained consistent across studies that are based on material from different years (cf. [23]). Finally, the present samples originate from a total of 160 locations, sampled over a period of four months, a design that largely prevents dominance of local or single effects. All this suggests that the snapshot view of the present study is not an exception.
Disruptive ecological selection has been identified as a causal root of population divergence also in other lake fishes, including, for example, species pairs of sticklebacks or coregonids, diverging in ecomorphological traits along the benthic-limnetic axis [8,49,72,73]. In contrast with the comparatively frequent examples for bimodal population divergence, the importance of disruptive ecological selection in complex species flock formation has only rarely been demonstrated [17,68]. Among the available examples are two bird radiations, Geospiza finches and Loxia crossbills, where disruptive selection acts on beak size [14,15]. In fishes, the adaptive landscapes of the Cyprinodon pupfish radiation [17] show patterns similar to those detected in the sharpfin radiation, though landscape topographies are more complex in sharpfins. This complexity largely fits the presence of five phenotypic groups in sharpfins, compared with only two ecological specialists and one generalist in the pupfish example [17,18].
The natural occurrence of interbreeding between sailfin morphospecies, largely separated by divergent ecological selection, provides a plausible explanation for both the complex adaptive landscapes and the porous population boundaries observed. Substantial but incomplete rspb.royalsocietypublishing.org Proc. R. Soc. B 283: 20152342 reproductive isolation in the sharpfin flock is apparent from the weak, but significant population genetic structure (electronic supplementary material, table S2). Assortative mating has been reported from most of Lake Matano's Telmatherina species, including three of the five sharpfins (T. sarasinorum, T. wahjui, T. sp. 'thicklip'; see [21,22,27,24]). Together with significant differences in habitat use (figure 1; electronic supplementary material, figure S3), genetic population structure and assortative mating provide strong arguments to reject the hypothesis of only one hypervariable species, and are in agreement with very similar patterns reported from its sister clade, Lake Matano's 'roundfins' [25]. In both Telmatherina clades, roundfins and sharpfins [24], divergence apparently follows ecological selection pressure under sympatric conditions, leading to the consensus view that the 'Matano flock' probably represents an early stage of adaptive radiation [1,2,5].
However, evidence for early stages of species flock formation appears surprising in the light of recent estimations for the age of divergence in sailfin silversides. Stelbrink et al. [28] applied different molecular clock approaches, and argued that the initial divergence of Telmatherina within Lake Matano, leading to the formation of sharpfins and roundfins, dates back approx. 1.9 Myr. This is in general agreement with the estimated geological age of the lake [20], and implies ages of approx. 1 Myr for the onset of divergence within both sharpfins and roundfins [28].
Three potential explanations for resolving this conflict appear moderately plausible. First, environmental instabilities may repeatedly shift the balance between disruptive selection and homogenizing gene flow, and may even lead to a breakdown of the barriers to gene flow and the reamalgamation of incipient species (e.g. [74,75]). However, we consider a 'speciation reversal' scenario rather unlikely in the case of extraordinarily deep Lake Matano, which is considered long-term stable [20,76]. There are so far no indications for significant earlier disturbances, or increased turbidity or other changes that might plausibly trigger the breakdown of reproductive barriers, despite the lake nowadays being clearly affected by human activities [77].
Second, intralake divergence may have proceeded, for unknown reasons, generally at rather low speed. This would be surprising in a species flock characterized by pronounced sexual and ecological selection pressure [78], and rates of speciation reported from other fish radiations [79]. As a third alternative, speciation processes might remain incomplete in Telmatherina of this lake. This would imply that sharpfins and roundfins remain in a stage of the speciation continuum [2,5] where disruptive ecological selection is strong enough to maintain a stable phenotypic diversity, but fails to complete the speciation process. Superficially similar cases have been reported from walking-stick insects or haplochromine cichlids (reviewed in [67]), but we are not aware of other adaptive radiations where speciation likewise remained incomplete about a million years after the onset of speciation processes.
Ethics. The study was carried out under research permits from the Kementerian Negara Riset dan Teknology (RISTEK), Indonesia ( permit no. 992/FRP/SM/V/2010) and in cooperation with the Indonesian Institute of Sciences (LIPI). No endangered or protected species were involved. Fishes were sacrificed with MS-222. All procedures followed the Guidelines for the Use of Fishes in Research from the American Fisheries Society and the legal requirements of Indonesia and Germany.