Skip to main content
  • Research article
  • Open access
  • Published:

Live fast, diversify non-adaptively: evolutionary diversification of exceptionally short-lived annual killifishes



Adaptive radiations are triggered by ecological opportunity – the access to novel niche domains with abundant available resources that facilitate the formation of new ecologically divergent species. Therefore, as new species saturate niche space, clades experience a diversity-dependent slowdown of diversification over time. At the other extreme of the radiation continuum, non-adaptively radiating lineages undergo diversification with minimal niche differentiation when ‘spatial opportunity’ (i.e. areas with suitable ‘ancestral’ ecological conditions) is available. Traditionally, most research has focused on adaptive radiations, while empirical studies on non-adaptive radiations remain lagging behind. A prolific clade of African fish with extremely short lifespan (Nothobranchius killifish), show the key evolutionary features of a candidate non-adaptive radiation – primarily allopatric species with minimal niche and phenotypic divergence. Here, we test the hypothesis that Nothobranchius killifish have non-adaptively diversified. We employ phylogenetic modelling to investigate the tempo and mode of macroevolutionary diversification of these organisms.


Nothobranchius diversification has proceeded with minor niche differentiation and minimal morphological disparity among allopatric species. Additionally, we failed to identify evidence for a role of body size or biogeography in influencing diversification rates. Diversification has been homogeneous within this genus, with the only hotspot of species-richness not resulting from rapid diversification. However, species in sympatry show higher disparity, which may have been caused by character displacement among coexisting species.


Nothobranchius killifish have proliferated following the tempo and mode of a non-adaptive radiation. Our study confirms that this exceptionally short-lived group have diversified with minimal divergent niche adaptation, while one group of coexisting species seems to have facilitated spatial overlap among these taxa via the evolution of ecological character displacement.


Ecological opportunity – the instance whereby novel environments with abundant resources become available – is a widespread source of divergent natural selection that can trigger the proliferation of an ancestor into multiple species adapted to different regions of the niche space [1,2,3,4]. Therefore, lineage diversification via this process, known as adaptive radiation [2], is expected to leave a phylogenetic signature consisting of an ‘early burst’ of species accumulation via niche-filling, followed by a slowdown in diversification rate as newly emerging species saturate ecological opportunity [5,6,7,8,9]. Consequently, adaptive diversification is diversity-dependent, as saturation of niche space reaches the ecological capacity of a given environment [10,11,12].

Evolutionary radiations are not exclusively adaptive [13], however, and can take place when species diversify with retention of the ancestral niche [14]. During this process of non-adaptive radiation [15], episodes of speciation are not the result of divergent natural selection and thus, newly emerging species remain ecologically and phenotypically similar [16]. Importantly, it has been suggested that in contrast to adaptive radiations, diversification via non-adaptive radiation is more likely to take place when newly emerging species radiate across geographically non-overlapping areas (allopatry), which can accommodate species with fundamentally similar niche demands [15, 17, 18]. The availability of such areas equipped with suitable ‘ancestral’ ecological conditions is what we refer to as ‘spatial opportunity’. Given that this form of radiation is bound by the availability of spatial opportunity, rather than by environments with an abundance of diverse niche space, the process of species accumulation is not expected to always exhibit a slow-down over time. Instead, an arithmetic accumulation of species is expected when incipient species encounter new spatial opportunity to continue clade expansion producing a phylogenetic late-burst [18]. Alternatively, if an abundance of spatial opportunity is available and exploited early in the clades’ history, it will produce an early burst, akin to an adaptive radiation [8]. Consequently, depending on the availability of vacant ancestral niches, non-adaptive radiations may diversify under a range of radiation trajectories.

Empirical investigations quantifying the tempo and mode of adaptive diversification have been historically popular and conducted across an array of continental and island radiations [6, 19,20,21,22,23]. The fundamental principles underlying non-adaptive radiations are well-studied ([24] and references therein). In contrast, the phylogenetic characterisation of the tempo and mode of non-adaptive radiations has emerged more recently and remains a pending task in our understanding of the role of adaptation during diversification [8, 18, 25, 26]. The extent to which niche and phenotypic disparity are preserved across prolific clades, and the impact of coexistence among some of the species to exert a triggering phenomenon such as ecological character displacement remain poorly known. Body size is frequently used to investigate both phenotypic and niche evolution across adaptive and non-adaptive radiations due to its central role in ecological and evolutionary processes and evolutionary correlation to other character traits [27,28,29,30].

The Nothobranchius genus of African killifish consists of over 70 freshwater species adapted to annually desiccating savannah pools, and which survive desiccation by diapausing embryos during the dry season specifically in alluvial vertisol soil [31,32,33,34]. As a result, Nothobranchius species have evolved extremely short lifespans, often limited to a few weeks [35, 36]. The genus shows a history of vicariance-driven diversification, with the mode of dispersal between pools still unclear [37, 38]. Nothobranchius species (71 [39]) have predominantly radiated allopatrically, currently ranging from Sudan to South Africa along the East of Africa, including Zanzibar and Mafia Islands [37, 40]. This spatial structure is coupled with minimal ecological and morphological disparity among species [40,41,42]. An area of lowland Tanzania apparently deviates from this tendency, as a high number of species coexist over a relatively small area, forming the only hotspot of Nothobranchius species richness [37, 43]. This hotspot (Fig. 1) is hypothesised to have resulted from dispersal into the hotspot or high diversification, either via secondary contact after speciation in allopatry, or via sympatric speciation. Both hypotheses predict that the interspecific coexistence within this area is conditioned by the evolution of ecological character displacement among species, which could have taken place during sympatric speciation (i.e. consistent with a process of adaptive diversification) or from divergent selection pressures forcing species to diversify phenotypically during earlier stages of allopatric diversification prior to secondary contact.

Fig. 1
figure 1

Nothobranchius species richness map. Richness and distance scale bottom left, with the blue outline representing the species hotspot. The mapped distributions are based on original data and the map was created using ArcGIS v10.0

Here, we employ phylogenetic modelling to investigate the tempo and mode of macroevolutionary diversification of the genus Nothobranchius. These analyses are implemented to test the core hypothesis that the Nothobranchius clade has predominantly radiated via non-adaptive diversification. The tempo of diversification was investigated by determining whether diversification was heterogeneous through time and across the phylogeny. Non-adaptive diversification was tested by modelling the evolution of body size to elucidate the evolutionary process that produced the pattern of minimal phenotypic disparity. We investigate whether any diversification heterogeneity across the phylogenetic tree is the result of species’ body size, biogeography, or an unmeasured trait. Finally, by investigating the effect of biogeography on diversification we explore why a hotspot has arisen and whether these species have shifted along the ‘radiation continuum’.


Data collection

To quantify diversification in phenotype and the role of the phenotype in diversification rates, we gathered a dataset of body size for 48 Nothobranchius species, using maximum standard length (SLmax) (a linear measure from the anterior tip down to the midlateral posterior of the hypural plate). Body size represents the primary descriptor for niche space in fishes as it integrates the most significant information on ecological characteristics of particular species [28, 44]. Measures of body size were obtained from the data collected directly by one of us (MR), and collected from multiple literature sources (Additional file 1). Prior to all analyses, body size data were transformed into natural logarithm to meet requirements of data distribution of variance homogeneity (Additional file 2). Secondly, we created a new dataset consisting of all occurrence records known for all Nothobranchius species present in the phylogenetic tree (see below). These data were obtained from extensive field records collected over the course of 10 years by one of us (MR, [45]), as well as from coordinates recorded by overlapping the maps published by Wildekamp [46] with a geographic information system and a set of secondary resources (Fig. 1; Additional files 3 and 4).

Phylogenetic tree

