Variable paralog expression underlies phenotype variation

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

Human faces are variable; we look different from one another. Craniofacial disorders further increase facial variation. To understand craniofacial variation and how it can be buffered, we analyzed the zebrafish mef2ca mutant. When this transcription factor encoding gene is mutated, zebrafish develop dramatically variable craniofacial phenotypes. Years of selective breeding for low and high penetrance of mutant phenotypes produced strains that are either resilient or sensitive to the mef2ca mutation. Here, we compared gene expression between these strains, which revealed that selective breeding enriched for high and low mef2ca paralog expression in the low- and high-penetrance strains, respectively. We found that mef2ca paralog expression is variable in unselected wild-type zebrafish, motivating the hypothesis that heritable variation in paralog expression underlies mutant phenotype severity and variation. In support, mutagenizing the mef2ca paralogs, mef2aa, mef2b, mef2cb, and mef2d demonstrated modular buffering by paralogs. Specifically, some paralogs buffer severity while others buffer variability. We present a novel, mechanistic model for phenotypic variation where variable, vestigial paralog expression buffers development. These studies are a major step forward in understanding the mechanisms of facial variation, including how some genetically resilient individuals can overcome a deleterious mutation.

Editor's evaluation

In this elegant genetic study, Bailon-Zambrano et al. draw on classical genetic concepts to address the clinically pertinent question of how genetic variants in the same gene can yield wildly different phenotypes in different individuals. Specifically, this work makes significant contributions to our understanding of how a mutant phenotype may be modified by the expression levels of paralogues of the mutant gene.

https://doi.org/10.7554/eLife.79247.sa0

Introduction

Human faces are variable

Human craniofacial variation allows us – and even computers – to identify specific individuals. A recent study indicated that human facial structures are more variable than other human anatomical features and that evolution selected for increased variation in craniofacial-associated genomic regions (Sheehan and Nachman, 2014). While human genomic variation contributes to craniofacial phenotypic variation (Liu et al., 2012; Paternoster et al., 2012), facial variation persists among more genetically homogeneous populations such as Finnish people (Sheehan and Nachman, 2014) and even monozygotic twins (Kohn, 1991). Therefore, even among individuals with similar genomes, craniofacial development may remain noisy or sensitive to minor fluctuations, compared with other developmental processes.

Phenotype severity and variation can be measured in different ways. Penetrance is the frequency of a phenotype associated with a genotype. Expressivity is different degrees of the same phenotype associated with a genotype. Here, we use both penetrance and average expressivity to measure phenotype severity. We also use the distribution in expressivity to measure phenotypic variation. Furthermore, there are two types of variation. Among-individual variation can be quantified by measuring differences between individuals. Within-individual variation is measured by quantifying departures from symmetry on left versus right sides of an individual known as fluctuating asymmetry (Valen, 1962). It is unknown if among- and within-individual variation are products of the same biological mechanisms (Hallgrimsson et al., 2019). There is evidence to support that these two types of variation are associated (Breuker et al., 2006; Santos et al., 2005); however, other studies indicate that they are independent (Debat et al., 2009; Breno et al., 2011). Further work is needed to resolve this question.

Variability in human craniofacial disease

In the 1940s, Waddington observed that mutant organisms are often associated with increased phenotypic variation compared to wild types (Waddington, 1942), and more recent studies support this finding (Fish, 2016; Hallgrímsson et al., 2006; Scharloo, 1991; Parsons et al., 2008). Human genetic craniofacial disease phenotypes also appear more variable than normal human facial phenotypes (Ward et al., 2000). For example, facial clefting among monozygotic twins can be incompletely penetrant or unilateral (Trotman et al., 1993; Young et al., 2021), demonstrating both among- and within-individual variation. How the same genetic disease allele can have devastating consequences in some individuals, while other resilient individuals can overcome the deleterious mutation is not well understood (Chen et al., 2016).

One human genetic disease that presents variable craniofacial phenotypes is MEF2C haploinsufficiency syndrome. Patients heterozygous for mutations affecting the transcription factor encoding gene MEF2C show variable facial dysmorphologies (Le Meur et al., 2010; Zweier et al., 2010; Gordon et al., 2018; Tonk et al., 2011). To our knowledge, no MEF2C homozygous mutant patients have been identified, likely due to embryonic lethality. Additionally, craniofacial asymmetry is documented in some patients with this disorder (Nowakowska et al., 2010; Novara et al., 2013). Thus, both among- and within-individual variation are present in MEF2C haploinsufficiency syndrome patients. Numerous MEF2C mutant alleles cause this disorder. Therefore, a model system in which the same allele can be studied in many different individuals is needed.

Variable buffering by gene family members might underlie among-individual variation

During evolution, whole genome duplications produced multiple mef2 genes, or paralogs, in vertebrate genomes. Thus far, no mef2 paralogs are associated with craniofacial development besides mef2c. When whole genome duplications occur, the most likely outcome is the loss of one of the duplicate genes through the accumulation of deleterious mutations and eventual nonfunctionalization (Takahata and Maruyama, 1979; Watterson, 1983). However, sometimes mutations occur in gene regulatory elements partitioning gene expression domains among the new duplicates. This preserves both copies and is called subfunctionalization (Force et al., 1999; Ohno, 2013). Subfunctionalized duplicates may experience only partial loss of expression subdomains (Force et al., 1999), retaining relics of their ancestral expression pattern even if the gene is no longer required for the original function.

Compensation, or the ability of paralogs or gene family members to make up for the loss of a gene due to an overlap in function, has long been recognized as a source of genetic robustness (Wagner, 2000; Kirschner and Gerhart, 1998; Gu et al., 2003; Conant and Wagner, 2004; Kamath et al., 2003). Robustness is defined as the ability of a biological system to overcome genetic or environmental perturbation (Diss et al., 2014; de Visser et al., 2003). Buffering lessens the impact of perturbation on a system. Although buffering by paralogs has long been proposed, whether paralogs retain buffering capacity after subfunctionalization has not been sufficiently addressed.

We hypothesize that vestigial expression can buffer against loss of another paralog and that variation in vestigial expression underlies phenotypic variation. This hypothesis is supported by work in systems ranging from yeast to human cells demonstrating that paralogs contribute to genetic robustness (Diss et al., 2014; De Kegel and Ryan, 2019; Dandage and Landry, 2019; Li et al., 2010; Hsiao and Vitkup, 2008; Wagner, 2008; Dean et al., 2008). However, whether paralog expression variation underlies phenotypic variation has not been directly tested. Finally, whether variation in paralogous compensation is subject to selection is unknown.

We previously found evidence to support this hypothesis in zebrafish. mef2ca single mutants produce dramatically variable craniofacial phenotypes (Nichols et al., 2016). mef2cb is the most closely related mef2ca paralog. However, mef2cb is not required for craniofacial development, and homozygous mutants are viable and indistinguishable from wild types (Hinits et al., 2012). However, we reported that ventral cartilage defects associated with mef2ca mutation become more severe when we removed a single functional copy of mef2cb from mef2ca homozygous mutants (De Laurier et al., 2014). Our vestigial buffering hypothesis predicts that although mef2cb is no longer overtly required for craniofacial development, remaining mef2cb expression in neural crest cells might partially substitute for mef2ca loss. There is further evidence from mouse mutants that Mef2 paralogs can functionally substitute for one another (Majidi et al., 2019). Because mutations in none of the other four zebrafish paralogous mef2 genes have been reported, their function is unknown. These paralogs have human orthologs, MEF2A, MEF2B, and MEF2D. None of these genes have been associated with human craniofacial development or disease.

We developed a system for understanding craniofacial development and variability

One prominent phenotype in zebrafish mef2ca homozygous mutants is expansion of the opercle, a bone supporting the gill flap. This phenotype is remarkably variable (De Laurier et al., 2014; Kimmel et al., 2003); mef2ca homozygous mutants have many opercle shapes, and some even develop wild-type-looking opercle bones; the phenotype is incompletely penetrant (Nichols et al., 2016). With some genes, mutants encoding a premature termination codon (PTC) upregulate compensating genes through transcriptional adaptation, explaining why some mutant alleles do not produce a phenotype (Rossi et al., 2015; El-Brolosy et al., 2019). Transcriptional adaptation does not contribute to incomplete penetrance in our system (Sucharov et al., 2019). Specifically, we pharmacologically inhibited transcriptional adaptation (by inhibiting RNA decay) and did not observe changes in penetrance.

We demonstrated that the factors underlying opercle variation are heritable through selective breeding, which shifted the penetrance of the expanded bone phenotype to generate strains with consistently low or high penetrance of this phenotype (Nichols et al., 2016). Unlike traditional mutagenesis modifier screens (St Johnston, 2002; Kile and Hilton, 2005), our selective breeding paradigm likely enriched for standing genetic variation without the need for further mutagenesis. Thus, we likely bred naturally occurring genetic modifiers in our background to fixation. Similar approaches have proven successful in other systems (Vu et al., 2015; Gasch et al., 2016). mef2ca is pleiotropic, and although we only selected on opercle bone penetrance, nearly all mef2ca-associated phenotypes increased or decreased penetrance in the high- and low-penetrance strains, respectively (Sucharov et al., 2019). The only phenotype remaining fully (100%) penetrant in the low- and high-penetrance strains is a shortened symplectic cartilage, a linear, rod-shaped cartilage that functions as a jaw support structure. In unselected lines, mef2ca-associated phenotypes are only found in mef2ca homozygous mutants (Miller et al., 2007). However, in high-penetrance heterozygotes, we observed the shortened symplectic cartilage, suggesting this phenotype is highly sensitive to mef2ca loss (Sucharov et al., 2019). Penetrance is a binary measurement (a shortened symplectic is present or not), and we have not yet examined expressivity (to what extent a mutant symplectic is shortened). We do not know if the average length of the shortened symplectic or the distribution of symplectic length is different between low- and high-penetrance strains.

Here, we capitalize on the strengths of our zebrafish system to address the mechanisms underlying craniofacial variation. We quantified expressivity by taking linear measurements of the symplectic cartilage, allowing us to measure severity, among-individual variation, and within-individual variation in the zebrafish craniofacial skeleton. Combining these measurements with penetrance scoring, we compared phenotype severity and variation between selectively bred strains. We found that increased severity is associated with increased variation. In the high-penetrance strain, we found a more severe symplectic phenotype, and more among-individual variation but not within-individual variation compared to the low-penetrance strain.

