Big fish, little divergence: phylogeography of Lake Tanganyika's giant cichlid, Boulengerochromis microlepis.

The largely endemic cichlid species flocks of the East African Great Lakes are among the prime examples for explosive speciation and adaptive radiation. Speciation rates differ among cichlid lineages, and the propensity to radiate has been linked to intrinsic and extrinsic factors such as sexual selection and ecological opportunity. Remarkably, only one cichlid tribe-the Boulengerochromini-comprises just a single species, Boulengerochromis microlepis, a predominantly piscivorous endemic of Lake Tanganyika and the world's largest cichlid. While the lineage diverged from its closest relatives at the onset of the Lake Tanganyika radiation >8 MYA, mitochondrial control region sequences collected in this study dated the most recent common ancestor of B. microlepis to ~60-110 KYA. There was no evidence of phylogeographic structure in the lake-wide sample. Patterns of genetic diversity and demographic analyses were consistent with slow and steady population growth throughout the reconstructed timescale. Additionally, the shallow divergence within the species may be related to a possibly large variance in reproductive success in this highly fecund species. Trophic niche space restriction by sympatric piscivores, lack of geographic structure, low potential for sexual selection arising from the monogamous mating system and extinction may have contributed to keeping the lineage monotypic.

Introduction

The largely endemic cichlid species flocks of the East African Great Lakes represent one of the most outstanding examples of explosive speciation and adaptive radiation (Fryer & Iles, 1972; Salzburger,2009; Sturmbauer et al., 2011). It is one characteristic of adaptive radiation that not all lineages that seeded these radiations actually diversified and if they did, to a similar extent. Some lineages are species-rich whereas others have produced only a few species. Such differences in species richness among subgroups of an endemic species assemblage have been linked to a different potential to radiate due to preadaptive phenomena such as the possession of particular lineage-specific key innovations (e.g., Mayr, 1963; Liem, 1973; Sanderson & Donoghue, 1994). Thus, the particular anatomy of the pharyngeal apophysis providing a second set of jaws and a highly specialized reproductive behavior (Liem, 1973; Crapon de Crapona, 1986) have been proposed as key innovations underlying the evolutionary success of cichlids as compared to other sympatric fish lineages. However, not all cichlid lineages that colonized emerging lakes diversified to a similar extent. Rather, some lineages are species-rich whereas others have produced a few or even a single species, despite similar age. These differences in species richness have been attributed to differences in ecological characteristics, in particular habitat preferences of the different cichlid lineages. In Lake Malawi, for example, highly stenotopic rock-dwelling cichlid species tend to be philopatric and gene flow between geographically close populations separated by habitat discontinuities in the form of sandy stretches of shoreline, river estuaries, or deep water is low (e.g., van Oppen et al., 1997; Markert et al., 1999; Rico & Turner, 2002; Pereyra et al., 2004). Cichlids occupying shallow sandy habitats show less population structure (Taylor & Verheyen, 2001; Pereyra et al., 2004; Anseeuw et al., 2011), and benthopelagic species show hardly any or no significant genetic structure at all on a lake-wide scale (Shaw et al., 2000; Genner et al., 2008, 2010b). For Lake Tanganyika, thus far, population genetic data are only available for rock-dwelling cichlid species and taxa preferring the intermediate habitat at the rock–sand interface in shallow water. Nevertheless, similar trends have been observed in Lake Tanganyika: Highly stenotopic rock-dwellers show high levels of geographic structuring, whereas less stenotopic taxa and taxa inhabiting the intermediate habitat are less structured (Meyer et al., 1996; Taylor et al., 2001; Duftner et al., 2006; Koblmüller et al., 2007, 2009,2011; Koch et al., 2007; Sefc et al., 2007; Takahashi et al., 2009; Wagner & McCune, 2009; Nevado et al.,2009, 2013; Van Steenberge et al., 2014). The overall pattern suggests that dispersal capacity is modulated by adaptation to particular habitat types and ecological niches during the radiation process (Sturmbauer,1998). Assuming that geographic population structuring is an important prerequisite for allopatric speciation (Turelli et al., 2001), lineages composed of stenotopic rock-dwellers should be particularly species-rich, whereas lineages comprising mobile benthopelagic taxa should be comparatively species poor. Indeed, in each of the three East African Great Lakes, Tanganyika, Malawi, and Victoria, patterns of cichlid diversity appear to follow this model (Turner, 1996; Seehausen et al., 1997; Koblmüller et al., 2008). Nonetheless, that a particular lineage did not diversify at all, is quite unusual. This is the case for the endemic Lake Tanganyika cichlid tribe Boulengerochromini, which only comprises the species Boulengerochromis microlepis, one of the few benthopelagic top-predators of the lake.