Macroevolutionary phylogenetic analyses of diversification were performed on a time-calibrated molecular phylogenetic tree that includes 49 of the 71 documented Nothobranchius species from Dorn et al. [37] (Additional file 5; [39]). The tree was built based on the analysis of five nuclear markers and one mitochondrial marker; fossil-based calibration was not possible due to the lack of a Nothobranchius fossil record, but used secondary calibration from a fossil-calibrated tree of spiny-rayed fish [37, 47]. For lineage diversification analyses and biogeographic analysis we used the full phylogeny (n = 49). For body size and trait-dependent diversification analyses, we used a phylogenetic tree of the 48 species for which these data were available.

Lineage diversification analyses

Tempo of diversification – the relationship between speciation and extinction rate – was qualitatively analysed using a lineage-through-time (LTT) plot using the package ‘ape’ [48] in R version 3.3.1 [49]. The LTT is plotted against a pure-birth null model displayed as confidence intervals ranging from 50 to 99%, to identify significant pulses of lineage accumulation. Taxa missing from the incomplete phylogeny (n = 49) are not corrected for in LTT analysis. The gamma (γ) statistic was calculated using phylogenetic node distances to detect whether the tempo of diversification is significantly different from the diversification tempo of a pure-birth null model, and can be used to quantify whether diversification underwent an early burst or late burst [50]. A property of the ultrametric tree is all tips (i.e. extant taxa) are of equal distance to the root of the tree (i.e. common ancestor), thus the bifurcation events that occur indicate tempo changes through time [51]. The γ statistic was produced using the R package ‘ape’ [48], to test for significance a critical value was calculated by running a Monte Carlo Constant Rate (MCCR) analysis using R package ‘laser’ [52]. MCCR analysis accounts for species absent from the phylogenetic tree, as an increase in type I error and negative bias arises when using an incomplete phylogeny [53]. Bias-correction used a proposed clade size of 71 Nothobranchius species. MCCR analysis was run for 10,000 iterations.

Hidden Markov model based maximum-likelihood analysis of clade diversification was conducted using R package ‘DDD’ [8]. Four models derived from the birth-death process were fitted to analyse whether diversification rate is constant or a function of species richness. Two diversity-independent models were fitted: the Yule model (pure-birth model), defined by a non-zero speciation rate and zero extinction rate, and a constant rate birth-death (crBD) model with non-zero speciation and extinction rates. In addition, two diversity-dependent models were fitted: diversity-dependent linear (DDL + E), and diversity-dependent exponential (DDE + E) models, that model speciation rate as a linear and exponential function of number of species in the phylogeny at any point in time, respectively, while accounting for a non-zero extinction rate [8]. Goodness-of-fit was established by ranking model on the bias-corrected Akaike Information Criterion (AICc) using maximum likelihood values showing the best-fitting model with minimum information loss [54, 55]. Similar to the MCCR test, clade diversification analysis incorporates taxa absent from the phylogeny. In the case of Nothobranchius, the phylogenetic tree assemblage of 49 is missing 22 species from the proposed complete clade of 71 [39]. All models were re-run to check whether model selection is consistent when the proportion of missing species changes, to account for the number of divergent populations that may be potentially described as separate species, as the conservative number of described species may bias diversification towards diversity-dependence. Parametric bootstrap likelihood ratio tests – at a significance value (α) of 0.05 – were used to ensure the best-fitting models chosen using AICc were reliable [56]. The best-fitting diversity-dependent model (DDL + E or DDE + E) from maximum likelihood analysis was compared with the constant rate birth-death model. The likelihood ratio test was carried out for 1000 bootstrap iteration at both 22 and 50 missing taxa.

To detect whether diversification rate is heterogeneous through time, we analysed the phylogenetic tree for rate shifts using Bayesian Analysis of Macroevolutionary Mixtures (BAMM) version 2.5 [57]. Phylogeny-specific priors for the number of diversification shifts, speciation and extinction parameters were estimated using the ‘BAMMtools’ package in R [58]. The priors are estimated based on the root age of the phylogenetic tree and using a pure-birth model. These estimates are subsequently used by BAMM to calculate a lineage-specific posterior distribution of speciation and extinction rates, which are then used to model diversification shifts. BAMM uses a reversible jump Markov Chain Monte Carlo (MCMC) approach, which we ran for 10,000,000 iterations to ensure convergence on four metropolis-coupled MCMC chains (deltaT = 0.1, swapPeriod = 1000), to detect diversification rate-shifts. BAMM incorporates species missing from the phylogeny, by assuming 22 species are missing and that phylogenetic species are sampled at random [59]. The R package ‘coda’ [60] was used to ensure MCMC convergence (effective sample size [ESS] for likelihood and number of shifts > 500) for posterior distribution accuracy after a 10% burnin. Model goodness-of-fit was delineated using Bayes factors (BF) from BAMM analysis [61], for a non-zero shift model to be accepted over the null hypothesis it must have a BF over 3 [62].

Phenotypic diversification analyses

Diversification of phenotypic traits was examined by trait disparity-through-time (DTT) analysis using the R package ‘geiger’ [63, 64]. This analysis uses body size of extant Nothobranchius species to reconstruct ancestral body size values and model disparity between species. DTT is carried out by the average squared Euclidean distance between species on a Euclidean (multivariant) space [19]. The null model of trait evolution is Brownian motion (BM), a stochastic evolution model of constant variance, produced by averaging 10,000 BM iterations [65]. Confidence intervals for DTT were calculated using the rank envelope test as this shows better type I and II error rate than the original pointwise envelope test [66]. Morphological disparity index (MDI) quantified the disparity from the reconstructed evolution of body size from extant species and the median trait values under Brownian evolution simulations. Whereby a positive MDI represents a greater disparity than Brownian expectation, conversely a negative MDI represents less disparity than Brownian expectation [19, 67]. Ancestral state reconstruction of body size was carried out using a maximum likelihood BM model of trait evolution and plotted onto the phylogenetic tree using the R package ‘phytools’ [65, 68, 69].

The quantification of body size evolution was carried out by fitting four models of trait evolution: Brownian motion (BM), Ornstein-Uhlenbeck (OU), Early-Burst (EB), and the Delta model using the R package ‘geiger’ [21, 63, 64]. BM models trait evolution between species as a correlation to the time since divergence as the trait values evolve under a random pattern with variance determined by the scaling parameter (σ) [65]. OU models trait evolution under Brownian motion with an additional stabilizing selection term, driving towards an optimum value (θ) at a rate of adaptation (α) [70, 71]. The EB model of evolution is the exponential decrease in trait evolution over time, indicative of the early divergence in adaptive radiations [21]. The Delta model of evolution detects the changes in tempo of trait evolution [72]. The delta (δ) statistic infers whether trait evolution happens early in the clades’ history and then experienced a slow-down (δ < 1), analogous to EB; or whether the trait evolution has occurred in the relative recent past (δ > 1), suggestive of adaptation to new fitness optima. When δ = 1 it indicates constant tempo through the lineage, equal to the tempo described by the BM model. All models assume gradual trait evolution, not testing for punctuated equilibrium [73]. Goodness-of-fit was determined using AICc.

We conducted an a posteriori approach to determine whether body size evolution is under a single- or multiple-optima OU model; testing whether morphology has converged around one or more fitness optima. This was carried out using the R package ‘SURFACE’ [74] and ‘ℓ1ou’ [75]. SURFACE determines the number and location of selection regimes (OU optima) across the phylogeny. Successive regimes are modelled across the tree using a stepwise approach until best-fit model using AICc, then regimes are collapsed until the most parsimonious regimes and shifts are chosen using AICc [74]. This stepwise AICc method has limitations in its false positive detections rate for optima shifts [75, 76]. The ℓ1ou uses a phylogenetic lasso method to detect and locate OU shifts, resolving to some extent the limitations of the SURFACE method by using a phylogenetic Bayesian information criterion (pBIC) [75]. To quantify support for each shift detected using ℓ1ou, bootstrapping was run for 1000 iterations. These methods are used comparatively as SURFACE will detect small magnitude OU shifts in a few or single species, while ℓ1ou will provide a more conservative estimate of shifts. Adaptive peaks were plotted onto a phenogram of reconstructed ancestral lineage’s body size using maximum likelihood BM in the R package ‘phytools’ [65, 68, 69].