What factors were selected upon from the original genetic background to increase or decrease penetrance? Work in nematodes demonstrates that expression variation in otherwise redundant genes contributes to variable penetrance in mutants (Raj et al., 2010). Therefore, we hypothesized that one or more of the five zebrafish mef2ca paralogs might be redundant with mef2ca and ameliorate the phenotype in the low-penetrance line or reciprocally intensify the phenotype in the high-penetrance line. We compared the expression of the mef2 paralogs between selectively bred strains and found that many of the paralogs are more highly expressed in the low-penetrance strain compared to the high-penetrance strain. Furthermore, we uncovered standing variation in paralog expression in unselected strains. To determine if the increased paralog expression we discovered in the low-penetrance strain buffers the mef2ca mutant phenotype, we mutagenized the mef2 paralogs. Double mutant analyses indicate that the different paralogs modularly buffer different aspects of the pleiotropic mef2ca mutant phenotype: some affect severity, while others affect variation and some affect both. These findings demonstrate that heritable, variable paralog expression is a major factor affecting phenotype severity and among-individual phenotypic variation but does not contribute to within-individual variation.

Results

Craniofacial phenotypes are more variable in mef2ca mutants compared with wild types

The phenotypic variation often associated with human genetic diseases might be partially due to individuals inheriting different mutant alleles retaining different levels of functional activity, expression domains, and/or transcriptional levels. A strength of the zebrafish system is that variation can be studied in many different individuals with the same mutant allele. Wild-type mef2ca encodes an N-terminal MADS box (MCM1, agamous, deficiens, and SRF) and an adjacent MEF2 domain (Chen et al., 2017; Figure 1A). These highly conserved domains mediate dimerization, DNA binding, and co-factor interactions (Potthoff and Olson, 2007). The mef2ca^b1086 mutant allele arose from a forward genetic screen and produces a PTC just downstream of the MADS box. This mutant allele variably displays phenotypes including ectopic bone near the opercle, interhyal joint fusion, dysmorphic ceratohyal, reduced Meckel’s cartilage, jaw joint fusion, and shortened symplectic (Figure 1B and C). Skeletal preparations from two full-sibling individuals illustrate both among- and within-individual variation easily observed in mef2ca mutants (Figure 1B). Comparing the two mutant individuals (B3 and B4) demonstrates among-individual variation associated with the opercle bone phenotype; one individual (B3) has phenotypically wild-type opercles, while the other individual has bilateral mutant phenotypes (B4). Meanwhile, within-individual (left-right) symplectic cartilage phenotype variation is present in one of the mef2ca^b1086 animals (B4); in this animal, the symplectic cartilage is longer on one side than the other.

Figure 1 with 1 supplement see all

Download asset Open asset

*mef2ca* mutant craniofacial phenotypes are more variable than wild types.

(A) Schematic of *mef2ca* exonic structure. The *mef2ca^b1086* mutant allele used in this study and regions encoding known functional domains are annotated. (**B1–B4**) Zebrafish heterozygous for *mef2ca* were pairwise intercrossed, and 6 days post fertilization (dpf), larvae were stained with Alcian blue and Alizarin red to label cartilage and bone. The individuals were then genotyped, flat mounted, and imaged. Two examples (upper and lower) are provided for each genotype. The following craniofacial skeletal elements are indicated in a wild-type individual: opercle bone (op), branchiostegal ray (br), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Indicated phenotypes associated with *mef2ca* mutants include: ectopic bone (arrowheads), interhyal and jaw-joint fusions (^), reduced mc (double arrowhead), and a shortened sy (red arrows). Dashed outline indicates symplectic cartilage. Scale bar: 50 μm. (C) The penetrance of *mef2ca* mutant-associated phenotypes observed in 6 dpf homozygous mutant larvae is indicated. (D) Schematic indicating how the symplectic cartilage length was measured in this study. (E) Symplectic cartilage length was measured from 6 dpf wild type or homozygous mutant larvae. The p-value from a Welch’s t-test is indicated (****≤0.0001). F-test value testing for significant differences in variation between genotypes is indicated. (F) Symplectic cartilage length on left and right sides of 6 dpf zebrafish was measured to determine fluctuating asymmetry, or the absolute difference between left and right, for wild type or mutant larvae. The p-value from a Welch’s t-test is indicated (****≤0.0001). For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. For E and F, n=22 for wild types and 44 for mutants.

Figure 1—source data 1 Penetrance data and raw symplectic cartilage length measurements.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig1-data1-v2.xlsx
Download elife-79247-fig1-data1-v2.xlsx

To quantify the variation associated with mef2ca mutants, we first scored penetrance of the various mef2ca phenotypes in an unselected AB background (Figure 1C). Based on penetrance, the symplectic cartilage is the craniofacial structure most sensitive to mef2ca loss (Sucharov et al., 2019). Moreover, we found that symplectic length correlates with opercle bone area (Figure 1—figure supplement 1), suggesting the shortened symplectic phenotype is a good proxy for the entire craniofacial complex. We developed a quantitative phenotyping assay where the symplectic length can be rapidly measured in many animals (Figure 1D). We used the total linear symplectic length (left side plus right side) to examine expressivity (Figure 1E). The average symplectic length was significantly shorter in mutants compared to wild types. The among-individual variation was significantly greater in mutants compared to wild types. When we examined within-individual variation by determining the absolute value of the difference between left- and right-symplectic cartilage length, we found that mutants have more within-individual variation than the wild type (Figure 1F).

Selective breeding reveals that severity and among-individual variation, but not within-individual variation, is heritable and segregates with penetrance

Another strength of the zebrafish system is that penetrance can be altered by selective breeding to better understand the influence of genetic background on severity and variation. An ongoing, long-term selective-breeding experiment in our laboratory derived two strains of zebrafish with consistently low and high penetrance of mef2ca-associated phenotypes in mef2ca^b1086 homozygous mutants (Nichols et al., 2016; Brooks and Nichols, 2017; Figure 2A). Our previous work examined the phenotypes present in homozygous mutant fish from these strains (Sucharov et al., 2019). Here, we closely examined mef2ca wild types (mef2ca^+/+) from these strains. We were surprised to observe that occasionally some mef2ca-associated phenotypes, like shortened symplectic cartilages (19%, n=16) and fused interhyal joints (13%, n=16), were observed in mef2ca^+/+ individuals from the high-penetrance strain (Figure 2B). In contrast, we never observed these phenotypes in mef2ca^+/+ individuals from the low-penetrance strain. Importantly, these phenotypes are not likely due to general developmental delay because other skeletal structures like pharyngeal teeth are unaffected. Thus, the phenotypes we discovered in high-penetrance mef2ca^+/+ individuals are specific to developmental processes associated with mef2ca function. Quantifying the severity of the shortened symplectic cartilage phenotype, we found significant differences between low- and high-penetrance strains for both mef2ca^+/+ as well as mef2ca^+/-, but homozygous mutants were statistically the same by this measure in the two strains (Figure 2C).

Figure 2

Download asset Open asset

Selective breeding affects *mef2ca*-associated phenotype severity and variation in *mef2ca* wild types, heterozygotes, and homozygous mutants.

(A) Selective breeding pedigree illustrating ectopic bone phenotype penetrance inheritance. Dashed boxes indicate families used in this study. Six generations of full-sibling inbreeding produced the animals used here, a more complete pedigree extending back >10 generations can be found in our previous publication (Sucharov et al., 2019). (B) Alcian blue- and Alizarin red-stained animals from the low- and high-penetrance strains were genotyped, and *mef2ca* homozygous wild types were flat mounted and imaged. The following craniofacial skeletal elements are indicated in a wild-type individual from the low-penetrance strain: opercle bone (op), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Phenotypes normally associated with *mef2ca* homozygous mutants are present in some wild types from the low-penetrance strain including: ih joint fusions (^) and shortened sy (red arrows). Bars indicating sy length are presented to illustrate the shortened symplectic phenotype present in some high-penetrance wild types but not low-penetrance wild types. A stage-appropriate complement of ankylosed of pharyngeal teeth (pt) are present, and normal sized op bones are present in the individual with shortened sy, indicating the phenotypes we discovered in high-penetrance *mef2ca^+/+* are not due to general delay. Dashed outline indicates symplectic cartilage. Scale bar: 50 μm (C) Symplectic cartilage length was measured from 6 days post fertilization (dpf) larvae from wild types, heterozygotes, and homozygous mutants from both the low- and high-penetrance strains. p-Values from a Dunnet’s T3 test are indicated (**≤0.01). (D) The coefficient of variation for symplectic length in all three genotypes from both strains was plotted (E) Table listing F-test values testing for significant differences in variation between strains comparing the same genotype. (F) Symplectic cartilage length on left and right sides of 6 dpf zebrafish was measured to determine fluctuating asymmetry or the absolute difference between left and right for all three genotypes from both strains. (G) Table listing the Dunnett’s test for significant differences in fluctuating asymmetry between all three genotypes from both strains. For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values.

Figure 2—source data 1 Symplectic cartilage length measurements from selectively bred strains.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig2-data1-v2.xlsx
Download elife-79247-fig2-data1-v2.xlsx

When we examined variation in symplectic cartilage length, we compared like genotypes between strains and observed that among-individual variation was greater in the high-penetrance strain for homozygous wild types and heterozygotes, and there was no significant difference in variation between strains for homozygous mutants (Figure 2D and E). Thus, the genetic background producing more severe phenotypes also exhibited more among-individual variation, even in wild types. Within-individual variation was not significantly different between backgrounds (Figure 2F and G).

Six mef2 paralogs in the zebrafish genome share highly conserved amino acid sequences

We sought to explore potential mechanisms underlying the differences in severity and variation between strains. Work from budding yeast demonstrates that gene duplication contributes to genetic robustness against null mutations and that the probability of compensation by gene duplicates is correlated with sequence similarity among the gene duplicates (Gu et al., 2003). Therefore, we hypothesized that mef2ca duplicates with high sequence similarity would modulate severity and variation. There are five mef2ca paralogs in the zebrafish genome. To determine which are the most similar to mef2ca and therefore more likely to compensate for its loss, we used Clustal Omega to compare the amino acid sequences of mef2 paralogs. While the C-terminal domain is divergent among different mef2 genes (Molkentin et al., 1996; Figure 3—figure supplement 1), the zebrafish mef2 paralogs each encode a MADS box and MEF2 domain which are remarkably similar in sequence across all paralogs. We found that the paralog with the most sequence similarity to mef2ca is mef2cb (Figure 3A and B). Together, these genes are the co-orthologs of mammalian Mef2c. mef2aa and mef2ab, co-orthologs of mammalian Mef2a, also share similar sequences to each other and the mef2c pair. By these analyses, mef2d and mef2b sequences are more divergent. There is only one zebrafish ortholog for Mef2d and Mef2b, presumably the other co-ortholog for each was lost. These results and interpretations are consistent with previous analyses of the evolutionary history of the mef2 family (Chen et al., 2017; Wu et al., 2011a; Martindale et al., 2004; He et al., 2019; Wu et al., 2011b). Thus, these duplicates are strong candidates for mechanistically underlying the differences we observe between strains.

Figure 3 with 1 supplement see all

Download asset Open asset

Zebrafish *mef2* paralogs encode highly conserved MADS box (MCM1, agamous, deficiens, and SRF) and MEF2 domains.