Previous studies have shown that environmental influences, mainly climatically driven lake level fluctuations (Cohen et al., 1997; Scholz et al., 2007; McGlue et al., 2008), have synchronized population divergence and patterns of past population size changes of stenotopic rock-dwelling cichlids within and across Lake Tanganyika and Malawi (Sturmbauer et al., 2001; Genner et al., 2010a; Koblmüller et al.,2011; Nevado et al., 2013). Whether population histories of highly mobile benthopelagic or truly pelagic species have been equally affected by large lake level fluctuations remains largely unknown.

Lake Tanganyika’s giant cichlid, B. microlepis, the only member of the tribe Boulengerochromini (Takahashi, 2003), is distantly related to all other cichlid species of this lake, so that it is assumed that its ancestor colonized the lake as one of the seeding lineages of the radiation (Salzburger et al., 2002). The phylogenetic relationships of B. microlepis are still tentative and a matter of discussion. Molecular phylogenetic analyses based on mtDNA (Klett & Meyer, 2002; Koblmüller et al., 2005) placed B. microlepis as closely related to the tribe Tilapiini (sensu Dunz & Schliewen, 2013), which does not occur in the lake itself but has a widespread distribution from the Congo basin and Lake Malawi southwards to South Africa. Nuclear data, on the other hand, imply that B. microlepis is part of the lacustrine assemblage, representing either its most basal lineage or being most closely related to the mainly bathypelagic genera Bathybates, Hemibates, and Trematocara (Nishida, 1997; Muschick et al., 2012; Dunz & Schliewen, 2013; Friedman et al., 2013). However, independent of the genetic marker used to infer phylogenetic relationships and the approaches employed to estimate divergence times, the Boulengerochromini are an old lineage that diverged from its closest relatives at least ~8 MYA (Genner et al.,2007; Koblmüller et al., 2008; Schwarzer et al., 2009; Friedman et al., 2013). With a total length exceeding 80 cm (Matthes, 1961; personal observation by SK), B. microlepis is considered the world’s largest cichlid fish species (Coulter, 1991). It is found over a rather wide depth range from the very shallow water down to a depth of 150 m in all kinds of habitat—though only infrequently over rocky substrate—preying mainly upon various fish species, but also on crabs, shrimps, molluscs, and insect larvae (Kawabata & Mihigo,1982; Coulter, 1991; Bayona, 1991a; Sturmbauer et al., 2008). B. microlepis is a substrate breeding, probably semelparous species, breeding in rather shallow sandy or intermediate habitat, with both male and female guarding their numerous (up to 12,000) fry (Poll, 1956; Kuwamura, 1986; Konings, 1998; Büscher, 2009). Occasionally, large numbers of breeding pairs can be observed in very shallow water (personal observation by SK). As a highly valued and pricy staple fish, B. microlepis is heavily targeted by local fishermen. One of its popular local names, “English fish,” reflects its popularity among well-to-do, typically English-speaking, people. Previous catch statistics for beach-seining indicate that the species used to be very abundant and constituted a remarkably high portion of the total yield (Bayona, 1991b). As beach-seining has been prohibited over most parts of Lake Tanganyika, the species is nowadays mainly caught by hook and line, often in very deep water (personal observation by SK). Currently, it does not seem to be threatened by overfishing as catch statistics for the very southern part do not indicate a population decline over the last decade (Sinyinza, unpublished data).

To date there have been no studies of spatial or temporal population genetic structure of B. microlepis or any other large predatory and highly mobile Lake Tanganyika cichlid fish species. Such studies may not only increase our understanding about factors and processes shaping intraspecific diversity, but might also help to identify potentially segregated fish stocks and thus provide important insights for conservation and fisheries management. Here, we characterize the genetic diversity of B. microlepis and reconstruct its phylogeographic structure and past population size trajectories based on DNA sequences of the most variable part of the mitochondrial control region. The findings are discussed in the light of the hydrologic history of Lake Tanganyika and the biological characteristics of the species.

Materials and methods