To determine whether Nothobranchius body size evolution displayed significant rate shifts throughout the phylogeny and whether this would mirror any diversification shifts we ran BAMM for trait evolution [57]. Priors for rate of trait evolution and expected number of shifts were calculated using the ‘BAMMtools’ R package [58]. The trait priors are estimated based on the root age of the phylogenetic tree and maximum likelihood under a Brownian motion model. BAMM calculates posterior distributions to model trait evolution. BAMM was implemented under the same procedure as diversification analysis (see above). MCMC convergence (likelihood and number of shifts ESS > 500, burnin 10%) was tested using the R package ‘coda’ [60]. BF were calculated for model selection [61].

Modelling trait-dependent diversification

In order to determine whether the rate at which lineages diversify is dependent on that lineage’s body size we ran trait-dependent diversification analysis using; a tip rate correlation technique termed inverse of equal-splits with simulated null model (ES-sim) [77], and maximum likelihood approach Quantitative State Speciation Speciation and Extinction (QuaSSE) [78]. These were used comparatively as for a phylogeny of ~ 50 species ES-sim provides a better Type I error rate but QuaSSE has higher statistical power [77]. ES-sim measures the tip-specific speciation rate and simulates a null distribution under a given model of trait evolution in order to test for significance. We ran ES-sim for both Brownian and OU null distributions by running 10,000 simulations, using the ‘essim’ code in R [77]. Spearman’s correlation was used to determine a significant monotonic (i.e. linear and sigmoidal) trait-dependent diversification relationship using body size data of extant species. QuaSSE is a state-dependent speciation-extinction model based on BiSSE [79] using the R package ‘diversitree’ [80]. Analysis using a birth-death process was performed using the phylogenetic tree and body size data to model: linear-, sigmoidal- and modal-diversification, for both stochastic trait evolution, and directional trait evolution [78]. Assuming random sampling, missing species are corrected for to avoid speciation-extinction bias for an incomplete phylogeny [59]. QuaSSE analysis calculates maximum likelihood for model selection, from which AIC was calculated to delineate model goodness-of-fit.

Biogeographic analysis

Of particular interest is a region of high species density in Tanzania. The boundary of the hotspot is delineated using 1o latitude × 1o longitude grid cells and includes all cells with greater than 10% of Nothobranchius species in the phylogeny (Fig. 1). To test whether species’ biogeography influences diversification we tested for biogeographic-dependent diversification using a geographic state and hidden state speciation-extinction model (GeoHiSSE) analysis, using R package ‘HiSSE’ [81,82,83]. This model uses a three-state Markov Chain, to allow for species that live in distinct biogeographical area (A or B), or occupy both areas (AB). Parameters of the model are speciation rate in each area (sA, sB) as well as both areas (sAB); the extinction rate of each area (xA, xB), and dispersal between areas (dA = A➔B, dB = B➔A). Missing species are corrected for using the same technique as QuaSSE [59]. We fit a seven parameter model without hidden states (GeoSSE), an equal model with hidden states (GeoHiSSE), and then two character-independent null models with hidden states (CID GeoSSE and CID GeoHiSSE). The null models have equal complexity to ensure for proper hypothesis testing and that the importance of geography in the formation of the hotspot [83]. Best-fit models were chosen using AICc, while parameter estimates are taken from model-averaging using Akaike weights (AICw). We do not attempt to delineate between range contraction and widespread species extinction in one geographic region due to insufficient phylogeny size.


Tempo and mode of lineage diversification

The accumulation of lineages is constant through time. Diversification is consistent with the pure-birth null model through the clades’ phylogenetic history (Fig. 2). A non-significant pulse in the late Pleiocene-early Pleistocene (~ 2.5–2 Mya), possibly represents a decreased rate of diversification. The slightly negative γ value (γ = − 0.145) indicates that the branching events of the reconstructed phylogenetic tree possesses a higher density of internal nodes marginally closer to the root of the tree than the tree tip compared to the pure-birth null model. The MCCR analysis found a critical value to reject a pure-birth model of − 2.298, consequently the γ statistic is highly non-significant (p = 0.731, 95% CI one-tailed test), rejecting an early- or late-burst of diversification.

Fig. 2
figure 2

Top, lineage-through-time (LTT) accumulation curve for Nothobranchius radiation (black line). Unit of time (x-axis) is in millions of years before present (Mya). Coloured area represents confidence intervals for pure-birth model. Confidence intervals range from 50% (0.5) to 99% (0.99). Bottom, disparity-through-time (DTT) for body size (solid line), against a Brownian motion null model (dotted line). Area within dashed lines is the confidence intervals at 95%. Time before present in in millions of years before present. Inset in the DTT plot is the Brownian motion reconstruction of body size across each lineage with subclade A labelled

Lineage diversification analysis elucidated that the Yule model was the best-fit for a proposed clade size of 71 and 99 species (Table 1). Constant rate birth-death also explained the data closely, while diversity-dependent models showed markedly less support (ΔAICc> 2), suggesting lineage tempo has not decreased with time. Comparison of the likelihood ratio tests from parametric bootstrapping was in accordance with the maximum likelihood analysis by rejecting the diversity-dependent (DDL + E) in favour of the constant rate birth-death model for both clade sizes (Additional file 6a, b, Additional file 7a and b). Therefore, all diversification analyses are congruent indicating no significant slow-down in diversification through time.

Table 1 Diversity-dependent diversification maximum likelihood analysis using bias-corrected Akaike Information Criterion (AICc) for model selection

BAMM analysis corroborates the constant rate through time with constant rate across lineages. After ensuring the MCMC had converged (likelihood ESS = 1461.15, number of shifts ESS = 1398.48), model comparison using BF showed that the null model with no shifts in diversification rate was the best-fitting model, a single shift model is the next best-fitting (BF = 1.38). The most frequent model selected was the null model of homogeneous diversification (frequency = 0.85), this is due to the majority of BAMM runs not detecting any significant rate shifts (Additional file 8). Posterior distribution detected a model which had a diversification shift on an ancestral lineage (Additional file 8), though this model was substantially less frequent (frequency = 0.1).

Body size evolution

The observed minimal morphological disparity through time is consistent with a non-adaptive diversification during the clades’ history. No significant disparity (95% CI) was detected in the clades’ history, although relatively high disparity is shown in the recent past (~ 1 Mya to present) indicating recent phenotypic divergence. The DTT plot (Fig. 2) showed consistency with the BM null model throughout most of the lineage’s history. Disparity of body size within clades is equal to or slightly higher than between clades (MDI = 0.113), but deviation from the null is not significant (p = 0.184). The positive MDI value is characteristic of species slightly overlapping in the body size morphospace. Reconstruction of species ancestral body size under BM on the phylogeny shows homogenous body size evolution with all species being morphologically similar (Fig. 2). Although a subclade (hereafter referred to as subclade A, see Fig. 2) with a common ancestor ~ 5 Mya had divergent body size evolution in the recent past. Other species that show extreme trait values diverged from a common ancestor in the distant past (N. ocellatus and N. thierryi).

