Recurrent changes in cichlid dentition linked to climate-driven lake-level ﬂ uctuations

. Climate-induced habitat change has often been invoked as an important driver of speciation and evolutionary radiation in cichlid ﬁ shes, yet studies linking morphological change directly to long-term environmental ﬂ uctuations are scarce. Here, we track changes through time in the oral dentition of Oreochromis hunteri , the endemic and only indigenous ﬁ sh species inhabiting the East African crater Lake Chala (Kenya/ Tanzania), in relation to climate-driven lake-level ﬂ uctuations during the last 25,000 yr. Fossil teeth were recovered from six periods representing alternating high and low lake-level stands associated with late-Pleistocene and Holocene trends in rainfall and drought. To construct a reference framework for morphological variation in fossil assemblages, we ﬁ rst analyzed the variability in oral tooth morphology within the modern-day ﬁ sh population. These analyses established that, like other cichlids, O. hunteri gradually replace their bicuspid/tricuspid oral dentition with more unicuspid teeth as they grow. In the fossil assemblages, we found systematic and recurrent changes in the relative abundance of these oral tooth types, with a higher proportion of unicuspids during low-stands. Moreover, O. hunteri living during low-stands systematically developed unicuspid dentition at a smaller body size, compared with conspeci ﬁ cs living during high-stands. Considering that low-stand conditions created a sizable area of oxygenated soft-bottom habitat that is lacking during high-stands, we propose that the associated change in oral dentition re ﬂ ects the ﬁ shes ’ exploitation of this new habitat either for feeding or for reproduction. The recurrent nature of the observed shifts provides evidence for the ability of O. hunteri to systematically adapt to local habitat change, and strongly suggests that morphological change in oral dentition promoted its long-term population persistence in an aquatic ecosystem presenting the dual challenge of being both sensitive to climate change (creating constantly shifting selective pressures) and isolated from surrounding populations (limiting gene ﬂ ow). In Lake Chala, ancient climate-driven lake-level ﬂ uctuations did not directly promote speciation in O. hunteri , but the population bottlenecks they may have caused potentially contributed to its divergence from sister species Oreochromis jipe .