Fin clips were taken from 88 individuals of B. microlepis collected with gill net or harpoon from 13 localities in Lake Tanganyika, or obtained at local fish markets in Bujumbura and Mpulungu or from artisanal fishermen on the lake, during several field trips between 2001 and 2013 (Fig. 1; Supplementary Table 1), and preserved in 96% ethanol. Whole genomic DNA was extracted following a standard Chelex protocol (Walsh et al., 1991). The most variable part of the mitochondrial control region was amplified and sequenced according to the protocols described in Koblmüller et al. (2011) and Duftner et al. (2005), respectively. The primers used for PCR and chain termination sequencing were L-Pro-F_Tropheus (Koblmüller et al., 2011) and TDK-D (Lee et al.,1995). DNA fragments were purified with Sephadex™ G-50 (Amersham Biosciences) and visualized on an ABI 3130×l capillary sequencer (Applied Biosystems). Sequences were aligned by eye in Mega 5.1 (Tamura et al., 2011). The length of the final alignment (including gaps) was 358 bp. Sequences are deposited in GenBank under the accession numbers KJ438258–KJ438286.

Fig. 1

a Map of Lake Tanganyika showing the sampling localities. Numbers in parentheses refer to sample sizes. The three deep-water basins are indicated by gray shading. b Median-joining (MJ) network of B. microlepis haplotypes. Circle sizes are proportional to haplotype frequency and connecting lines are proportional to mutation events between haplotypes. Different colors refer to different sampling localities as shown in a

Phylogenetic relationships among haplotypes were visualized by a full median-joining (MJ) network (Bandelt et al., 1999), by setting the weighted genetic distance parameter ε to 90 [10 (default character weight) × 9 (the observed maximum number of pairwise differences)] and activating the “MJ square” option, followed by maximum parsimony post-processing (Polzin & Daneschmand, 2003) as implemented in Network (version 4.6; available at www.fluxus-engineering.com/sharenet.htm). Haplotype (Hd) and nucleotide diversity (p) were calculated in DnaSP 5.10 (Librado & Rozas, 2009). Spatial genetic structure of B. microlepis was analyzed by spatial analysis of molecular variance using SAMOVA 1.0 (Dupanloup et al., 2002), which defines groups of populations that are maximally differentiated from each other in a geographically homogeneous environment by maximizing the proportion of genetic variance (FCT; Wright,1978) among K groups. Analyses for K = 2–7 were conducted with 100 simulated annealing runs. For each K, the configuration with the largest FCT values after all 100 independent simulated annealing processes (5,000 iterations) was retained as the best spatially genetic clustering. To test for signals of past population expansion, we calculated a mismatch distribution as well as Fu’s Fs (Fu, 1997) and Tajima’s D (Tajima,1989) in Arlequin 3.5.1.2 (Excoffier & Lischer, 2010). The fit of the observed mismatch distribution to the expectations based on growth parameter estimates was evaluated by the sum of squared differences (SSD) and the raggedness index (rg). Furthermore, past population size trajectories were inferred by means of a Bayesian coalescent approach [Gaussian Markov random field (GMRF) skyride tree prior; Minin et al., 2008] as implemented in BEAST 1.8.0 (Drummond & Rambaut, 2007). We employed the model of molecular evolution selected by the Bayesian information criterion (BIC) in jModelTest 0.1 (Posada, 2008), assuming a strict molecular clock and a substitution rate of 0.0324 and alternatively 0.057 per site per MY (Genner et al., 2007, 2010a; Koblmüller et al., 2009). Two independent MCMC runs of 1 million generations each were conducted, sampling every 1000th step with a burn-in of the first 10% of sampled generations. Verification of effective sample sizes (ESS > 200 for all parameters), trace of MCMC runs, and visualization of past demographic changes were done in Tracer 1.5 (Rambaut & Drummond, 2009).

Results

In total 29 haplotypes were detected in 88 individuals. High haplotype (Hd) diversity (Hd = 0.871) contrasts with low nucleotide diversity (π) and little genetic divergence between haplotypes (π = 0.00951; maximum number of pairwise differences between haplotypes = 9). The MJ network did not indicate geographic structure. Haplotypes are shared between geographically distant localities, sometimes even from opposite ends of the lake (e.g., Bujumbura, Burundi, and Chimba, Zambia) (Fig. 1b). SAMOVA identified four distinct clusters as the grouping with the highest FCT (FCT = 0.37894, P < 0.001) (Supplementary Fig. 1). However, other solutions (with K ranging from 2 to 7) achieved similar FCT values, and all solutions consisted of one large group spanning a lakewide distribution, and variable numbers of small and geographically restricted groups. Importantly, the scenario of K = 1, i.e., one single panmictic population, cannot be tested in the SAMOVA framework. Based on the structure of the haplotype network and the failure of SAMOVA to group individuals geographically, we conclude that B. microlepis is connected by high levels of gene flow across Lake Tanganyika. This does not preclude the possibility of population genetic differentiation, which might be detectable with appropriate sampling.