The OU model had the best-fit in continuous trait analysis, with a strong pull (α = 0.272). The Delta model of evolution comprehensively explains the pattern of the data, while BM showed marginal explanation, and EB was substantially incongruent (Table 2). The selection of OU over BM suggests there is a stabilizing selection pressure present around to one or more adaptive optima and falsifying a random walk mode of trait evolution. The greater likelihood of the Delta model of trait evolution over BM and EB – with a δ value more than 1 – suggests an adaptive proliferation late in the clades’ history. Stabilizing selection under an OU model would drive species to overlap in the morphospace around body size fitness optima. The two approaches to detect multiple OU optima identified different numbers of shifts. Both detected a shift on the terminal branch to N. ocellatus, with good bootstrap support from ℓ1ou (0.709). SURFACE detected this shift but also detected an additional four shifts on the adaptive landscape (Fig. 3). The two adaptive optima ℓ1ou detected were θA = 40.25 mm and θB = 115.79 mm, whereas SURFACE detected four adaptive optima within the body size range (θ1 = 41.26 mm, θ2 = 95.58 mm, θ4 = 59.74 mm, and θ5 = 30.57 mm) along with an optimum distant from all body sizes (θ3 = 4.76 mm). SURFACE found convergence between two species (N. ocellatus and N. orthonotus). Similar body size across many species is evidence of a stabilizing selection pressure, proposed in the OU model of trait evolution. Divergence was detected at three nodes of the phylogenetic tree all within subclade A using surface, while ℓ1ou did not detect any shifts in this subclade (Additional file 9).

Table 2 Continuous trait evolution analysis, using bias-corrected Akaike Information Criterion for model selection: using Brownian motion, Ornstein-Uhlenbeck, Early-Burst, and Delta models
Fig. 3
figure 3

Phenogram of trait evolution of each lineage reconstructed under a Brownian motion model, with shading around lineages the 95% CI. Fitness optima (θ) under a multi-optima OU model shown on y-axis. θA and θB are optima from ℓ1ou, and θ1, θ2, θ3, θ4, and θ5 are optima from SURFACE

After ensuring MCMC convergence (likelihood ESS = 1801.00, number of shifts ESS = 1475.34), BAMM analysis of trait evolution found that the null model of no rate shifts of body size evolution was the best-fitting, as all other models had BF < 3 (Additional file 10). Although the posterior distribution found several models which include a single rate shift, the most congruent is the null model with trait evolution not undergoing significant shifts. The analysis detected four shifts all occurring in subclade A (Additional file 10).

Trait-dependent diversification

Trait evolution analysis showed that OU is the best-fit model for Nothobranchius body size with a relatively strong pull (Table 2), thus the OU null model for ES-sim is the most appropriate to test. The OU null model showed diversification is not dependent on body size when testing for a monotonic relationship (P = 0.963), with minimum correlation between speciation rate and trait (Spearman’s rho = − 0.014; Additional file 11). Under a Brownian null model, diversification also showed insignificant trait-dependent diversification showing an equal relationship (Spearman’s rho = − 0.014, P = 0.957). QuaSSE analysis distinguished the modal model of diversification under stochastic trait evolution to best support the phylogenetic diversification with regards to body size (Table 3). Identifying body size to significantly affect the diversification of a clade (p0.05), inferring that species with intermediate body size have the highest diversification rate. Even though all other models display significance they show poorer support relative to the modal model (ΔAIC> 0). The models of trait-dependent diversification under the directional function deteriorated the goodness-of-fit of all models, indicating stochastic trait evolution.

Table 3 Quantitative State Speciation and Extinction (QuaSSE) analysis for body size dependency on diversification. Linear, Sigmoidal and Modal models were used

Biogeography of diversification

Our analyses comparing diversification patterns between the hotspot in lowland Tanzania and the surrounding areas (Fig. 1) that host low species diversities revealed the null, geography-independent model fitted the data best (Table 4). The single hidden state character-independent model (CID GeoSSE) was the best-fit, while the GeoSSE and CID GeoHiSSE models had the next best fit (ΔAICc> 9). These models can be interpreted to highlight the importance of non-biogeographic characters in diversification rate heterogeneity, and thus particular biogeographic regions are not driving diversification. Parameter estimates attained from model-averaging using AICw found that widespread species residing both in and out of the hotspot have the highest diversification rate (0.675 Sp/My) compared to species in the hotspot (0.265 Sp/My) and outside the hotspot (0.410 Sp/My).

Table 4 Geographic and hidden state speciation-extinction (GeoHiSSE) analysis for biogeographic region effect on diversification


Our study provides supportive evidence for the hypothesis that the radiation of Nothobranchius killifish, a clade of exceptionally short-lived vertebrates, has proliferated into over 70 species with minimal niche diversification among them, thus relying on the availability of ‘spatial opportunity’ for species accumulations over time. Interestingly, diversification rates are lower in the one area in which multiple Nothobranchius species coexist (i.e., a species-richness hotspots), suggesting a potential slowdown of the speed of speciation events as species accumulate in the same geographic area. We argue that this link reinforces the view that diversification mediated by availability of ‘spatial opportunity’ is likely to be the dominant engine for Nothobranchius proliferation, relative to diversification mediated by niche divergence (i.e., adaptive radiation). Lineage diversification is coupled with body size evolution, which was also found to be homogeneous through time, with no rapid body size evolution at any point within the genus’ phylogenetic history. Body size disparity is minimal as all species except one are evolving under one adaptive regime.

Tempo and mode of the Nothobranchius radiation

The niche-filling nature of adaptive radiations predicts a diversity-dependent signature of species accumulation in phylogenies. Whereas, diversification via non-adaptive radiation has more relaxed predictions about the diversification dynamics, which can range from early to late bursts of species accumulation. As we suggest above, as long as ‘spatial opportunity’ is available, species will encounter the conditions to radiate allopatrically [15, 18]. For example, the tempo of two non-adaptive radiations – Plethodon salamanders [8, 25] and Rattus rats [26] – have shown an early burst of species accumulation with a diversity-dependent decline, while another non-adaptive clade – Phymaturus lizards [18] – showed a late burst. Our study reveals a constant, diversity-independent rate of diversification over time without evidence for a slow-down. In addition, when species nested in the phylogeny access abundant spatial opportunity they may diversify faster than the sister lineages. We find no evidence of diversification rate heterogeneity, which again supports gradual access to suitable niche space for diversification in allopatry. However, the power to detect rate heterogeneity in a phylogeny of this size is low [84], especially when only few species may be diversifying under a different macroevolutionary regime [85, 86]. The savannah pools replicate the insular biogeographic conditions of archipelago systems and thus, impose dispersal limitations, and so diversification is due to extrinsic factors that promote allopatry (i.e. vicariance or flooding events, [37]). The widespread geographic range of Nothobranchius species indicates that there are environments with viable niches across East Africa. If dispersal was autonomous, as in terrestrial non-adaptive radiations [8, 25, 26], an early burst may take place as spatial opportunity could be accessed early in the clades’ history with a diversity-dependent decline. We argue that non-adaptive radiations cannot be determined based on the tempo of diversification alone; instead, criteria should include niche conservatism and phenotypic stasis.

Nothobranchius species have fundamentally diversified in allopatry, a characteristic of non-adaptive radiations, as confinement to micro-lacustrine environments is unconducive to sympatric speciation [15, 37, 38, 87, 88]. Stochastic environmental events that produce allopatric populations will likely have promoted diversification. East African major climatic fluctuations have taken place during Nothobranchius’ history (~ 10 Mya to present) [89]. The establishment of Africa’s savannah in the late-Miocene (8–9 Mya) would have produced an abundance of spatial opportunity for the killifish, which develop via embryonic diapause [90,91,92] and can cope with seasonal droughts. Environmental fluctuations from 5 Mya to 1 Mya would have produced vicariance and may explain the presence of several sister species arising around this time [93]. These environmental drivers of diversification remain reasonable mechanisms for alternative dating of Nothobranchius diversification, starting either approximately ~ 14 Mya [94] or, as suggested by Furness et al. [95] even at around 50 Mya. Once allopatry is established speciation without ecological divergence potentially occurs via intense genetic drift in found populations without gene flow, driving reproductive isolation [38]. Alternatively, Nothobranchius exhibit a diverse array of karyotypes, therefore chromosomal speciation may cause speciation and may prevent hybridization upon secondary contact [96]; this process is detected in another non-adaptive radiation [26].