(A) Neighbor joining tree generated by Clustal Omega Multiple Sequence Alignment tool depicts the evolutionary relationships between the different zebrafish *mef2* paralogs. The distance values (branch length) are indicated, which represent the evolutionary distance between the individual amino acid sequences and a consensus sequence. (B) *mef2*-encoded protein sequence alignment reveals high conservation of MADS box (green) and MEF2 (yellow) domains among all six paralogs. These domains are responsible for DNA binding, dimerization, and cofactor interactions. Transcript IDs used for alignment using the HHalign algorithm are listed in the Materials and methods.

Closely related mef2 paralogs share similar temporospatial expression dynamics

Our previous work, and that of others, demonstrates when (20–30 hpf) (hours post fertilization) and where (cranial neural crest cells) mef2ca functions during craniofacial development (De Laurier et al., 2014; Sucharov et al., 2019; Miller et al., 2007). We therefore hypothesized that mef2ca duplicates with temporospatial expression dynamics similar to mef2ca would modulate severity and variation. To examine paralog temporal expression dynamics, we performed RT-quantitative PCR (qPCR) (quantitative polymerase chain reaction) on wild-type zebrafish heads from 20 to 30 hpf (Figure 4A). We found that closely related paralogs share similar gene expression dynamics. For example, mef2ca and mef2cb have both early (~24 hpf) and late (~28 hpf) expression peaks. mef2aa and mef2ab share an early expression peak. mef2d and mef2b have a single late and early expression peak, respectively.

Figure 4 with 1 supplement see all

Download asset Open asset

Wild-type gene expression studies reveal *mef2* paralog expression dynamics, strain-specific expression levels, and standing paralog expression variation.

(A) Wild-type head expression of all *mef2* paralogs was quantified by RT-quantitative PCR (qPCR) at 1 hr intervals from 20 to 30 hpf. Expression of each paralog was normalized to *rps18*. Error bars are SD. (B) Expression of each paralog in cranial neural crest cells at 24 hpf was determined by single cell RNA-sequencing on sorted cells. Seurat-based clustering subdivided the cells into four populations as described (Mitchell et al., 2021): anterior arches (aa), frontonasal (fn), posterior arches (pa), and a satellite population containing melanocyte lineage cells (m). (C) Wild-type head expression of all *mef2* paralogs was quantified by RT-qPCR at 24 hpf to compare paralog expression levels between the low- and high-penetrance strains. Expression of each paralog was normalized to *rps18*. Asterisks indicate significant difference (*≤0.05 and **≤0.01). Error bars are SD. (D) Wild-type head expression of all *mef2* paralogs was quantified by RT-qPCR at 24 hpf to compare paralog expression levels between two families from the unselected AB strain. Expression of each paralog was normalized to *rps18*. Asterisks indicate significant difference (*≤0.05, ***≤0.001, and ****≤0.0001). Error bars are SEM.

Figure 4—source data 1 qPCR raw values.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig4-data1-v2.xlsx
Download elife-79247-fig4-data1-v2.xlsx
Figure 4—source data 2 qPCR raw values.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig4-data2-v2.xlsx
Download elife-79247-fig4-data2-v2.xlsx
Figure 4—source data 3 qPCR raw values.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig4-data3-v2.xlsx
Download elife-79247-fig4-data3-v2.xlsx
Figure 4—source data 4 qPCR raw values.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig4-data4-v2.xlsx
Download elife-79247-fig4-data4-v2.xlsx

To examine paralog expression in the cells that give rise to the craniofacial skeleton, we examined our previously published single cell RNA-sequencing dataset to monitor mef2 paralog expression in isolated 24 hpf wild-type cranial neural crest cells (Mitchell et al., 2021). While this method does not capture all the cells orchestrating craniofacial development, it does allow us to specifically assay cranial neural crest cells, including those residing in the anterior arches, which are the precursors of the cells which form the skeletal elements affected by mef2ca loss. We found that mef2ca has the strongest expression of all the paralogs across different populations of cranial neural crest cells and that expression is strongest in the anterior arches (Figure 4B). The closely related paralog mef2cb is the next highest expressed, while other paralogs are more weakly expressed in these craniofacial progenitor cells at this stage.

mef2 paralogs are differentially expressed in the low- versus high-penetrance strains, and paralog expression varies in an unselected background

Data showing that paralog transcriptional adaptation does not account for the phenotypic differences between the low- and high-penetrance strains (Sucharov et al., 2019) do not rule out the possibility that selective breeding changed paralog expression between strains in mef2ca wild types. To examine this possibility, we used RT-qPCR to compare mef2 paralog expression between wild types from the high- and low-penetrance strains. Strikingly, we found that mef2aa, mef2ab, mef2b, mef2ca, and mef2d were all significantly differentially expressed in heads from wild types from the high-penetrance strain compared with heads from wild types from the low-penetrance strain (Figure 4C). These findings strongly suggest that one outcome of selective breeding for low- and high-mef2ca phenotype penetrance is generally increased and decreased expression, respectively, of the mef2 paralogs. We do not see overall increases in transcription in the low-penetrance strain compared with the high-penetrance strain; housekeeping genes are not significantly upregulated in the low-penetrance strain (Figure 4—figure supplement 1).

The differences in paralog expression between selectively bred strains suggest that we might be selecting upon standing variation in paralog expression present in the unselected strain. To test this, we examined paralog gene expression in heads from offspring of two separate families of unselected AB wild types. We found significant differences in mef2ab, mef2b, mef2cb, and mef2d expression between these families. Of note, the direction of paralog expression differences between families was consistent across paralogs. That is, family 1 had generally higher mef2 paralog expression compared with family 2. These data support that paralog expression variation is present in unselected lines and that our selective breeding for penetrance enriched for these expression differences.

mef2cb buffers against mef2ca loss

We hypothesized that the paralog expression differences that we discovered between strains underlie the differences in severity and variation between the low- and high-penetrance strains. To test this hypothesis, we studied mef2ca mutants in the context of mef2 paralog mutations. For mef2cb, we crossed a previously generated allele (De Laurier et al., 2014) into the mef2ca low-penetrance strain. We then maintained this stock by outcrossing to unselected AB. mef2cb is the paralog with the highest degree of sequence similarity to mef2ca and is the second highest expressed mef2 paralog in cranial neural crest cells behind mef2ca (Figures 3 and 4). We confirmed that mef2cb homozygous mutants do not have an overt craniofacial phenotype and are homozygous viable (De Laurier et al., 2014; Figure 5A, B and C), indicating that mef2cb is not required for zebrafish craniofacial development in an otherwise wild-type background. However, when we removed one functional copy of mef2ca from mef2cb homozygous mutants, mef2ca mutant-associated phenotypes developed (Figure 5B and C), phenocopying the high-penetrance strain (Sucharov et al., 2019). We find further evidence for phenocopy; mef2ca homozygous mutant phenotype penetrance increased when we removed one copy of mef2cb (Figure 5B and D). Removing both copies of mef2cb from mef2ca homozygous mutants produces severe, nonspecific defects which make larvae impossible to meaningfully study (De Laurier et al., 2014) although their craniofacial skeletons are severely affected (Figure 5—figure supplement 1). Measuring the symplectic cartilage length further demonstrates that when mef2cb is fully functional, development is buffered against partial loss of mef2ca; there is no difference in symplectic length between wild types and mef2ca heterozygotes (wild type versus mef2ca^+/-;mef2cb^+/+). In contrast, when mef2cb is disabled, development is sensitive to partial loss of mef2ca (wild type versus mef2ca^+/-;mef2cb^-/-) (Figure 5E). Removing copies of mef2cb from mef2ca wild types does not significantly change symplectic cartilage length but does significantly increase symplectic cartilage variation (Figure 5F). Thus, even in the mef2ca wild-type context, this paralog buffers against phenotypic variation. We conclude that mef2cb buffers against mef2ca-associated phenotype severity, and among-individual variation, but not within-individual variation (Figure 5D, E and Figure 5—figure supplement 2). However, our gene expression study in high- and low-penetrance strains indicated that several mef2 paralogs are differentially expressed between strains. Therefore, we examined how other paralogs might also buffer against mef2ca loss.

Figure 5 with 2 supplements see all

Download asset Open asset

*mef2cb* function buffers against *mef2ca* loss.

(A) Schematic of *mef2cb* exonic structure, mutant allele used in this study, and regions encoding proposed functional domains are annotated. (B) Zebrafish heterozygous for both *mef2ca^b1086* and *mef2cb* were pairwise intercrossed. 6 days post fertilization (dpf) larvae were stained with Alcian blue and Alizarin red to label cartilage and bone. Stained larvae were genotyped, flat mounted, and imaged. The following craniofacial skeletal elements are indicated in a wild-type individual: opercle bone (op), branchiostegal ray (br), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Indicated phenotypes associated with *mef2ca* mutants include: ectopic bone (arrowheads), interhyal and jaw-joint fusions (^), dysmorphic ch (arrows), reduced mc (double arrowhead), and a shortened sy (red arrows). Dashed outline indicates symplectic cartilage. Scale bar: 50 μm (C and D) The penetrance of *mef2ca* mutant-associated phenotypes observed in 6 dpf larvae is indicated. Asterisk indicates significant difference in penetrance between the indicated genotypes by Fishers exact test. (E) Symplectic cartilage length was measured from 6 dpf larvae from the indicated genotypes. Asterisks indicate significant differences in symplectic length. The p-values from a Dunnet’s T3 test are indicated (*≤0.05, ***≤0.001, and ****≤0.0001) (F) Table listing F-test values for significant differences in variation between genotypes. For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. N’s for all analyses are indicated in C and D.

Figure 5—source data 1 Symplectic length cartilage measurements.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig5-data1-v2.xlsx
Download elife-79247-fig5-data1-v2.xlsx

mef2d buffers against mef2ca homozygous mutant severity but not variability

We detected mef2d transcripts in the anterior arch population of 24 hpf wild-type cranial neural crest cells, and mef2d expression is significantly lower in the high-penetrance strain compared with the low-penetrance strain (Figure 4C and D). To explore the function of this paralog, we generated a mef2d mutant allele (Figure 6A). For mef2d and all subsequent paralog functional experiments (mef2b and mef2aa), we generated new mutant alleles using (Clustered regularly interspaced short palindromic repeats) CRISPR/Cas9 mutagenesis in the low-penetrance strain. We then maintained these stocks by outcrossing to unselected AB. Homozygous mef2d mutants did not develop any overt skeletal phenotypes in an otherwise wild type, unselected background (Figure 6B). Intercrosses between heterozygotes produced homozygous mutant adults at the expected Mendelian frequency, indicating that mef2d is not required for craniofacial development or viability in laboratory conditions. Similarly, Mef2d mutant mice are viable and display no overt phenotypic abnormalities in standard laboratory conditions (Arnold et al., 2007; Kim et al., 2008). To test our hypothesis that reduced levels of mef2d expression associated with selective breeding for high mef2ca penetrance contribute to severity, we examined offspring from mef2ca;mef2d double heterozygous parents (Figure 6B). When we examined the penetrance of all mef2ca mutant-associated skeletal phenotypes in mef2ca homozygous mutants, we found that the frequency of ventral cartilage defects was significantly increased by removing a single functional copy of mef2d and further increased in double homozygous mutants (Figure 6C). The penetrance of other mef2ca-associated phenotypes was not significantly changed when mef2d function was removed. While symplectic length was not affected in mef2d single mutants, we found that removing both copies of mef2d from mef2ca homozygous mutants significantly shortened symplectic cartilage length. We did not observe significant changes in among-individual variation when we removed mef2d function from mef2ca homozygotes, but symplectic cartilage length is more variable in mef2ca heterozygotes when mef2d is mutated (Figure 6E).