The fit of the observed mismatch distribution to the expectations based on growth parameter estimates, with nonsignificant SSD and rg values (Fig. 2a), and a significantly negative Fu’s Fs (Fs = −16.162, P < 0.001) are consistent with recent populationgrowth. Tajima’s D was negative, as expected for recent population growth, but not significantly different from zero (D = −0.946, P = 0.176; but note that the power to detect population expansion is considerably lower for Tajima’s D as compared to Fu’s Fs; Ramos-Onsins & Rozas, 2002). Congruently, a GMRF skyride indicated slight continuous population growth over the last approximately 35–60 KY (dependingon the substitution rate assumed) (Fig. 2b). However, the width of the 95% HPD intervals does not reject a scenario with constant population size. The time to the most recent common ancestor (MRCA) was estimated as 63.349–111.447 KYA (respective 95% HPDs: 35.679–94.256, 62.769–165.821), and present female effective population size as 1,203,169–2,110,173(respective 95% HPDs: 413,400–4,166,021; 725,041–7,306,560), both depending on the assumed substitution rate.

Fig. 2

Demographic history of B. microlepis. a Mismatch distribution. Black columns represent the observed frequency of pairwise differences. Gray lines refer to the expected distribution based on parameter estimates and their 95% confidence limits simulated under a model of population growth. Sum of squared differences (SSD) and raggedness index (rg) and their respective P values are given to describe the fit of the observed mismatch distribution to the expectations based on growth parameter estimates. b GMRF skyride plot of past population size trajectories. The skyride plot shows the product of female effective population size (fN) and mutation rate (μ) through time, assuming substitution rates of 5.7 and 3.24% per site per MY (Genner et al., 2007, 2010a; Koblmüller et al., 2009). The thick black line represents the median values; the thin gray lines denote 95% highest posterior density (HPD) intervals

Discussion

Analysis of the most variable region of the mitochondrial control region revealed shallow population divergence and a lack of phylogeographic structure in B. microlepis, with several haplotypes shared between individuals from very distant parts of the lake. Mitochondrial divergence within B. microlepis lies within the range of divergence found in regional samples of stenotopic littoral cichlids for which the same mitochondrial DNA fragment was analyzed. For example, 2.5% maximum divergence in B. microlepis from across the entire lake compare to 2.2-4.2% (mean 3.5%) maximum divergence in Variabilichromis moorii, Eretmodus cyanostictus, Ophthalmotilapia ventralis, Neolamprologus caudopunctatus, and Perissodus microlepis, all sampled within southern Zambia (Duftner et al., 2006; Koblmüller et al., 2007, 2009; Sefc et al., 2007). In contrast, maximum divergence in a lake-wide sample of Neolamprologus pulcher/brichardi amounted to 8.7% (Duftner et al., 2007). Our sampling is strongly biased toward the very south of Lake Tanganyika, with fewer individuals and locations sampled in the middle and northern parts of the lake. Therefore, tests for population structure suffered from reduced power, and it is possible that deviations from panmixis remained undetected in this study.

Boulengerochromis microlepis is one of the toppredators in Lake Tanganyika. As such, B. microlepis is highly mobile and not restricted to a particular type of habitat or depth and apparently disperses over long distances without habitat-imposed restrictions. Similar to B. microlepis, phylogeographic structure is also lacking in benthopelagic cichlids in the genera Rhamphochromis and Diplotaxodon from Lake Malawi (Shaw et al., 2000; Genner et al., 2008,2010b). Both genera comprise medium-sized to large offshore predators that feed in the water column, with large fish being almost exclusively piscivorous (Turner et al., 2002). One Rhamphochromis species studied in detail, R. longiceps, uses lagoons and satellite lakes as nursery areas, which might promote population differentiation if fish remained close to their breeding grounds. Yet, no geographic population structure has been detected, which has been explained by an apparent lack of homing behavior, and dispersal mainly determined by prey availability and opportunistic selection of suitable nursery grounds (Genner et al., 2008). Similarly, four of the eight Diplotaxodon species studied so far show no spatial population subdivision, whereas the other four species show slight but significant spatial genetic differences among breeding grounds, indicating natal homing to breeding grounds in these species (Genner et al., 2010b). Fish in our sample were mainly caught by hook and line (by artisanal fishermen) and are therefore unlikely to comprise a large portion of breeding and fry-guarding individuals, which purportedly do not feed at all (Poll,1956; Kuwamura, 1986; Fohrman, 1994) and hence cannot be caught with bait. Consequently, our sample does not allow inferences on natal homing, as these large and mobile fish may have been collected far from their breeding grounds.