Phenotypic evolution

Phenotypic evolution is the result of the selection pressures from biotic and abiotic interactions. In view of that fact, the stasis or disparity of a species’ body size will show how it is interacting with ecological space (niche) [44]. Body size disparity was consistent with variance expected under BM evolution for all of the genus’ history. Phylogenetic niche conservatism across the geographic range resulted in an OU process, with low disparity a result of evolution towards two certain adaptive optima, while additional optima may exist on the adaptive landscape. The close proximity of three peaks from SURFACE may be variance around a single peak, around which 47 species located, or the result of different diet at each peak, as size is a factor in diet composition [97,98,99]. Indeed, a common coexistence of three species often matches the three peaks in body size [45]. We acknowledge, however, that the reliability of the OU model from extant taxa data is questionable as over complex models are selected with high type I error in small trees [76, 100]; furthermore ancestral reconstruction under BM when trait is evolving under an OU model can produce misleading disparity [101]. A single species (N. ocellatus) was found to have largely divergent body size, evolving into predatory species that preys on smaller Nothobranchius. This evolutionary transition has not produced further divergence and N. ocellatus possesses relative wide geographic distribution. Interestingly, in another clade of annual killifish from southern Neotropics, Austrolebias, the same evolutionary transition has led to a group of five predatory species. Whether they all are the result of a single evolutionary transition with consequent diversification into allopatric species [102] or evolved repeatedly [103, 104], is unresolved. Apart from the rare evolution of gigantism, there is relatively low morphological and ecological disparity in killifish; with similar niche space occupied by African and South American genera [42, 105]. Thus, it can be predicted that killifish’s dependency on a specific microhabitat, due to their life history adaptation of annualism, has constrained their ability to exploit ecologically different niche space.

Interestingly, the spatial overlap among the species in the hotspot that our study has identified is hypothesised to have promoted divergence in body size via ecological character displacement to reduce competition for niche space in the spatially restrictive pools, a process that has been identified in Nothobranchius species under certain conditions [98]. Seasonal pools in costal Tanzania, in particular, support coexistence of species that are spatially separated between very shallow marginal habitats and deeper water, while their feeding morphology also suggests some level of trophic differentiation [42]. Most small-bodied species from coastal Tanzania were not included in the Nothobranchius phylogeny [37], hindering insights into the evolutionary transitions to their functional niche. Frequent flooding in coastal regions of Tanzania would prove a viable mode of dispersal from which sympatry is established from secondary contact [106]. Therefore, once secondary contact is established, divergent natural selection could drive body size diversification. As divergent evolution of body size is only strongly supported in one species and not repeatedly found in allopatric species, it is unlikely that the divergent morphology of sympatric species occurs prior to secondary contact.

The phenotypic and biogeographic nature of diversification

The hypothesis that variation in body size influences variation in diversification rates is rejected by our analyses, although results provided contradictory outcomes. Speciation was found to not show a linear or sigmoidal relationship with body size with ES-sim, while QuaSSE analysis found both models as well as a modal model to have a significant relationship with diversification. We accept the results of ES-sim, due to its lower type I error [70], and as body size is predicted to be unimportant when diversification is dependent on abiotic factors (i.e. flooding). GeoHiSSE parameter estimates showed that diversification rate within the hotspot is lower than in the surrounding regions, finding diversification has not the caused the high density of species [107]. The higher species richness is potentially linked to the frequent flooding which increases dispersal into the hotspot but negates divergence through gene flow when populations are reconnected compared to inland Africa where flooding, and connection between suitable habitats is less frequent [32]. The East African rift system uplift caused increased rainfall in inland Tanzania, with coastal regions near Dar es Salaam frequented by flooding [106, 108]. Hidden states allow for adequate testing of the GeoSSE model, however, although diversification rate heterogeneity is found to be minimal the inability to model continuous characters with hidden states leaves the QuaSSE model vulnerable to a high type I error rate, as well as the limited power of SSE models on small phylogenies [109,110,111,112].


Our study investigates the macroevolution of a non-adaptive radiation model system – the African Nothobranchius killifish. The genus has diversified at a constant rate through time due to gradual exploitation of spatial (rather than ecological) opportunity. The evolution of embryo diapause in this continental genus, driven by natural selection arising from the desiccation of their habitats every year, has facilitated their diversification across eastern Africa. However, we suggest that this life history trait has constrained Nothobranchius and other killifish adaptively, most likely due to the substrate dependency of alluvial vertisol soils. Therefore, variance in phenotypic evolution has been relatively constrained as ancestral niche space is conserved from diversification via spatial opportunity, and thus phenotypic disparity is minimal. The exception of a hotspot of species richness which displays higher body size diversity indicates that under certain conditions Nothobranchius may be able to diverge phenotypically. Collectively, our study contributes to the accumulating evidence about a process of evolution that results in the (sometimes) prolific diversification of lineages in the absence of niche-filling, and thus, in the absence of adaptive niche diversification.



Akaike Information Criterion


Bias-corrected AIC


Bayesian Analysis of Macroevolutionary Mixtures


Bayes factors


Brownian motion


constant rate birth-death

DDE + E:

diversity-dependent exponential

DDL + E:

diversity-dependent linear


disparity through time




effective sample size


equal-splits with simulated null model


Geographic State Speciation-Extinction


lineage through-time


Monte Carlo Constant Rate


Markov Chain Monte Carlo


morphological disparity index


million years ago




Quantitative State Speciation Speciation and Extinction

SLmax :

maximum standard length


