Skip to main content

Ecomorphological divergence and habitat lability in the context of robust patterns of modularity in the cichlid feeding apparatus



Adaptive radiations are characterized by extreme and/or iterative phenotypic divergence; however, such variation does not accumulate evenly across an organism. Instead, it is often partitioned into sub-units, or modules, which can differentially respond to selection. While it is recognized that changing the pattern of modularity or the strength of covariation (integration) can influence the range or rate of morphological evolution, the relationship between shape variation and covariation remains unclear. For example, it is possible that rapid phenotypic change requires concomitant changes to the underlying covariance structure. Alternatively, repeated shifts between phenotypic states may be facilitated by a conserved covariance structure. Distinguishing between these scenarios will contribute to a better understanding of the factors that shape biodiversity. Here, we explore these questions using a diverse Lake Malawi cichlid species complex, Tropheops, that appears to partition habitat by depth.


We construct a phylogeny of Tropheops populations and use 3D geometric morphometrics to assess the shape of four bones involved in feeding (mandible, pharyngeal jaw, maxilla, pre-maxilla) in populations that inhabit deep versus shallow habitats. We next test numerous modularity hypotheses to understand whether fish at different depths are characterized by conserved or divergent patterns of modularity. We further examine rates of morphological evolution and disparity between habitats and among modules. Finally, we raise a single Tropheops species in environments mimicking deep or shallow habitats to discover whether plasticity can replicate the pattern of morphology, disparity, or modularity observed in natural populations.


Our data support the hypothesis that conserved patterns of modularity permit the evolution of divergent morphologies and may facilitate the repeated transitions between habitats. In addition, we find the lab-reared populations replicate many trends in the natural populations, which suggests that plasticity may be an important force in initiating depth transitions, priming the feeding apparatus for evolutionary change.


Characterizing the pattern and magnitude of covariation among traits has been a central theme of evolutionary biology for more than 200 years [1, 2]. However, it was not until Olson and Miller [3] that the understanding of trait covariation was formalized into a statistical framework. They suggested that both developmental and functional interactions result in the observed correlations between traits that they called morphological integration. Since Olson and Miller, much effort has gone into characterizing the patterns and strength of trait correlations at various scales – from between populations to across phyla (e.g., [4, 5]). These studies often characterize the types of correlation that exist using two interdependent terms: modularity and integration. Suites of traits that appear more correlated with each other than to other traits are termed modules. The number and identity of modules across a structure is generally referred to as the pattern of modularity, whereas the strength of correlation among traits within a module is termed integration. ‘Tinkering’ (sensu [6]) with both the pattern of modularity and the magnitude of integration provides a means to alter phenotypic variation in a way that may impact how a population can respond to selection [7, 8]. For example, an anatomical structure may be more able to respond (i.e., more evolvable) if the direction of selection aligns with the major axis of covariation [9]. In addition, by parsing an organism into discrete anatomical units, each module can become a separate target for natural selection. If regions of an organism can develop and evolve independently, then this could permit an increase in morphological diversity (i.e., disparity [10]), open up unique or unoccupied niches [11], and influence the rate of evolution [7, 12].

Empirical and simulation studies suggest three scenarios for potential relationships between modularity and shape (Fig. 1): 1) Differences in the pattern of modularity or magnitude of integration are associated with a corresponding change in shape. 2) A similar pattern of modularity or magnitude of integration is associated with differences in shape among populations. 3) Differences in the pattern of modularity or magnitude of integration are associated with no concomitant change in shape. Whereas scenario 1 is well supported by the literature [13,14,15,16], scenarios 2 and 3 appear to be a less common occurrence [17,18,19]. While all these scenarios are supported by empirical evidence, each scenario can send morphological evolution on a drastically different trajectory. As outlined above, these divergent trajectories may result in rate or disparity differences, and likely have implications for evolvability (the ability to generate adaptive phenotypic variation [20]), depending on how variation is partitioned among modules.

Fig. 1
figure 1

Three possible scenarios depicting the relationship between modularity and organismal morphology. Points on mandible schematics reflect example landmark coordinates, and arrows denote covariation between landmarks. Red colors reflect a pattern of covariation based on two mandibular modules, blue color reflects three mandibular modules. Numbers refer to descriptions noted in the main text

Cichlid fish from the East African Rift Valley provide an opportunity to examine these different scenarios of concordant or discordant changes in modularity and morphology. This lineage displays many convergent phenotypes that have evolved to exploit similar habitats and trophic niches in different lakes [17, 21]. It has been proposed that similarity in the pattern of modularity across lakes may have facilitated the iterative nature of cichlid evolution [22]. In fact, selection may act on the degree of covariation between traits, given modularity itself can evolve (reviewed by [23]), and there appears to be little overlap between loci that control shape of individual traits and those that control covariation between traits [18, 23,24,25]. Lake Malawi boasts the greatest taxonomic and morphological diversity of cichlids of any other African Rift Valley lake [26]. Taxa frequently partition their habitat by depth. Position along a depth gradient correlates with large differences in light, temperature, and oxygen that lead to differences in diet, predators, physiology, and sensory systems [27]. Previous studies have documented substantial morphological differences among such eco-types, especially with respect to the craniofacial skeleton [17, 22, 28]. However, there appears to be a general conservation in integration across broad taxonomic levels and feeding morphologies [29], and only minor differences in patterns of modularity [18, 30].

Ecomorphological divergence can occur via genetic mechanisms or via phenotypic plasticity. While the former may facilitate adaptation over many generations, the latter permits populations to respond to variation in environmental conditions within a generation by remodeling their morphology. Many studies have documented how populations exposed to different diets can adaptively remodel their trophic morphology to permit more efficient feeding behaviors [31,32,33]; however, the extent to which levels of integration or patterns of modularity can be influenced by plasticity remains an open question [but see [34] for some examples]. As Lake Malawi cichlids are known to partition habitat by depth, a feature that has led to broad differences in trophic morphology due to divergent feeding behaviors [17], phenotypic plasticity could be a means to facilitate depth transitions allowing adaptive morphological changes that could later become canalized. The mechanisms that would underlie these plastic morphological changes could center on changes to the pattern of modularity and/or the level of integration. If plasticity can facilitate a change in modularity and/or integration, this would change how variation is partitioned and accumulated in different regions of the feeding apparatus. Alternatively, if plasticity retains a particular pattern of modularity or level of integration, this may facilitate rapid morphological change along a ‘line of least resistance’ in the feeding apparatus of a population via genetic assimilation in order to adapt to a new habitat [35, 36].

Here we use the Lake Malawi cichlid genus, Tropheops, to more explicitly explore the relationship between morphology, modularity and evolution. Tropheops are one of the most rapidly evolving clades of cichlids, with speciation rates estimated to be as high as 1 per 1000 years [37]. Given this rapid diversification, we use amplified fragment length polymorphism (AFLP) to construct a phylogenetic tree of Tropheops populations to examine the frequency of transitions between depths and to more formally test hypotheses concerning morphological evolution. Furthermore, we focus on characterizing the feeding apparatus of Tropheops (mandible, lower pharyngeal jaw, maxilla, pre-maxilla), as members of this species complex occupy a spectrum of depths from shallow sediment-free conditions that require a more robust feeding apparatus to pluck attached filamentous algae from rocks, to deep sediment-rich habitats that require a more gracile feeding apparatus to sift through sediment on and between the rocks [28, 38]. First we compare the patterns of craniofacial modularity and within-module integration between Tropheops in shallow and deep habitats. We then examine rates of morphological evolution and disparity both between depths, and among modules of a given skeletal element. If Tropheops undergo multiple transitions between deep and shallow habitats, we expect these transitions to be accompanied by morphological change, due to differences in feeding behavior, and increases in the rate of morphological evolution in functionally important modules (i.e., muscle attachment sites). Similarly, we examine differences in disparity to assess whether specific depths facilitate the exploration of novel regions of morphospace, indicative of greater dietary range or more relaxed selection. Finally, we examine the role of plasticity in Tropheops morphological evolution by attempting to experimentally recapitulate the trends in morphology, modularity, and disparity observed in the natural populations.

We predict that Tropheops species exhibit a general conservation in the patterning of their craniofacial modules, and that this attribute is associated with the repeated evolution of ‘shallow’ and ‘deep’ ecomorphs. While we expect shallow and deep populations to exhibit similar rates of morphological evolution and disparity overall, we predict the fastest rates of morphological evolution should arise in those modules with the greatest functional importance. We also predict that we can mirror those patterns observed in natural populations via experimentally inducing plastic differences in morphology by raising Tropheops in conditions that mimic shallow or deep environments. We expect that plastic change in morphology will occur without concomitant changes in the pattern of modularity. Taken together, this would support that assertion that the cichlid craniofacial skeleton is characterized by robust patterns of modularity, despite differences in functional demands and morphology, and that this attribute may facilitate the rapid colonization of divergent habitats.


Tropheops phylogenetic tree