Despite the old evolutionary age of the B. microlepis-lineage, which dates back to the very onset of the Lake Tanganyika radiation, its mtDNA genealogy coalesces in the very recent past. The high haplotype diversity and low nucleotide diversity observed in B. microlepis, combined with the old age of the lineage, are typical for species that experienced population growth following a period of low effective population size (Grant & Bowen, 1998). Indeed, our data do not reject a scenario of population growth in the recent past, although the signature of expansion is not entirely compelling (Fig. 2). A possible scenario includes a decline in numbers during the East African megadrought period at ~135–75 KYA, when the lake level dropped by up to 435 m below the present level (which was still not enough to separate the three subbasins into three distinct lakes; McGlue et al., 2008) and the inhabitable lake area was reduced considerably. The subsequent rise of the lake level brought an expansion of the available habitat, which may have triggered a steady rise of B. microlepis stock numbers.

Additionally, large variance in reproductive success and stochasticity of juvenile mortality might have contributed to the low levels of mtDNA divergence (Hedgecock, 1994). Recent theoretical work has demonstrated that large variance in reproductive success produces shallow mitochondrial genealogies in species with high reproductive potential (e.g., Hedrick, 2005; Sargsyan & Wakeley, 2008; Hoban et al., 2013). Empirical support for this hypothesis comes from numerous studies on highly fecund taxa that showed unexpectedly shallow phylogenies and small effective population sizes (e.g., Li & Hedgecock, 1998; Arnason, 2004; Christie et al., 2010; Hedgecock & Pudovkin, 2011). B. microlepis is a highly fecund substrate breeding cichlid that produces up to 12,000 eggs per brood. Parents defend their offspring for at least 9 months. At that time, the offspring has grown to an average length of about 15 cm and has attained a size too big to be preyed upon by most other predators in the lake, except the large lates perches (Latidae) (Konings, 1998; Büscher,2009). The fry are certainly most vulnerable shortly after becoming mobile, when they hover above the nest in a dense cloud. Usually, predators rapidly decimate the brood to a considerable extent (Konings,1998). Large clutches and the long guarding phase predispose B. microlepis to considerable variance in reproductive success. Given the high mobility of this species, haplotypes of particularly successfully reproducing mothers could spread quickly and periodically replace current common ancestors.

Another puzzle, given the age of the lineage, is its failure to seed more than a single species. The other lineages of predatory benthopelagic cichlids in Lakes Malawi and Tanganyika are typically species poor as well, but nevertheless radiated into at few distinct taxa (Turner, 1996; Genner et al., 2008; Koblmüller et al.,2008). A meta-analysis of lacustrine cichlid assemblages suggested that the major factors that predispose cichlids to adaptive radiations are ecological opportunity and intrinsic lineage-specific traits related to sexual selection (Wagner et al., 2012). Compared to the benthopelagic, polygynous mouthbrooders in Lakes Tanganyika and Malawi, sexual selection has probably been weak in B. microlepis, judging by its monogamous mating system and lack of sexual dimorphism. Moreover, the community structure in Lake Tanganyika may also have played a significant role in preventing a radiation of the Boulengerochromini. In Lake Tanganyika, numerous piscivorous fish species from other (cichlid and non-cichlid) lineages coexist with B. microlepis, including the smaller piscivores in the genera Lepidiolamprologus and Bathybates as well as Hemibates stenosoma, and several endemic species of clariid catfish and lates perches that grow considerably larger than B. microlepis. Hence, available niche space of B. microlepis, as determined by body size, is probably restricted by heterospecific trophic competition. The extant piscivorous cichlid lineages colonized Lake Tanganyika roughly simultaneously, shortly after the establishment of deep-water conditions (Koblmüller et al.,2008; no data are available for late perches and clariid catfish). Sexual selection in the mouthbrooding piscivores and geographic isolation among the littoral substrate breeders may have given these lineages a head start in speciation, such that the Boulengerochromini ancestor quickly found the neighboring niches occupied. Contrary to restricted diversification, extinction may have eliminated traces of past speciation in the Boulengerochromini. Unfortunately, no means exist to reconstruct historic diversity beyond the coalescence of the extant mitochondrial haplotypes.

To conclude, we show that the mitochondrial genealogy of B. microlepis is very shallow despite the old age of the lineage and that there is no indication of phylogeographic structuring. The observed patterns of genetic diversity and divergence may be the product of recent population growth and large variance in reproductive success. Low potential for sexual selection and limited trophic niche space may have hindered speciation of the Boulengerochromini. Despite their old age, the Boulengerochromini are the only known cichlid tribe comprising only a single species, the giant cichlid from Lake Tanganyika.