Species per million years


  1. Simpson GG. The major features of evolution. New York: Columbia University Press; 1953.

    Google Scholar 

  2. Schluter D. The ecology of adaptive radiation. Oxford: Oxford University Press; 2000.

    Google Scholar 

  3. Losos JB. Adaptive radiation, ecological opportunity, and evolutionary determinism. Am Nat. 2010;175:623–39.

    Article  PubMed  Google Scholar 

  4. Yoder JB, Clancey E, Des Roches S, Eastman JM, Gentry L, Godsoe W, Hagey TJ, Jochimsen D, Oswald BP, Robertson J, Sarver BAJ, Schenk JJ, Spear SF, Harmon LJ. Ecological opportunity and the origin of adaptive radiation. J Evol Biol. 2010;23:1581–96.

    Article  CAS  PubMed  Google Scholar 

  5. Rabosky DL, Lovette IJ. Explosive evolutionary radiations: decreasing speciation or increasing extinction through time? Evolution. 2008;62:1866–75.

    Article  PubMed  Google Scholar 

  6. Rabosky DL, Lovette IJ. Density-dependent diversification in north American wood warblers. Proc R Soc B. 2008;275:2363–71.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Phillimore AB, Price TD. Density-Dependent Cladogenesis in Birds. PLoS Biol. 2008;6:e71.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  8. Etienne RS, Haegeman B, Stadler T, Aze T, Pearson PN, Purvis A, Phillimore AB. Diversity-dependence brings molecular phylogenies closer to agreement with the fossil record. Proc R Soc B. 2012;279:1300–9.

    Article  PubMed  Google Scholar 

  9. Price TD, Hooper DM, Buchanan CD, Johansson US, Tietze DT, Alström P, Olsson U, Ghosh-Harihar M, Ishtiaq F, Gupta SK, Martens J, Harr B, Singh P, Mohan D. Niche filling slows the diversification of Himilayan songbirds. Nature. 2014;509:222–5.

    Article  CAS  PubMed  Google Scholar 

  10. Rabosky DL. Ecological limits and diversification rate: alternative paradigms to explain the variation in species richness among clade and regions. Ecol Lett. 2009;12:735–43.

    Article  PubMed  Google Scholar 

  11. Rabosky DL, Glor RE. Equilibrium speciation dynamics in model adaptive radiation of island lizards. PNAS. 2010;107:22178–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Ezard THG, Purvis A. Environmental changes define ecological limits to species richness and reveal the mode of macroevolutionary competition. Ecol Lett. 2016;19:899–906.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Simões M, Breitkreuz L, Alvarado M, Baca S, Cooper JC, Heins L, Herzog K, Lieberman BS. The evolving theory of evolutionary radiations. Trends Ecol Evol. 2016;31:27–34.

    Article  PubMed  Google Scholar 

  14. Wiens JJ, Graham CH. Niche conservatism: integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol Syst. 2005;36:519–39.

    Article  Google Scholar 

  15. Rundell RJ, Price TD. Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends Ecol Evol. 2009;24:394–9.

    Article  PubMed  Google Scholar 

  16. Wiens JJ. Speciation and ecology revisited: phylogenetic niche conservatism and the origin of species. Evolution. 2004;58:193–7.

    Article  PubMed  Google Scholar 

  17. Gittenberger E. What about non-adaptive radiation? Biol J Linn Soc. 1991;43:263–72.

    Article  Google Scholar 

  18. Reaney AM, Saldarriage-Córdoba M, Pincheira-Donoso D. Macroevolutionary diversification with limited niche disparity in a species-rich lineage of cold-climate lizards. BMC Evol Biol. 2018;18:16.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Harmon LJ, Schulter JA, Larson A, Tempo LJB. Mode of evolutionary radiation in Iguanian lizards. Science. 2003;301:961–4.

    Article  CAS  PubMed  Google Scholar 

  20. Seehausen O. African cichlid fish: a model system in adaptive radiation research. Proc R Soc B. 2006;273:1987–98.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Harmon LJ, Losos JB, Davies TJ, Gillespie RG, Gittleman JL, Jennings WB, Kozak KH, McPeek MA, Moreno-Roark F, Near TJ, Purvis A, Ricklefs RE, Schluter D, Schulter JA, Seehausen O, Sidlauskas BL, Torres-Carvajal O, Weir JT, Mooers AØ. Early bursts of body size and shape evolution are rare in comparative data. Evolution. 2010;64:2385–96.

    PubMed  Google Scholar 

  22. Burbrink FT, Ruane S, Pyron RA. When are adaptive radiation replicated in areas? Ecological opportunity and unexceptional diversification in west Indian dipsadine snakes (Colubridae: Alsophiini). J Biogeogr. 2012;39:465–75.

    Article  Google Scholar 

  23. Pincheira-Donoso D, Harvey LP, Ruta M. What defines an adaptive radiation? Macroevolutionary diversification dynamics of an exceptionally species-rich continental lizard radiation. BMC Evol Biol. 2015;15:153.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Wilke T, Benke M, Albrecht C, Bichain JM. The neglected side of the coin: non-adaptive radiations in spring snails (Bythinella spp.). In: Glaubrecht M, editor. Evolution in action. Berlin: Springer; 2010. p. 551–78.

    Chapter  Google Scholar 

  25. Kozak KH, Weisrock DW, Larson A. Rapid lineage accumulation in a non-adaptive radiation: phylogenetic analysis of diversification rates in eastern north American woodland salamanders (Plethodontidae: Plethodon). Proc R Soc B. 2006;273:539–46.

    Article  CAS  PubMed  Google Scholar 

  26. Rowe KC, Aplin KP, Baverstock PR, Moritz C. Recent and rapid speciation with limited morphological disparity in the genus Rattus. Syst Biol. 2011;60:188–203.

    Article  PubMed  Google Scholar 

  27. LaBarbera M. Analyzing body size as a factor in ecology and evolution. Annu Rev Ecol Syst. 1989;20:97–117.

    Article  Google Scholar 

  28. Peters RH. The ecological implication of body size. Cambridge: Cambridge University Press; 1983.

    Book  Google Scholar 

  29. Uyeda JC, Hansen TF, Arnold SJ, Pienaar J. The million-year wait for macroevolutionary bursts. PNAS. 2011;108:15908–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Rabosky DL, Santini F, Eastman J, Smith SA, Sidlauskas B, Chang J, Alfaro ME. Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation. Nat Commun. 2013;4(1958).

  31. Wourms JP. The developmental biology of annual fishes. III. Pre-embryonic and embryonic diapause of variable duration in the eggs of annual fishes. J Exp Zool. 1972;182:389–414.

    Article  CAS  PubMed  Google Scholar 

  32. Watters BR. The ecology and distribution of Nothobranchius fishes. J Am Killifish Assoc. 2009;42:37–76.

    Google Scholar 

  33. Dorn A, Ng’oma E, Janko K, Reichwald K, Polačik M, Platzer M, Cellerino A, Reichard M. Phylogeny, genetic variability and colour polymorphism of an emerging animal model: the short-lived annual Nothobranchius fishes from southern Mozambique. Mol Phylogenet Evol. 2011;61:739–49.

    Article  CAS  PubMed  Google Scholar 

  34. Bartáková V, Reichard M, Janko K, Polačik M, Blažek R, Reichwald K, Cellerino A, Bryja J. Strong population genetic structuring in an annual fish, Nothobranchius furzeri, suggests multiple savannah refugia in southern Mozambique. BMC Evol Biol. 2013;13:196.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Valdesalici S, Cellerino A. Extremely short lifespan in the annual fish Nothobranchius furzeri. Proc R Soc B. 2003;270:S189–91.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Vrtílek M, Žák J, Polačik M, Blažek R, Reichard M. Longitudinal demographic study of wild populations of African annual killifish. Sci Rep. 2018;8:4774.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Dorn A, Musilova Z, Platzer M, Reichwald K, Cellerino A. The strange case of east African annual fishes: aridification correlates with diversification for a savannah aquatic group? BMC Evol Biol. 2014;14:210.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Bartáková V, Reichard M, Blažek R, Polačik M, Bryja J. Terrestrial fishes: rivers are barriers to gene flow in annual fishes from the African savanna. J Biogeogr. 2015;42:1832–44.

    Article  Google Scholar 

  39. Froese R, Pauly D. FishBase. 2016. Accessed Oct 2016.

  40. Wildekamp RH. A world of killies: atlas of the oviparous cyprindontiform fishes of the world. 4th ed. Elyria: American killifish association; 2004.

    Google Scholar 

  41. Reichard M. The evolutionary ecology of African annual fishes. In: Berois N, García G, de Sá R, editors. Annual fishes: life history strategy, diversity, and evolution. Boca Raton: CRC Press; 2015. p. 133–58.

    Chapter  Google Scholar 

  42. Costa WJEM. Comparative morphology, phylogeny, and classification of African seasonal killifishes of the tribe Nothobranchiini (Cyprinodontiformes: Aplocheilidae). Zoo J Linn Soc. 2018;184:115–35.

  43. Seegers L. Aqualog Killifishes of the World: Old World Killis II. Mörfelden-Walldorf: Aquaristik; 1997.

    Google Scholar 

  44. Pincheira-Donoso D, Tregenza T, Butlin RK, Hodgson D. Sexes and species as rival units of niche saturation during community assembly. Glob Ecol Biogeogr. 2018;27:593–603.

    Article  Google Scholar 

  45. Reichard M, Janáč M, Polačik M, Blažek R, Vrtílek M. Community assembly in Nothobranchius annual fishes: nested patterns, environmental niche and biogeographic history. Ecol Evol. 2017;7:2294–306.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Wildekamp RH. A world of killies: atlas of the oviparous cyprindontiform fishes of the world. 1st ed. Elyria: American Killifish Association; 1994.

    Google Scholar 

  47. Near TJ, Dornburg A, Eytan RI, Keck BP, Smith WL, Kuhn KL, Moore JA, Price SA, Burbrink FT, Friedman M, Wainwright PC. Phylogeny and tempo of diversification in the superradiation of spiny-raying fishes. PNAS. 2013;109:13698–703.

    Article  Google Scholar 

  48. Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;63:289–90.

    Article  Google Scholar 

  49. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2017.

    Google Scholar 

  50. Pybus OG, Harvey PH. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc R Soc B. 2000;267:2267–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Pybus OG, Rambaut A, Holmes EC, Harvey PH. New inferences from tree shape: numbers of missing taxa and population growth rates. Syst Biol. 2002;51:881–8.

    Article  CAS  PubMed  Google Scholar 

  52. Rabosky DL. LASER: a maximum likelihood toolkit for detecting temporal shifts in diversification rates from molecular phylogenies. Evol Bioinforma. 2006;2:247–50.

    Article  Google Scholar 

  53. Cusimano N, Renner SS. Slowdowns in diversification rates from real phylogenies may not be real. Syst Biol. 2010;59:458–64.

    Article  PubMed  Google Scholar 

  54. Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control. 1974;19:716–23.

    Article  Google Scholar 

  55. Burnham KP, Anderson DR. Model Selection and Multimodel inference: a practical information-theoretic approach. New York: Springer; 2002.

    Google Scholar 

  56. Etienne RS, Pigot AL, Phillimore AB. How reliably can we infer diversity-dependent diversification from phylogenies? Methods Ecol Evol. 2016;7:1092–9.

    Article  Google Scholar 

  57. Rabosky DL. Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS One. 2014;9:e89543.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  58. Rabosky DL, Grundler M, Anderson C, Title P, Shi JJ, Brown JW, Huang H, Larson JG. BAMMtools: an R package for the analysis of evolutionary dynamics on phylogenetic trees. Methods Ecol Evol. 2014;5:701–7.

    Article  Google Scholar 

  59. FitzJohn RG, Maddison WP, Otto SP. Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies. Syst Biol. 2009;58:595–611.

    Article  PubMed  Google Scholar 

  60. Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6:7–11.

    Google Scholar 

  61. Shi JJ, Rabosky DL. Speciation dynamics during the global radiation of extant bats. Evolution. 2015;69:1528–45.

    Article  PubMed  Google Scholar 

  62. Kass RE, Raftery AE. Bayes factors. J Am Stat Assoc. 1995;90:773–95.

    Article  Google Scholar 

  63. Harmon LJ, Weir JT, Brock CD, Glor RE, Challenger W. GEIGER: investigating evolutionary radiation. Bioinformatics. 2008;24:129–31.

    Article  CAS  PubMed  Google Scholar 

  64. Pennell MW, Eastman JM, Slater GJ, Brown JW, Uyeda JC, FitzJohn RG, Alfaro ME, Harmon LJ. Geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics. 2014;30:2216–8.

    Article  CAS  PubMed  Google Scholar 

  65. Felsenstein J. Phylogenies and the Comparative method. Am Nat. 1985;125:1–15.

    Article  Google Scholar 

  66. Murrell DJ. A global envelope test to detect non-random bursts of trait evolution. Methods Ecol Evol. 2018;9:1739–48.

    Article  Google Scholar 

  67. Slater GJ, Price SA, Santini F, Alfaro ME. Diversity versus disparity and the radiation of modern cetaceans. Proc R Soc B. 2010;277:3097–104.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Revell LJ. Two new graphical methods for mapping trait evolution on phylogenies. Methods Ecol Evol. 2013;4:754–9.

    Article  Google Scholar 

  69. Revell LJ. Phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012;3:217–23.

    Article  Google Scholar 

  70. Hansen TF. Stabilizing Selection and the Comparative analysis of adaptation. Evolution. 1997;51:1341–51.

    Article  PubMed  Google Scholar 

  71. Butler MA, King AA. Phylogenetic Comparative analysis: a modeling approach for adaptive evolution. Am Nat. 2004;164:683–95.

    Article  PubMed  Google Scholar 

  72. Pagel M. Inferring the historical patterns of biological evolution. Nature. 1999;401:877–84.

    Article  CAS  PubMed  Google Scholar 

  73. Gould SJ, Eldredge N. Punctuated equilibrium comes of age. Nature. 1993;366:223–7.

    Article  CAS  PubMed  Google Scholar 

  74. Ingram T, Mahler DL. SURFACE: detecting convergent evolution from comparative data by fitting Ornstein-Uhlenbeck models with stepwise Akaike information criterion. Methods Ecol Evol. 2013;4:416–25.

    Article  Google Scholar 

  75. Khabbazian M, Kriebel R, Rohe K, Ané C. Fast and accurate detection of evolutionary shifts in Ornstein-Uhlenbeck models. Methods Ecol Evol. 2016;7:811–24.

    Article  Google Scholar 

  76. Ho LST, Ané C. Intrinsic inference difficulties for trait evolution with Ornstein-Uhlenbeck models. Methods Ecol Evol. 2014;5:1133–46.

    Article  Google Scholar 

  77. Harvey MG, Rabosky DL. Continuous traits and speciation rates: alternatives to state-dependent diversification models. Methods Ecol Evol. 2017;9:984–93.

    Article  Google Scholar 

  78. FitzJohn RG. Quantitative traits and diversification. Syst Biol. 2010;59:619–33.

    Article  PubMed  Google Scholar 

  79. Maddison WP, Midford PE, Otto SP. Estimating a binary character’s effect on speciation an extinction. Syst Biol. 2007;56:701–10.

    Article  PubMed  Google Scholar 

  80. FitzJohn RG. Diversitree: comparative phylogenetic analyses of diversification in R. Methods Ecol Evol. 2012;3:1084–92.

    Article  Google Scholar 

  81. Goldberg EE, Lancaster LT, Ree RH. Phylogenetic inference of reciprocal effects between geographic range evolution and diversification. Syst Biol. 2011;60:451–65.

    Article  PubMed  Google Scholar 

  82. Beaulieu JM, O’Meara BC. Detecting hidden diversification shifts in models of trait-dependent speciation and extinction. Syst Biol. 2016;65:583–601.

    Article  PubMed  Google Scholar 

  83. Caetano DS, O’Meara BC, Beaulieu JM. In press hidden state models improve state-dependent diversification approaches, including biogeographic models. Evolution.

  84. Kodandaramaiah U, Murali G. What affects power to estimate speciation rate shifts? PeerJ. 2018;6:e5495.

    Article  PubMed  PubMed Central  Google Scholar 

  85. Moore BR, Höhna S, May MR, Rannala B, Huelsenbeck JP. Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures. PNAS. 2016;113:9569–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. Rabosky DL, Mitchell JS, Chang J. Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models. Syst Biol. 2017;66:477–98.

    Article  PubMed  PubMed Central  Google Scholar 

  87. Kisel Y, Barraclough TG. Speciation has a spatial scale that depends on levels of gene flow. Am Nat. 2010;175:316–34.

    Article  PubMed  Google Scholar 

  88. Seehausen O, Wagner CE. Speciation in freshwater fishes. Annu Rev Ecol Evol Syst. 2014;45:621–51.

    Article  Google Scholar 

  89. Ivory SJ, Blome MW, King JW, McGlue MM, Cole JE, Cohen AS. Environmental change explains cichlid adaptive radiation at Lake Malawi over the past 1.2 million years. PNAS. 2016;113:11895–900.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Murphy WJ, Collier GE. A molecular phylogeny for aplocheiloid fishes (Atherinomorpha, Cyprinodontiformes): the role of vicariance and the origins of annualism. Mol Biol Evol. 1997;14:790–9.

    Article  CAS  PubMed  Google Scholar 

  91. deMenocal PB. African climate change and faunal evolution during the Pliocene-Pleistocene. Earth Plant Sci Lett. 2004;220:3–24.

    Article  CAS  Google Scholar 

  92. Trauth MH, Maslin MA, Deino A, Strecker MR. Late Cenozoic moisture history of East Africa. Science. 2005;309:2051–3.

    Article  CAS  PubMed  Google Scholar 

  93. Trauth MH, Larrasoaña JC, Mudelsee M. Trends, rhythms and events in Plio-Pleistocene African climate. Quat Sci Rev. 2009;28:399–411.

    Article  Google Scholar 

  94. Helmstetter AJ, Papadopulos AST, Igea J, Van Dooren TJM, Leroi AM, Savolainen V. Viviparity stimulations diversification in an order of fish. Nat Commun. 2016;7:11271.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Furness AI, Reznick DN, Springer MS, Meredith RW. Convergent evolution of alternative developmental trajectories associated with diapause in African and South American killifish. 2015;282:20142189.

  96. Krysanov E, Demidova T, Nagy B. Divergent karyotypes of the annual killifish genus Nothobranchius (Cyprinodontiformes, Nothobranchiidae). Comp Cytogenet. 2016;10:439–45.

    Article  PubMed  PubMed Central  Google Scholar 

  97. Polačik M, Reichard M. Diet overlap among three sympatric African annual killifish species (Nothobranchius spp.) from Mozambique. J Fish Biol 2010;77:754–768.

  98. Polačik M, Harrod C, Blažek R, Reichard M. Trophic niche partitioning in communities of African annual fish: evidence from stable isotopes. Hydrobiologia. 2014;721:99–106.

    Article  CAS  Google Scholar 

  99. Cellerino A, Valenzano DR, Reichard M. From the bush to the bench: the annual Nothobranchius fishes as a new model system in biology. Biol Rev. 2016;91:511–33.

    Article  PubMed  Google Scholar 

  100. Cooper N, Thomas GH, Venditti C, Meade A, Freckleton RP. A cautionary note on the use of Ornstein Uhlenbeck models in macroevolutionary studies. Biol J Linn Soc. 2016;118:64–77.

    Article  Google Scholar 

  101. Mahler DL, Weber MG, Wagner CE, Ingram T. Pattern and process in the Comparative study of convergent evolution. Am Nat. 2017;190:S13–28.

    Article  PubMed  Google Scholar 

  102. Costa WJEM. Parallel evolution in ichthyophagous annual killifishes of South America and Africa. Cybium. 2011;35:39–46.

    Google Scholar 

  103. García G, Wlasiuk G, Lessa EP. High levels of mitochondrial cytochrome b divergence and phylogenetic relationships in the annual killifishes of the genus Cynolebias (Cyprinodontifomes, Rivulidae). Zool J Linnean Soc. 2000;129:93–110.

    Article  Google Scholar 

  104. Van Dooren TJM, Thomassen HA, Smit F, Helmstetter AJ, Savolainen V. A case for sympatric speciation by cannibalism in South-American annual killifish (Austrolebias). 2018; doi: https//

  105. Furness AI. The evolution of an annual life cycle in killifish: adaptation to ephemeral aquatic environments through embryonic diapause. Biol Rev. 2015;91:796–812.

    Article  PubMed  Google Scholar 

  106. Kalala AM, Msanya BM, Amuri NA, Semoka JM. Pedological characterization of some typical alluvial soils of Kilombero District, Tanzania. Am J Agric For. 2017;5(1):11.

    Google Scholar 

  107. Huang D, Goldberg EE, Chou LM, Roy K. The origin and evolution of coral species richness in a marine biodiversity hotspot. Evolution. 2018;72:288–302.

    Article  PubMed  Google Scholar 

  108. Sepulchre P, Ramstein G, Fluteau F, Schuster M, Tiercelin JJ, Brunet M. Tectonic uplift and East Africa Aridification. Science. 2006;313:1419–23.

    Article  CAS  PubMed  Google Scholar 

  109. Davis MP, Midford PE, Maddison W. Exploring power and parameter estimation of the BiSSE method for analysing species diversification. BMC Evol Biol. 2013;13:38.

    Article  PubMed  PubMed Central  Google Scholar 

  110. Maddison WP, FitzJohn RG. The unsolved challenge to phylogenetic correlation tests for categorical characters. Syst Biol. 2015;64:127–36.

    Article  PubMed  Google Scholar 

  111. Rabosky DL, Goldberg EE. Model inadequacy and mistaken inferences of trait-dependent speciation. Syst Biol. 2015;64:340–55.

    Article  CAS  PubMed  Google Scholar 

  112. Feldman A, Sabath N, Pyron RA, Mayrose I, Meiri S. Body sizes and diversification rates of lizards, snakes, amphisbaenians and the tuatara. Glob Ecol Biogeogr. 2016;25:187–97.

    Article  Google Scholar 