Our Bayesian tree constructed using AFLP markers exhibited monophyletic groupings for the Tropheops and Maylandia species complexes within the mbuna (Figure S1; Table S1). Nodal support for many clades was typically low (i.e., posterior probabilities < 0.75). Low nodal support is not unusual for Malawi cichlids given their rapid radiation, and high hybridization rates [39]. Notably, support for monophyly of Tropheops ‘species’ is rare, and in many cases species are conspicuously polyphyletic. Exceptions to this include species from isolated/island populations such as Chinyamwezi and Chinyamkwazi. Groupings by locality are also rare. Different species from the same locality (e.g., Mazinzi Reef, “MZ” in Figure S1) are widely distributed across the phylogeny. However, replicate individuals from the same locality and species almost always cluster together, indicating the robustness of the genetic data. There are only three Tropheops taxa that are not monophyletic: T. microstoma from Domwe Island, T. lilac from Thumbi West, and T. sp “zebra mumbo” from Mumbo Island. Individuals from these species/populations group within a larger clade that appears to be more reticulate.

We also examined a number of tree statistics output from MrBayes. The MCMC run produced a high effective sample size (ESS), > 500 for all parameters, indicating that the trace contained few correlated samples, and represented the posterior distribution well. We also report a potential scale reduction factor (PSRF) value close to 1 for our convergence diagnostic, indicating that we have a good sample from the posterior probability distribution. We then discarded 25% of the initial trees for burn-in and randomly sampled 1000 trees from the Bayesian posterior distribution (BPD). We used the BPD trees for comparative methods that assessed differences in rates, disparity, and modularity.

Tropheops morphology

When species were grouped based on their occupation of either shallow or deep habitats, bone shapes exhibited some overlap, but generally occupied distinct regions of principal component (PC) morphospace (Fig. 2; Table S2, S3). Previous work has shown that shallow water species typically possess short, stout jaws to forage on clean, filamentous strands of algae, whereas deep water species are generally characterized by longer jaws to comb and shift loose material from sediment-covered substrate [28]. Our morphometric analyses confirm and extend this trend. We discuss results from the mandible and lower pharyngeal jaw in the main text, maxilla and pre-maxilla results can be found in the supplementary information.

Fig. 2
figure 2

Morphospace occupation for natural Tropheops feeding bones from different depth regimes. Wireframe models reflect morphology at the extreme of a given axis. a, mandible PC1–3; b, lower pharyngeal jaw PC1–3. Red, individuals from shallow populations; black, individuals from deep populations


The first PC axis for the mandible characterized the height of the ascending arm and mandible length (20.20% of the variation). PC2 explained differences in the position of the ascending arm (anteriorly to posteriorly projected), the depth of the dentary, and the width of the mandible (10.99% of the variation). PC3 represented change in the lateral compression, or depth, of the mandible (9.61% of the variation). Mandible morphology differed based on occupation of deep or shallow habitats (pMANOVA; F = 16.6, P < 0.01), however much of this separation was confined to PC1 (Table S3).

Lower pharyngeal jaw

The first PC axis for the pharyngeal jaw reflected change in wing length, tooth-plate width, and jaw length (28.48% of the variation). PC2 explained differences in the depth of the jaw base and keel (10.32% of the variation). PC3 represented change in the wing height and the concave to convex curvature of the tooth plate (9.01% of the variation). Pharyngeal jaw morphology separated based on habitat depth (pMANOVA; F = 12.8, P < 0.01), with much of this separation confined to PC1 (Table S3).

Comparative methods

Our phylogenetic tree was constructed using multiple Tropheops ‘species’ from different localities and across a wide range of depth regimes to characterize transitions between deep and shallow habitats. We used stochastic character mapping (SIMMAP) on our 1000 BPD trees to quantitatively assess transition rates between deep and shallow habitats [40, 41]. We found 11.3 transitions between deep and shallow habitats on average, with more changes occurring in the direction of shallow to deep habitats (7.9 transitions), rather than deep to shallow habitats (3.4 transitions).

Modularity and integration

We next examined the pattern of modularity and the strength of integration within modules in species from deep or shallow habitats to understand whether differences in feeding morphology are linked to divergent covariation patterns. We used the software package EMMLi [42] to assess the fit of different modularity hypotheses (Fig. 3; Table S4) in a phylogenetic framework with our 1000 BPD trees and recorded the frequency at which each modularity hypothesis was supported from populations in each habitat. Given modularity analyses can be biased by differences in sample size between populations, we also tested the sensitivity of our modularity models to uneven sample size as our shallow population had almost double that of the deep population. Evidence for retention of a single pattern of modularity between depths would suggest that modularity is robust to change despite the opposing functional demands associated with these different habitats. It would also provide another example of different morphologies evolving within the context of a conserved covariation structure (e.g., Fig. 1 [17, 18];). On the other hand, evidence for differences in the pattern of modularity between habitats would suggest that breaking of the covariation structure was involved in the evolution of these divergent eco-morphologies.

Fig. 3
figure 3

Partitioning schematics for competing modularity hypotheses. Colors reflect module partitions to be assessed by EMMLi. Letters correspond to the partitioning scheme (Table S4). The best fitting modularity hypothesis is illustrated by *. a, Mandible; b, lower pharyngeal jaw

We found that the pattern of modularity was largely consistent between Tropheops from deep or shallow environments (Table 1; Table S5). In both groups, the best-fitting partitions were based on functional units including the attachment sites for muscles or ligaments, tooth bearing regions, and joints. This indicates that patterns of modularity within craniofacial bones (1) arise due to functional demands, (2) are robust to differences in habitat and foraging mode, and (3) are not associated with morphological divergence. While we found the pattern of modularity is similar between habitats, we found that the level of integration within modules differs between habitats, which suggests that eco-morphological divergence is associated with changes in the strength of covariation within modules.

Table 1 EMMLi output highlighting best supported modularity model(s) for the mandible and pharyngeal jaws in both natural and experimental populations. Additional competing models are presented if within two AICc units of the best supported model


For both shallow and deep datasets, a six module pattern was supported at the highest frequency for shallow individuals and ~ 30% of runs for deep individuals (Table 1, S5; Fig. 4a). This module-partitioning hypothesis breaks the mandible into a tooth-bearing region, lateral line canal, retro-articular process, quadrate-articular joint, ascending arm of the articular, and articular excurvation (Fig. 3a). This six-module pattern is highly similar to that reported by Parsons et al. [18] from a 2D landmark dataset in closely related mbuna species. In each habitat, EMMLi returned distinct values of integration (ρ) for both within- and among module comparisons for all six modules. Estimated values of ρ were highest for the lateral line canal and ascending arm modules (modules B and E respectively), lowest for the tooth-bearing module (module A), and intermediate for the quadrate-articular joint, articular excurvation, and retro-articular modules (modules C, D, and F respectively). Notably, Tropheops from deeper habitats exhibited relatively higher values of ρ (i.e., high levels of integration) within modules compared to species from shallow habitats (Fig. 5c; Table S6).

Fig. 4
figure 4

Support for competing modularity hypotheses between habitats as determined by EMMLi. a, mandible; b, lower pharyngeal jaw. The ‘Shallow Sampling’ plot reflects the frequency of module model selection derived from a sample size that matches the deep population. Colors used in the module model frequency plots are placed as background colors on module partition schematics that they represent

Fig. 5
figure 5

Violin plots depicting parameter values output from the rate of morphological evolution, disparity, and integration analyses across habitats for each module. Anatomical schematics illustrate the location of the module being tested in red. The depicted modules, and their associated letters, are based on the best-fitting modularity hypothesis determined by EMMLi. See Fig. 3 and Table S4 for full range of modularity hypotheses we tested. a-c, Mandible; d-f lower pharyngeal jaw. Red, shallow habitat; Black, deep habitat

Tropheops from the deep habitat exhibited some variation in which modularity hypothesis was best supported. The five- and six-module models exhibited AICc scores within four units of each other in the deep sample, indicating overall similarity in model fit (Table 1). As a result, we cannot conclusively say which modularity model best-fits our mandible data. The five-module hypothesis differs from the six-module hypothesis in that it unites the articular excurvation and lateral line canal modules. The overall pattern is therefore fairly similar between the two hypotheses, and the ability of EMMLi to detect a difference between the five- and six-module hypotheses may be hindered by a small sample size in the deep taxa. Indeed, the five-module hypothesis was selected in ~ 50% of the models in taxa from the deep habitat, but when we assessed the sensitivity of our data to low sample size by sub-sampling the shallow taxa such that they mimic the sample size of the deep taxa, we find additional competing modularity hypotheses have increased levels of support (Fig. 4a).

Lower pharyngeal jaw