Figure 6

Download asset Open asset

*mef2d* function buffers against *mef2ca* loss.

(A) Schematic of *mef2d* exonic structure, mutant allele used in this study, and regions encoding proposed functional domains are annotated. (B) Zebrafish heterozygous for both *mef2ca ^b1086* and *mef2d* were pairwise intercrossed. 6 days post fertilization (dpf) larvae were stained with Alcian blue and Alizarin red to label cartilage and bone. Stained larvae were genotyped, flat mounted, and imaged. The following craniofacial skeletal elements are indicated in a wild-type individual: opercle bone (op), branchiostegal ray (br), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Indicated phenotypes associated with *mef2ca* mutants include: ectopic bone (arrowheads), interhyal and jaw-joint fusions (^), dysmorphic ch (arrows), reduced mc (double arrowhead), and a shortened sy (red arrows). Dashed outline indicates symplectic cartilage. Scale bar: 50 μm (C) The penetrance of *mef2ca* mutant-associated phenotypes observed in 6 dpf larvae is indicated. Asterisk denotes significant difference in penetrance between the indicated genotypes by Fishers exact test. (D) Symplectic cartilage length was measured from 6 dpf larvae from the indicated genotypes, and asterisk indicates significant differences in symplectic length (*≤0.05, **≤0.01, and ***≤0.001). (E) Table listing F-test values testing for significant differences in variation between genotypes. For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. N’s for all analyses are indicated in C.

Figure 6—source data 1 Symplectic cartilage length measurements.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig6-data1-v2.xlsx
Download elife-79247-fig6-data1-v2.xlsx

mef2b buffers against mef2ca homozygous mutant variation but not severity

mef2b is the most divergent mef2ca paralog by amino acid sequence and is only minimally expressed in cranial neural crest cells at 24 hpf (Figures 3 and 4). However, we did observe significantly higher mef2b expression in the low-penetrance strain compared with the high-penetrance strain (Figure 4D). We generated a mef2b mutant allele to test for a role for this gene in craniofacial development (Figure 7A). Homozygous mef2b mutants did not exhibit any overt skeletal phenotypes (Figure 7B) and intercrosses between animals heterozygous for this allele produced homozygous mutant adults at the expected Mendelian frequency. When we removed functional copies of mef2b from mef2ca homozygous mutants, we did neither observe any significant changes in mef2ca mutant-associated phenotype penetrance (Figure 7C) nor did we detect further reductions in symplectic cartilage length when functional copies of mef2b were removed from mef2ca homozygous mutants (Figure 7D). These findings indicate that mef2b function does not affect mef2ca mutant severity. However, we do observe a modest increase in among-individual variation in mef2ca homozygous mutants when mef2b is disabled (Figure 7E).

Figure 7

Download asset Open asset

*mef2b* function buffers against *mef2ca* loss.

(A) Schematic of *mef2b* exonic structure, mutant allele used in this study, and regions encoding proposed functional domains are annotated. (B) Zebrafish heterozygous for both *mef2ca ^b1086* and *mef2b* were pairwise intercrossed. 6 days post fertilization (dpf) larvae were stained with Alcian blue and Alizarin red to label cartilage and bone. Stained larvae were genotyped, flat mounted, and imaged. The following craniofacial skeletal elements are indicated in a wild-type individual: opercle bone (op), branchiostegal ray (br), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Indicated phenotypes associated with *mef2ca* mutants include: ectopic bone (arrowheads), interhyal and jaw-joint fusions (^), dysmorphic ch (arrows), reduced mc (double arrowhead), and a shortened sy (red arrows). Scale bar: 50 μm. (C) The penetrance of *mef2ca* mutant-associated phenotypes observed in 6 dpf larvae is indicated. (D) Symplectic cartilage length was measured from 6 dpf larvae from the indicated genotypes, and asterisk indicates significant differences in symplectic length (***≤0.001 and ****≤0.0001). (E) Table listing F-test values testing for significant differences in variation between genotypes. For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. N’s for all analyses are indicated in C.

Figure 7—source data 1 Symplectic cartilage length measurements.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig7-data1-v2.xlsx
Download elife-79247-fig7-data1-v2.xlsx

mef2aa buffers against mef2ca heterozygous phenotypes but not mef2ca homozygous mutant phenotypes

mef2aa is minimally expressed in cranial neural crest cells (Figure 4). However, mef2aa expression differs between strains (Figure 4C). We generated a mef2aa mutant (Figure 8A) that does not develop any overt skeletal phenotypes when homozygous (Figure 8B). When we removed mef2aa from mef2ca homozygous mutants, we did not observe any significant changes in mef2ca-associated phenotype penetrance (Figure 8C), symplectic cartilage length (Figure 8D), or variation (Figure 8E). However, when we removed one functional copy of mef2ca from mef2aa homozygous mutants, mef2ca-associated phenotypes developed with low frequency (Figure 8B). Specifically, some mef2ca heterozygous animals developed nubbins (small lumps of cartilage) and shortened symplectic cartilages; phenotypes traditionally only associated with homozygous mutants (Piotrowski et al., 1996). Neither of these phenotypes are ever seen in unselected mef2ca heterozygotes but do develop in mef2ca heterozygotes from the high-penetrance strain (Sucharov et al., 2019). Therefore, disabling mef2aa partially phenocopies the high-penetrance strain.

Figure 8

Download asset Open asset

*mef2aa* function buffers against *mef2ca* partial loss.

(A) Schematic of *mef2aa* exonic structure, mutant allele used in this study, and regions encoding proposed functional domains are annotated. (B) Zebrafish heterozygous for both *mef2ca ^b1086* and *mef2aa* were pairwise intercrossed. 6 days post fertilization (dpf) larvae were stained with Alcian blue and Alizarin red to label cartilage and bone. Stained larvae were genotyped, flat mounted, and imaged. The following craniofacial skeletal elements are indicated in a wild type: opercle bone (op), branchiostegal ray (br), Meckel’s (mc), ceratohyal (ch), symplectic (sy) cartilages, interhyal (ih), and jaw (jj) joints. Indicated phenotypes associated with *mef2ca* mutants include: cartilage nubbin fused to the mc symphysis (brown arrowhead), ectopic bone (black arrowheads), ih and jj fusions (^), dysmorphic ch (arrows), reduced mc (double arrowhead), and a shortened sy (red arrows). Scale bar: 50 μm. (C) The penetrance of *mef2ca* mutant-associated phenotypes observed in 6 dpf larvae is indicated. (D) Symplectic cartilage length was measured from 6 dpf larvae from the indicated genotypes, and asterisk indicates significant differences in symplectic length (****≤0.0001). (E) Table listing F-test values testing for significant differences in variation between genotypes. For box and whisker plots, the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. N’s for all analyses are indicated in C.

Figure 8—source data 1 Symplectic cartilage length measurements.: https://cdn.elifesciences.org/articles/79247/elife-79247-fig8-data1-v2.xlsx
Download elife-79247-fig8-data1-v2.xlsx

Discussion

Vestigial paralog expression may provide developmental robustness

Vertebrate mef2 functions downstream of endothelin signaling in the developing craniofacial skeleton (Miller et al., 2007; Verzi et al., 2007). The endothelin pathway was subfunctionalized following whole genome duplications in vertebrates (Square et al., 2020). Thus, mef2 genes may have been subfunctionalized following genome duplications. Comparing mice and zebrafish further supports mef2 subfunctionalization. In mice, Mef2c is required for both heart and craniofacial development (Verzi et al., 2007; Lin et al., 1997). In zebrafish, both co-orthologs (mef2ca and mef2cb) function redundantly in the heart (Hinits et al., 2012), while craniofacial function has been subfunctionalized to just mef2ca (Figure 9). However, it is possible that the ancestral craniofacial function and expression pattern of mef2cb are partially retained following subfunctionalization, even though it is no longer required for this function. We propose that while duplicated genes can evolve new expression domains and functions, vestiges of their original expression pattern remain and can buffer against loss of another paralog. Consistently, with the exception of mef2ca, the mef2 paralog mutants we analyzed here do not have an overt single mutant skeletal phenotype but do modify mef2ca mutant phenotypes. In our system, selective breeding likely fixed existing paralog expression variation. Thus, in the low-penetrance strain, selection may replicate paralog ancestral expression, restoring near full redundancy (Figure 9).

Figure 9

Download asset Open asset

Model for *mef2* gene duplicate evolution and replication of ancestral, redundant expression via selective breeding.

In our model, an ancestral *mef2* gene existed with distinct expression domains. In this example, ancestral *mef2c* had expression in the head and the heart. Following whole genome duplication, initially the regulatory and coding sequences would be identical. Through divergence and natural selection, factors controlling gene regulation would acquire mutations that dampen, but do not eliminate, some expression domains. In this example, *mef2cb* expression is dampened in the head. These gene expression changes result in subfunctionalization but with vestigial retention of the original expression. *mef2cb* is no longer required for craniofacial development, yet traces of the original craniofacial expression remain. The amount of vestigial expression is variable among individuals, resulting in variable buffering against *mef2ca* mutant phenotypes across a population. Our selective breeding selected on this variable expression resulting in higher paralog expression in the low-penetrance strain compared with the high-penetrance strain. Thus, *mef2* paralog expression in the low-penetrance strain resembles the ancestral condition following whole genome duplication, when both copies’ expression profiles were highly similar, and the genes were redundant for craniofacial development.

Our experiments support a model where we selected upon pre-existing alleles that either amplify or dampen existing vestigial paralog expression variation, depending on the direction of selection. We found that expression of several paralogs was affected by selection. Previously however, we reported that the rapid response to selection in our system indicates that relatively few heritable factors contribute to penetrance (Nichols et al., 2016; Sucharov et al., 2019). To account for this paradox, we propose that a few shared genetic factors regulate paralog expression. By selecting for penetrance, expression of several paralogs was increased or decreased via selection for these shared genetic factors. In support, we found evidence of shared gene regulation in our experiments examining mef2 gene expression dynamics. The expression dynamics of paralogs closely related in sequence and phylogeny are similar, suggesting shared regulatory modules. In this model, few factors are inherited that control the expression of many paralogs. These mef2 paralog regulatory factors might be either shared mef2 enhancers (in cis) or might be in shared upstream transcriptional regulators (in trans) that simultaneously regulate several mef2 paralogs. Future studies will discern between these two, not mutually exclusive, possibilities.