Download references


We thank to R. Etienne for help with using DDD parametric bootstrapping technique, P. Assheton for assistance in using an alternative statistical test in ES-sim analysis, and F. Grattarola for assisting with the mapping of species-richness. We would also like to thank two anonymous reviewers, whose comments greatly improved the paper.


This project was not funded. MR was funded through Czech Science Foundation grant 16-00291S – the funding body played no role in the design of the study, in the collection, analysis, and interpretation of data, or in writing the manuscript.

Availability of data and materials

The datasets used in this research are included within the article and in additional supporting files.

Author information

Authors and Affiliations



Original idea: JWL, DP-D, design of study: JWL, DP-D, data collection: MR, JWL, data analyses: JWL, wrote the manuscript: JWL, DP-D, MR. All authors read and approved the manuscript.

Corresponding author

Correspondence to Daniel Pincheira-Donoso.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

List of references from which the body size data is taken. (DOCX 15 kb)

Additional file 2:

Body size data of 48 Nothobranchius species used in this study. (PDF 35 kb)

Additional file 3:

List of references from which locality data is taken. (PDF 4 kb)

Additional file 4:

Geographic locality data for 49 Nothobranchius species corresponding to the species in the phylogeny. (DOCX 18 kb)

Additional file 5:

Nothobranchius phylogeny (n = 49) with colours representing geographic region; blue: endemic to hotspot, red: widespread, green: non-endemic to hotspot. (DOCX 16 kb)