We found the pharyngeal jaw best fit a five-module hypothesis for both shallow and deep Tropheops members (Table 1, S5; Fig. 4b). EMMLi found support for partitioning into modules that include the tooth plate, pharyngeal wings, and attachment sites for the pharyngoclithralis internus, pharyngoclithralis externus, and pharyngohyoideus muscles (Fig. 3b). While there was some evidence for a two-module hypothesis in shallow taxa, reflecting partitions that divide the anterior and posterior portions of the lower pharyngeal jaw, this pattern arose in less than 10% of cases. The data may be somewhat sensitive to sample size, as additional modularity hypotheses gain support when we sub-sample the shallow taxa (Fig. 4b). Despite a small sample size in the deep populations, a five-module hypothesis was selected in 100% of the models. In both habitats, estimated values of ρ were typically highest for the three modules defined by muscle attachment sites (modules C, D, and E), and lowest for the tooth plate and pharyngeal wing modules (modules A and B respectively). Similar to the mandible, we found that the relative differences among module ρ were generally similar between deep and shallow species, but that Tropheops from deeper habitats exhibited higher values of ρ (Fig. 5f; Table S6).

Morphological disparity and rates of evolution

Morphological disparity and rates of morphological evolution represent different measures of evolutionary potential (i.e., evolvability). Evidence for differences in disparity and/or rates between shallow- and deep-water habitats would suggest that evolvability of the feeding apparatus is influenced by foraging environment.

Morphological disparity

We detected no statistical difference in disparity for any feeding bones between Tropheops from shallow versus deep environments (Table S6). Similarly, we found no evidence for a difference in disparity between habitats when we compared subsets of landmarks defined by the best fitting modularity hypothesis suggested by EMMLi (Fig. 5b, e; Table S6). Thus, with respect to magnitudes of morphological variation, evolvability appears to be similar across shallow and deep foraging habitats. In spite of a lack of statistical difference, we note that for the mandible, disparity was consistently higher in the shallow population across modules (Fig. 5b). Distinct trends in the lower pharyngeal jaw were more difficult to determine, as disparity in shallow populations was higher in two of the five modules (Fig. 5e).

Rates of morphological evolution

Evolutionary rates were also compared across foraging habitats as well as between modules within bones (Table S6). For the former, we found no difference in rates of morphological evolution for craniofacial bones between species from shallow or deep environments. Further, no differences in rates were observed between depths when bones were partitioned into modules defined by EMMLi (Fig. 5a, d; Table S6). When the entire bone is considered, the results from both our disparity and rates analyses suggest that evolvability of the feeding apparatus is the same across foraging habitats. Again, despite the lack of statistical significance, we note deep populations typically exhibited faster rates of morphological evolution in both the mandible and lower pharyngeal jaw (Fig. 5a, d).

Alternatively, we observed statistically significant differences among several mandible and lower pharyngeal jaw modules (Table S7). In general, those modules representing bony processes or muscle attachment sites that have a direct association with feeding biomechanics evolved the most rapidly. For example, the ascending arm of the mandible (module E) is evolving more than 1.5x faster than the quadrate-articular joint (module C), (ascending arm, σ2 = 1.09 × 10− 5 (95% CI = 9.81 × 10− 6, 1.33 × 10− 5); quadrate-articular, σ2 = 6.94 × 10− 6 (95% CI = 6.30 × 10− 6, 8.55 × 10− 6); p < 0.01). Similarly, the pharyngeal jaw muscles attaching to the posterior surface were evolving almost three times as fast as muscles attaching to the anterior keel (pharyngoclithralis internus muscle attachment site (module C), σ2 = 2.08 × 10− 5 (95% CI = 1.88 × 10− 5, 2.56 × 10− 5); pharyngohyoideus muscle attachment site (module E), σ2 = 7.07 × 10− 6 (95% CI = 6.94 × 10− 6, 9.47 × 10− 6); p < 0.001). These data suggest that certain regions of the feeding apparatus exhibit greater evolutionary potential than others.

Experimental recapitulation of shallow-deep water environments

Deep- and shallow-water foraging environments broadly mimic the benthic-limnetic eco-morphological axis that characterizes multiple cichlid radiations [17, 43, 44]. Tropheops from shallow water habitats generally possess benthic feeding morphologies in order to more efficiently forage on attached filamentous algae, whereas species that live at depth tend to possess more gracile, limnetic morphologies in order to suck and sift loose material from sediment covered rocks [28, 30]. Since the benthic-limnetic foraging axis can be re-created in the lab, we sought to test whether patterns of morphological divergence, modularity, and disparity observed across natural populations of Tropheops could be replicated within a single species reared under alternate foraging environments in the lab. If lab-reared Tropheops mimic the divergence, modularity, or disparity results found in the evolutionary sample, this would be consistent with a role for phenotypic plasticity in response to alternate kinematic demands in influencing patterns of evolutionary change. Alternatively, if patterns in lab-reared Tropheops do not match those from natural populations, this would suggest a larger role for genetic divergence driving evolution in this lineage.

Mandible morphology

The first PC axis for the mandible characterized differences in the position of the ascending arm (anteriorly to posteriorly projected), the length of the coronoid process, and the length of the mandible (13.48% of the variation). PC2 explained differences in the depth of the mandible, and RA length (11.19% of the variation). PC3 represented change in the height of the ascending arm and the height of the dentigerous portion of the mandible (9.33% of the variation). Unlike the natural populations, mandible morphology did not differ based on benthic or limnetic treatments (MANOVA; F = 1.73, P = 0.16), and there was substantial overlap in morphospace (Fig. 6a; Table S3).

Fig. 6
figure 6

Morphospace occupation for experimental Tropheops feeding bones mimicking different depth regimes. Wireframe models reflect morphology at the extreme of a given axis. a, mandible PC1–3; b, pharyngeal jaw PC1–3. Red, benthic-shallow members; black, limnetic-deep members

Pharyngeal jaw morphology

The first PC axis for the pharyngeal jaw reflected change in the width of the tooth plate and depth of the keel (16.14% of the variation). PC2 explained differences in jaw length (14.24% of the variation). PC3 represented change in the depth of the jaw base (11.6% of the variation). Pharyngeal jaw morphology separated based on feeding treatment (MANOVA; F = 6.31, P < 0.01), with much of this separation confined to PC2 and PC3 (Fig. 6b; Table S3).

Modularity, integration, and disparity

While patterns of shape divergence between treatments were generally not the same as what was observed among natural populations (Table S8), patterns of modularity were highly similar. In particular, the experimental populations exhibited the same six-module hypothesis for the mandible and five-module hypothesis for the lower pharyngeal jaw as observed in the natural populations (Table 1, S9). We also found that a single pattern of modularity was favored regardless of whether the experimental population was subjected to an environment that mimicked a deep or shallow habitat. Thus, patterns of modularity within craniofacial bones appear to be highly robust among Tropheops.

We also found that the magnitude of within-module integration (ρ) differed between modules within a bone, and between modules across treatments (Table 2). Within the mandible, the greatest values of ρ were observed for the ascending arm in the deep/limnetic habitat, which is similar to what was observed in natural populations, although a larger difference in ρ was observed between treatments for this module in the experimental population. Another similarity between experimental and natural populations is that the lowest values of ρ were observed for the retro-articular and tooth bearing modules in lab-reared Tropheops. Across treatments, the experimental population also mirrored trends observed in natural populations insofar as deep/limnetic animals exhibiting generally higher levels of ρ relative to shallow/benthic animals.

Table 2 Comparison of within-module disparity and integration between depths for the mandible and pharyngeal jaw in the experimental populations

Finally, we found no differences in disparity between Tropheops individuals raised in environments that mimic deep or shallow habitats (Table 2). Similarly, we found no differences in disparity when we partition our landmarks into those modules defined by EMMLi and then compare disparity between depths (Table 2). These trends mirror those observed in the natural Tropheops populations (i.e., Fig. 5, Table S6). For example, the tooth-bearing module exhibited the highest disparity across treatments/environment for both the mandible and lower pharyngeal jaw.

In general, while differences were noted, trends in shape, disparity, and modularity were broadly similar between lab-reared and natural populations of Tropheops. These observations are consistent with the hypothesis that plasticity plays an important role in facilitating evolutionary divergence in response to selection for distinct foraging habitats (e.g., [45]).


Many species, two ecomorphs

Lake Malawi hosts the greatest number of cichlid species of any of the East African Rift Valley lakes, and within this species-flocks, Tropheops exhibit one of the highest speciation rates [37]. Following the construction of an AFLP phylogenetic tree using a subset of Tropheops from the southern portion of the lake, we find that the evolution of this species complex may be characterized by multiple transitions between deep and shallow environments (approximately 11 transitions). Tropheops residing in more shallow habitats exhibited a more robust feeding apparatus, while those Tropheops members residing in deeper habitats exhibited a more slender and gracile feeding apparatus [38]. The specific differences in morphology are predicted to influence feeding performance and biomechanics [46, 47], and include the area for muscle or ligament attachment sites (i.e., increased maxillary shank area for A1 attachment, increase pharyngeal jaw keel for pharyngohyoideus attachment), as well as the length of bony process that would have the effect of changing the mechanical advantage (i.e., the ascending arm of the mandible, the wings of the pharyngeal jaw). These data formalize and extend previously published trends (e.g., [28, 30]), and suggest the existence of two general Tropheops ecomorphs within Lake Malawi.

Two ecomorphs, one covariance structure