INTRODUCTION
The ability of species to adapt their ecology and behavior to changing environments is a driving force of speciation (Thompson et al. 2018). Especially when environmental change gives rise to new, unexploited habitats (i.e., when it creates ecological opportunity), rapid adaptation provides ample opportunities for populations to diverge into separate species (Schluter 2000). Such ecological speciation is prominent in the cichlid fishes (Cichlidae) of the African Great Lakes, where it has rapidly produced highly diverse species flocks (Fryer and Iles 1972). In these lakes, lake-level fluctuations induced by Quaternary climate change have caused phases of extreme contraction (Lake Malawi), fragmentation (Lake Tanganyika), or even complete desiccation (Lake Victoria; reviewed by Salzburger et al. 2014). Alternatingly creating and eliminating certain types of freshwater habitat, these fluctuations have been invoked as an important driver of evolutionary radiation in African cichlids (Rossiter 1995, Danley et al. 2012, Ivory et al. 2016. Intrinsic biological factors also influence whether cichlids radiate, as colonization of new lakes and habitats does not automatically result in diversification and radiation (Wagner et al. 2012). One group of cichlids whose members rarely radiate is the paraphyletic species group formerly referred to as tilapia (Brawand et al. 2014). Compared to the often stenotopic and haplochromine cichlids, the tilapiine cichlids, including the genus Oreochromis, display a less specialized but highly versatile morphology, making them more tolerant to habitat change and less prone to extinction (Ribbink 1990). One trait that has enabled cichlids in general to cope with new environments is the large flexibility of their trophic apparatus, which consists of both oral and pharyngeal jaws (Muschick et al. 2011, van Rijssel et al. 2015. Oreochromines are generally herbivorous mouth-brooders and hence have a trophic morphology reflecting trade-offs between multiple functions (Trewavas 1983). Nevertheless, subtle differences in jaw and tooth morphology still permit them to specialize in certain ecological niches, especially during periods of food scarcity.
Morphological variation in dentition among Oreochromis populations has been directly linked to differences in habitat (Ndiwa et al. 2016), but evidence of morphological shifts within a single population in response to local habitat change is rare, especially over the substantial periods typically required for speciation. Study of fossils has already proven informative in this respect for other teleost fishes (Purnell et al. 2007, Bellwood et al. 2014), but cichlids are plagued by a scarcity of continuous fossil deposits, the typically disarticulated nature of their remains, and difficulty to assign individual fossils to specific lineages. These problems have hampered detailed study of long-term patterns of morphological change in a coupled ecological-evolutionary context (Murray 2001). Cichlid teeth and scales preserved in the sediments of existing African lakes hold the best hope for a continuous record (Cohen 2012). Although such records are now tentatively being explored (Reinthal et al. 2011, Muschick et al. 2018, diverse methodological and confounding factors have thus far prevented detailed study of environment-phenotype relations in individual cichlid species lineages (Dieleman et al. 2015).
One African lake of particular interest in this context is Lake Chala, an isolated crater lake on the border between Kenya and Tanzania near Mt. Kilimanjaro (Fig. 1b,c). This small (4.2 km 2 ) but~90-m deep lake with permanently stratified lower water column (Buckles et al. 2014, Wolff et al. 2014) has a finely laminated sediment record indicating long-term continuity of lacustrine conditions (Verschuren et al. 2009, Moernaut et al. 2010, and it contains abundant fossil teeth, bones, and scales of cichlid fishes throughout the last 25,000 yr (Fig. 1a). Currently, Lake Chala harbors the endemic Oreochromis hunteri (Chala tilapia) as only indigenous fish species, but following a series of introductions in the 1980s, this species now co-occurs with three nonnative cichlids, namely another but less welldefined Oreochromis species originally attributed to O. korogwe; the tilapiine Coptodon rendalli; and a much smaller haplochromine species of the Astatotilapia bloyeti complex (Dadzie et al. 1988, Seegers et al. 2003, Dieleman et al. 2019, Moser et al. 2018. Lake Chala is the focus of a multifaceted research program aiming to reconstruct the longterm environmental and climate history of equatorial East Africa. These studies revealed that the lake experienced pronounced fluctuations during the last 25,000 yr (Moernaut et al. 2010;Fig. 1a), creating a succession of high-stand and lowstand conditions sufficiently contrasting to have had an impact on the diversity and nature of aquatic habitats available to O. hunteri (Fig. 1d). Assuming that the indigenous tilapiine would over time have adapted to these habitat changes  (bones, scales) in the sediment record of Lake Chala, in relation to the reconstructed sequence of past lake-level changes during the last 25,000 yr (Moernaut et al. 2010). The horizontal gray bars delineate the six time intervals studied here, representing four ancient high-stand episodes and two ancient low-stand episodes. (b) Map of equatorial East Africa with location of the principal great lakes Victoria, Tanganyika, and Malawi, and with Lake Chala indicated as a star. (c) Bathymetry of Lake Chala as recorded in March 2003 with indication of the coring site (star), depth contours at 10-meter intervals (except the deepest contour at 94 m), and the crater-rim catchment boundary (dashed line). The dark gray shading shows the approximate extent of lake bottom area that remained anoxic throughout the past 25,000 yr. The light gray shading highlights the area of gently sloping lake bottom at the base of the steep underwater rock face nearshore that was likely at least seasonally oxygenated during the ancient low-stands and hence provided additional habitat for Oreochromis hunteri at those times. (d) Schematic representation of lacustrine habitat available to O. hunteri during high-stand (left) and low-stand (right) episodes, the former situation also characterizing Lake Chala today. The water column is divided into the approximate to ensure optimal exploitation of available food and reproductive resources, we hypothesized that this process would have induced morphological changes in the oral dentition of O. hunteri, which might be traceable through time by morphometric analysis of fossil teeth preserved in the sediment record. Lake Chala offers a uniquely controlled setting for such a study, because of its isolation from other freshwater systems, its continuous sedimentary archive containing abundant fish fossils, its well-studied environmental past, and the known recent timing of the arrival of other fish species. However, despite being the type species of its genus, currently only very general morphological and ecological descriptions of O. hunteri are available (Trewavas 1983). Reliable interpretation of its fossil record hence necessitated the construction of a quantitative reference frame for observed variation in oral tooth morphology. Thus, we first quantify the intraspecific variability in oral tooth morphology among O. hunteri in the modernday population, followed by analyzing variability in the morphology of fossil oral teeth deposited during six periods with contrasting habitat conditions over the last 25,000 yr. Finally, we assess whether variation in tooth shape observed through time can be linked to independently documented environmental changes, and whether the observed change can be attributed to the sorting of genotypes under natural selection or to phenotypic plasticity.

Study system
The water level of modern-day Lake Chala is high: Steep, rocky shorelines all around quickly drop down to~55-60 m depth, from where soft sediments slope more gently toward a maximum depth of~90 m in the middle of the lake (Moernaut et al. 2010; Fig. 1d). Monthly monitoring of dissolved oxygen throughout the water column indicates that during the long dry season of July-September the well-oxygenated zone extends to a depth of~45 m (Wolff et al. 2011, Buckles et al. 2014. For most of the year, only the uppermost 20-30 m is oxygenated, thereby restricting fish habitat to the steep-sloping rocky crater walls and the open-water (limnetic) zone. The rocks are covered by epilithic algae and interspersed with small sandy patches, providing a relatively continuous food source throughout the year. By contrast, limnetic productivity is largely restricted to a phytoplankton bloom in the dry season (Buckles et al. 2014, van Bree et al. 2018; during the rest of the year, low phytoplankton productivity implies poor open-water feeding conditions. The depositional center of Lake Chala has been continuously anoxic throughout the last 25,000 yr (Verschuren et al. 2009), allowing for excellent preservation of organic fossils (Meyer et al. 2018). However, the major lake-level declines (up to~40 m, Moernaut et al. 2010;Fig. 1a) which occurred during prolonged periods of regionally drier climate (Verschuren et al. 2009) must have created at least seasonally oxygenated conditions in soft-sediment areas of lake bottom at the base of the underwater rock face, which are currently too deep to be supplied with oxygen ( Fig. 1d). Considering the low primary productivity of Lake Chala overall, availability of this additional benthic habitat during ancient low-stands is likely to have affected the ecology of O. hunteri.

Modern-day study material
To study the extant population of O. hunteri, we purchased 88 specimens from fishermen, fresh-caught at monthly intervals between January 2014 and September 2015. Fish were photographed in lateral view, and their standard length (SL; distance in cm between the tip of the snout and the posterior end of the last vertebra) depth zones with year-round oxygenation (0-20 m depth, light blue), seasonal oxygenation (20-45 m depth, light-to-dark blue gradient), and permanent anoxia (>45 m depth, dark blue). During high-stand episodes, O. hunteri is restricted to rocky benthic habitat nearshore and the limnetic habitat offshore; food resources associated with soft-bottom benthic habitat are accessible only during low-stand episodes. Fish and substrate drawings by Tim Goethals, combined with a seismic-reflection profile (cf. Moernaut et al. 2010) showing a representative cross section of Lake Chala bathymetry.
( Fig. 1. continued) ❖ www.esajournals.org was measured to generate a body size-frequency distribution of the population (Fig. 2c). The heads of 12 specimens with SL values covering the full range of non-juvenile body sizes (here ≥16 cm, based on fisherman's catches) were frozen and transported to Belgium. Obtaining insight into variation of tooth morphology across the largest possible range of body sizes is important because individual cichlid fish develop multiple sets of teeth during their lives (polyphyodonty; Fryer and Iles 1972), all of which may end up in the fossil record. Upon arrival in the laboratory, the fish were preserved in 10% formaldehyde for later tooth extraction. For this purpose, the oral jaws (left dentary, DTL; left premaxilla, PML) were dissected, cleared, and stained using a protocol adapted from Taylor and Van Dyke (1985); for details, see Dieleman et al. (2015). Oreochromis species typically possess many oral teeth (between~350 and~950 in O. hunteri), implanted in the jaws in multiple rows (Trewavas 1983). O. hunteri generally has bicuspid teeth in the front row (FR; the outer series in official terminology; Barel et al. 1975) and tricuspid teeth in the second row (SR) and inner row (IR) of the jaw, together the so-called inner series. This large morphological disparity among the oral teeth of a single fish implies that establishing differences in average tooth morphology between specimens or populations (either in space or in time) requires measuring large numbers of teeth (Dieleman et al. 2015). For morphometric analysis of the modern-day population, we extracted 60 teeth from each specimen. This number was determined by iteratively reducing the sample included in morphometric analysis of two test specimens (Dieleman et al. 2015) to find the lowest number that did not affect the total combined level of morphological variation represented in the dataset. These 60 teeth were comprised of 10, 5, and 15 teeth extracted randomly from respectively the FR, SR, and IR of both the DTL and the PML of each analyzed specimen.
Teeth were photographed at 609 magnification under a binocular microscope. The curvature of longer teeth hampered standardized orientation of tooth crowns. Therefore, the crowns were clipped off and photographed again in a standardized manner. To optimize digitization, the contrast of photographs was enhanced using Photoshop CS6, and PML teeth were mirrored for direct comparison with DTL teeth. The edited photographs were bundled per jaw bone and per specimen into .tps files using tpsUtil version 1.58 (Rohlf 2010a).

Fossil study material
In 2003 and 2005, the CHALLACEA project recovered a 20.82 m long composite sequence of finely laminated lacustrine muds from the center of Lake Chala, deposited continuously during the past 25,000 yr (Verschuren et al. 2009). For this study, we first assessed the temporal distribution of fish fossils (teeth, scales, and bones) over this period by determining their abundance in 1026 contiguous 2-cm intervals throughout the core sequence, starting at 30 cm depth. For each interval, 2 g of untreated wet mud was rinsed on a 250-lm mesh sieve, and the retained residue was scanned under a binocular microscope at 25-409 magnification. These counts were visualized in R v.3.3.2 (R Development Core Team 2016) using the package rioja v.0.9.15 (Juggins 2017). Next, we used the lake-level reconstruction based on seismic stratigraphy (Moernaut et al. 2010) to select six time intervals of betweeñ 500 and~1500 yr duration, from alternating lake high-stands and low-stands and which promised sufficient fossil teeth for comprehensive morphometric analyses (Fig. 1a). Given the considerable volume of sediment required for extracting sufficient fossil teeth, the exact position and duration of selected intervals also depended on availability of adequate core material (Table 1). Therefore, samples for this study were not extracted from the already partly depleted master core sequence, but from overlapping duplicate core sections. The continuous lamination of Lake Chala sediments permits unambiguous cross-correlation of these duplicate sections to the master sequence.
Extraction of teeth was done by stirring sediments for 30 min in H 2 O, after which the solution was sieved into two size fractions (150-250 lm and >250 lm) that were both searched for fossils under a binocular microscope. Only oral teeth were retrieved, representing 27% of all fossil teeth encountered (Table 1). Pharyngeal teeth were counted but not analyzed further. The collected teeth were photographed at 609 magnification, clipped, re-photographed, and treated following the same procedure as adopted for modern-day teeth. Photographs of the strongly asymmetric bicuspid teeth were mirrored for direct comparisons, but the quasi-symmetric tricuspid and unicuspid teeth were not. Fossil teeth were grouped in .tps files per time interval.

Analysis of tooth size and shape
Data on tooth size and shape were obtained following methods described in Dieleman et al. (2015). Tooth length and narrowest width in micrometer were measured using ImageJ v. 1.49 m (Rasband 1997; Fig. 2a). Measurements were log-transformed (base 10) and further analyzed in R. The neck of fossil teeth was often broken, preventing measurement of total length, although their tops were usually sufficiently intact for inclusion in the morphometric dataset. Therefore, we examined to what extent enameloid (or more precisely, acrodin) width (EW, the Euclidean distance between landmarks 1 and 3 in Fig. 2a) can serve as a proxy for tooth length (TL) and body size (SL). Specifically, we correlated these measures from modern-day specimens for each tooth type separately and tested Table 1. Overview of the six fossil teeth assemblages (indicated with Roman numerals I-VI) extracted from sediment samples that represent alternating lake high-stand (H) and low-stand (L) periods during the last 25,000 yr, with indication of their composite core depth, sample wet weight, and the numbers of oral teeth that were respectively recovered and included in the morphometric dataset. the least-square correlations with Spearman's rank tests. Finally, we compared the median EW of fossil teeth, per tooth type, from high-stands and low-stands using the two-sample Wilcoxon tests. Non-parametric tests were used throughout, due to violation of the assumptions for parametric tests. The shape of tooth crowns was quantified with semi-landmark analysis including three landmarks (LMs) and two curves containing 20 semilandmarks each, anchored to the LM (Dieleman et al. 2015;Fig. 2a). Tooth shape was digitized in tpsDig2 v.2.17 (Rohlf 2010b) and imported in IMP6 (Sheets 2008) for Procrustes superimposition and perpendicular projection of semi-landmarks. Teeth for which one or more landmarks could not be assigned unambiguously were excluded, resulting in a dataset of 671 modernday teeth and 886 fossil teeth, which was further analyzed in R. These morphometric data were examined with non-metric multidimensional scaling (NMDS) to reconstruct patterns of morphospace occupation (i.e., the distribution of morphological trait variation in ordination space) in two dimensions, using functions in the packages vegan v.2.4.3 (Oksanen et al. 2017) and MASS v.7.3.45 (Venables and Ripley 2002). The resulting stress value (8.9) is low according to the criteria of Kruskal (1964) and Clarke (1993), indicating no risk for false interpretations.
As mentioned, tooth-shape variation within the modern-day O. hunteri population was analyzed separately to develop a comparative framework for morphological variation in the fossil teeth. Although the oral teeth of cichlids are generally described as unicuspid, bicuspid, or tricuspid (Trewavas 1983;Fig. 2b), the assignment of individual teeth to one of these tooth types is often somewhat ambiguous. Hence, we examined the robustness of our results under two quantitative classification methods, one based on the angles between cusps and one based on the clustering of tooth shape in morphospace using a Gaussian mixture model. As the two methods yielded highly similar results (Appendix S1: Figs. S1, S2), we here only present results based on inter-cusp angles. In this method, the number of cusps is determined by counting the number of crown incisions with an angle <135°. Tricuspids possess two such incisions, bicuspids one, and unicuspids none.
Considering the unavoidable difference in sampling methods, our collections of modern-day and fossil teeth cannot be considered fully equivalent, precluding detailed quantitative comparisons. However, as both collections derive from the same species (Dieleman et al. 2015(Dieleman et al. , 2019, a general comparison of the total morphospace occupied by the collections of fossil and modern-day teeth, and their respective overall average morphology, seemed appropriate. Subsequently, we statistically compared the morphospace occupation of fossil teeth from the two low-stand and four high-stand episodes, first by pooling tooth assemblages per type of lake phase and then pairwise between each of the six fossil assemblages, with Bonferroni correction. Assumptions for multivariate parametric tests, that is, multivariate normality and homogeneity of the variance-covariance, were not met; consequently, we used permutation tests as implemented in the package Hotelling v.1.0.4 (Curran 2017). Using the intercusp angle method, teeth were then assigned to one of the three major tooth types, allowing examination of the relative abundances of these tooth types in the six fossil assemblages. The significance of differences in tooth-type abundances between low-stands and high-stands was tested by comparing the observed abundances with those obtained for 10,000 randomly sampled datasets of the same sample size.

Inferring the mechanism of morphological change
Observed differences in tooth morphology between low-stand and high-stand episodes may be caused by genotype sorting driven by natural selection or by phenotypic plasticity. To differentiate between these potential drivers, we examined the variation in morphospace occupation among fossil assemblages, based on the expectation that selection intensity would be different between low-stands and high-stands, and hence produce a systematic difference in the amount of morphological variation. Alternatively, if observed morphological change is mainly driven by phenotypic plasticity, we expected to find similar levels of morphological variation during lowstands and high-stands. For this analysis, we pooled all specimens in a high-stand and a lowstand group, and compared the morphospace occupation between both lake states with the betadisper function of vegan. The multivariate homogeneity of variances between both was then tested with an ANOVA on the Euclidean distances of group members to the group centroid after principal coordinate analysis.

Modern-day reference framework
The SL of the 88 O. hunteri specimens collected in Lake Chala ranged from 16 to 30 cm, with 81% of all specimens measuring between 23 and 29 cm (Fig. 2c). Modern-day teeth become both longer and wider as body size increases, and because the increase in width is proportionally larger, they become more robust (Fig. 3). Enameloid width of all three tooth types correlates strongly with tooth length (unicuspids, r = 0.85, P < 0.001; bicuspids, r = 0.82, P < 0.001; tricuspids, r = 0.90, P < 0.001; Appendix S1: Fig. S3a) and with body size (unicuspids, r = 0.55, P < 0.001; bicuspids, r = 0.71, P < 0.001; tricuspids, r = 0.74, P < 0.001; Appendix S1: Fig. S3b). Therefore, we can use EW in fossils as a proxy of tooth length, and the EW distributions of fossil fish assemblages as an estimate of median body size.
The morphospace occupation of all 671 modern-day teeth is plotted in Fig. 4a, with a gradient of symbol color reflecting body size of the source specimens. Non-metric multidimensional scaling axis 1 discriminates between unior bicuspids (negative values) and tricuspids (positive values; Fig. 4b), whereas NMDS axis 2 reflects the distinction between unicuspids (negative values) and bi-or tricuspids (positive values; Fig. 4c). Tooth-shape changes related to growth are evident along both NMDS axes: The predominantly bi-and tricuspid teeth of small specimens (green symbols) are gradually replaced by unicuspid teeth in larger specimens (red symbols). The same morphospace with teeth classified by tooth type based on intercusp angles is visualized in Fig. 5a. Despite some overlap, teeth assigned to each of the three tooth types mostly occupy distinct locations in morphospace. The FR consists predominantly of bicuspid teeth, whereas the SR and IR consist mostly of tricuspid teeth (Fig. 5b), especially in small individuals. As fish grow, their dentition becomes increasingly unicuspid (Fig. 5c).

Fossil teeth assemblages
Enameloid width distributions show that the recovered fossil oral teeth (n = 886) are significantly smaller, on average, than our reference collection of modern-day O. hunteri teeth ( Fig. 6a; W = 114,790, P < 0.001). Grouping fossil teeth assemblages according to lake level reveals that teeth deposited during a high-stand episode (subtotal n = 360) are significantly smaller, on average, than those deposited during a lowstand (subtotal n = 526; W = 77, 434, P < 0.001). However, in all three tooth types separately the individual assemblages display a marked decline in median tooth size from older (VI) to younger (I) age since deposition (Appendix S1: Fig. S4), indicating that the former result may be biased by the abundant material from assemblage III (Table 1). Also, although most samples are skewed toward larger teeth, this skewness is enhanced in older samples (Appendix S1: Fig. S4), suggesting that post-depositional corrosion results in a disproportionate loss of smaller teeth over long time periods.
The morphospace occupation of fossil teeth, as represented by six fossil assemblages covering the past 25,000 yr, overlaps almost completely with that of the modern-day population of O. hunteri (Fig. 6b). When the fossil assemblages are grouped according to lake level, a significant difference is observed between the average shape of oral teeth from O. hunteri living during high-stands and those living during low-stands (T 2 = 85.25, P < 0.001). Moreover, when comparing the two low-stand and four high-stand assemblages separately, this shift in tooth morphology is shown to be recurrent through time (Fig. 7a). The difference is mainly situated along NMDS2, with on average more negative values for "lowstand" teeth. Since the variation in fossil tooth  shape explained by the two NMDS axes is identical to that in modern-day teeth (Fig. 4b, c), on average more negative NMDS2 values similarly indicate a greater proportion of unicuspid teeth in ancient O. hunteri which lived during lowstands. As might be expected, the average tooth morphology of modern-day O. hunteri, which experience high-stand conditions (Fig. 1a), plotted closest to the centroids of fossil teeth assemblages deposited during past high-stands (Fig. 7a), of which the most recent (from episode I) are between~500 and~1000 yr old (Table 1). Multivariate permutation tests yield significant differences in five out of eight Bonferroni-corrected pairwise comparisons of the individual high-and low-stand assemblages (Table 2). Classifying fossil teeth based on inter-cusp angles reveals a somewhat larger overlap between tooth types in fossils (Fig. 7b) than in modern-day teeth (Fig. 5a). Consistent with our observations in NMDS morphospace (Fig. 7a), the proportion of unicuspid teeth is found to be markedly higher during the two ancient low-stand episodes than during high-stand episodes (Fig. 7c), with the average increase in unicuspid teeth during low-stands (~10%) being significantly larger (P < 0.01) than expected by chance. Finally, the variance in morphospace occupation of fossil teeth from ancient high-stands is significantly larger than that of those from low-stands (F = 5.7742, DF = 1; 888, P = 0.016; Fig. 8), indicating that overall morphological variation is reduced during low-stands.

DISCUSSION
The oral dentition of modern-day Oreochromis hunteri Our data show that during their lifetime, modern-day O. hunteri in Lake Chala undergo significant changes in oral dentition. First, the teeth become gradually wider and more robust with increasing body size (Fig. 3). Secondly, larger fishes develop a greater proportion of unicuspid teeth (Figs. 4,5). Specifically, small fishes have bicuspid FR teeth followed by tricuspid SR and IR teeth, whereas very large fishes possess an almost entirely unicuspid oral dentition. Individuals of intermediate body size have an intermediate dentition: They largely retain the bicuspid/ tricuspid arrangement, but the major/middle tooth cusp becomes more dominant relative to the side cusp(s). The change toward a predominantly unicuspid dentition in very large specimens is thus gradual rather than abrupt (Fig. 4).
This observed change to unicuspid FR teeth and, eventually, an (almost) entirely unicuspid oral dentition is known to occur in large individuals of several other Oreochromis species, mainly in sexually mature males (Trewavas 1983). Most cichlids display dietary shifts as they grow (Njiru et al. 2008), and these are often accompanied by changes in different aspects of their trophic morphology (Streelman et al. 2007). This study is, to our knowledge, the first to document this ontogenetic change in cichlid oral dentition with geometric morphometrics, and it reveals that the changes occur gradually, at least in O. hunteri. In other cichlid taxa, a unicuspid oral dentition is often linked to feeding at a higher trophic level (Fryer and Iles 1972), but oreochromines are rarely entirely piscivorous or insectivorous. Most Oreochromis species, including modern-day O. hunteri in Lake Chala, are herbivorous (Trewavas 1983). We therefore concur with Trewavas (1983) that the shift to unicuspid dentition in O. hunteri is likely associated instead with sexual maturation and enhanced territorial behavior linked to breeding.

Trends in fossil teeth abundance and preservation
The presence of fish remains throughout the studied sediment record (Fig. 1a) demonstrates that fish have been continuously present in Lake Chala over at least the past 25,000 yr. Lack of a systematic contrast in the total abundance of fish fossils between high-stand and low-stand episodes further suggests that lake-level fluctuations had a negligible influence on post-mortem transport of fish remains to the mid-lake coring site. Nevertheless, fossil abundances vary considerably over the last 25,000 yr, showing a distinct increase in abundance~10,500 yr ago, that is, shortly after the Pleistocene-Holocene transition.  Fig. 7. (a) Morphospace occupation plot of the shape variation in oral tooth crowns among the six fossil teeth assemblages from Lake Chala, as reconstructed by non-metric multidimensional scaling (NMDS) and with separate colors for each of the four ancient high-stand episodes (green triangles) and two ancient low-stand episodes (yellow and ochre triangles). Also indicated are their respective centroid scores (filled stars), in relation to the centroid score of tooth-shape variation in the modern-day Oreochromis hunteri population (red-filled star). (b) The same morphospace occupation plot of all fossil teeth as in (a) but with colors representing the three principal tooth types. (c) Proportions of the three tooth types in each of the six fossil assemblages from alternating highstand (I, II, IV, and VI) and low-stand (III and V) episodes. Low-stand assemblages (gray bars) display a significantly larger proportion of unicuspid teeth.
( Fig. 7. continued) Along this same transition, a marked increase in the organic carbon content of the sediments is observed (Blaauw et al. 2011), presumably reflecting increased aquatic primary production (Meyer et al. 2018) associated with regional climate warming around that time (Loomis et al. 2017). This enhanced productivity may in turn have increased the lake's carrying capacity to sustain a larger fish population. No evidence of a pronounced demographic expansion within the last 25,000 yr was found in genetic analyses of modern-day O. hunteri (Dieleman et al. 2019); however, the sample size of that study may have been too small to pick up such a signature. In principle, the observed trend in fossil abundance down-core could also partly be due to taphonomic effects. Chemical corrosion by organic acids in sediments makes teeth progressively more brittle over time, which may not only decrease the overall abundance of teeth in older samples, but also progressively reduce the proportion of smaller teeth, as they are likely to corrode faster. However, the relatively modest changes in mean tooth size with age across the six assemblages (Appendix S1: Fig. S4) compared to the abrupt, order-of-magnitude increase in the abundance of fish fossils~10,500 BP (Fig. 1a) indicate that the latter almost certainly reflects a true increase in the O. hunteri population at that time.
Our data on EW indicate that fossil teeth are somewhat smaller, on average, than those of modern-day fishes (Fig. 6a). The main reason for this difference is probably that the smallest specimens in our modern-day survey are 16 cm, due to the prescribed large mesh size of fishing nets and the minimum size of fish deemed suitable for selling, whereas teeth of smaller fish also contribute to the fossil record. On the other hand, these smaller teeth may be more susceptible to corrosion after long-term burial, and/or our 150lm mesh sieve may not have retained all of those present in the sediment. Also, cichlids are polyphyodont (Fryer and Iles 1972), and each tooth has an average life span of around 100 d (Tuisku and Hildebrand 1994). Assuming that the life span of O. hunteri is similar to that of other Oreochromis species (up to~10 yr; Froese and Pauly 2017), each large specimen contributes up to three dozen sets of teeth to the sediment record during its life, further increasing the proportion of small teeth. Thus, the smaller median size of fossil teeth conforms to expectation. In the end, whereas the minimum tooth size in our modern-day reference collection is mainly controlled by the mesh size of fish nets, that in fossil assemblages is controlled by the sieves used for sediment processing. Modern-day and fossil assemblages are hence not exactly equivalent as to tooth size; however, their morphospace occupation is nevertheless highly similar (Figs. 6b, 7a). Fig. 8. Principal coordinate analysis plot showing that overall morphological variation is significantly larger among fossil teeth deposited during ancient high-stand episodes (black circles; n = 363) than among fossil teeth deposited during low-stand episodes (red triangles; n = 527), as measured by the Euclidean distances of individual teeth to their respective group centroid (solid black and red circles). Notes: Direct pairwise comparison between the morphology of fossil and modern-day tooth assemblages is not appropriate, because the two datasets are not fully equivalent (see Methods and Discussion). * Significance at the 5% level.

Temporal trends in oral tooth shape
The near-complete congruence between fossil and modern-day oral teeth in morphospace (Fig. 6b) confirms pre-1980s historical accounts (G€ unther 1889(G€ unther , Dadzie et al. 1988) and genetic evidence (Dieleman et al. 2019) that O. hunteri has been the only fish species inhabiting Lake Chala before the recent introductions, and this for at least the last 25,000 yr. The overall variability in tooth shape in the living population has thus been remarkably stable over time, implying that our examination of modern-day O. hunteri provides a robust framework for interpretation of the fossil assemblages representing ancient O. hunteri populations. The similarity of modern-day teeth with those of ancient O. hunteri living under comparable high-stand conditions as recently as~500-1000 yr ago (Fig. 7a) also implies that any change in oral tooth morphology that may have occurred since the 1980s, due to several decades of co-habitation and, hence, possible interference competition with the introduced fish species (character displacement; Brown and Wilson 1956), is modest compared with the mean difference in tooth morphology between ancient high-stand and low-stand populations.
Focusing on the fossil tooth assemblages, our data revealed significant and recurrent changes in the morphospace occupation of oral teeth from fish populations living at different times in the past (Fig. 7), synchronous with the independently documented lake-level changes (Moernaut et al. 2010;Fig. 1a). Specifically, during past episodes of low lake level, an average increase of 10% in the relative abundance of unicuspid teeth is observed (Fig. 7c). The differences in oral dentition between past high-stand and low-stand populations of O. hunteri in Lake Chala are overall similar to the ontogenetic changes in oral dentition which occur in the modern-day population. However, averaged per tooth type the mean size of teeth from low-stand assemblages is not systematically larger than that of teeth from high-stand assemblages (Appendix S1: Fig. S4), arguing against a larger mean body size of low-stand O. hunteri populations. Rather, our data suggest that during low-stand episodes O. hunteri shifted toward a unicuspid dentition earlier in life, likely triggered by a change in available habitat. For example, the unicuspid dentition could have been advantageous to prey on the invertebrate macrobenthos that likely developed during low-stands in the then oxygenated and accessible soft-bottom sediments at the foot of the steep underwater rock face (Fig. 1d). However, insectivorous oreochromines are rare in all of the diverse aquatic habitats they occupy, so the probability that O. hunteri would have acquired that habit when presented the opportunity may be rather low. Alternatively, an earlier shift to unicuspid dentition might indicate that low-stand fishes reached sexual maturity at a smaller body size (i.e., earlier age) than those living during high-stands. Population-level intraspecific variability in body size and age at maturity is common in Oreochromis species, and strongly linked to environmental variation (Duponchelle and Panfili 1998).
Besides marked trends in the proportion of unicuspid teeth potentially reflecting changes in lifehistory traits, we also documented significantly smaller overall variation in the oral tooth morphology of fishes living during low-stand periods compared to high-stand periods (Fig. 8). This finding suggests that natural selection either reduced the amount of standing variation in the oral dentition of low-stand populations, for example, by increased competition or a reduction of total population size, or alternatively that diversifying selection (as opposed to simply relaxing selection) may have occurred during high-stands. Our present data and earlier genetic analysis (Dieleman et al. 2019) indicate that such selection did not cause population divergence in O. hunteri, contrasting with the rapid phenotypic and genomic differentiation and diversification being reported in haplochromine cichlids (Malinsky et al. 2015, Moser et al. 2018. Apparently, the recurrent morphological changes documented here have enabled the endemic population of O. hunteri to cope with successive events of local environmental change without promoting speciation. In this and possibly other oreochromine cichlids, long-term population separation and the genetic drift associated with any population bottlenecks during this time may be more probable mechanisms for the formation of new species than adaptation to ecologically contrasting habitats. The operation of these mechanisms is exemplified in Lake Chala, where long-term isolation of the local fish population from neighboring populations of Oreochromis jipe (at least 25,000 yr but possibly >200,000 yr; Dieleman et al. 2019) eventually led to the evolution of endemic O. hunteri. Thus, even if a sequence of ancient lake-level fluctuations did not directly promote speciation in O. hunteri, they may still have been instrumental in creating genetic divergence from its more widespread sister species O. jipe.
Modification of the dental apparatus in response to climate-driven habitat change has frequently been hypothesized to be an important step in cichlid radiation (Kocher 2004, Danley et al. 2012, Ivory et al. 2016). However, evidence supporting this hypothesis is thus far mainly derived from phenotype-environment correlations that do not necessarily imply causality. In this study, we exploited the unique setting of Lake Chala's sediment record and a reference framework of tooth morphology in modern-day O. hunteri to link morphological trends in fossil teeth assemblages directly to the independently reconstructed evidence for long-term habitat change. Specifically, our analyses revealed systematic and recurrent changes in O. hunteri tooth morphology in response to a sequence of ancient lake-level fluctuations. Our data further suggest that selection played an important role in the observed morphological change, notwithstanding the substantial level of phenotypic plasticity that has been documented for the dental apparatus of oreochromine and other cichlids (Ribbink 1990, van Rijssel et al. 2015, Ndiwa et al. 2016). Currently, we cannot exclude the possibility that phenotypic plasticity played a role in providing the initial variation on which selection subsequently acted. In fact, the long-term persistence of O. hunteri in Lake Chala, and other oreochromines inhabiting fluctuating aquatic ecosystems (Ribbink 1990), may result from plasticity-driven adaptation, that is, the ability to plastically modify morphological features and subsequently refine them through selection. Such adaptation may play an important role in speciation (Thompson et al. 2018) and in evolutionary radiation (West-Eberhard 2005, Pfennig et al. 2010, Muschick et al. 2011, Laland et al. 2014; however, in Oreochromis species it is perhaps primarily a powerful mechanism to cope with recurrent environmental fluctuation. Our data, certainly, provide evidence for the ability of O. hunteri to rapidly and systematically adapt to local habitat change, and to repeatedly do so over very long periods of time. This adaptive response increased its chance of long-term population persistence in an aquatic ecosystem presenting the dual challenge of being both sensitive to climate change (producing constantly shifting selective pressures) and hydrographically isolated from neighboring populations (limiting gene flow). Indeed, this endemic fish population survived at least four major climate-driven environmental upheavals (potential population bottlenecks) within the currently documented part of lake history. In a broader context, the results of this study support the often invoked but rarely documented hypothesis that past climate-driven lake-level fluctuations have stimulated evolutionary radiation in the cichlid communities of East Africa's large rift lakes.

ACKNOWLEDGMENTS
This study was carried out under Memorandum of Understanding A14/TT/0923 between Ghent University and the National Museums of Kenya (NMK), and institutional affiliation of Dirk Verschuren with NMK. We thank Caxton Oluseno and the fishermen of Lake Chala for field assistance, and Maarten Van Daele for assistance in cross-correlating duplicate core sections to the dated master core sequence. Walter Salzburger, Jos Snoeks, George Turner, Ann Huysseune, Dries Bonte, and Dominique Adriaens provided many insightful comments which improved this paper. This work was financed by Ghent University Research Council through Collaborative Research Activity "Deep-CHALLA" (2013-2018; Jorunn Dieleman and Dirk Verschuren) and by ANR-JCJC-EVOLINK (Bert Van Bocxlaer). The lake-sediment cores used in this study were collected under permit 13/001/11C of the Kenyan Ministry of Education, Science and Technology, in the framework of the European Science Foundation's Euro-CLIMATE project CHALLACEA (2005)(2006)(2007)(2008).