Evidence for genetic assimilation

It is fascinating that a phenotype originally associated with a homozygous loss of function mutation can appear in homozygous wild types after selective breeding for high penetrance. We previously observed a similar phenomenon with heterozygotes in the high-penetrance strain and proposed that this is a form of genetic assimilation (Sucharov et al., 2019; Waddington, 1953). In Waddington’s genetic assimilation experiments (Waddington, 1942; Waddington, 1953; Waddington, 1959; Waddington, 1952; Waddington, 1956; Waddington, 1957), he studied phenotypes that were originally only present in perturbed conditions. Following selective breeding, the phenotypes arose without the perturbation. Similarly, we studied phenotypes that were originally only present in perturbed conditions, mef2ca mutants. Following selective breeding, the phenotypes arose without the perturbation, in wild types. We offer a mechanistic explanation for our findings. Analyzing craniofacial skeletons in the context of the mef2ca mutation revealed cryptic paralog expression variation. Specifically, in wild types, paralog expression variation does not impact phenotype. However, in mutants, the consequences of variable expression levels are unveiled, allowing us to select upon them. Because the paralogs seem to have shared regulatory systems, we enriched for variants producing lower expression of several paralogs, including mef2ca, by selective breeding for high penetrance. This pan-paralog expression decrease led the high-penetrance strain to be both extremely sensitive to coding mutations in mef2ca and to manifest mef2ca mutant-associated phenotypes in mef2ca wild types from this strain. The latter is likely due to decreased mef2ca expression in high-penetrance mef2ca wild types.

Decoupling buffering mechanisms

We found that paralogs modularly buffer the mef2ca homozygous mutants. For example, disabling mef2cb affects most mef2ca-related craniofacial phenotypes, while mef2d mutations primarily affect penetrance of ventral cartilage phenotypes, and only variation is affected when mef2b function is removed. Although we were unable to mutagenize mef2ab, despite significant effort, this is not likely to affect our general conclusions. Similar to work in other systems (Gu et al., 2003), the paralog with the highest expression (mef2ca) is the one associated with a single mutant phenotype, and the paralog with the most sequence similarity (mef2cb) is the one playing the largest buffering role. However, others have reported that partial gene expression overlap is also predictive of buffering capacity (Kafri et al., 2005; Kafri et al., 2006), in line with our results that even distantly related paralogs with more dissimilar expression profiles (mef2d) can modify penetrance in our system. Different types of compensation events have previously been classified as either active or passive (Diss et al., 2014). In our system, compensation is likely passive, but subject to selection, which can change robustness across generations.

Determining whether the specific portion of the phenotype that each paralog buffers is related to that paralog’s wild-type expression pattern would expand our understanding of modularity. Unfortunately, in situ gene expression protocols were not sensitive enough to detect the low expression of the individual paralogs. Moreover, the quantitative phenotyping in this study was limited to the symplectic cartilage, as a proxy for the whole system (Figure 1—figure supplement 1). However, quantitative phenotyping of other parts of the developing craniofacial complex (Sasaki et al., 2013; Kimmel et al., 2015) might reveal more information about modular buffering of severity and variation by paralogs. For example, penetrance of ceratohyal cartilage defects was affected by mef2d, and therefore quantitative variation of this structure, which would be difficult to measure, might also be affected by loss of mef2d.

Our experiments also decouple the mechanisms that buffer among- and within-individual variation. Examining selectively bred strains and the different paralog mutants demonstrates that buffering in our system only regulates among-individual variation, not within-individual variation (Figure 5—figure supplement 2). These two types of variation are buffered by different mechanisms in our system. Our results shed some light on the long-standing debate surrounding similarities and differences between mechanisms buffering these two types of variation (Hallgrimsson et al., 2019).

These studies advance our understanding of craniofacial variation. We propose that cryptic variation in vestigial paralog expression is the noise underlying variable craniofacial development. Whole genome duplications producing paralogs are not zebrafish specific or specific to this gene family. In fact, previous work demonstrates that paralogs contribute to robustness in diverse systems (Gu et al., 2003; Kamath et al., 2003; De Kegel and Ryan, 2019). Thus, even human craniofacial variation might be due to cryptic paralog expression variation and may explain how some genetically resilient humans (Chen et al., 2016) can overcome a deleterious mutant allele.

Materials and methods

Zebrafish strains and husbandry

Request a detailed protocol

All fish were maintained and staged according to established protocols (Kimmel et al., 1995; Westerfield, 1993). Selective breeding was performed as previously described (Nichols et al., 2016; Sucharov et al., 2019; Brooks and Nichols, 2017). Briefly, because the mef2ca mutation is lethal when homozygous, heterozygous full-sibling parents were intercrossed, and their homozygous mutant offspring were scored for penetrance when determining what family was to be propagated for the next generation. The mef2ca^b1086 and mef2cb^fh288 mutant alleles have been previously described (Hinits et al., 2012; Miller et al., 2007) and were maintained by outcrossing to the unselected AB background. When we first introduced the mef2cb^fh288 allele to mef2ca^b1086, we crossed it to the low-penetrance mef2ca^b1086 selectively bred strain and subsequently maintained these double heterozygotes by outcrossing to the unselected AB background.

CRISPR/Cas9-induced mutant alleles

Request a detailed protocol

We generated germline mutant alleles using CRISPR/Cas9 mutagenesis (Hwang et al., 2013) with modifications as described (Mitchell et al., 2021). Briefly, we designed sgRNAs(Single-guide RNA) within or just downstream of the MADS or MEF2 domain. The XbaI-digested pT3TS-nCas9n plasmid (Addgene plasmid #46757) was used as a template to transcribe Cas9 mRNA with the T3 mMESSAGE kit (Invitrogen). We transcribed sgRNAs (see table below) from PCR-generated templates using the MEGAscript T7 Kit (Thermo Fisher Scientific). One cell-stage embryos were injected with a mix of 200 ng/µl Cas9 mRNA and 50 ng/µl of each gene-specific sgRNA. Injected embryos were raised, and founders identified by amplifying the genomic region containing the sgRNA site and identifying banding size shifts indicating insertions and/or deletions (see genotyping assay table for primers). All new paralog mutants were originally generated in the low-penetrance strain and were subsequently maintained on the AB background for at least three generations. The following sgRNAs were used: mef2d^co3013: 5’-GGACAAATACCGGAAGAGCG-3’; mef2b^co2012: 5’-CACGAGAGCCGCACTAACAC-3’; mef2aa^co3017: 5’-TCATGGACGACCGTTTCGGC-3’.

We generated six independent sgRNAs for mef2ab, and none of them mutagenized this locus. Precise sequences of mutant alleles are indicated below:

Gene	Mutant insertion (underlined) and deletion (italicised) sequence
mef2d ^co3013	Exon3: AATACCGGAAGATCGAGGAGCTGGATATCCTC
mef2b ^co3012	Exon3: AACCTCACGAGAGCCGCACTAACAC
mef2aa^co3017	Exon5: TCATGCCCCTGGACGACCGTTTCGGCAAA

Cartilage and bone staining and imaging

Request a detailed protocol

Fixed animals were stained with Alcian blue and Alizarin red as described previously (Brooks and Nichols, 2017; Walker and Kimmel, 2007). Alcian blue- and Alizarin red-stained 6 days post fertilization (dpf) skeletons were dissected and flat mounted for Nomarski imaging on a Leica DMi8 inverted microscope equipped with a Leica DMC2900 as previously described (Nichols et al., 2013).

Phenotype scoring

Request a detailed protocol

For penetrance scoring, 6 dpf Alcian blue- and Alizarin red-stained skeletons were genotyped then scored for the proportion of animals with a given genotype that exhibit a particular phenotype. In the interest of strong rigor and reproducibility, phenotypes were scored by three observers blinded to genotype. All three agreed with the number of animals in each phenotypic class, indicating that phenotype penetrance can be reproducibly identified by different observers. For symplectic cartilage measurements, whole mount Alcian blue- and Alizarin red-stained skeletons were imaged under a transmitted light dissecting scope. Images were captured with Zeiss ZEN software. This software was then used to measure the linear distance in microns from the posterior most point of the interhyal cartilage and the distal tip of the symplectic cartilage similar to a previous study (Talbot et al., 2016). For each individual, the total symplectic length was calculated by summing the measurements from the left and right side. These same measurements were also used to calculate the absolute value of the left minus the right symplectic cartilage lengths to determine developmental instability for each individual. All raw phenotype data are presented in source data tables.

Reverse transcription-quantitative polymerase chain reaction

Request a detailed protocol

Gene expression studies were performed as previously described (Sucharov et al., 2019). For the time course study from unselected AB and for comparing low- and high-penetrance wild types, live individual 24 hpf embryos from each strain had their heads removed. Decapitated bodies were genotyped to identify homozygous wild types. Heads from five to six identified homozygous wild types were pooled, and total RNA was extracted with TRI Reagent. cDNA was prepared with Superscript III from Invitrogen. qPCR experiments utilized a real-time PCR StepOnePlus system from Applied Biosystems and SYBR green. A standard curve was generated from serially diluted (1:2:10) cDNA pools, and primers with a slope of –3.3 ±0.3 were accepted. The relative quantity of target cDNA was calculated using Applied Biosystems StepOne V.2.0 software and the comparative Ct method. After surveying the expression of many housekeeping genes at multiple stages, we determined that rps18 expression was the most consistent across stages, genotypes, and strains. Target gene expression in all experiments was normalized to rps18. Reactions were performed in technical triplicate, and the results represent two to six biological replicates. The following primers were used: rps18 FW, 5’-CTGAACAGACAGAAGGACATAA-3’ and rps18 REV 5’-AGCCTCTCCAGATCTTCTC-3’, mef2ca FW, 5’-GTCCAGAATCCGAGGACAAATA-3’ and mef2ca REV 5’-GAGACAGGCATGTCGTAGTTAG-3’, mef2cb FW, 5’-AGTACGCCAGCACAGATA-3’ and mef2cb REV 5’-AGCCATTTAGACCCTTCTTTC –3’, mef2aa FW, 5’-CCACGAGAGCAGAACCAACTC-3’ and mef2aa REV 5’-GTCCATGAGGGGACTGTGAC-3’, mef2ab FW, 5’-AACCTCACGAGAGCAGAACC-3’ and mef2ab REV 5’-AGGACATATGAGGCGTCTGG-3’, mef2b FW, 5’-CCGATATGGACAAAGTGCTG-3’ and mef2b REV 5’-CCAATCCCAATCCTTTCCTT-3’, mef2d FW, 5’-TTCCAGTATGCCAGCACTGA-3’ and mef2d REV 5’-CGAATCACGGTGCTCTTTCT-3’. All qPCR numerical data and statistical analyses are reported in supplementary data table.

scRNA-seq analysis

Request a detailed protocol

We analyzed our published data set (Mitchell et al., 2021) from 24 hpf wild-type AB zebrafish cranial neural crest cells for mef2 paralog expression. Seurat’s ‘VlnPlot’ function was used to generate mef2 paralog plots indicating the distribution of cells expressing each paralog across each of the four clusters we described in the previous publication. The raw, feature-barcode matrix for this dataset can be accessed from the GEO database (accession number GSE163826).

Genotyping assays

Request a detailed protocol

mef2cab1086 was genotyped by KASP as previously described (Brooks and Nichols, 2017). PCR-based genotyping assays are as follows:

Gene^allele	Primers	Enzyme	wt	mut
mef2d^co3013	Two separate reactions are run per sample. First reaction: Fw Full, Rv Full, and Fw wt Second reaction: Fw Full, Rv Full, and Rev mut Fw Full: 5’-AAGAAAGGCTTTAACGGTTGC-3’ Rv Full: 5’-AAGAGAAGGACGGAGGTTAGA-3’ Fw wt: 5’-GACAAATACCGGAAGAGCGA-3’ Rv mut: 5’-AGAGGATATCCAGCTCCTCTTC-3’	None	173 bp and 98 bp	121 bp and 173 bp
mef2aa^co3017	Fw: 5’-TTGACCCAACGGTTTACAGA-3’ Rv: 5’-CACAAAGCCAAGCAAAAACA-3’	NlaIII	291 bp and 145 bp	Uncut 442 bp
mef2cb ^fh288	Fw: 5'-TCCCTGCTTCTCTCTAGGTGACATTTACATCG-3 ‘ Rv: 5'-TCGTGTGGCTCGTTGTACTC-3’	TaqaI	190 bp and 10 bp	Uncut 200 bp
mef2b ^co3012	Fw: 5’-CGAGATCGCTCTCATCATCTT-3’ Rv: 5’-GACATACTGGAGGTATACAGACCAAA-3’	AciI	112 bp and 38 bp	Uncut 139 bp

Sequence alignments

Request a detailed protocol

We used the Clustal Omega tool by EMBL-EBI for multiple sequence alignment of the different paralog gene products. We obtained the following transcripts and protein products from ENSEMBL for alignments: mef2aa-206 ENSDART00000171594.2, mef2ab-201 ENSDART00000173414.2, mef2b-202 ENSDART00000166300.3, mef2ca-202 ENSDART00000099134.5, mef2cb-207 ENSDART00000183585.1, and mef2d-203 ENSDART00000132589.2ENSDART00000132589.2.

Statistical analyses

Request a detailed protocol

Penetrance scores were compared using Fisher’s exact test to determine significance. All scoring data and exact p-values are reported in source data tables. Studies with symplectic cartilage measurements are presented as box and whisker plots, and the box extends from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median, and the bars are minimum and maximum values. We used a Welch’s t-test or Dunnet’s T3 test to compare total symplectic cartilage (left plus right sides) between genotypes. We used F-test to test for significant differences in variation between genotypes. For developmental instability, the absolute value of the difference between left and right symplectic cartilage lengths were grouped by genotype, and Welch’s t-test, Brown-Forsythe, or Welch’s ANOVA tests were used to determine significant differences in left-right asymmetry between genotypes. Power analyses were used to determine the number of animals to be examined for each experiment. All statistical analyses and exact p-values are reported in source data.

Data availability

All raw data are provided in supplementary data table. Sequencing dataset have been deposited in GEO. The raw, feature-barcode matrix can be accessed from the GEO database (accession number GSE163826).

The following previously published data sets were used

1. Mitchell, et al
(2021) NCBI Gene Expression Omnibus
ID GSE163826. The alx3 gene shapes the zebrafish neurocranium by regulating frontonasal neural crest cell differentiation timing.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163826

References

1. Arnold MA
2. Kim Y
3. Czubryt MP
4. Phan D
5. McAnally J
6. Qi X
7. Shelton JM
8. Richardson JA
9. Bassel-Duby R
10. Olson EN
(2007) Mef2C transcription factor controls chondrocyte hypertrophy and bone development
Developmental Cell 12:377–389.

https://doi.org/10.1016/j.devcel.2007.02.004
- PubMed
- Google Scholar
(2011) No relationship between canalization and developmental stability of the skull in a natural population of Mastomys natalensis (Rodentia: Muridae)
Biological Journal of the Linnean Society 104:207–216.

https://doi.org/10.1111/j.1095-8312.2011.01702.x
- Google Scholar
(2006) A single basis for developmental buffering of Drosophila wing shape
PLOS ONE 1:e7.

https://doi.org/10.1371/journal.pone.0000007
- PubMed
- Google Scholar
1. Brooks EP
2. Nichols JT
(2017) Shifting zebrafish lethal skeletal mutant penetrance by progeny testing
Journal of Visualized Experiments 1:56200.

https://doi.org/10.3791/56200
- PubMed
- Google Scholar
1. Chen R
2. Shi L
3. Hakenberg J
4. Naughton B
5. Sklar P
6. Zhang J
7. Zhou H
8. Tian L
9. Prakash O
10. Lemire M
11. Sleiman P
12. Cheng WY
13. Chen W
14. Shah H
15. Shen Y
16. Fromer M
17. Omberg L
18. Deardorff MA
19. Zackai E
20. Bobe JR
21. Levin E
22. Hudson TJ
23. Groop L
24. Wang J
25. Hakonarson H
26. Wojcicki A
27. Diaz GA
28. Edelmann L
29. Schadt EE
30. Friend SH
(2016) Analysis of 589,306 genomes identifies individuals resilient to severe Mendelian childhood diseases
Nature Biotechnology 34:531–538.

https://doi.org/10.1038/nbt.3514
- PubMed
- Google Scholar
1. Chen X
2. Gao B
3. Ponnusamy M
4. Lin Z
5. Liu J
(2017) Mef2 signaling and human diseases
Oncotarget 8:112152–112165.

https://doi.org/10.18632/oncotarget.22899
- PubMed
- Google Scholar
1. Conant GC
2. Wagner A
(2004) Duplicate genes and robustness to transient gene knock-downs in Caenorhabditis elegans
Proceedings. Biological Sciences 271:89–96.

https://doi.org/10.1098/rspb.2003.2560
- PubMed
- Google Scholar
1. Dandage R
2. Landry CR
(2019) Paralog dependency indirectly affects the robustness of human cells
Molecular Systems Biology 15:e8871.

https://doi.org/10.15252/msb.20198871
- PubMed
- Google Scholar
1. Dean EJ
2. Davis JC
3. Davis RW
4. Petrov DA
(2008) Pervasive and persistent redundancy among duplicated genes in yeast
PLOS Genetics 4:e1000113.

https://doi.org/10.1371/journal.pgen.1000113
- PubMed
- Google Scholar
(2009) Plasticity, canalization, and developmental stability of the Drosophila wing: joint effects of mutations and developmental temperature
Evolution; International Journal of Organic Evolution 63:2864–2876.

https://doi.org/10.1111/j.1558-5646.2009.00774.x
- PubMed
- Google Scholar
1. De Kegel B
2. Ryan CJ
(2019) Paralog buffering contributes to the variable essentiality of genes in cancer cell lines
PLOS Genetics 15:e1008466.

https://doi.org/10.1371/journal.pgen.1008466
- PubMed
- Google Scholar
1. De Laurier A
2. Huycke TR
3. Nichols JT
4. Swartz ME
5. Larsen A
6. Walker C
7. Dowd J
8. Pan L
9. Moens CB
10. Kimmel CB
(2014) Role of mef2ca in developmental buffering of the zebrafish larval hyoid dermal skeleton
Developmental Biology 385:189–199.

https://doi.org/10.1016/j.ydbio.2013.11.016
- PubMed
- Google Scholar
1. de Visser J
2. Hermisson J
3. Wagner GP
4. Ancel Meyers L
5. Bagheri-Chaichian H
6. Blanchard JL
7. Chao L
8. Cheverud JM
9. Elena SF
10. Fontana W
11. Gibson G
12. Hansen TF
13. Krakauer D
14. Lewontin RC
15. Ofria C
16. Rice SH
17. von Dassow G
18. Wagner A
19. Whitlock MC
(2003) Perspective: evolution and detection of genetic robustness
Evolution; International Journal of Organic Evolution 57:1959–1972.

https://doi.org/10.1111/j.0014-3820.2003.tb00377.x
- PubMed
- Google Scholar
1. Diss G
2. Ascencio D
3. DeLuna A
4. Landry CR
(2014) Molecular mechanisms of paralogous compensation and the robustness of cellular networks
Journal of Experimental Zoology. Part B, Molecular and Developmental Evolution 322:488–499.

https://doi.org/10.1002/jez.b.22555
- PubMed
- Google Scholar
1. El-Brolosy MA
2. Kontarakis Z
3. Rossi A
4. Kuenne C
5. Günther S
6. Fukuda N
7. Kikhi K
8. Boezio GLM
9. Takacs CM
10. Lai SL
11. Fukuda R
12. Gerri C
13. Giraldez AJ
14. Stainier DYR
(2019) Genetic compensation triggered by mutant mRNA degradation
Nature 568:193–197.

https://doi.org/10.1038/s41586-019-1064-z
- PubMed
- Google Scholar
1. Fish JL
(2016) Developmental mechanisms underlying variation in craniofacial disease and evolution
Developmental Biology 415:188–197.

https://doi.org/10.1016/j.ydbio.2015.12.019
- PubMed
- Google Scholar
1. Force A
2. Lynch M
3. Pickett FB
4. Amores A
5. Yan YL
6. Postlethwait J
(1999) Preservation of duplicate genes by complementary, degenerative mutations
Genetics 151:1531–1545.

https://doi.org/10.1093/genetics/151.4.1531
- PubMed
- Google Scholar
(2016) The power of natural variation for model organism biology
Trends in Genetics 32:147–154.

https://doi.org/10.1016/j.tig.2015.12.003
- PubMed
- Google Scholar
1. Gordon CT
2. Tessier A
3. Demir Z
4. Goldenberg A
5. Oufadem M
6. Voisin N
7. Pingault V
8. Bienvenu T
9. Lyonnet S
10. de Pontual L
11. Amiel J
(2018) The association of severe encephalopathy and question mark ear is highly suggestive of loss of MEF2C function
Clinical Genetics 93:356–359.

https://doi.org/10.1111/cge.13046
- PubMed
- Google Scholar
1. Gu Z
2. Steinmetz LM
3. Gu X
4. Scharfe C
5. Davis RW
6. Li WH
(2003) Role of duplicate genes in genetic robustness against null mutations
Nature 421:63–66.

https://doi.org/10.1038/nature01198
- PubMed
- Google Scholar
(2006) The brachymorph mouse and the developmental-genetic basis for canalization and morphological integration
Evolution & Development 8:61–73.

https://doi.org/10.1111/j.1525-142X.2006.05075.x
- PubMed
- Google Scholar
1. Hallgrimsson B
2. Green RM
3. Katz DC
4. Fish JL
5. Bernier FP
6. Roseman CC
7. Young NM
8. Cheverud JM
9. Marcucio RS
(2019) The developmental-genetics of canalization
Seminars in Cell & Developmental Biology 88:67–79.

https://doi.org/10.1016/j.semcdb.2018.05.019
- PubMed
- Google Scholar
1. He M
2. Zhou D
3. Ding NZ
4. Teng CB
5. Yan XC
6. Liang Y
(2019) Common carp MEF2 genes: evolution and expression
Genes 10:588.

https://doi.org/10.3390/genes10080588
- Google Scholar
1. Hinits Y
2. Pan L
3. Walker C
4. Dowd J
5. Moens CB
6. Hughes SM
(2012) Zebrafish mef2ca and mef2cb are essential for both first and second heart field cardiomyocyte differentiation
Developmental Biology 369:199–210.

https://doi.org/10.1016/j.ydbio.2012.06.019
- PubMed
- Google Scholar
1. Hsiao TL
2. Vitkup D
(2008) Role of duplicate genes in robustness against deleterious human mutations
PLOS Genetics 4:e1000014.

https://doi.org/10.1371/journal.pgen.1000014
- PubMed
- Google Scholar
1. Hwang WY
2. Fu Y
3. Reyon D
4. Maeder ML
5. Tsai SQ
6. Sander JD
7. Peterson RT
8. Yeh JRJ
9. Joung JK
(2013) Efficient genome editing in zebrafish using a CRISPR-Cas system
Nature Biotechnology 31:227–229.

https://doi.org/10.1038/nbt.2501
- PubMed
- Google Scholar
(2005) Transcription control reprogramming in genetic backup circuits
Nature Genetics 37:295–299.

https://doi.org/10.1038/ng1523
- PubMed
- Google Scholar
1. Kafri R
2. Levy M
3. Pilpel Y
(2006) The regulatory utilization of genetic redundancy through responsive backup circuits
PNAS 103:11653–11658.

https://doi.org/10.1073/pnas.0604883103
- PubMed
- Google Scholar
1. Kamath RS
2. Fraser AG
3. Dong Y
4. Poulin G
5. Durbin R
6. Gotta M
7. Kanapin A
8. Le Bot N
9. Moreno S
10. Sohrmann M
11. Welchman DP
12. Zipperlen P
13. Ahringer J
(2003) Systematic functional analysis of the Caenorhabditis elegans genome using RNAi
Nature 421:231–237.

https://doi.org/10.1038/nature01278
- PubMed
- Google Scholar
1. Kile BT
2. Hilton DJ
(2005) The art and design of genetic screens: mouse
Nature Reviews. Genetics 6:557–567.

https://doi.org/10.1038/nrg1636
- PubMed
- Google Scholar
1. Kim Y
2. Phan D
3. van Rooij E
4. Wang DZ
5. McAnally J
6. Qi X
7. Richardson JA
8. Hill JA
9. Bassel-Duby R
10. Olson EN
(2008) The MEF2D transcription factor mediates stress-dependent cardiac remodeling in mice
The Journal of Clinical Investigation 118:124–132.

https://doi.org/10.1172/JCI33255
- PubMed
- Google Scholar
(1995) Stages of embryonic development of the zebrafish
Developmental Dynamics 203:253–310.

https://doi.org/10.1002/aja.1002030302
- PubMed
- Google Scholar
1. Kimmel CB
2. Ullmann B
3. Walker M
4. Miller CT
5. Crump JG
(2003) Endothelin 1-mediated regulation of pharyngeal bone development in zebrafish
Development 130:1339–1351.

https://doi.org/10.1242/dev.00338
- PubMed
- Google Scholar
(2015) Patterns of variation and covariation in the shapes of mandibular bones of juvenile salmonids in the genus Oncorhynchus
Evolution & Development 17:302–314.

https://doi.org/10.1111/ede.12135
- PubMed
- Google Scholar
1. Kirschner M
2. Gerhart J
(1998) Evolvability
PNAS 95:8420–8427.

https://doi.org/10.1073/pnas.95.15.8420
- PubMed
- Google Scholar
1. Kohn LAP
(1991) The role of genetics in craniofacial morphology and growth
Annual Review of Anthropology 20:261–278.

https://doi.org/10.1146/annurev.an.20.100191.001401
- Google Scholar
1. Le Meur N
2. Holder-Espinasse M
3. Jaillard S
4. Goldenberg A
5. Joriot S
6. Amati-Bonneau P
7. Guichet A
8. Barth M
9. Charollais A
10. Journel H
11. Auvin S
12. Boucher C
13. Kerckaert JP
14. David V
15. Manouvrier-Hanu S
16. Saugier-Veber P
17. Frébourg T
18. Dubourg C
19. Andrieux J
20. Bonneau D
(2010) Mef2C haploinsufficiency caused by either microdeletion of the 5q14.3 region or mutation is responsible for severe mental retardation with stereotypic movements, epilepsy and/or cerebral malformations
Journal of Medical Genetics 47:22–29.

https://doi.org/10.1136/jmg.2009.069732
- PubMed
- Google Scholar
1. Li J
2. Yuan Z
3. Zhang Z
(2010) The cellular robustness by genetic redundancy in budding yeast
PLOS Genetics 6:e1001187.

https://doi.org/10.1371/journal.pgen.1001187
- PubMed
- Google Scholar
1. Lin Q
2. Schwarz J
3. Bucana C
4. Olson EN
(1997) Control of mouse cardiac morphogenesis and myogenesis by transcription factor MEF2C
Science 276:1404–1407.

https://doi.org/10.1126/science.276.5317.1404
- PubMed
- Google Scholar
1. Liu F
2. van der Lijn F
3. Schurmann C
4. Zhu G
5. Chakravarty MM
6. Hysi PG
7. Wollstein A
8. Lao O
9. de Bruijne M
10. Ikram MA
11. van der Lugt A
12. Rivadeneira F
13. Uitterlinden AG
14. Hofman A
15. Niessen WJ
16. Homuth G
17. de Zubicaray G
18. McMahon KL
19. Thompson PM
20. Daboul A
21. Puls R
22. Hegenscheid K
23. Bevan L
24. Pausova Z
25. Medland SE
26. Montgomery GW
27. Wright MJ
28. Wicking C
29. Boehringer S
30. Spector TD
31. Paus T
32. Martin NG
33. Biffar R
34. Kayser M
(2012) A genome-wide association study identifies five loci influencing facial morphology in Europeans
PLOS Genetics 8:e1002932.

https://doi.org/10.1371/journal.pgen.1002932
- PubMed
- Google Scholar
1. Majidi SP
2. Reddy NC
3. Moore MJ
4. Chen H
5. Yamada T
6. Andzelm MM
7. Cherry TJ
8. Hu LS
9. Greenberg ME
10. Bonni A
(2019) Chromatin environment and cellular context specify compensatory activity of paralogous MEF2 transcription factors
Cell Reports 29:2001–2015.

https://doi.org/10.1016/j.celrep.2019.10.033
- PubMed
- Google Scholar
(2004) Investigating the origins of triploblasty: “ mesodermal ” gene expression in a diploblastic animal, the sea anemone Nematostella vectensis (phylum, Cnidaria; class, Anthozoa)
Development 131:2463–2474.

https://doi.org/10.1242/dev.01119
- PubMed
- Google Scholar
1. Miller CT
2. Swartz ME
3. Khuu PA
4. Walker MB
5. Eberhart JK
6. Kimmel CB
(2007) Mef2ca is required in cranial neural crest to effect endothelin1 signaling in zebrafish
Developmental Biology 308:144–157.

https://doi.org/10.1016/j.ydbio.2007.05.018
- PubMed
- Google Scholar
(2021) The Alx3 gene shapes the zebrafish neurocranium by regulating frontonasal neural crest cell differentiation timing
Development 148:dev197483.

https://doi.org/10.1242/dev.197483
- PubMed
- Google Scholar
(1996) Mutational analysis of the DNA binding, dimerization, and transcriptional activation domains of MEF2C
Molecular and Cellular Biology 16:2627–2636.

https://doi.org/10.1128/MCB.16.6.2627
- PubMed
- Google Scholar
1. Nichols JT
2. Pan L
3. Moens CB
4. Kimmel CB
(2013) Barx1 represses joints and promotes cartilage in the craniofacial skeleton
Development 140:2765–2775.

https://doi.org/10.1242/dev.090639
- PubMed
- Google Scholar
(2016) Ligament versus bone cell identity in the zebrafish hyoid skeleton is regulated by mef2ca
Development 143:4430–4440.

https://doi.org/10.1242/dev.141036
- PubMed
- Google Scholar
1. Novara F
2. Rizzo A
3. Bedini G
4. Girgenti V
5. Esposito S
6. Pantaleoni C
7. Ciccone R
8. Sciacca FL
9. Achille V
10. Della Mina E
11. Gana S
12. Zuffardi O
13. Estienne M
(2013) Mef2C deletions and mutations versus duplications: a clinical comparison
European Journal of Medical Genetics 56:260–265.

https://doi.org/10.1016/j.ejmg.2013.01.011
- PubMed
- Google Scholar
1. Nowakowska BA
2. Obersztyn E
3. Szymańska K
4. Bekiesińska-Figatowska M
5. Xia Z
6. Ricks CB
7. Bocian E
8. Stockton DW
9. Szczałuba K
10. Nawara M
11. Patel A
12. Scott DA
13. Cheung SW
14. Bohan TP
15. Stankiewicz P
(2010) Severe mental retardation, seizures, and hypotonia due to deletions of MEF2C
American Journal of Medical Genetics. Part B, Neuropsychiatric Genetics 153B:1042–1051.

https://doi.org/10.1002/ajmg.b.31071
- PubMed
- Google Scholar
Book
1. Ohno S
(2013) Evolution by Gene Duplication
New York: Springer Science & Business Media.

https://doi.org/10.1007/978-3-642-86659-3
- Google Scholar
(2008) Phenotypic variability and craniofacial dysmorphology: increased shape variance in a mouse model for cleft lip
Journal of Anatomy 212:135–143.

https://doi.org/10.1111/j.1469-7580.2007.00845.x
- PubMed
- Google Scholar
1. Paternoster L
2. Zhurov AI
3. Toma AM
4. Kemp JP
5. St Pourcain B
6. Timpson NJ
7. McMahon G
8. McArdle W
9. Ring SM
10. Smith GD
11. Richmond S
12. Evans DM
(2012) Genome-Wide association study of three-dimensional facial morphology identifies a variant in PAX3 associated with nasion position
American Journal of Human Genetics 90:478–485.

https://doi.org/10.1016/j.ajhg.2011.12.021
- PubMed
- Google Scholar
(1996) Jaw and branchial arch mutants in zebrafish II: anterior arches and cartilage differentiation
Development 123:345–356.

https://doi.org/10.1242/dev.123.1.345
- PubMed
- Google Scholar
1. Potthoff MJ
2. Olson EN
(2007) Mef2: a central regulator of diverse developmental programs
Development 134:4131–4140.

https://doi.org/10.1242/dev.008367
- PubMed
- Google Scholar
(2010) Variability in gene expression underlies incomplete penetrance
Nature 463:913–918.

https://doi.org/10.1038/nature08781
- PubMed
- Google Scholar
1. Rossi A
2. Kontarakis Z
3. Gerri C
4. Nolte H
5. Hölper S
6. Krüger M
7. Stainier DYR
(2015) Genetic compensation induced by deleterious mutations but not gene knockdowns
Nature 524:230–233.

https://doi.org/10.1038/nature14580
- PubMed
- Google Scholar
(2005) Genetics and geometry of canalization and developmental stability in Drosophila subobscura
BMC Evolutionary Biology 5:1–28.

https://doi.org/10.1186/1471-2148-5-7
- PubMed
- Google Scholar
(2013) Edn1 and HAND2 interact in early regulation of pharyngeal arch outgrowth during zebrafish development
PLOS ONE 8:e67522.

https://doi.org/10.1371/journal.pone.0067522
- PubMed
- Google Scholar
1. Scharloo W
(1991) Canalization: genetic and developmental aspects
Annual Review of Ecology and Systematics 22:65–93.

https://doi.org/10.1146/annurev.es.22.110191.000433
- Google Scholar
1. Sheehan MJ
2. Nachman MW
(2014) Morphological and population genomic evidence that human faces have evolved to signal individual identity
Nature Communications 5:1–10.

https://doi.org/10.1038/ncomms5800
- PubMed
- Google Scholar
1. Square TA
2. Jandzik D
3. Massey JL
4. Romášek M
5. Stein HP
6. Hansen AW
7. Purkayastha A
8. Cattell MV
9. Medeiros DM
(2020) Evolution of the endothelin pathway drove neural crest cell diversification
Nature 585:563–568.

https://doi.org/10.1038/s41586-020-2720-z
- PubMed
- Google Scholar
1. St Johnston D
(2002) The art and design of genetic screens: Drosophila melanogaster
Nature Reviews. Genetics 3:176–188.

https://doi.org/10.1038/nrg751
- PubMed
- Google Scholar
1. Sucharov J
2. Ray K
3. Brooks EP
4. Nichols JT
(2019) Selective breeding modifies mef2ca mutant incomplete penetrance by tuning the opposing Notch pathway
PLOS Genetics 15:e1008507.

https://doi.org/10.1371/journal.pgen.1008507
- PubMed
- Google Scholar
1. Takahata N
2. Maruyama T
(1979) Polymorphism and loss of duplicate gene expression: a theoretical study with application of tetraploid fish
PNAS 76:4521–4525.

https://doi.org/10.1073/pnas.76.9.4521
- PubMed
- Google Scholar
(2016) Pharyngeal morphogenesis requires fras1-itga8-dependent epithelial-mesenchymal interaction
Developmental Biology 416:136–148.

https://doi.org/10.1016/j.ydbio.2016.05.035
- PubMed
- Google Scholar
1. Tonk V
2. Kyhm JH
3. Gibson CE
4. Wilson GN
(2011) Interstitial deletion 5q14.3q21.3 with MEF2C haploinsufficiency and mild phenotype: when more is less
American Journal of Medical Genetics. Part A 155A:1437–1441.

https://doi.org/10.1002/ajmg.a.34012
- PubMed
- Google Scholar
(1993) Analyses of craniofacial and dental morphology in monozygotic twins discordant for cleft lip and unilateral cleft lip and palate
The Angle Orthodontist 63:135–139.

https://doi.org/10.1043/0003-3219(1993)063<0135:AOCADM>2.0.CO;2
- PubMed
- Google Scholar
1. Valen LV
(1962) A study of fluctuating asymmetry
Evolution; International Journal of Organic Evolution 16:125–142.

https://doi.org/10.1111/j.1558-5646.1962.tb03206.x
- Google Scholar
1. Verzi MP
2. Agarwal P
3. Brown C
4. McCulley DJ
5. Schwarz JJ
6. Black BL
(2007) The transcription factor MEF2C is required for craniofacial development
Developmental Cell 12:645–652.

https://doi.org/10.1016/j.devcel.2007.03.007
- PubMed
- Google Scholar
1. Vu V
2. Verster AJ
3. Schertzberg M
4. Chuluunbaatar T
5. Spensley M
6. Pajkic D
7. Hart GT
8. Moffat J
9. Fraser AG
(2015) Natural variation in gene expression modulates the severity of mutant phenotypes
Cell 162:391–402.

https://doi.org/10.1016/j.cell.2015.06.037
- PubMed
- Google Scholar
1. Waddington CH
(1942) Canalization of development and the inheritance of acquired characters
Nature 150:563–565.

https://doi.org/10.1038/150563a0
- Google Scholar
1. Waddington CH
(1952) Selection of the genetic basis for an acquired character
Nature 169:625–626.

https://doi.org/10.1038/169625b0
- PubMed
- Google Scholar
1. Waddington CH
(1953) Genetic assimilation of an acquired character
Evolution; International Journal of Organic Evolution 7:118–126.

https://doi.org/10.1111/j.1558-5646.1953.tb00070.x
- Google Scholar
1. Waddington CH
(1956) Genetic assimilation of the bithorax phenotype
Evolution; International Journal of Organic Evolution 10:1.

https://doi.org/10.2307/2406091
- Google Scholar
1. Waddington CH
(1957) The genetic basis of the ‘assimilated bithorax’ stock
Journal of Genetics 55:241–245.

https://doi.org/10.1007/BF02981639
- Google Scholar
1. Waddington CH
(1959) Canalization of development and genetic assimilation of acquired characters
Nature 183:1654–1655.

https://doi.org/10.1038/1831654a0
- PubMed
- Google Scholar
1. Wagner A
(2000) Robustness against mutations in genetic networks of yeast
Nature Genetics 24:355–361.

https://doi.org/10.1038/74174
- PubMed
- Google Scholar
1. Wagner A
(2008) Gene duplications, robustness and evolutionary innovations
BioEssays : News and Reviews in Molecular, Cellular and Developmental Biology 30:367–373.

https://doi.org/10.1002/bies.20728
- PubMed
- Google Scholar
1. Walker MB
2. Kimmel CB
(2007) A two-color acid-free cartilage and bone stain for zebrafish larvae
Biotechnic & Histochemistry 82:23–28.

https://doi.org/10.1080/10520290701333558
- PubMed
- Google Scholar
(2000) Quantitative approach to identifying abnormal variation in the human face exemplified by a study of 278 individuals with five craniofacial syndromes
American Journal of Medical Genetics 91:8–17.

https://doi.org/10.1002/(sici)1096-8628(20000306)91:1<8::aid-ajmg2>3.0.co;2-#
- PubMed
- Google Scholar
1. Watterson GA
(1983) On the time for gene silencing at duplicate loci
Genetics 105:745–766.

https://doi.org/10.1093/genetics/105.3.745
- PubMed
- Google Scholar
Book
1. Westerfield M
(1993)
The Zebrafish Book: A Guide for the Laboratory Use of Zebrafish

Eugene: University of Oregon Press.
- Google Scholar
1. Wu W
2. de Folter S
3. Shen X
4. Zhang W
5. Tao S
(2011a) Vertebrate paralogous MEF2 genes: origin, conservation, and evolution
PLOS ONE 6:e17334.

https://doi.org/10.1371/journal.pone.0017334
- PubMed
- Google Scholar
1. Wu W
2. Huang X
3. Cheng J
4. Li Z
5. de Folter S
6. Huang Z
7. Jiang X
8. Pang H
9. Tao S
(2011b) Conservation and evolution in and among SRF- and MEF2-type MADS domains and their binding sites
Molecular Biology and Evolution 28:501–511.

https://doi.org/10.1093/molbev/msq214
- PubMed
- Google Scholar
1. Young JI
2. Slifer S
3. Hecht JT
4. Blanton SH
(2021) Dna methylation variation is identified in monozygotic twins discordant for non-syndromic cleft lip and palate
Frontiers in Cell and Developmental Biology 9:656865.

https://doi.org/10.3389/fcell.2021.656865
- PubMed
- Google Scholar
1. Zweier M
2. Gregor A
3. Zweier C
4. Engels H
5. Sticht H
6. Wohlleber E
7. Bijlsma EK
8. Holder SE
9. Zenker M
10. Rossier E
11. Grasshoff U
12. Johnson DS
13. Robertson L
14. Firth HV
15. Ekici AB
16. Reis A
17. Rauch A
18. Cornelia Kraus
(2010) Mutations in MEF2C from the 5q14.3q15 microdeletion syndrome region are a frequent cause of severe mental retardation and diminish MeCP2 and CDKL5 expression
Human Mutation 31:722–733.

https://doi.org/10.1002/humu.21253
- PubMed
- Google Scholar

Article and author information

Author details

Raisa Bailon-Zambrano

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Formal analysis, Investigation, Writing – review and editing

Contributed equally with
Juliana Sucharov

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5848-3952
Juliana Sucharov

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Conceptualization, Formal analysis, Investigation, Writing – review and editing

Contributed equally with
Raisa Bailon-Zambrano

Competing interests
No competing interests declared
Abigail Mumme-Monheit

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Investigation, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0090-1418
Matthew Murry

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Investigation

Competing interests
No competing interests declared
Amanda Stenzel

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Investigation

Competing interests
No competing interests declared
Anthony T Pulvino

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Investigation

Competing interests
No competing interests declared
Jennyfer M Mitchell

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Investigation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4222-5235
Kathryn L Colborn

Department of Surgery, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Methodology

Competing interests
No competing interests declared
James T Nichols

Department of Craniofacial Biology, University of Colorado Anschutz Medical Campus, Aurora, United States

Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing

For correspondence
JAMES.NICHOLS@UCDENVER.EDU

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-7263-1704

Funding

National Institute of Dental and Craniofacial Research (R01 DE029193)

James T Nichols

National Science Foundation (Graduate Research Fellowships Program 201569)

Raisa Bailon-Zambrano

National Institute of Dental and Craniofacial Research (F32 DE029995)

Jennyfer M Mitchell

National Science Foundation (Graduate Research Fellowships Program)

Abigail Mumme-Monheit

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Charles Kimmel for careful reading of this manuscript, members of the department of craniofacial biology for insightful discussions, Austin Tillery for contributions in the early stages of this study, and the University of Colorado zebrafish care staff.

Ethics

All of our work with zebrafish has been approved by the University of Colorado Institutional Animal Care and Use Committee (IACUC), Protocol # 00188. Animals were euthanized by hypothermic shock followed by 1.5% sodium hypochlorite.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.