While we observed conspicuous differences in morphology between shallow and deep habitats, we found little evidence for differences in patterns of modularity within any bone examined. These results are consistent with previous studies [18], and suggest that patterns of cichlid craniofacial modularity can be robust to changes in habitat and functional demands. They also suggest that morphological diversification in Tropheops does not require ‘tinkering’ of an underlying pattern of covariation to produce adaptive phenotypic change. Indeed, recent evidence suggests there may be a finite amount of modularity patterns possible in the teleost skull [19], indicating a general conservation of covariation patterning could be common among taxa [7].

One possible evolutionary consequence of a conserved pattern of modularity is the ability of a lineage to undergo multiple transitions between discrete environments (e.g., [48]). Results from our stochastic character mapping supports this hypothesis, and suggests that transitions are fairly common and can occur in either direction between deep and shallow habitats. Seehausen [27] theorized that repeated transitions between these two habitats across cichlid lineages may occur due to ecological character displacement. In this model, a cichlid population would undergo intense local resource competition driving disruptive selection from shallow habitats toward deeper habitats. Subsequently, reproductive isolation would occur via depth-induced changes in signaling phenotypes (e.g., coloration due to visual sensitivities at different depths [49]). As resource competition grows in the deeper population, a transition back to the shallows may occur, and so on. We hypothesize that this iterative transitioning between habitats is aided by the conservation of craniofacial modularity. Felice et al. [50] used the metaphor of a fly in a tube to describe how trait covariation can influence morphological evolution, and demonstrated how a population can move around in morphospace, exhibiting different trait values, under a given covariational structure. Changing the pattern of covariation will, 1) take time, as this would likely involve the accumulation of multiple mutations over a period of generations, and 2) will change the range of morphologies that can evolve. Retaining a single pattern of covariation between habitats could allow morphological evolution to occur more rapidly (if the direction of selection aligns with the pattern of covariation [9]), but only if the limits of the morphospace includes those adaptive peaks for each habitat.

Alternatively, it is possible that a seemingly conserved pattern of modularity is a secondary consequence of shape being more evolvable than modularity. For instance, it is possible that Tropheops members have repeatedly and rapidly transitioned between shallow and deep environments due to fluctuating water levels in Lake Malawi. In this scenario, divergent patterns of modularity may not be able to evolve in such short time intervals between habitat transitions [51].

Evolvability of the feeding apparatus is consistent across habitats but distinct among anatomical modules

Rates of morphological evolution and magnitudes of disparity represent different metrics of evolutionary potential (e.g., evolvability), and we find that both measures are similar between habitats, indicating that neither deep nor shallow habitats have the capacity to constrain or facilitate morphological evolution. However, when comparing modules within bones, notable differences in evolutionary potential arise. Specifically, modules differ with respect to both rates and magnitudes of evolutionary change. In general, those regions of the anatomy more directly associated with muscle and ligamentous input evolve at a higher rate than modules with no obvious functional role. This trend is consistent with the adaptive value of these traits with respect to feeding kinematics. Indeed, the trend in functionally relevant components of anatomical structures exhibiting fast rates and high evolvability has been noted in other teleosts groups (i.e., Electric fishes [52, 53], cichlids, pomacentrids, centrarchids, and labrids [54]), and may reflect a general trend across organisms [55, 56]. In terms of disparity, it is notable that the tooth-plate module exhibited the highest levels for both the mandible and lower pharyngeal jaw. Given the ability of dentition to respond plastically to shifts in diet [57, 58], this trend suggests an important role for phenotypic plasticity in promoting disparity in the cichlid feeding apparatus.

Evolvability is not associated with the magnitude of phenotypic integration

We also document differences in the magnitude of phenotypic integration (ρ) between modules within bones. Differences in the strength of integration are predicted to have the capacity to either constrain or facilitate morphological evolution based on the direction of selection [59]. However, there is as yet no consensus for how (or whether) phenotypic integration influences evolutionary potential. Empirical studies looking at the association between evolutionary rates and integration have produced mixed results, ranging from no association [7], to positive [60] or negative associations between the two [10]. Simulation and theoretical studies have demonstrated that while morphological disparity can be governed by the strength of integration, more uncertainty surrounds the association between rates of morphological evolution and integration [7, 50]. Our results are more in line with these recent simulations. Specifically, they show little to no relationship between phenotypic integration and evolutionary rates or disparity. The only exception to this general trend is observed within the mandible (Fig. 5a-c), where there appears to be a negative correlation between disparity and integration; however, this association is largely driven by the tooth-bearing module which exhibits the highest level of disparity and lowest integration. Additionally, when we compare module parameters between habitats we find that, on average, taxa from the deeper habitats exhibit stronger within-module integration and faster rates of morphological evolution. Trends in disparity are more variable among feeding bones for a given module; shallow taxa typically exhibit greater mandible disparity, while deep taxa typically exhibit greater lower pharyngeal jaw disparity.

Plasticity plays a large role in influencing patterns of evolutionary change

By experimentally subjecting a shallow dwelling Tropheops species to environments that mimic either a shallow or deep habitat, we were able to replicate patterns detected across the Tropheops species complex. In other words, Tropheops species appear capable of craniofacial remodeling that mirrors variation observed at the species level. In terms of morphology, similar aspects of shape are being affected, despite the fact that they sometimes load on different PC axes (e.g., height of the ascending arm loads on PC1 for natural populations and PC3 in experimental populations). We find it especially notable that lab-reared populations of a single species, forced to feed using alternate strategies, possess the same pattern of modularity described for the larger subsample of the Tropheops species complex. These data suggest that modularity within craniofacial bones appears robust between natural and laboratory environments, and between foraging habitats and behaviors. Finally, many similarities are also observed in terms of within-module levels of integration and disparity. The conservation of these morphological patterns illustrates how plasticity in the feeding apparatus can be responsive to change within a specific covariation structure that promotes the accumulation of morphological variation in areas that would have direct biomechanical implications for feeding. If plasticity can replicate the patterns we observe in the natural population, this suggests that the environment has the capacity to play a large role in shaping Tropheops phenotypic evolution, at least during early stages of divergence. Empirical studies have demonstrated that phenotypic plasticity facilitates adaptive phenotypic diversification [61] and can vary in its effects among taxa [62]. Thus, plasticity may help to expedite transitions between depths allowing them to be ‘morphologically primed’ for diversification.

While plasticity can assist in this morphological divergence, it is likely that the occupation of a different habitat following a transition would also require some genetic and developmental evolution. We find evidence for this based on the larger magnitude of divergence between the shallow and deep members from the natural population relative to the experimental populations. It appears that plasticity cannot match the range of morphologies present in the natural populations, at least under the conditions we provided. Previous studies have also shown differences in the allele frequency of ptch1, a gene involved in directing the development of bones involved in the feeding system, among Tropheops members from different depth regimes [63, 64]. These studies suggest that genetic differentiation could underlie some of the differences observed in the natural populations of Tropheops, and plasticity is not solely responsible for differences between the natural populations. Indeed, previous work has demonstrated variability in the ability of closely-related cichlids to plastically remodel the craniofacial skeleton in response to alternate feeding regimes, whereby ecological specialists exhibited little response, while generalists exhibited a marked, and predicable, morphological response [62]. As such, in more ecologically generalist taxa such as Tropheops there may be a benefit to retaining plasticity in order to facilitate rapid colonization of new niches, and the transition between niches, rather than becoming genetically canalized [65]. By possessing a feeding system that is readily able to respond to biomechanical stimulus and remodel itself, Tropheops can avoid intense competition with both congeners and other ecological specialists. As a consequence, this may allow Tropheops to exploit a much greater spectrum of depth regimes. Thus, plasticity appears an important first step toward an evolutionary response [45].


Cichlid fish from the East African rift lakes are characterized by their extraordinary taxonomic and phenotypic diversity. Among the lakes, depth is considered a major axis of trophic niche partitioning among various cichlid populations [28, 66]. Habitat partitioning is the first component to many adaptive radiations [67], and the wealth of freshwater reefs, shallow bays, and isolated island environments provide ample ecological opportunity for cichlids to diversify along a habitat (i.e., depth) gradient [27]. Here we construct a phylogenetic tree focused primarily on a single species complex Tropheops, known to occupy multiple depth regimes in Lake Malawi, to examine how depth has shaped morphological evolution in the feeding apparatus. We found that selection drives feeding bone morphology toward different shapes based on their occupation of either shallow or deep environments. Despite differences in morphology we found no difference in disparity, rates of morphological evolution, or the pattern of modularity between Tropheops members residing at different depths. This indicates Tropheops exhibit a conserved pattern of modularity between depths, despite differences in morphology, a feature that may have facilitated the rapid depth transitions and subsequent colonizations. However, we do find differences in phenotypic integration and rates of morphological evolution among our modules, and these differences appear are most pronounced in those modules with functional roles (i.e., muscle attachment sites). Functionally salient modules are candidates for the greatest amount of morphological change given they must remodel in response to the differing biomechanical conditions experienced between the two depth regimes. Taken together, our data represent a rare example of morphological divergence coupled with module pattern conservation, indicating that multiple adaptive phenotypes can arise from a single pattern of covariation.


Specimen collection

We collected 58 individuals from across the Tropheops species complex from 12 different localities in the southern part of Lake Malawi during two field trips in 1996 and 2001. Our Tropheops sampling strategy obtained taxa from a wide depth gradient (deep, n = 24; shallow, n = 34). We skeletonized our specimens and then extracted four bones critical for food capture and processing (mandible, pharyngeal jaw, maxilla, and pre-maxilla). To gain 3D models of our bones of interest we used micro computed tomography (μCT) and scanned all specimens at 20–25 μm resolution at 90 kV and 75 μA. All μCT scans were performed using an X-Tek HMXST 225 μCT scanner (Nikon Corporation). We segmented the bones using Mimics v19 (Materialise NV), and then exported the 3D models for morphometric analysis.

We gathered a total of 136 individuals from multiple Lake Malawi cichlid genera along the southeast arm of Lake Malawi for phylogenetic analysis. Our sampling strategy focused on the subgenus Tropheops tropheops, and included 90 T. tropheops individuals from ~ 20 subspecies. Based on detailed transects performed by Ribbink et al. [38], we were able to comprehensively sample the Tropheops species complex from the southern part of Lake Malawi, implying that our results should provide a good estimate of the habitat transition rate. Many members of the T. tropheops species complex lack formal description, therefore species identification and nomenclature follow Ribbink et al. [38] and Konings [68]. We added a number of Lake Malawi cichlids to our sample from outside the T. tropheops complex including 18 Maylandia from four species and 17 non-mbuna from five genera. We also included two Tropheus duboisi individuals from Lake Tanganyika, which were purchased through the aquarium trade, bringing the total to 138 individuals (Table S1). We collected tissue samples from live fish via pectoral fin clips of each specimen and stored them in 95% EtOH prior to DNA extraction. Following DNA extraction, all specimens were skeletonized and stored. Whenever possible, individuals used in the phylogenetic analysis were also used as part of the morphological assessment.

Phylogenetic tree construction

We extracted genomic DNA from the fin clips of 1–3 individuals of each taxon by phenol-chloroform extraction. We used Amplified Fragment Length Polymorphisms (AFLP) to generate a character matrix for phylogenetic analysis. AFLP is a DNA fingerprinting technique that characterizes thousands of restriction polymorphisms spread throughout the genome [69]. Genomic DNA was first double-digested using two restriction enzymes EcoRI and MseI. Double stranded adapters are then ligated onto the overhanging ends of the fragments. An initial PCR reaction amplifies a subset of fragments that match adapter primers containing an additional nucleotide (EcoRI-A and MseI-C). The product of this first amplification is then used as the template for a further 18 different amplifications performed with primers containing an additional 2-nucleotide extension (E-ACA, M-CAA, M-CAG, M-CTA, M-CTT; E-ACC, M-CAA, M-CAC, M-CAT, M-CTA; E-ACT, M-CAG, M-CAT, M-CTA, M-CTG, M-CTT; E-AGC, M-CAG, M-CAT, M-CTA, M-CTG, M-CTT). For detailed restriction-ligation and PCR protocol information see [69, 70].

Fragments were separated using a Beckman Coulter CEQ 8000 capillary sequencer. Peaks were scored using a quartic model with a slope threshold of 2.0% and relative peak height of 5.0%. Bands were scored as present/absent using Beckman Coulter’s Fragment Analysis Module. The presence of each fragment was confirmed manually. Fragments between 80 and 500 bp in size were binned (1 nucleotide bin width) using Beckman Coulter’s AFLP Analysis Software. The binary output was imported to an Excel spreadsheet and formatted for MrBayes v3.2.6.

We performed a Bayesian phylogenetic analysis of our AFLP data in MrBayes v3.2.6 [71] running on the Cipres computing environment [72]. Our complete dataset contained 7953 binary characters scored for 138 individuals, with 90 individuals coming from the Tropheops species complex. We used the F81-like restriction site model intended for analyzing restriction fragment data with the ‘noabsencesites’ coding bias (characters that are absent in all taxa are unobservable) correction parameter [73]. While this binary state model implemented in MrBayes is a simplistic model of the actual genetic processes operating on the evolution of AFLP markers, methods which provide a more accurate representation of AFLP marker evolution, are computationally demanding for larger datasets [74, 75]. We set the Dirichlet prior for the state frequencies to 0.73, 0.27 matching the empirical 0/1 frequencies in the dataset and selected gamma variation across sites. We set our outgroup taxa as Tropheus duboisi from Lake Tanganikya, and forced monophyly for all ingroup Lake Malawi taxa. We then placed two topology constraints on our tree. We placed constraints to force monophyly in the node leading to all non-mbuna (mostly sand dwellers) and mbuna (rock dwellers), and then forced monophyly at the node leading to all Mbuna. To estimate the posterior distribution of our model parameters and calculate tree topology we performed two runs of a Markov chain Monte Carlo (MCMC) analysis for 2 million generations with ten chains (one cold, nine hot) at a chain temperature of 0.1 and sampled every 200 generations. We assessed stationary distributions based on the potential scale reduction factor (PSRF) and effective sample size (ESS) statistics reported in MrBayes [76]. We discarded the first 25% trees for burn-in, based on the stationary distributions reviewed in Tracer v1.7 [77].

We then randomly sampled 1000 post burn-in AFLP trees from the Bayesian posterior distribution (BPD) to account for phylogenetic uncertainty in all future comparative method analyses. These 1000 BPD trees were then converted into ultrametric trees. To force the AFLP trees to become ultrametric we used a penalized likelihood method contained in the chronos function from the R package ape [78]. Chronos uses a penalized likelihood method to convert branch lengths from the number of substitution per site to time. We set one calibration point at the split between Rhamphochromis and the rest of the Lake Malawi radiation using soft bounds of 1–2 million years, corresponding to the age of Lake Malawi [79].

To assess the transition rate between deep and shallow habitats we applied stochastic character maps (SIMMAP) to each of our 1000 BPD trees [40, 41] using the R package phytools v.0.6–44 [80]. We then calculated the average number of transitions between habitats, and examined the direction of these transitions to ascertain whether shallow to deep, or deep to shallow transitions occurred more frequently.

Geometric morphometric analysis

We placed landmarks (LMs) across all four bones to best characterize functionally and developmentally relevant aspects of shape change occurring in the Tropheops species complex (Fig. 7; Figure S2; maxilla: 4 fixed LMs, 40 semi-LMs; pre-maxilla: 5 fixed LMs, 20 semi-LMs; mandible: 14 fixed LMs, 95 semi-LMs; pharyngeal jaw: 10 fixed LMs, 35 semi-LMs). We present results for the maxilla and pre-maxilla bones in the supporting information (available online), and the mandible and lower pharyngeal jaw results in the main text, as these two bones have been the focus of many previous studies into trophic evolution and plasticity in cichlids. All landmark data was collected using Landmark v3.6 [81] and processed using the geomorph package v3.0.7 [82] in R v3.5.1 [83]. We used the digit.curves function in geomorph to resample and array our semi-landmarks equidistantly along a curve between two fixed landmarks [84]. Following the placement of landmarks we performed a Procrustes superimposition to remove the effects of size, translation, and rotation such that the landmark configurations are in register [85]. To investigate the effects of allometry on our shape data we performed a Procrustes ANOVA between centroid size and shape for each bone. We found a significant effect of allometry for two bones (mandible: r2 = 0.050; P < 0.001; lower pharyngeal jaw: r2 = 0.053; P = 0.006). To remove the effects of allometry on all of our bones, we performed a regression of shape on geometric centroid size to generate a landmark data set based on residuals. Following allometric correction, we calculated mean shape configurations for each operational taxonomic unit (OTU) for use with future comparative methods, with each OTU represented by approximately three individuals (deep, OTU n = 8; shallow OTU n = 14). We then conducted a principal components (PC) analysis on all landmark configurations to reduce the shape data into a series of orthogonal axes that best explained the major variation in bone shapes among Tropheops taxa. We extracted five PCs from each bone, as subsequent PC axes accounted for < 5% of the variation in shape (Table S2). We constructed morphospaces for each bone from the first three PCs, and colored each individual based on their depth assignment, either shallow or deep, to illustrate the morphological variation present in the sample.

Fig. 7
figure 7

Landmarking scheme for Tropheops feeding bones. a, mandible lateral view; b, mandible dorsal view; c, pharyngeal jaw dorsal view; d, pharyngeal jaw ventral view. Red circles, fixed landmarks; blue landmarks, semi-landmark curve positions

Modularity, integration, and disparity


We used EMMLi (evaluating modularity with maximum likelihood) to simultaneously test multiple modularity hypotheses based on different landmark partitions in our Tropheops taxa (Fig. 3), and then assess the strength of covariation within those modules [42]. The EMMLi method compares a suite of modularity hypotheses using AICc to determine the best fitting modularity hypothesis. EMMLi extracts a vector based correlation matrix with size determined by the number of landmarks. This symmetric matrix provides a single number for covariance between landmarks, as opposed to one for each dimension, and are considered more biologically meaningful as each landmark is considered a separate trait [42]. We computed two vector correlation matrices for each bone; one matrix calculated using Tropheops taxa from the deep habitat and another matrix for those from a shallow habitat. EMMLi tests the fit of different partition hypotheses by allowing modules to either constrain or vary the correlation coefficient (ρ) within and among different modules. To correct for the evolutionary associations among our taxa, we used the phylogenetic independent contrasts (PICs) of shape in the EMMLi analyses [86, 87]. PICs of shape were generated from the 1000 post burn-in AFLP trees and used to derive a distribution of modularity analyses and ρ values. We used these 1000 EMMLi analyses to determine how frequently, and how strongly (based on AICc score), competing modularity hypotheses were selected across both habitats. Given the discrepancy in sample size between deep and shallow taxa (deep, OTU n = 8; shallow OTU n = 14), we randomly sampled individuals from our shallow population without replacement to mimic the sample size of the deep population 1000 times. To determine how sensitive our results were to sample size differences we then re-ran our EMMLi analysis on the ‘sampled’ shallow population and compared this with results from the original shallow population. The frequency of incorrectly selected modularity models will give an indication into how sensitive the data are to sample size.

Module partition hypotheses varied from more simplistic models whereby modules are split into anterior-posterior or dorsal-ventral landmark partitions, to more complex models that attempt to capture regions of functional importance, such as sites of muscle or ligament attachment, or regions where bones articulate (see Table S4 for all partition hypotheses). We performed the EMMLi analysis on two different datasets with the intention of investigating any similarity or differences in the best-supported partitions. The first analysis contained all Tropheops individuals sampled from the lake (natural population), and a second analysis contained a single species of Tropheops from a feeding trial experiment (experimental population, see the ‘Dietary plasticity’ section below). Differences in partitions between depths would suggest habitat might have differentially influenced the covariation structure of our bones, which would have implications for evolvability and differences in function or/and development.

Comparative methods

As we often had multiple individuals for each Tropheops ‘species’, we calculated a mean landmark configuration and PC score for each species. We performed all subsequent comparative methods on these mean data. For all comparative methods we pruned any tips from the tree where we lacked trait data, which left a total of 22 Tropheops taxa: 14 from the shallow environment and 8 from the deep. All comparative methods were performed using a random sample of 1000 post burn-in AFLP trees to account for phylogenetic uncertainty and generate distributions of output parameters.

Phylogenetic MANOVA

We used a phylogenetic multivariate analysis of variance (pMANOVA) to test for differences in bone shape between Tropheops from shallow and deep habitats. The test statistic for the pMANOVA was calculated using the first five PC scores from each bone and compared to a null distribution of PC scores. We generated a null distribution by simulating 1000 new dependent variables based on the rate matrix obtained from the phylogenetic tree. The pMANOVA was conducted using the R package Geiger v.2.0.6 with the aov.phylo function [88].


We compared the disparity of Tropheops taxa residing in either shallow or deep environments using Procrustes variance (PV). PV measures the distribution of taxa around the mean shape for a given habitat grouping. We used the morphol.disparity function in the R package geomorph to conduct the PV analysis [82]. During the disparity calculation for the mean Tropheops shapes we used landmark configurations that retained their allometric shape components and added in centroid size as a variable in the model. Phylogeny was accounted for in the model following the use of the procD.pgls function in the R package geomorph. Disparity is therefore calculated using the residuals obtained from a phylogenetic least squares estimation of coefficients, instead of the ordinary least squares estimation used in the standard PV analysis on the experimental Tropheops bone shapes. We also repeated the disparity analysis using subsets of landmarks defined by the best fitting modularity hypothesis determined by EMMLi and compared module-specific disparity between habitats. Differences in disparity between habitats would suggest some type of morphological constraint or morphological opportunity is operating on one habitat relative to the other, or would imply some differences exist in developmental robustness.

Rates of morphological evolution

We used the function compare.evol.rates in geomorph to compare the rate of morphological evolution between habitats using a complete landmark dataset for each bone, and a landmark dataset subset by the best-supported modularity hypotheses defined by EMMLi [89]. Note that we are limited to assuming a Brownian Motion (BM) model of evolution for this analysis, however the results returned should be conservative relative to a model of evolution that more closely fits our data.

We used the function compare.multi.evol.rates in the R package geomorph to compare the rate of morphological evolution among different modules in our mandible and lower pharyngeal jaw data sets under a BM model of evolution. Significance is determined via simulating tip data under BM and comparing this to the empirical data [90].

Dietary plasticity

To assess whether trends in the pattern of modularity and morphological disparity observed in the natural populations can be replicated via phenotypic plasticity, we conducted a feeding trial experiment to replicate the conditions experienced in deep or shallow environments using a single species, Tropheops ‘red cheek’ (TRC). While TRC naturally resides in more shallow habitats, previous studies have demonstrated that this ectomorph can exhibit a plastic phenotype [32]. We divided our experimental population of TRC into two tanks at 1 month old and subjected them two diet treatments that differed in the biomechanical demand placed on the feeding system to mimic the shallow-deep habitats. We raised 19 individuals in a shallow environment habitat whereby individuals consumed puréed algae flake food mixed with 1.5% food-grade agar spread over two lava rocks fed once per day. TRCs in the shallow environment must use a ‘tooth scraping’ feeding mode to separate the algae from the rocks, a feeding behavior that is considered more biomechanically challenging. We also raised 17 individuals in an environment mimicking a deep habitat. TRCs in the deep environment were provided with the same amount of food but ground down and sprinkled into the tank and fed once per day. This forced the TRCs to exhibit a greater degree of suction and ram feeding behaviors to consume their food, a feeding behavior that is considered less biomechanically challenging [32]. The feeding trial continued for 5 months on a 12 h light-dark cycle in water conditions mimicking that of Lake Malawi. Upon completion, all animals were euthanized by overdose with tricaine methanesulfonate (Aquatic Ecosystems Inc.), fixed overnight at room temperature in 4% paraformaldehyde (Sigma) in 1x phosphate-buffered saline, and stored in 70% ethanol. The Institutional Animal Care and Use Committee (IACUC) at the University of Massachusetts Amherst and the University of New Hampshire approved all protocols.

We μCT scanned and landmarked all experimental fish as described in the ‘Specimen collection’ and ‘Geometric morphometric analysis’ sections above. We placed landmarks for all individuals using the same schemes used in the evolutionary analysis and again tested for allometry in our bones of interest. We found a significant effect of allometry on our shape data (mandible: r2 = 0.158; P < 0.001; lower pharyngeal jaw: r2 = 0.207; P < 0.001). We removed the allometric component of shape as described above to generate a landmark data set based on residuals. Finally, we replicated the EMMLi, disparity, and MANOVA analyses described above in the experimental populations to assess differences between those individuals raised in deep versus shallow conditions. Note that we used non-phylogenetic versions of the MANOVA and disparity analyses as we are examining the plastic response of a single Tropheops species.

Availability of data and materials

All morphometric data and phylogenetic trees are contained in the supplementary information.



Amplified Fragment Length Polymorphism


Akaike Information Criterion corrected;


Brownian Motion


Bayesian Posterior Distribution


Evaluating Modularity with Maximum Likelihood


Effective Sample Size


Institutional Animal Care and Use Committee




Multivariate Analysis of Variance


Markov Chain Monte Carlo


Operational Taxonomic Unit


Phylogenetic Independent Contrasts


Principal Component


Principal Component Analysis


Procrustes Variance


Potential Scale Reduction Factor


Tropheops ‘red cheek


Micro-Computed Tomography


  1. Darwin CR. On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life. London: John Murray; 1859.

    Google Scholar 

  2. Rolian C, Willmore KE. Morphological integration at 50: patterns and processes of integration in biological anthropology. Evol Biol. 2009;36:1–4.

    Article  Google Scholar 

  3. Olson EC, Miller RL. Morphological integration. Chicago: University of Chicago Press; 1958.

    Google Scholar 

  4. Goswami A. Morphological integration in the carnivoran skull. Evolution. 2006;60:169–83.

    PubMed  Google Scholar 

  5. Goswami A, Polly PD. The influence of modularity on cranial morphological disparity in Carnivora and Primates (Mammalia). PLoS One. 2010;5:e9517.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Jacob F. Evolution and tinkering. Science. 1977;196:1161–6.

    CAS  PubMed  Google Scholar 

  7. Goswami A, Smaers JB, Soligo C, Polly PD. The macroevolutionary consequences of phenotypic integration : from development to deep time. Philos Trans R Soc B Biol Sci. 2014;369:20130254.

    CAS  Google Scholar 

  8. Klingenberg CP. Studying morphological integration and modularity at multiple levels: concepts and analysis. Philos Trans R Soc Lond Ser B Biol Sci. 2014;369:20130249.

    Article  Google Scholar 

  9. Klingenberg CP. Evolution and development of shape: integrating quantitative approaches. Nat Rev Genet. 2010;11:623–35.

    CAS  PubMed  Google Scholar 

  10. Felice RN, Goswami A. Developmental origins of mosaic evolution in the avian cranium. Proc Natl Acad Sci. 2018;115:555–60.

    CAS  PubMed  Google Scholar 

  11. Wagner GP, Pavlicev M, Cheverud JM. The road to modularity. Nat Rev Genet. 2007;8:921–31.

    CAS  PubMed  Google Scholar 

  12. Conith AJ, Meagher MA, Dumont ER. The influence of climatic variability on morphological integration, evolutionary rates, and disparity in the Carnivora. Am Nat. 2018;191:704–15.

    Article  PubMed  Google Scholar 

  13. Young NM, Hallgrímsson B. Serial homology and the evolution of mammalian limb covariation structure. Evolution. 2005;59:2691–704.

    PubMed  Google Scholar 

  14. Sanger TJ, Mahler DL, Abzhanov A, Losos JB. Roles for modularity and constraint in the evolution of cranial diversity among Anolis lizards. Evolution. 2012;66:1525–42.

    Article  PubMed  Google Scholar 

  15. Larouche O, Zelditch ML, Cloutier R. Modularity promotes morphological divergence in ray-finned fishes. Sci Rep. 2018;8:7278.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Smith AJ, Nelson-Maney N, Parsons KJ, James Cooper W, Craig AR. Body shape evolution in sunfishes: divergent paths to accelerated rates of speciation in the Centrarchidae. Evol Biol. 2015;42:283–95.

    Article  Google Scholar 

  17. Cooper WJ, Parsons K, McIntyre A, Kern B, McGee-Moore A, Albertson RC. Bentho-pelagic divergence of cichlid feeding architecture was prodigious and consistent during multiple adaptive radiations within African rift-lakes. PLoS One. 2010;5:e9551.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Parsons KJ, Márquez E, Albertson RC. Constraint and opportunity: the genetic basis and evolution of modularity in the cichlid mandible. Am Nat. 2012;179:64–78.

    Article  PubMed  Google Scholar 

  19. Parsons KJ, Son YH, Crespel A, Thambithurai D, Killen S, Harris MP, et al. Conserved but flexible modularity in the zebrafish skull: implications for craniofacial evolvability. Proc R Soc B Biol Sci. 2018;285:20172671.

    Google Scholar 

  20. Payne JL, Wagner A. The causes of evolvability and their evolution. Nat Rev Genet. 2019;20:24–38.

    Article  CAS  PubMed  Google Scholar 

  21. Fryer G, Iles TD. The cichlid fishes of the Great Lakes of Africa. Their biology and evolution. Edinburgh: Oliver and Boyd; 1972.

    Google Scholar 

  22. Parsons KJ, Cooper WJ, Albertson RC. Modularity of the Oral jaws is linked to repeated changes in the craniofacial shape of African cichlids. Int J Evol Biol. 2011;2011:1–10.

    Article  Google Scholar 

  23. Hu Y, Parsons KJ, Albertson RC. Evolvability of the cichlid jaw: new tools provide insights into the genetic basis of phenotypic integration. Evol Biol. 2014;41:145–53.

    Article  Google Scholar 

  24. Cooper WJ, Wernle J, Mann K, Albertson RC. Functional and genetic integration in the skulls of Lake Malawi cichlids. Evol Biol. 2011;38:316–34.

    Google Scholar 

  25. Albertson RC, Powder KE, Hu Y, Coyle KP, Roberts RB, Parsons KJ. Genetic basis of continuous variation in the levels and modular inheritance of pigmentation in cichlid fishes. Mol Ecol. 2014;23:5135–50.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Turner GF, Seehausen O, Knight ME, Allender CJ, Robinson RL. How many species of cichlid fishes are there in African lakes? Mol Ecol. 2001;10:793–806.

    CAS  PubMed  Google Scholar 

  27. Seehausen O. Process and pattern in cichlid radiations – inferences for understanding unusually high rates of evolutionary diversification. New Phytol. 2015;207:304–12.

    PubMed  Google Scholar 

  28. Albertson RC. Morphological divergence predicts habitat partitioning in a Lake Malawi cichlid species complex. Copeia. 2008;2008:689–98.

    Article  Google Scholar 

  29. Parsons KJ, Son YH, Albertson RC. Hybridization promotes Evolvability in African cichlids: connections between Transgressive segregation and phenotypic integration. Evol Biol. 2011;38:306–15.

    Google Scholar 

  30. Albertson RC, Pauers MJ. Morphological disparity in ecologically diverse versus constrained lineages of Lake Malaŵi rock-dwelling cichlids. Hydrobiologia. 2018.

  31. Anderson PSL, Renaud S, Rayfield EJ. Adaptive plasticity in the mouse mandible. BMC Evol Biol. 2014;14:85.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Parsons KJ, Concannon M, Navon D, Wang J, Ea I, Groveas K, et al. Foraging environment determines the genetic architecture and evolutionary potential of trophic morphology in cichlid fishes. Mol Ecol. 2016;25:6012–23.

    PubMed  Google Scholar 

  33. Parsons KJ, Rigg A, Conith AJ, Kitchener AC, Harris S, Zhu H. Skull morphology diverges between urban and rural populations of red foxes mirroring patterns of domestication and macroevolution. Proc R Soc B Biol Sci. 2020;287:20200763.

    Article  CAS  Google Scholar 

  34. Randau M, Sanfelice D, Goswami A. Shifts in cranial integration associated with ecological specialization in pinnipeds (Mammalia, Carnivora). R Soc Open Sci. 2020;6:190201.

    Article  Google Scholar 

  35. Klingenberg CP. Morphological integration and developmental modularity. Annu Rev Ecol Evol Syst. 2008;39:115–32.

    Google Scholar 

  36. Waddington CH. Canalization of development and the inheritance of acquired characters. Nature. 1942;150:563–5.

    Google Scholar 

  37. Won Y-J, Sivasundar A, Wang Y, Hey J. On the origin of Lake Malawi cichlid species: a population genetic analysis of divergence. Proc Natl Acad Sci. 2005;102(suppl 1):6581–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Ribbink AJ, Marsh BA, Marsh AC, Ribbink AC, Sharp BJ. A preliminary survey of the cichlid fishes of rocky habitats in Lake Malawi. South African J Zool. 1983;18:149–309.

    Google Scholar 

  39. Malinsky M, Svardal H, Tyers AM, Miska EA, Genner MJ, Turner GF, et al. Whole-genome sequences of Malawi cichlids reveal multiple radiations interconnected by gene flow. Nat Ecol Evol. 2018;2:1940–55.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Huelsenbeck JP, Nielsen R, Bollback JP. Stochastic mapping of morphological characters. Syst Biol. 2003;52:131–58.

    Article  PubMed  Google Scholar 

  41. Bollback JP. SIMMAP: stochastic character mapping of discrete traits on phylogenies. BMC Bioinformatics. 2006;7:88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Goswami A, Finarelli JA. EMMLi: a maximum likelihood approach to the analysis of modularity. Evolution. 2016;70:1622–37.

    Article  PubMed  Google Scholar 

  43. Moser FN, Van Rijssel JC, Mwaiko S, Meier JI, Ngatunga B, Seehausen O. The onset of ecological diversification 50 years after colonization of a crater lake by haplochromine cichlid fishes. Proc R Soc B Biol Sci. 2018;285:20180171.

    Google Scholar 

  44. Malinsky M, Challis RJ, Tyers AM, Schiffels S, Terai Y, Ngatunga BP, et al. Genomic islands of speciation separate cichlid ecomorphs in an east African crater lake. Science. 2015;350:1493–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. West-Eberhard MJ. Developmental plasticity and evolution. New York: Oxford University Press; 2003.

    Google Scholar 

  46. Westneat MW. Feeding mechanics of teleost fishes (Labridae; Perciformes): A test of four-bar linkage models. J Morphol. 1990;205:269–95.

    Article  PubMed  Google Scholar 

  47. Conith AJ, Lam DT, Albertson RC. Muscle-induced loading as an important source of variation in craniofacial skeletal shape. Genesis. 2019;57:e23263.

    PubMed  Google Scholar 

  48. Cooper WJ, Westneat MW. Form and function of damselfish skulls: rapid and repeated evolution into a limited number of trophic niches. BMC Evol Biol. 2009;9.

  49. Seehausen O, Terai Y, Magalhaes IS, Carleton KL, Mrosso HDJ, Miyagi R, et al. Speciation through sensory drive in cichlid fish. Nature. 2008;455:620–6.

    CAS  PubMed  Google Scholar 

  50. Felice RN, Randau M, Goswami A. A fly in a tube : macroevolutionary expectations for integrated phenotypes. Evolution. 2018;72:2580–94.

    PubMed  PubMed Central  Google Scholar 

  51. Ivory S, Cohen AS, Ivory SJ, Blome MW, King JW, Mcglue MM, et al. Environmental change explains cichlid adaptive radiation at Lake Malawi over the past 1.2 million years. Proc Natl Acad Sci. 2016;113:11895–900.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. Evans KM, Vidal-García M, Tagliacollo VA, Taylor SJ, Fenolio DB. Bony patchwork: mosaic patterns of evolution in the skull of electric fishes (Apteronotidae: Gymnotiformes). Integr Comp Biol. 2019;59:420–31.

    Article  PubMed  Google Scholar 

  53. Evans KM, Kim LY, Schubert BA, Albert JS. Ecomorphology of Neotropical electric fishes: an integrative approach to testing the relationships between form, function, and trophic ecology. Integr Org Biol. 2019;1.

  54. Holzman R, Collar DC. Price S a, Hulsey CD, Thomson RC, wainwright PC. Biomechanical trade-offs bias rates of evolution in the feeding apparatus of fishes. Proc R Soc B Biol Sci. 2012;279:1287–92.

    Article  Google Scholar 

  55. Evans KM, Waltz BT, Tagliacollo VA, Sidlauskas BL, Albert JS. Fluctuations in evolutionary integration allow for big brains and disparate faces. Sci Rep. 2017;7:40431.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Muñoz MM, Hu Y, Anderson PSL, Patek SN. Strong biomechanical relationships bias the tempo and mode of morphological evolution. Elife. 2018;7:e37621.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Gunter HM, Fan S, Xiong F, Franchini P, Fruciano C, Meyer A. Shaping development through mechanical strain: the transcriptional basis of diet-induced phenotypic plasticity in a cichlid fish. Mol Ecol. 2013;22:4516–31.

    PubMed  Google Scholar 

  58. Muschick M, Barluenga M, Salzburger W, Meyer A. Adaptive phenotypic plasticity in the Midas cichlid fish pharyngeal jaw and its relevance in adaptive radiation. BMC Evol Biol. 2011;11:116.

    PubMed  PubMed Central  Google Scholar 

  59. Hallgrímsson B, Jamniczky H, Young NM, Rolian C, Parsons TE, Boughner JC, et al. Deciphering the palimpsest: studying the relationship between morphological integration and phenotypic Covariation. Evol Biol. 2009;36:355–76.

    PubMed  PubMed Central  Google Scholar 

  60. Hu Y, Ghigliotti L, Vacchi M, Pisano E, Detrich HW, Albertson RC. Evolution in an extreme environment: developmental biases and phenotypic integration in the adaptive radiation of antarctic notothenioids. BMC Evol Biol. 2016;16:142.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Rajkov J, Weber AA-T, Salzburger W, Egger B. Adaptive phenotypic plasticity contributes to divergence between lake and river populations of an east African cichlid fish. Ecol Evol. 2018;8:7323–33.

    PubMed  PubMed Central  Google Scholar 

  62. Parsons KJ, Taylor AT, Powder KE, Albertson RC. Wnt signalling underlies the evolution of new phenotypes and craniofacial variability in Lake Malawi cichlids. Nat Commun. 2014;5:1–11.

    Article  CAS  Google Scholar 

  63. Hu Y, Albertson RC. Hedgehog signaling mediates adaptive variation in a dynamic functional system in the cichlid feeding apparatus. Proc Natl Acad Sci. 2014;11:8530–4.

    Google Scholar 

  64. Roberts RB, Hu Y, Albertson RC, Kocher TD. Craniofacial divergence and ongoing adaptation via the hedgehog pathway. Proc Natl Acad Sci. 2011;108:13194–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. Wund MA, Baker JA, Clancy B, Golub JL, Foster SA. A test of the “flexible stem” model of evolution: ancestral plasticity, genetic accommodation, and morphological divergence in the threespine stickleback radiation. Am Nat. 2008;172:449–62.

    PubMed  Google Scholar 

  66. Parnell NF, Streelman JT. The macroecology of rapid evolutionary radiation. Proc R Soc B Biol Sci. 2011;278:2486–94.

    Google Scholar 

  67. Streelman TJ, Danley PD. The stages of vertebrate evolutionary radiation. Trends Ecol Evol. 2003;18:126–31.

    Article  Google Scholar 

  68. Konings A. Malawi Cichlids in Their Natural Habitat. 4th edition. El Paso: Cichlid Press; 2007.

  69. Vos P, Hogers R, Bleeker M, Reijans M, van de Lee T, Hornes M, et al. AFLP: a new technique for DNA fingerprinting. Nucleic Acids Res. 1995;23:4407–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Albertson RC, Markert JA, Danley PD, Kocher TD. Phylogeny of a rapidly evolving clade: the cichlid fishes of Lake Malawi, East Africa. Proc Natl Acad Sci. 1999;96:5107–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61:539–42.

    PubMed  PubMed Central  Google Scholar 

  72. Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: Proceedings of the Gateway Computing Environments Workshop (GCE). New Orleans; 2010. p. 1–8.

  73. Felsenstein J. Phylogenies from restriction sites: a maximum-likelihood approach. Evolution. 1992;46:159–73.

    PubMed  Google Scholar 

  74. Luo R, Hipp A, Larget B. A Bayesian Model of AFLP Marker Evolution and Phylogenetic Inference. Stat Appl Genet Mol Biol. 2007;6:article 11.

    Google Scholar 

  75. Koopman WJM, Wissemann V, De Cock K, Van Huylenbroeck J, De Riek J, Sabatino GJH, et al. AFLP markers as a tool to reconstruct complex relationships: a case study in Rosa (Rosaceae). Am J Bot. 2008;95:353–66.

    CAS  PubMed  Google Scholar 

  76. Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed Phylogenetics and dating with confidence. PLoS Biol. 2006;4:e88.

    PubMed  PubMed Central  Google Scholar 

  77. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian Phylogenetics using tracer 1.7. Syst Biol. 2018;67:901–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. Paradis E, Claude J, Strimmer K. APE: analyses of Phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90.

    Article  CAS  PubMed  Google Scholar 

  79. Sturmbauer C, Salzburger W, Duftner N, Schelly R, Koblmüller S. Molecular Phylogenetics and evolution evolutionary history of the Lake Tanganyika cichlid tribe Lamprologini ( Teleostei : Perciformes ) derived from mitochondrial and nuclear DNA data. Mol Phylogenet Evol. 2010;57:266–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  Google Scholar 

  81. Wiley DF. Landmark Editor 3.0. 2006.

    Google Scholar 

  82. Adams DC, Otárola-Castillo E. Geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol Evol. 2013;4:393–9.

    Article  Google Scholar 

  83. R Core Team. R: A Language and Environment for Statistical Computing. 2018.

    Google Scholar 

  84. Bookstein FL. Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape. Med Image Anal. 1997;6:225–43.

    Google Scholar 

  85. Rohlf FJ. On applications of geometric morphometrics to studies of ontogeny and phylogeny. Syst Biol. 1998;47:147–67.

    CAS  PubMed  Google Scholar 

  86. Felsenstein J. Phylogenies and the comparative method. Am Nat. 1985;125:1–15 Accessed 29 Mar 2016.

    Google Scholar 

  87. Bardua C, Wilkinson M, Gower DJ, Sherratt E, Goswami A. Morphological evolution and modularity of the caecilian skull. BMC Evol Biol. 2019;19:30.

    Article  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  Google Scholar 

  89. Adams DC. Quantifying and comparing phylogenetic evolutionary rates for shape and other high-dimensional phenotypic data. Syst Biol. 2014;63:166–77.

    Article  PubMed  Google Scholar 

  90. Denton JSS, Adams DC. A new phylogenetic test for comparing multiple high-dimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lanternfishes (Myctophiformes; Myctophidae). Evolution. 2015;69:2425–40.

    PubMed  Google Scholar 

Download references


We thank Aaron Dobberfuhl, Mathius Msafiri, Alex Pollen, Suzy Renn, and Caroly Shumway for assistance in the field and Rodney Rohde (Texas State University – San Marcos) for access to the Beckman Coulter CEQ8000. We also wish to acknowledge members of the Albertson lab for comments on early versions of this manuscript.


AJC and RCA were supported by a National Institutes of Health award #R01DE026446. TDK was supported by a National Science Foundation award #9905127. The funding body was not involved in the design of the study and collection, analysis, and interpretation of data, and in writing the manuscript.

Author information

Authors and Affiliations



AJC devised study, performed morphometric analysis, built phylogenetic tree, and wrote the manuscript. MRK collected specimens, extracted DNA, gathered sequence data, and wrote the manuscript. TDK collected specimens, and wrote the manuscript. RCA devised study, collected specimens, and wrote the manuscript. All authors have read and approved the final version of the manuscript.

Corresponding author

Correspondence to Andrew J. Conith.

Ethics declarations

Ethics approval and consent to participate

All work was performed in compliance with the Institutional Animal Care and Use Committee at UMass Amherst (#2018–0094 to RCA) and University of New Hampshire (#990603 to TDK). Collection permits were issued though the University of Malawi and the Malawi government.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

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

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Conith, A.J., Kidd, M.R., Kocher, T.D. et al. Ecomorphological divergence and habitat lability in the context of robust patterns of modularity in the cichlid feeding apparatus. BMC Evol Biol 20, 95 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: