- Open Access
Developmental basis of phenotypic integration in two Lake Malawi cichlids
EvoDevovolume 7, Article number: 3 (2016)
Cichlid fishes from the Rift Lakes of East Africa have undergone the most spectacular adaptive radiations in vertebrate history. Eco-morphological adaptations in lakes Victoria, Malawi and Tanganyika have resulted in a vast array of skull shapes and sizes, yet primary axes of morphological variation are conserved in all three radiations, prominently including the size of the preorbital region of the skull. This conserved pattern suggests that development may constrain the trajectories of cichlid head morphological evolution.
Here, we (1) present a comparative analysis of adult head morphology in two sand-dweller cichlids from Lake Malawi with preorbital size differences representative of the main axis of variation among the three lakes and (2) analyze the ontogeny of shape and size differences by focusing on known developmental modules throughout the head. We find that (1) developmental differences between the two species correlate with known developmental modules; (2) differences in embryonic cartilage development result in phenotypically integrated changes among all bones derived from a single cartilage, while differences in dermal bone development tend to influence isolated regions within a bone; and lastly (3) species-specific morphologies appear in the embryo as subtle differences, which become progressively amplified throughout ontogeny. We propose that this amplification takes place at skeletal growth zones, the locations and shapes of which are patterned during embryogenesis.
This study is the most anatomically comprehensive analysis of the developmental differences underlying cichlid skull evolution in the Rift Lakes of East Africa. The scale of our analysis reveals previously unnoticed correlations between developmental modules and patterns of phenotypic integration. We propose that the primary axes of morphological variation among East African cichlid adaptive radiations are constrained by the hierarchical modularity of the teleost head skeleton.
Dawkins (1988) introduced the term “evolvability” to describe the capacity of a lineage to generate heritable phenotypic variation (see also Alberch 1991) [1, 2]. The manner and degree to which developmental mechanisms determine evolvability is one of the major questions of the extended evolutionary synthesis [3–5]. The evolution of the head has been described as the key innovation that allowed the further diversification and evolutionary success of the vertebrate lineage . However, our understanding of the developmental factors that have enabled the massive morphological and functional diversification of the head remains limited. Despite considerable progress in identifying cranial developmental modules and their underlying gene regulatory networks [7–13], little is known about how these modular units structure evolvability.
Some of the best-characterized cranial developmental modules are the pharyngeal arches, an array of serially homologous, bilateral structures that segment the head ventrally . Their derivatives form the entire facial skeleton, which in ray-finned fishes include the ethmoid, upper and lower jaws, suspensorium, opercular series and gill arches. The modular nature of these arches along the antero-posterior axis has been demonstrated by the segment-specific expression and functions of Hox genes in patterning arch identity [15–17]. A second set of genes including, but not limited to, jag1b, nkx3.2 and hand2 subdivides individual pharyngeal arches into dorsal, intermediate and ventral domains, respectively, that give rise to correspondingly located skeletal elements with specific shapes and sizes . The overlap of these two orthogonal patterning systems defines semiautonomous developmental modules through which individual skeletal elements could evolve independently from each other. Changes affecting higher-level modules may affect several adjacent skeletal elements, resulting in phenotypically integrated changes in shape and size, while, in contrast, changes in lower level modules may affect individual bones. In addition, skeletal modularity depends to some extent on the mode of ossification, e.g., individual cartilages may subdivide into multiple distinct bones, while dermal bones often result from the fusion of independent ossification centers . The hierarchical modularity of the vertebrate head may provide useful clues to understanding how patterns of head evolution have been shaped by the developmental constraints imposed by this organization.
Rapid diversification in cranial form has repeatedly generated extensive ecological diversity in multiple lineages of cichlid fishes, which makes them ideal for exploring the developmental properties associated with evolvability [20–24]. Cichlids from the Great East African Rift Lakes (Victoria, Malawi and Tanganyika) have produced some of the most spectacular adaptive radiations in vertebrate history. Their cranial skeletons have evolved extremely quickly, and this has facilitated the invasion of an extraordinary diversity of trophic niches [20–24]. Young et al.  and Cooper et al.  showed that the primary axes of evolutionary divergence in all three lakes were strongly associated with the size of the preorbital region, suggesting that developmental constraints may have played an important role in shaping their skull evolution.
The adaptive radiation of the Lake Malawi cichlids is currently at a particularly interesting and useful stage for evolutionary developmental biologists. Although it is still very young, a tremendous range of skull shapes has already arisen. The evolution of skull diversity has outpaced the evolution of absolute barriers to interbreeding, such that many species with extremely different cranial morphologies will interbreed in captivity. This dramatically facilitates using both developmental and genetic mapping approaches to study cichlid head evolvability [25–29].
There are two broad categories of environments in Lake Malawi that are accessible to cichlids, and early in their evolutionary history, this lineage split into two sister clades, each of which then diversified further within one of these realms. The names of these lineages reflect their primary habitat use: the mbuna, or rock-dwellers, and the sand-dwellers . The major evolutionary axis for both of these clades mirrors the pattern of morphological divergence seen in all three lakes (i.e., diversification in preorbital size and jaw length; [22, 24]). Genetic and developmental investigations of Lake Malawi cichlid skull evolution have primarily focused on isolated modules, such as the jaw skeleton in mbuna species [25, 26, 31–36]. In order to test whether the hierarchical modularity of the cichlid head has constrained its evolution, we conducted the most anatomically complete comparative analysis of cichlid head skeleton development to date. By comparing the cranial development of two sand-dweller species with extreme differences in skull morphology both to each other and to the offspring of a hybrid cross between them, we were able to gain insight into the ontogenetic basis of their divergent morphologies. Since the parent species exhibit strong differences in preorbital size and jaw length, this study has strong relevance to determining the developmental factors that have shaped the cichlid radiations within the Great East African Rift Lakes.
We examined the morphological divergence and craniofacial development of Copadichromis azureus (CA), which is an omnivorous cichlid with small jaws ; and Dimidiochromis compressiceps (DC), also known as the “Malawi eye-biter,” which is a predator of small fishes (mainly juveniles of Utaka or other shoaling species) with a large mouth and fast biting jaws . We find that morphological differences in adult head shape originate from a few centers that correlate with known developmental modules. Ontogenetic divergence in skeletal shape and/or size is detected for many elements at embryonic and larval stages. As expected, phenotypically integrated differences in adult bones derived from a common cartilage precursor seem to originate from differences in the development of the precursor cartilage itself. In contrast, differences in adult dermal bone morphology seem more independent from each other, which correlates with their more autonomous development. Our observations also indicate that large differences in adult morphology start as subtle module-specific differences in skeletal shape, which are amplified at subsequent stages through differential growth. This reveals a previously unsuspected role for skeletal growth zones in the morphological divergence of cichlids, which will form the focus of our future developmental and genetic mapping studies of cichlid head evolvability.
Results and discussion
Morphological differences between the heads of D. compressiceps, C. azureus and their hybrid progeny
Head size and shape differences
Differences in head length (HL), head height (HH) and head width (HW) are expressed as a percentage of trunk length (TL), which is a better normalizing factor than standard length (SL) here as it excludes HL (Fig. 1). Both species differ significantly in HL and HW (Tukey’s test of honest significant difference (HSD) p < 0.001; Fig. 1a, b, d, e; Table 1), where the longer (63 % TL), narrower DC head (14 % TL) contrasts with the shorter (45 % TL), wider head of CA (20 % TL). The HL of the F1 hybrids (49 % TL) is closer to that of CA, which suggests dominance of CA alleles for this trait, while intermediate F1 HW values (17 % TL) suggest codominance of parental alleles (Fig. 1c–e). In contrast, HH does not differ significantly between CA (45 % TL), DC (47 % TL) and F1 hybrids (44 % TL; Fig. 1f; Table 1).
In order to determine which aspects of skull anatomy contribute to the observed morphological differences, we quantified shape variation between CA and DC using a morphometric analysis where axes of significant divergence were identified by canonical variate analysis (CVA). Eighteen anatomical landmarks were used on lateral views (Fig. 2b, d) and six landmarks on dorsal views (Fig. 2h, j). Lateral landmarks were largely based on  and selected to represent major aspects of the functional morphology of fish feeding. The CVA returned two axes of significant morphological difference between CA and DC for the lateral landmark data (Axis 1: Wilk’s λ = 0.0004, χ 2 = 193.4773, df = 64, p = 6.17057e−15; Axis 2: Wilk’s λ = 0.0490, χ 2 = 73.8914, df = 31, p = 2.31768e−05), which also distinguished each of them from their F1 hybrid progeny, with F1 lateral shape being intermediate between CA and DC (Fig. 2f).
Because of the similar HH between CA, DC and their F1 hybrids (Fig. 1f), landmarks best approximating head height (landmarks 11, 14; Fig. 1d) were set as Bookstein landmarks (asterisks) in order to visualize more precisely the deformation described by CV1 on all other landmarks (Fig. 2e). A clear pattern of anterior deformation is spatially restricted to all landmarks located anterior to the eye in the CA configuration (Fig. 2e, landmarks 1–5, 9, 12, 15–18), indicating a marked trend in anterior expansion of the preorbital region in DC in comparison with CA.
Dorsal head shape differences between CA and DC were also discriminated by CVA, with a single statistically significant axis (Wilk’s λ = 0.0311, χ 2 = 123.2334, df = 16, p < 2.22045e−16). Here, the antero-posterior (A-P) component of CV1 deformation vectors on preorbital landmarks 1, 2 and 6 of the CA configuration is indicative of an anteriorly elongated preorbital region in DC, while the lateral component of the CV1 deformation vectors on landmarks 2, 3, 5 and 6 indicates a narrower head in DC. As with the lateral measures of morphology (Fig. 2f), dorsal F1 shape was intermediate between CA and DC (Fig. 2l).
The results of these analyses confirm that the strong cranial differences in shape between CA and DC are largely restricted to differences in the size of the preorbital region. The skull shapes of CA and DC are therefore widely separated by the same shape axis that describes the primary component of cichlid cranial evolution that arose convergently among both the mbuna and sand-dwellers of Lake Malawi and among the cichlid radiations in lakes Victoria and Tanganyika as well.
To identify the basis of head shape differences, neurocrania of CA, DC and F1 hybrid adult specimens were dissected and examined in their lateral (Fig. 3), dorsal and ventral aspects (Fig. 4). Canonical variate analysis of lateral landmarks 9–12 and 16 found one axis of significant divergence (Wilk’s λ = 0.0371, χ 2 = 72.4799, df = 6, p = 1.26617e−13). Vectors reflecting CV1 variation around the CA landmark configuration supported a general difference in neurocranium aspect ratio (Fig. 3b), with a more pronounced deformation of preorbital landmarks 9, 12 and 16. This more detailed morphometric analysis reveals that expansion of the preorbital region in DC results from A-P elongation and reorientation of the vomer (VO) and lateral ethmoid (LE) in combination with a pronounced positional shift of the mesethmoid (ME; arrows and asterisks in Fig. 3h–j).
In lateral view, with the landmark configurations from all specimens backtransformed for a SL of 10 cm, preorbital length is 91.64 % longer in DC as measured by the distance between landmarks 9 and 16 (Fig. 3d; Table 1). The F1 values for this measurement are closer to those of DC than of CA. Directly posterior to the preorbital bones, above and below the eye, the frontal (FR) and parasphenoid (PS) bones also differ markedly in size along the A–P axis between CA and DC, and there is a pronounced difference in PS orientation as well (Fig. 3e, f, h, i).
In dorsal view, with the landmark configurations from all specimens backtransformed for a SL of 10 cm, the neurocranium is 46.25 % longer in DC (as measured between landmarks 1 and 4) and 67.47 % wider in CA (as measured between landmarks 3 and 5; Fig. 4a–c; Table 1). Dorsal and ventral measurements revealed that shape differences in the following bones contributed to the greater neurocranial length of DC: VO, LE, ME, FR and PS (Fig. 4d–o). Most of these bones are located in the preorbital region (VO, LE, ME; asterisks in Fig. 4g–i, m–o). In contrast, differences in neurocranial width were most pronounced in the postorbital region in DC, with noticeable medio-lateral (M-L) narrowing of bones ventral to the otic cartilage: sphenotic (SO) and pterotic (PO; Fig. 4d–o).
Dorsal and ventral measurements of F1 neurocrania suggest that DC alleles influencing preorbital and neurocranial length tend to be dominant over those of CA (Fig. 4b), while alleles of both species have roughly codominant influence on neurocranial width (Fig. 4c).
The two parental species display strong differences in jaw size and shape (Figs. 2, 5, 6 ). The upper jaws are composed of two bones: the tooth-bearing premaxilla (PMX) and the maxilla (MX; Fig. 5a). Canonical variate analysis of lateral landmarks 5 and 15–18 identified one significantly different CV axis (Wilk’s λ = 0.1530, χ 2 = 70.3949, df = 12, p = 2.7017e−10) (Fig. 5b, c). CV1 deformation vectors revealed three upper jaw shape differences, all of which were associated with PMX anatomy: (1) length of the ascending arm; (2) length of the dentigerous (tooth-bearing) arm; and (3) the angle between these two processes (anatomy sensu Barel et al. 1976) . PMX ascending arm length (distance between lateral landmarks 15 and 18) was 6.75 % longer in DC, but this difference was not significant (Tukey’s HSD; Fig. 5g; Table 1). PMX dentigerous arm length (as estimated by the distance between lateral landmarks 17 and 18) was significantly greater in DC (12.37 %; Tukey’s HSD; Fig. 5h; Table 1). A 10.4° reduction in the angle between these processes in DC was also significant (Tukey’s HSD; Fig. 5i; Table 1). Interestingly, F1 hybrid PMX process lengths and angles were very similar to those of DC (Fig. 5f–i; Table 1), which indicates a clear dominance of DC alleles over CA alleles in determining PMX morphology. Conversely, F1 tooth coverage of the PMX suggests dominance of CA alleles. Teeth cover the entire A-P extent of the PMX dentigerous arm in CA and the F1 progeny, but only cover roughly the anterior two-thirds of its length in DC (Fig. 5d–f).
The lower jaws are composed of two major bones: the tooth-bearing dentary (DE) and the articular (AR; Fig. 6a). Canonical variate analysis of lateral landmarks 1–5 identified a single axis of significant difference between CA and DC (Wilk’s λ = 0.0516, χ 2 = 111.1926, df = 12, p < 2.22045e−16). CV1 deformations of CA landmarks revealed three lower jaw shape differences: (1) lower jaw length (distance between landmarks 2 and 3); (2) DE thickness (distance between landmarks 1 and 2); and (3) AR primordial process length (distance between landmarks 4 and 5). There were significant differences in lower jaw length (27.4 % longer in DC) and DE depth (201.8 % thicker in DC), but not in AR length (only 0.35 % longer in DC; Fig. 6; Tukey’s HDS; Table 1). Lower jaw shape and size in F1 hybrids indicate rough allelic codominance for lower jaw length and dominance of CA alleles for DE depth.
Suspensorium and opercular series
The suspensorium articulates with the upper and lower jaws anteriorly, with the neurocranium posterio-dorsally and with the opercular series posteriorly. It contains the oral jaw biting muscles (the adductor mandibulae; see Fig. 1) and serves important functional roles in both feeding and ventilation. The suspensorium of DC is much elongated in comparison with CA (Figs. 1, 7). Canonical variate analysis of lateral landmarks 4–8, 13, 14 and 16 identified a single axis of significant difference between CA and DC (Wilk’s λ = 0.1168, χ 2 = 40.7961, df = 12 p = 5.29996e−05) (Fig. 7g). CV1 deformations of CA landmarks were strongest for landmarks 4, 6 and 7 (Fig. 7b), which indicates differences in the size of the quadrate (QA), metapterygoid (MP), symplectic (SY) and hyomandibular (HM) bones. There was a significant difference in posterior suspensorium depth (distance between landmarks 4 and 6) between CA and DC (26.3 % longer in DC; Tukey’s HSD; Fig. 7i; Table 1), but no significant difference in anterior suspensorium depth (distance between lateral landmarks 4 and 16) between CA and DC (only 7.50 % longer in DC; Tukey’s HSD; Fig. 7h; Table 1). F1 hybrid suspensorium depth was very similar to that of DC (Fig. 7i), which suggests a dominance of DC alleles.
Each suspensorium connects ventrally to a ceratohyal complex, which bears the ventral-most bones of the opercular series: the branchiostegal rays (Fig. 8). The ceratohyal complex was not included in our whole-head morphometric analysis due to its internal location, but its length was measured and found to be significantly longer (12.1 %; Tukey’s HSD; Table 1) in DC than in CA (Fig. 8d, e). F1 hybrid ceratohyal complex length recapitulated that of CA, indicating dominance of CA alleles on this character (Fig. 8e).
Gill skeleton and pharyngeal jaws
The gill skeleton determines the shape and size of the pharyngeal cavity, supports the gills, filters food particles and processes swallowed food. The overall differences in head size, length, width and diet between CA and DC led us to predict that the piscivorous DC would have a larger gill skeleton that would allow it to process large prey. Here again, we did not include the gill skeleton in our morphometric analysis due to its internal location. Instead, we measured the length of the first ceratobranchial bone (CB1) as a proxy for all CB elements (Fig. 9a–c) and found it to be 30.8 % longer in DC than in CA (Fig. 9d). Planktivorous fish also tend to have smaller, more tightly spaced gill rakers than piscivores, which function to capture small food particles before they exit the gills. CA individuals have more gill rakers on CB1 (14 versus 11.4; Fig. 9e; Table 1), which is shorter than in DC, effectively resulting in a significant difference in raker spacing between the two species (0.7 mm versus 1.2 mm; Tukey’s HSD; Fig. 9f; Table 1). F1 hybrid CB1 length, raker number (13.2) and raker spacing (0.8 mm) are more similar to CA (Fig. 9d–f; Table 1), which suggests dominance of CA alleles over DC alleles.
Cichlid fishes have dorsal and ventral pharyngeal tooth plates that are used for processing food and transporting it into the gut. The lower tooth plate is ossified over the CB5 element and fused along its midline (Fig. 9h–j). There were significant differences in lower tooth plate surface area, the angle formed by lower tooth plate arms and tooth spacing (Fig. 9g, k, l; Tukey’s HSD; Table 1): DC has 18 % smaller (surface area), narrower plates (19° smaller angle) with more widely spaced teeth (253 um in DC versus 187 um in CA) than CA, and F1 hybrid phenotypes suggest codominance of DC and CA alleles on lower tooth plate size and shape. However, F1 hybrid tooth spacing is transgressive—141 um (outside the range established by CA and DC), indicative of synergistic activities between CA and DC alleles. Similar differences in skeletal size and tooth density were observed in upper tooth plates (UTPs; Fig. 9m–o)—DC has 36.4 % smaller UTP3 + 4 elements with more widely spaced teeth (279 um in DC versus 247 um in CA; Fig. 9p, q; Table 1). In this case, CA alleles seem dominant over DC alleles in their effects on UTP3 + 4 size, while synergistic effects of DC and CA alleles on tooth spacing are suggested by a transgressive F1 hybrid phenotype (172 um).
Endochondral growth zones and dermal bone sutures
Skeletal shape and size differences must arise due to differential cartilage/bone growth and remodeling. To gain insight into species-specific differences in growth, we examined sutures, which are well-characterized regions of growth in dermal bones, as well as cartilage growth zones of endochondral bones [13, 19, 41]. We compared the shape and size of sutures and endochondral growth zones in skeletal elements of the neurocranium and suspensorium that differ between CA and DC: the VO, LE, ME and PS (Fig. 10a–d, g–j); the FR, PA and SOC (Fig. 10e, f); the PA, EO, PO and SO (Fig. 10k, l); the QA and MP (Fig. 10m–p); the ventral and dorsal HH and CH (Fig. 10q–s); and the CH and EH (Fig. 10t, u). Qualitative observations of suture organization indicated local differences between neurocranial preorbital bones, where suture patterns between the VO, ME and LE were less intricately folded in DC than in CA (Fig. 10c, d, i, j). These differences were not observed in neurocranial vault sutures (Fig. 10e, f).
Three presumptive zones of endochondral growth were observed in our skeletal preparations (Fig. 10i–l, n–u). Based on assumed organizational similarities with long-bone growth plates, we speculate that zone 1 contains resting and proliferative chondrocytes, as indicated by alcian blue staining. Zone 2 contains mature and hypertrophic chondrocytes ensheathed in bone matrix, as indicated by alizarin red staining associated with cellular cartilage. Zone 3 is made of bone and does not contain chondrocytes. Zones 1 and 2 were consistently larger in DC than in CA, at the expense of zone 3, in all elements examined—QA and MP (Fig. 10o, p) as well as CH and EH (Fig. 10t, u). In some instances, the distinction between cartilage- and dermal bone growth zones was blurred, as at the junctions between the dorsal and ventral HH and the CH (Fig. 10r, s), where endochondral growth zones were partly replaced by sutures, and this replacement was more advanced in CA than in DC. In the neurocranium, zones 1 and 2 were also observed at the junction between the ethmoid cartilage (EC) and the ME, as well as between the otic cartilage (OC) and the PA, EO and PO bones, where differences in zone 1–zone 2 junction shape were observed (Fig. 10i–l).
Developmental basis of divergent head skeleton morphologies
In order to identify the developmental events leading to the formation of the skeletal element shape and size differences observed between the heads of CA and DC, we fixed and stained specimens of each species at developmental stages chosen to demonstrate differences in embryonic cartilage patterning (5 dpf—1CFRE), larval cartilage growth (8 dpf—3CFRE) and/or bone development (15 dpf—9CFRE—11 mm) .
The spatially clustered morphological differences between the adult neurocrania of CA and DC in the preorbital and otic regions suggest that ontogenic differences affecting the cartilage precursors from which adult bones are derived may be responsible for these differences. For instance, differences in ethmoid cartilage (EC) embryonic patterning are predicted to affect the shape and size of all bones that derive from or form in close association with this cartilage. Bones derived from the EC are the ME and the LEs, while the VO is a dermal bone that differentiates directly ventral to and partially surrounds the EC (Fig. 11m–r). Bones derived from the embryonic otic capsule are the EO, SO, PO and EOC (Fig. 11k–p).
The EC first differentiates as a plate in between and extending anterior to the rod-like trabeculae (Fig. 11a–f, k, l). Measurements of ethmoid–trabecular size revealed slight differences in length and width at 8 dpf between CA and DC, but these were not significant (Fig. 11e–h, k, l; Tukey’s HSD; Table 2). Differences in neurocranial width also became noticeable by 8 dpf, but in an unexpected manner, as the head was significantly wider in DC (11.1 % wider at the epiphyseal bar and 5.4 % wider posteriorly; Tukey’s HSD; Table 2) than that in CA at this stage (Fig. 11i, j).
By 15 dpf the EC was 17.8 % longer and 11.3 % narrower in DC (Fig. 11s, t; Tukey’s HSD; Table 2), and the LE and ME had started to ossify (Fig. 11m–r; w–x). The VO was also ossified ventral to the EC. Differences in neurocranial width at this stage started to reflect the adult trend, with the CA neurocranium 4 % wider than that of DC (Fig. 11u, v). At this stage, the frontal bone makes up most of the neurocranial roof and develops as a membrane bone anterior and posterior to the epiphyseal bar cartilage (Fig. 11m–r). No difference in FR bone length was found at this stage (data not shown). The PS was also ossified, but showed no difference in size (data not shown, Fig. 11m, n, q, r). Five bones associated with the otic cartilage were ossified by 15 dpf: the PA, EO, SO, PO and EOC, but no differences were observed between these bones in CA versus DC (Fig. 11m, n).
Overall, a difference in EC length is detectable by 8 dpf and amplified by 15 dpf, supporting the hypothesis that this initial difference may underlie the evolutionary divergence of ME, LE and VO adult preorbital morphologies in CA and DC. In contrast, differences in otic cartilage morphology were not detected at the stages examined and must arise later as the bones grow and remodel.
Interspecific developmental differences of first and second arch derivatives
As with the neurocranium, the localized morphological differences observed between the adult jaw and suspensorial skeletons of CA and DC suggest that ontogenic differences in a few cartilage precursors may underlie shape and size differences in all of their derivatives. Namely, the concerted differences in QA, MP, SY and HM bone sizes between DC and CA may be initiated by size differences in their precursor cartilage: the palatoquadrate (PQ) and the hyosymplectic (HS; Fig. 12a, b), which are derived from the first and second pharyngeal arch dorsal modules, respectively. Strikingly, PQ size (surface area) was already 16.5 % larger in DC than in CA at 5 dpf, 40.3 % larger at 8 dpf and 47 % larger at 15 dpf (Fig. 12a–f, L; Tukey’s HSD; Table 2). HS was also larger in DC, although to a lesser extent—12.1 % at 5 dpf, 18.5 % at 8 dpf and 13.2 % at 15 dpf (Fig. 12a–f, n; Tukey’s HSD; Table 2). In contrast, Meckel’s (MC) and the ceratohyal (CH) were only 5.8 % and 9.7 % larger in DC at 15 dpf (Fig. 12a–f, k, m; Tukey’s HSD; Table 2). These results suggest that specific difference(s) in PQ and HS cartilage development underlie, at least in part, the differences in size and shape of their bone derivatives in adult DC and CA.
Unlike cartilage-derived bones, groups of dermal bones are not derived from common anlage, although they may be influenced by common developmental signals, such as endothelin-1 on the opercular series . The lower jaw skeleton represents an interesting case in which there is a cartilaginous lower jaw (Meckel’s cartilage) in the embryos and young larvae of teleosts, but where the two bones of the adult lower jaw (DE and retroarticular) do not result from endochondral ossification of this structure. They instead develop entirely or almost entirely from dermal ossifications that form around Meckel’s cartilage (there may be some limited endochondral ossification around the articular–quadrate joint), which is retained within the resulting bony envelope (Fig. 12e, f). Cranial dermal bones associated with the jaws stained with alcian/alizarin (Fig. 12e–j) revealed a subset of the differences at 15 dpf noted in upper and lower jaw dermal bones at adult stages: The PMX angle was 6.1° smaller and the DE already had 70.3 % greater depth in DC (Fig. 12o, s; Tukey’s HSD; Table 2). PMX dentigerous and ascending arms were 10.3 % and 11.56 % longer, respectively, in DC (Fig. 12p, q; Tukey’s HSD; Table 2), while DE length was not significantly different.
As exemplified by differences in PQ and HS formation, shape and size differences affecting cartilage bones in the adult pharyngeal skeleton originate in cartilage development during embryonic stages and thus result in morphologically integrated changes among all derived bones. This replicates the pattern seen with endochondral bone development in the preorbital region of the neurocranium, but the same is not true for dermal bone formation. Developmental changes that affect the shape and size of dermal ossifications have much more spatially restricted effects on the pharyngeal skeleton (e.g., the DE in the lower jaw).
Interspecific developmental differences of gill arch derivatives
Interspecific differences in adult gill arch morphology include: (1) CB length; (2) gill raker size and spacing; (3) pharyngeal plate size; and (4) tooth size and spacing. No significant differences in CB length were found at 5, 8 or 15 dpf (Fig. 13d, Tukey’s HSD; Table 2), while CB1 was 19.2 % larger (surface area) in DC at 15 dpf (Fig. 13e; Tukey’s HSD; Table 2), reflecting a difference in thickness. Gill raker buds were observed by 8 dpf in both species (Fig. 13f1, f2) and were already larger and more widely spaced in DC at this stage. By 15 dpf CA had 14 rakers on CB1, while DC had 10 (Fig. 13g1, g2, h), spaced 43.2 um and 63.8 um apart, respectively (Fig. 13i).
Similar to CB development, embryonic development of upper (UTP) and lower pharyngeal tooth plates (LTP) was remarkably similar between CA and DC (Fig. 13j1–k2, n1–o2) in terms of both shape and size, followed by a significant increase in size in DC by 15dpf. UTP and LTP plate sizes were larger in DC (24 % and 14 %, respectively) by 15 dpf (Fig. 13l1–m, p1–q; Tukey’s HSD; Table 2), which did not correlate with adult pharyngeal tooth plate sizes. In addition, the angle between the lateral arms of the LTP was 4.5° larger in DC (Fig. 13p1, p2, r), which suggests that widening of the pharyngeal cavity precedes its narrowing in DC. A similar number of equally spaced tooth buds were observed on the UTP (Fig. 13j1–l2) and LTP (Fig. 13n1–p2) at 5, 8 and 15 dpf in both species, which indicates that differences in tooth number and spacing appear at later stages.
Overall, our data indicate that unlike differences in neurocranial and anterior pharyngeal arch skeletal elements, shape and size differences in the 15 dpf gill arch skeleton do not correlate with adult morphology; DC larvae have thicker CB elements and larger tooth plates than CA larvae. While this morphology may better adapt DC juveniles for the capture of larger prey items, a shift in gill arch developmental trajectory must occur later in life to adopt the adult morphology. In contrast, gill raker number and spacing do correlate between adult and embryonic/larval stages.
Growth zone development
Because of the central role played by skeletal growth zones in post-embryonic skeletal development, we investigated the timing of their appearance during embryonic stages. The first cellular rearrangements indicative of growth zones were observed at 8 dpf in pharyngeal cartilages and were most obvious in larger, anterior elements: PQ, HS and CH (Fig. 14). By 8 dpf fields of chondrocytes flattened along semicircular patterns within these elements (Fig. 14d–g, l–o). This stereotypical pattern of cartilage development is characteristic of both growth plates in mammalian long bones and craniofacial cartilage development in larval zebrafish [41, 44]. The spatial layout of these arrays corresponded to those of future growth zones (Fig. 14a, b, h–k, p–s). At this stage, two obvious regions of coordinated chondrocyte flattening were observed in PQ, four in HS and three in the CH of both CA and DC (Fig. 14a). Chondrocyte hypertrophy and ossification were later observed (15 dpf) on one side of each of these arrays, as revealed by alizarin red-positive staining (Fig. 14a, b, h–k, p–s). Quantification of growth zone width in PQ, HS and CH revealed significant differences in width between CA and DC: Growth zone width was 29.0 % larger in PQ of DC than of CA and 9.3 % larger in HS of DC than of CA, while it was 5.9 % smaller in CH of DC (Fig. 14c; Tukey’s HSD; Table 2). While the significance of growth zone size remains to be determined, a reasonable interpretation is that wider growth zones may produce faster growth via a larger number of proliferative chondrocytes. In addition, the wider PQ and HS growth zones of DC correlate with the larger size of their derivative adult elements in DC.
Phenotypic integration and modular organization in the cichlid skull
Two cranial ossification mechanisms appear to produce separate patterns of skeletal element integration that underlie morphological differences in the crania of C. azureus and D. compressiceps. Endochondral ossification tends to produce integrated changes in multiple skeletal elements derived from, or whose development is intimately associated with, single precursor cartilages (Fig. 15a). The VO, ME and LE bones, for example, undergo coordinated shape and positional changes that directly correspond to differences in the development of the ethmoid cartilage. In contrast, dermal bones are prone to region-specific changes within individual skeletal elements (Fig. 15b). The PMX, for example, develops via dermal ossification, and the length of its dentigerous process varies independently of that of its ascending arm.
This opposition in modular development—division of a module into smaller ones (cartilage bone) or fusion of small modules into a larger one (dermal bone), provides distinct properties of phenotypic integration likely to influence the evolvability of skeletal morphology (Fig. 15a, b). In the head, the strategic placement of cartilage versus dermal bones provides distinct opportunities to evolve novel shapes. We speculate that cartilage bones give the opportunity to achieve integrated changes in the morphology of all their derivatives, which may be useful to create changes in overall head shape rapidly, such as length, height and/or width, but constrain possibilities for bone-specific shape evolution (Fig. 15a). In contrast, developmental changes in any of the distinct ossification centers that form individual dermal bones have the potential to produce novel individual bone shapes that may not be integrated with surrounding elements, but which may be useful where subtle changes in individual bone shape influence function, such as in jaw elements (Fig. 15b).
Phenotypic integration derived from the modular organization of complex structures has been a defining aspect of the developmental module concept since its formulation , and the mouse lower jaw skeleton is an established model to study the relationships between modularity and phenotypic integration [46–48]. However, the contrasting modular architecture of cartilage and dermal bones highlighted in our study provides a novel perspective on how development may have constrained the rapid diversification of craniofacial morphologies in East African Rift lake cichlids.
Role of skeletal growth zones in cichlid craniofacial evolution
Our findings also suggest that evolutionary changes in skeletal development at embryonic/larval stages become amplified by subsequent growth that generates profound differences in adult skeletal morphology, such as those we observe between DC and CA. A similar conclusion was recently drawn from the comparative analysis of shape ontogeny in six other Lake Malawi cichlids with distinct head morphologies . We speculate that this amplification occurs at spatially restricted growth zones, the properties of which are dictated by the type of bone they produce: dermal or cartilage bone. Cartilage bone growth occurs at endochondral growth zones, as exemplified for the QA and MP bones, and our results suggest that these zones influence skeletal shape and size through both their shape, which is patterned during embryonic and larval development, and their proliferative activity, (Fig. 15c, d). In addition, element-specific differences in cartilage proliferation may produce overall head shape differences. Cell proliferation assays will help determine how growth zones contribute to skeletal shape and size differences in the cichlid head.
In contrast, dermal bone outgrowth occurs through a variety of mechanisms in which deposition patterns of new mineralized matrix by osteoblasts differ. The best-characterized pattern of tetrapod dermal bone growth occurs at the edges of developing cranial bones. The profound shape changes occurring in the DC skull vault beyond embryonic development suggest an important role for differential growth at cranial sutures. Studies of the zebrafish operculum have revealed a similar pattern, where an increase in dermal bone size by osteoblasts leaves incremental growth bands in the bone matrix, while two other modes, “spur” and “veil” formation, produce new shapes with only minor increases in bone size . Many dermal bones differ in shape or size between DC and CA, the most spectacular being the dentary, and future studies will have to determine the location and mode(s) of bone outgrowth for each element.
Candidate developmental pathways underlying cichlid craniofacial evolution
The morphologically diverse cichlid species assemblies populating the Rift Lakes of East Africa have been described as natural mutant panels which can offer insight into the identity of genetic variants underlying adaptive phenotypic changes . This “evolutionary mutant model” approach complements the phenotypic analyses of induced mutations in conventional model organisms to explore phenotype to genotype connections . Although some loci without prior association with craniofacial development have been discovered in cichlids , polymorphisms mapped in cichlids and other natural populations are often associated with developmental pathways previously described in model organisms: Hedgehog, Wnt, Bmp signaling, etc. [26, 27, 52, 53]. The novelty of these natural variants is that they often reveal functionally important cis-regulatory elements [54, 55] and groups of loci that affect phenotypes in a quantitative manner. Furthermore, the phenotypic changes that they induce are functionally integrated into organisms in a way that does not impair their ability to survive. Ultimately, quantitative genetic studies of craniofacial traits in cichlids will fill a large knowledge gap in the structure and function of the gene regulatory networks (GRN) that control skull morphogenesis.
Our descriptive analysis of major trait differences between CA and DC and the ontogeny of these characters points toward several known developmental pathways, some of which act during embryonic patterning of the pharyngeal arches and palate region. For instance, the specific enlargement of cartilage elements derived from the dorsal first and second pharyngeal arches, the palatoquadrate and hyosymplectic, may result from divergence in the dorso-ventral pharyngeal patterning GRNs. As interactions between Edn1, Bmp and Wnt pathways pattern dorsal pharyngeal domains by restricting the ventral extant of jag1b expression [56–58], we speculate that this GRN may underlie expansion of dorsal elements by either ventrally stretching the jag1b expression domain, or stimulating increased proliferation of jag1b-expressing skeletal progenitors. Similarly, differences in the palate developmental program, which is regulated by Shh and Pdgf signaling, may underlie differences in ethmoid cartilage size between DC and CA and the preorbital region expansion of the cichlid neurocranium [59–61].
The F1 progeny examined here is the product of a mapping cross between CA and DC that we are currently using to identify: (1) the loci underlying craniofacial morphological divergence between these two species and (2) the developmental pathways underlying the adaptive radiations of East African Rift Lake Cichlids. This cross holds promise for contributing to our understanding of cichlid skull evolvability since the cranial differences between these species (Figs. 1, 2) are strongly associated with the most important axis of morphological evolution that has arisen in all Rift Lake cichlid radiations: divergence in preorbital size [22, 24]. We are currently examining the F2 progeny of this cross at multiple developmental stages. In combination with our detailed developmental descriptions of both parental species, this should allow us to gain insight into the genetic controls of bone development in multiple cranial regions and their connection to cichlid skull evolvability.
DC and CA brood stocks were purchased from pet stores and maintained separately at 28 °C in 60-gallon tanks. DC offspring were obtained from breeding tanks housing DC adults alone and CA offspring from those containing CA adults alone. Hybridization was achieved by crossing a DC male with a CA female in a 60-gallon tank containing a total of five CA females. Mouth-brooding females were separated from tank mates with a tank partition. Juvenile offspring were collected at 2 weeks post-fertilization and grown in 10-gallon tanks for 2–3 months, after which they were grown for a year in a 60-gallon tank. Juveniles were fed a mixture of live brine shrimp and crushed Hikari Cichlid Gold baby pellets, followed by baby pellets alone, mini pellets and finally medium pellets, depending on specimen age/size. All adults were fed Hikari Cichlid Gold medium pellets.
Forty-three adult specimens were used in this study: 12 DC individuals, 15 CA individuals and 16 F1 hybrids. Adult specimens of mixed sexes were collected at 12 months of age and older (CA: 5.7 cm < SL < 10.3 cm; DC: 9.3 < SL < 15.4; F1: 6.4 cm < SL < 11.0 cm), killed with tricaine as approved by the UCI IACUC and fixed for 3–7 days in 10 % formalin. Specimen heads were then skinned and skeletal staining was performed as described by . Internal skeletal elements (ceratohyal and gill arch complexes) were dissected in five individuals of each species/line. Bone type and abbreviations are indicated in Table 3.
Embryos and larvae were also killed with tricaine as approved by the UCI IACUC, at 5 days post-fertilization (dpf)/1 caudal fin ray elements (CFRE), 8 dpf/3 CFRE and 15 dpf/9 CFRE  and fixed in 4 % paraformaldehyde (PFA) for 3 days. The 15dpf/9 CFRE stage was chosen primarily because it precedes yolk depletion by approximately a day, making growth up to this stage independent of food intake. Specimens were cleared and stained as described in , and five specimens of each species were analyzed for each stage. Adult skeletal nomenclature follows .
Images of whole adult specimens were captured using a Panasonic GH1 camera with a Panasonic Lumix G 20 mm f/1.7 lens. Images of dissected skeletal elements from the neurocrania of adults, larvae and embryos were captured using an Olympus SZX12 dissecting microscope equipped with an Olympus DF PLAPO 1.2 × PF2 objective and an Olympus DP70 camera using DP-Controller software (Olympus). Images of flat-mounted embryos and larvae were captured with a Zeiss Axioplan 2 microscope equipped with Plan-NEOFLUAR objectives, and a Micropublisher 5.0 RTV camera using Volocity software (Improvision). ImageJ was used for all linear and surface area measurements. Landmark positions were also recorded as (x, y) coordinates in ImageJ.
Data transformations, geometric morphometrics and statistics
Individual differences in size, orientation and position were first minimized by performing a Procrustes superimposition of the landmarks of all adult specimens using the program CoordGen8. Procrustes transformations superimpose the landmarks of all specimens as much as possible and scale the landmark panels of each individual to the same centroid size. The software CVAGen8 was then used to perform a canonical variance (CV) analysis to identify canonical variate axes that best discriminated the three groups of landmark configurations: CA, DC and F1. Deformation vectors representing the variation identified by a given CV axis around a mean landmark configuration were also visualized in CVAGen8.
Linear measurements between landmarks in adult specimens or on dissected skeletal elements were plotted against individual standard length (SL) and backtransformed to a SL of 10 cm. Linear measurements in embryos and larvae were not transformed.
Data were graphed using ggplot2 (RStudio), and Tukey’s honest significant difference (HSD) tests were performed using TukeyHSD (RStudio).
gene regulatory network
honest significant difference
lower tooth plate
upper tooth plate
Alberch P. From genes to phenotypes—dynamic systems and evolvability. Genetica. 1991;84:5–11.
Dawkins R. The evolution of evolvability. In: Langton C, editors. First Los Alamos conference on artificial life, Los Alamos, NM. 1988.
Hendrikse JL, Parsons TE, Hallgrimsson B. Evolvability as the proper focus of evolutionary developmental biology. Evol Dev. 2007;9:393–401.
Klingenberg CP. Developmental constraints, modules and evolvability. In: Hallgrimsson B, Hall BK, editors. Variation: a central concept in biology. New York: Elsevier; 2005. p. 219–47.
Hallgrimsson B, Jamniczky H, Young NM, Rolian C, Parsons TE, Boughner JC, Marcucio RS. Deciphering the palimpsest: studying the relationship between morphological integration and phenotypic covariation. Evol Biol. 2009;36:355–76.
Gans C, Northcutt RG. Neural crest and the origin of vertebrates—a new head. Science. 1983;220:268–73.
Depew MJ, Compagnucci C. Tweaking the hinge and caps: testing a model of the organization of jaws. J Exp Zool Part B Mol Dev Evol. 2008;310B:315–35.
Hecht J, Stricker S, Wiecha U, Stiege A, Panopoulou G, Podsiadlowski L, Poustka AJ, Dieterich C, Ehrich S, Suvorova J, et al. Evolution of a core gene network for skeletogenesis in chordates. PLoS Genet. 2008;4:13.
Porto A, de Oliveira FB, Shirai LT, De Conto V, Marroig G. The evolution of modularity in the mammalian skull I: morphological integration patterns and magnitudes. Evol Biol. 2009;36:118–35.
Fish JL, Villmoare B, Koebernick K, Compagnucci C, Britanova O, Tarabykin V, Depew MJ. Satb2, modularity, and the evolvability of the vertebrate jaw. Evol Dev. 2011;13:549–64.
Mallarino R, Grant PR, Grant BR, Herrel A, Kuo WP, Abzhanov A. Two developmental modules establish 3D beak-shape variation in Darwin’s finches. Proc Natl Acad Sci USA. 2011;108:4057–62.
Kimmel CB. Skull developmental modularity: a view from a single bone - or two. J Appl Ichthyol. 2014;30:600–7.
Schilling TF, Le Pabic P. Neural crest cells in craniofacial skeletal development. In: Trainor P, editor. Neural crest cells; evolution, development and disease. Amsterdam: Elsevier; 2014. p. 127–51.
de Beer GR. The development of the vertebrate skull. London: Oxford University Press; 1937.
Rijli FM, Mark M, Lakkaraju S, Dierich A, Dollé P, Chambon P. A homeotic transformation is generated in the rostral branchial region of the head by disruption of Hoxa-2, which acts as a selector gene. Cell. 1993;75:1333–49.
Gendron-Maguire M, Mallo M, Zhang M, Gridley T. Hoxa-2 mutant mice exhibit homeotic transformation of skeletal elements derived from cranial neural crest. Cell. 1993;75:1317–31.
Le Pabic P, Scemama J-L, Stellwag EJ. Role of Hox PG2 genes in Nile tilapia pharyngeal arch specification: implications for gnathostome pharyngeal arch evolution. Evol Dev. 2010;12:45–60.
Clouthier DE, Garcia E, Schilling TF. Regulation of facial morphogenesis by endothelin signaling: insights from mice and fish. Am J Med Genet Part A. 2010;152A:2962–73.
Hall BK. Bones and cartilage; developmental and evolutionary skeletal biology. 2nd ed. Amsterdam: Elsevier; 2015.
Sturmbauer C. Explosive speciation in cichlid fishes of the African Great Lakes: a dynamic model of adaptive radiation. J Fish Biol. 1998;53:18–36.
Sturmbauer C, Baric S, Salzburger W, Ruber L, Verheyen E. Lake level fluctuations synchronize genetic divergences of cichlid fishes in African lakes. Mol Biol Evol. 2001;18:144–54.
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:A38–50.
Salzburger W, Mack T, Verheyen E, Meyer A. Out of Tanganyika: genesis, explosive speciation, key-innovations and phylogeography of the haplochromine cichlid fishes. BMC Evol Biol. 2005;5:15.
Young KA, Snoeks J, Seehausen O. Morphological diversity and the roles of contingency, chance and determinism in african cichlid radiations. PLoS ONE. 2009;4:8.
Hu YA, Albertson RC. Hedgehog signaling mediates adaptive variation in a dynamic functional system in the cichlid feeding apparatus. Proc Natl Acad Sci USA. 2014;111:8530–4.
Albertson RC, Streelman JT, Kocher TD, Yelick PC. Integration and evolution of the cichlid mandible: The molecular basis of alternate feeding strategies. Proc Natl Acad Sci USA. 2005;102:16287–92.
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:3629.
Powder KE, Cousin H, McLinden GP, Albertson RC. A nonsynonymous mutation in the transcriptional regulator lbh is associated with cichlid craniofacial adaptation and neural crest cell development. Mol Biol Evol. 2014;31:3113–24.
Powder KE, Milch K, Asselin G, Albertson RC. Constraint and diversification of developmental trajectories in cichlid facial morphologies. Evodevo. 2015;6:25.
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 USA. 1999;96:5107–10.
Parsons KJ, Marquez E, Albertson RC. Constraint and opportunity: the genetic basis and evolution of modularity in the cichlid mandible. Am Nat. 2012;179:64–78.
Hu YN, 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.
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.
Albertson RC, Cooper WJ, Mann KA. More than meets the eye: functionally salient changes in internal bone architecture accompany divergence in cichlid feeding mode. Int J Evol Biol. 2012;2012:538146.
Albertson RC, Kocher TD. Genetic and developmental basis of cichlid trophic diversity. Heredity. 2006;97:211–21.
Albertson RC, Streelman JT, Kocher TD. Genetic basis of adaptive shape differences in the cichlid head. J Hered. 2003;94:291–301.
Stauffer JR, Konings AF. Review of Copadichromis (Teleostei:Cichlidae) with the description of a new genus and six new species. Ichthyol Explor Freshw. 2006;17:9–42.
Fryer G, Iles TD. The cichlid fishes of the great lakes of Africa. Edinburgh: Oliver & Boyd Press; 1972.
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:17.
Barel CDN, Witte F, Van Oijen MJP. The shape of the skeletal elements in the head of a generalized Haplochromis species: Haplochromis elegans (Pisces, Cichlidae) with two examples of trophically correlated shape differences. Neth J Zool. 1976;26:163–265.
Kozhemyakina E, Lassar AB, Zelzer E. A pathway to bone: signaling molecules and transcription factors involved in chondrocyte development and maturation. Development. 2015;142:817–31.
Fujimura K, Okada N. Development of the embryo, larva and early juvenile of Nile tilapia Oreochromis niloticus (Pisces:Cichlidae). Developmental staging system. Dev Growth Differ. 2007;49:301–24.
Kimmel CB, Walker MB, Miller CT. Morphing the hyomandibular skeleton in development and evolution. J Exp Zool Part B Mol Dev Evol. 2007;308B:609–24.
LeClair EE, Mui SR, Huang A, Topczewska JM, Topczewski J. Craniofacial skeletal defects of adult zebrafish glypican 4 (knypek) mutants. Dev Dyn. 2009;238:2550–63.
Raff RA. The shape of life: genes, development, and the evolution of animal form. 1st ed. Chicago: University of Chicago Press; 1996. xxiii + 520 p.
Atchley WR, Hall BK. A model for development and evolution of complex morphological structures. Biol Rev Camb Philos Soc. 1991;66:101–57.
Klingenberg CP, Leamy LJ, Routman EJ, Cheverud JM. Genetic architecture of mandible shape in mice: effects of quantitative trait loci analyzed by geometric morphometrics. Genetics. 2001;157:785–802.
Klingenberg CP, Leamy LJ, Cheverud JM. Integration and modularity of quantitative trait locus effects on geometric shape in the mouse mandible. Genetics. 2004;166:1909–21.
Kimmel CB, DeLaurier A, Ullmann B, Dowd J, McFadden M. Modes of developmental outgrowth and shaping of a craniofacial bone in zebrafish. PLoS ONE. 2010;5:12.
Kocher TD. Adaptive evolution and explosive speciation: the cichlid fish model. Nat Rev Genet. 2004;5:288–98.
Albertson RC, Cresko W, Detrich HW III, Postlethwait JH. Evolutionary mutant models for human disease. Trends Genet. 2009;25:74–81.
Roberts RB, Hu Y, Albertson RC, Kocher TD. Craniofacial divergence and ongoing adaptation via the hedgehog pathway. Proc Natl Acad Sci USA. 2011;108:13194–9.
Cleves PA, Ellis NA, Jimenez MT, Nunez SM, Schluter D, Kingsley DM, Miller CT. Evolved tooth gain in sticklebacks is associated with a cis-regulatory allele of Bmp6 (vol 111, pg 13912, 2014). Proc Natl Acad Sci USA. 2014;111:18090.
Erickson PA, Cleves PA, Ellis NA, Schwalbach KT, Hart JC, Miller CT. A 190 base pair, TGF-beta responsive tooth and fin enhancer is required for stickleback Bmp6 expression. Developmental Biology. 2015;401:310–23.
Carroll SB. Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution. Cell. 2008;134:25–36.
Zuniga E, Stellabotte F, Crump JG. Jagged-Notch signaling ensures dorsal skeletal identity in the vertebrate face. Development. 2010;137:1843–52.
Alexander C, Zuniga E, Blitz IL, Wada N, Le Pabic P, Javidan Y, Zhang T, Cho KW, Crump JG, Schilling TF. Combinatorial roles for BMPs and endothelin 1 in patterning the dorsal-ventral axis of the craniofacial skeleton. Development. 2011;138:5135–46.
Alexander C, Piloto S, Le Pabic P, Schilling TF. Wnt signaling interacts with Bmp and Edn1 to regulate dorsal-ventral patterning and growth of the craniofacial skeleton. PLoS Genet. 2014;10:7.
Wada N, Javidan Y, Nelson S, Carney TJ, Kelsh RN, Schilling TF. Hedgehog signaling is required for cranial neural crest morphogenesis and chondrogenesis at the midline in the zebrafish skull. Development. 2005;132:3977–88.
Eberhart JK, Swartz ME, Crump JG, Kimmel CB. Early Hedgehog signaling from neural to oral epithelium organizes anterior craniofacial development. Development. 2006;133:1069–77.
Eberhart JK, He X, Swartz ME, Yan Y-L, Song H, Boling TC, Kunerth AK, Walker MB, Kimmel CB, Postlethwait JH. MicroRNA Mirn140 modulates Pdgf signaling during palatogenesis. Nat Genet. 2008;40:290–8.
Datovo A, Bockmann FA. Dorsolateral head muscles of the catfish families Nematogenyidae and Trichomycteridae (Siluriformes: Loricarioidei): comparative anatomy and phylogenetic analysis. Neotrop Ichthyol. 2010;8:193–246.
Walker MB, Kimmel CB. A two-color acid-free cartilage and bone stain for zebrafish larvae. Biotech Histochem. 2007;82:23–8.
Liem KF. Comparative functional anatomy of the Nandidae (Pisces: Teleostei). Fieldiana Zoology. 1970;56:1–164.
PL conceived, designed and performed the experiments, data quantifications and statistical analyses and drafted the manuscript. WJC helped to draft the manuscript. TFS participated in the design of the study and helped to draft the manuscript. All authors read and approved the final manuscript.
We thank Chris Rackauckas for his support with data analysis in RStudio. This work was supported by NIH/NIDCR F32DE21308 to PL as well as R01DE13828 and R01HD073182 to TFS.
The authors declare that they have no competing interests.