Additional file 6:

(A) The distribution of logarithms of likelihood ratio of constant rate Birth-Death (crBD) and (B) diversity-dependent diversification linear speciation and extinction (DDL + E) from parametric bootstrapping. Both models were run with 22 missing species and a significance value (α) set to 0.05. The blue arrow shows the logarithm of likelihood ratio for the significance value (2.99). The black arrow shows the logarithm of likelihood ratio for the data (0.12). (DOCX 31 kb)

Additional file 7:

(A) The distribution of logarithms of likelihood ratio of constant rate Birth-Death (crBD) and (B) diversity-dependent diversification linear speciation and extinction (DDL + E) from parametric bootstrapping. Both models were run with 50 missing species and a significance value (α) set to 0.05. The blue arrows is the logarithm ratio for the significance value (3.07). The black arrow shows the logarithm of likelihood ratio for the data (0.10). (PDF 27 kb)

Additional file 8:

Phylogenetic tree with speciation rate shown by colour gradient. Black dot depict significant shifts in diversification rate. Frequency (f) of the model selected show above respective model. (PDF 181 kb)

Additional file 9:

Left, ℓ1ou selective regimes, red line indicates divergence to a different selective regime, with asterisk and strength of shift on the node. Right, SURFACE selective regimes plotted regime shift are numbered on nodes, regimes are differentiate by colour, grey is divergence and red is convergence, (1) is the ancestral regime and thus is not shown on the tree. (PDF 182 kb)

Additional file 10:

Phylogenetic tree with rate of phenotypic evolution shown by colour gradient. Black dots show significant shifts in trait evolution and frequency (f) of model selected show above respective model. (PDF 22 kb)

Additional file 11:

Scatter plot of the logarithm of the speciation rate against the body size of Nothobranchius species. Speciation rate is measured using the ES-method. (PDF 36 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lambert, J.W., Reichard, M. & Pincheira-Donoso, D. Live fast, diversify non-adaptively: evolutionary diversification of exceptionally short-lived annual killifishes. BMC Evol Biol 19, 10 (2019).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: