Ligament versus bone cell identity in the zebrafish hyoid skeleton is regulated by mef2ca

When development goes awry in genetic mutants, the resulting phenotypes often vary substantially compared with the wild type (Waddington, 1942). Why some mutants have overt, dramatic phenotypes but their genetically similar mutant siblings appear wild type remains largely mysterious. We previously reported an extremely variable phenotype in the zebrafish mef2ca^b1086 mutant (DeLaurier et al., 2014). Mef2c encodes a MADS domain-containing transcription factor crucial for wild-type skeletal development in mice and zebrafish (Arnold et al., 2007; Miller et al., 2007; Verzi et al., 2007). The zebrafish mef2ca^b1086 allele, which contains an early stop codon and is thought to encode just the DNA-binding MADS domain, results in a notably variable ectopic bone phenotype involving the opercle, a bone supporting the gill cover (DeLaurier et al., 2014; Miller et al., 2007). Within this gain-of-bone phenotype, the recurrence of an ectopic bony strut led us to hypothesize that a cryptic developmental pathway present in wild types is revealed in mutants (DeLaurier et al., 2014). Therefore, even in the face of extreme variation, careful study reveals consistencies in the phenotype that can disclose important insights into the developmental mechanisms that underlie the phenotype and its variability.

The larval zebrafish head skeleton consists of a wide array of different cell types including bone-producing osteoblasts and ligament-generating cells (‘ligamentoblasts’). Osteoblasts and ligamentoblasts in the head largely originate from a field of neural crest-derived progenitors, which is exposed to signals that regulate the expression of transcription factors (Barske et al., 2016; Medeiros and Crump, 2012). In turn, the transcription factors execute cell fate choices promoting the correct cell type at the proper location. When transcription factor function is disrupted, as in genetic mutants, cells can choose an alternative cell fate such that the animal develops an excess of one cell type at the expense of another. For example, the transcription factor barx1 controls the choice to become a joint cell versus a cartilage cell in the zebrafish lower jaw, resulting in the gain of joint cells at the expense of cartilage cells in barx1 mutants (Nichols et al., 2013). Mef2c is involved in cell fate decisions in the mouse hematopoietic system (Schüler et al., 2008; Stehling-Sun et al., 2009), motivating the hypothesis that mef2ca functions either as a switch or upstream to a switch in the head skeleton. However, what cell types, if any, mef2ca^b1086 mutants lose when they gain osteoblasts has not yet been examined.

Whereas single gene mutations may allow cells to choose an alternative fate, mutant phenotypes are frequently variable and may even include phenotypically wild-type individuals among genotypic mutants, a phenomenon known as incomplete penetrance. It is likely that complex genetics, environmental variations and stochastic developmental noise contribute to the penetrance of mutant phenotypes. One model accounting for variable phenotypes that fall into two classes, affected versus unaffected, is the threshold character model of quantitative genetics (Falconer et al., 1996; Lynch and Walsh, 1998). In the threshold character model, a combination of genes and environment controls a normally unobservable continuous variable termed ‘liability’, which underlies an outwardly observable phenotype. The threshold model has been invoked to explain the inheritance and incomplete penetrance of human disease alleles (Dipple and McCabe, 2000a,b).

One liability proposed to underlie incomplete penetrance is transcriptional variability (Raj et al., 2010). A recently appreciated source of transcriptional variability that we explore in this study involves mobile repetitive DNA elements, transposons. Transposons can mediate transcriptional variation by interacting with host genes in cis to control transcription of nearby genes as well as variably altering the local epigenetic landscape (Feschotte, 2008; Rebollo et al., 2012; Slotkin and Martienssen, 2007; Whitelaw and Martin, 2001). Furthermore, stochastic variation in epigenetic silencing of transposons can account for phenotypic variation, even among genetically identical animals (Morgan et al., 1999; Whitelaw and Martin, 2001). The degree to which transposon-host interactions contribute to wild-type or mutant development remains largely unknown.

Here, we detect consistency among the variable mef2ca^b1086 mutant phenotypes. We find that the earliest manifestation of the ectopic bone phenotype is often the appearance of osteoblasts along a ‘track’ running perpendicular to the developing opercle. The cells normally in this track are those of a ligament, newly discovered in this study, which appear to be variably ‘fate switched’ to osteoblasts in the mef2ca^b1086 mutant. Fate-switching variability is heritable such that selective breeding can drive the penetrance of the phenotype up or down, and we have selectively bred two strains to high and low penetrance. We report differential expression of the mef2ca^b1086 mutant allele in our high and low penetrance strains and that overexpression of mef2ca^b1086 mutant transcripts in the low penetrance strain increases fate-switch penetrance whereas knockdown of mef2ca^b1086 translation in the high penetrance strain decreases fate-switch penetrance. Furthermore, differential mef2ca^b1086 expression is correlated with differential DNA methylation on a transposable element at the mef2ca locus. From these studies, we propose an epigenetic basis of the phenotypic variation: the expression of the deleterious, likely antimorphic mef2ca^b1086 allele is the liability underlying fate-switch penetrance, which is correlated with variable epigenetic transposon silencing.

The prominent opercle and branchiostegal ray bones in the larval zebrafish head skeleton (Fig. 1A) are variably expanded in mef2ca^b1086 mutants. The phenotype ranges from wild-type-like bones (Fig. 1B, mef2ca^b1086 mild) to large plate-like bony expansions (Fig. 1B, mef2ca^b1086 severe). One recurring mutant phenotype is a skeletal strut or bone bridge orthogonal to the direction of normal opercle growth (Fig. 1B, arrows), leading us to propose that mef2ca^b1086 mutants reveal a cryptic developmental pathway (DeLaurier et al., 2014). Furthermore, we observed that the first ectopic osteoblasts appear along a tissue track orthogonal to the developing opercle. These findings motivated us to assess carefully whether there are any regular patterns of ectopic bone formation in mutants. To this end, time-lapse recordings of mef2ca^b1086 mutant zebrafish carrying the sp7:EGFP osteoblast transgene were ‘skeletonized’ (see Materials and Methods) and aligned to observe more easily the location and direction of the protrusions of ectopic growth (Fig. 1C; Movie 1). In mef2ca^b1086 mutants, the nascent orthotopic opercle grows predominantly in a straight line in the ventral direction, as we have previously described for wild types (Kimmel et al., 2010). However, the skeleton representing the first ectopic osteoblasts in this time-lapse recording reveals their location in a protrusion 40 µm from the dorsal (joint) end of the opercle extending in an essentially straight line orthogonal to the direction of opercle growth (Fig. 1C, arrows). By the final frames, the skeletonized outline reveals the line of orthotopic osteoblasts in the dorsal-ventral axis with the ectopic outgrowth extending along a perpendicular track in the anterior-posterior axis. Thirteen skeletonized mef2ca^b1086 mutant time-lapse recordings were overlaid to uncover any consistencies in the location and orientation of ectopic bone in mef2ca^b1086 mutants (Fig. 1D). Given the extreme variation in late-stage phenotypes (Fig. 1B), we were surprised to discover that both the location and orientation of the earliest ectopic growth is very ordered. The location is consistently between 30 and 70 µm from the dorsal end of the bone along the major axis, and the orientation is almost exclusively orthogonal to the direction of primary opercle growth, thus defining the track of ectopic osteoblasts. These data support our hypothesis that a special population of cells residing in the track perpendicular to the opercle is competent to transform into bone in mef2ca^b1086 mutants.

A clue to the identity of the track of tissue came from observing skeletal mobility in non-anesthetized wild-type larvae during jaw opening (Movie 2). We observed the coordinated movement of the opercle and the adjacent ceratohyal cartilage (Fig. 1E), suggesting that, during the larval stage, some tether links them. Ligaments are cord-like structures that tether bone to bone or bone to cartilage, and we hypothesized that a ligament connects the anterior apex of the opercle to the posterior end of the ceratohyal. To our knowledge, there is no description of such a ligament in any fish.

To determine whether cells residing in the track between the opercle and ceratohyal have molecular characteristics associated with ligament identity, we developed methods of visualizing ligaments and examined them in wild-type animals. Scleraxis encodes a basic helix-loop-helix transcription factor associated with ligament identity in mouse, chick and zebrafish (Chen and Galloway, 2014; Cserjesi et al., 1995; Schweitzer et al., 2001). Expression of the zebrafish scleraxis a co-ortholog (scxa) has been reported in the developing pharyngeal skeleton (Chen and Galloway, 2014), although expression at the putative ligament attached to the opercle was not described. We found scxa mRNA expression along the track of tissue orthogonal and anterior to the developing opercle described above (Fig. 2A, arrow). A second marker, Thrombospondin 4b (Thbs4b) is an extracellular matrix protein required for tendon attachment in zebrafish (Subramanian and Schilling, 2014, 2015). We detected strong Thbs4b protein immunoreactivity in the track of tissue (Fig. 2B, arrow) extending orthogonal and anterior to the opercle (Fig. 2B, op) matching the observed track of ectopic bone forming cells in Fig. 1D. A third marker of the putative ligament is a member of the aristaless domain-containing family of transcription factors that is expressed in a subset of neural crest-derived craniofacial mesenchyme including early osteoblasts in mice (Hudson et al., 1998; Qu et al., 1997, 1998). In zebrafish, expression of one of these, alx4a, marks a subset of the mesenchyme that will develop into craniofacial structures in the head (Dee et al., 2013). Furthermore, alx4a regulatory sequences drive transgene expression in osteoblasts and fibroblast-like cells in the anterior pectoral fin (Nachtrab et al., 2013). We found that the alx4a transgene is expressed in elongated cells in the track of tissue extending orthogonal to the opercle (Fig. 2C, arrow), as well as more ventrally positioned rounded cells. Imaging the cartilage skeleton along with the alx4a reporter revealed that the alx4a track terminates at the ceratohyal cartilage as predicted for the ligament in question (Fig. S1), supporting the notion that alx4a is expressed in ligament cells. We further note that a subset of osteoblasts in the opercle where the putative ligament attaches to this bone also expresses alx4a (Fig. 2C, double label). These three labels (scxa, Thbs4b and alx4a) show that the track of tissue linking the anterior end of the opercle to the posterior end of the ceratohyal in larval zebrafish is a ligament, which we have named the operculohyoid ligament.

In mef2ca^b1086 mutant zebrafish, we found that whenever ectopic osteoblasts were detected by ectopic sp7 expression, ligament cells were correspondingly lost, revealed by decreased scxa expression (Fig. 2D, asterisk), Thrombospondin 4b antibody staining (Fig. 2E, asterisk) or alx4a transgene expression (Fig. 2F, asterisk). The complementary nature of the osteoblast and ligament markers in mutants suggest that bone cells are gained at the expense of ligament cells in mef2ca^b1086 mutants, motivating the hypothesis that cells fated to become ligamentoblasts in wild types can instead fate switch to become osteoblasts in mef2ca^b1086 mutants.

To test the fate switch hypothesis, that mef2ca^b1086 mutant progenitor cells adopt the osteoblast cell fate instead of the ligament cell fate, we monitored the differentiation of these two cell types by time-lapse analysis (Movie 3). In these experiments, we found that the first opercle cells initiate sp7 expression at approximately the same time, 58 h post-fertilization (hpf), and location in wild types and mutants (Fig. 3Aa,Ba, op). We observed that the opercle grows by addition of osteoblasts in the ventral direction in both wild types and mutants (Fig. 3Ab,Bb, op). However, the number of double-positive sp7;alx4a cells appeared to be reduced in the mutant, suggesting that the identity of the ligament connection site may already be lost at these early stages. Around 70 hpf, alx4a single-positive ligamentoblasts in developing wild-type animals were first detected along the track of tissue orthogonal and anterior to the op (Fig. 3Ac, arrows). Strikingly, in mef2ca^b1086 mutants, ectopic osteoblasts appeared at essentially the same time and place that the ligamentoblasts appeared in wild types (Fig. 3Bc, asterisk). Additional ectopic osteoblasts continued to appear in the mutant concurrent with the decreased appearance of ligament cells (compare Fig. 3Bd, asterisk, with 3Ad, arrow), although additional ectopic bone cells variably develop outside of the ligament track at later stages (Fig. 2D-F; Fig. 3Bd).

The appearance of osteoblasts in mutants at the time and place that ligamentoblasts appear in wild types led us to hypothesize that ectopic bone in mutants originates from progenitor cells that remain in position in the ligament track and differentiate into osteoblasts instead of ligamentoblasts. To test this hypothesis, we generated a transgenic line in which the pharyngeal arch ectomesenchyme progenitor field is labeled with nuclear-localized mCherry (fli1a:NLSmCherry) allowing observation of progenitor cells as they differentiate into osteoblasts (Movie 4). Based on these studies, we estimate that the progenitor field is similar in wild types and mutants with respect to size and number of nuclei (Fig. 4Aa,Ba) before any sp7-expressing osteoblasts were observed (Fig. 4Ca,Da). A few progenitor cells near the future opercle joint were the first to initiate the osteoblast program as indicated by sp7 transgene expression (Fig. 4Ab,Bb,Cb,Db, op). Osteoblasts were continuously added to the growing opercle by additional progenitor cells initiating sp7 expression (Fig. 4Ac,Bc,Cc,Dc). Following well-isolated progenitor cell nuclei in mutants (Fig. 4Bc,Bd, arrows) indicated that ectopic osteoblasts originate from the track of cells connecting the opercle and the ceratohyal. Alternatives to the transfating hypothesis are ectopic cell migration, or proliferation of the progenitor field. However, we detect no evidence of these, which would be observable in our time-lapse recordings with fli1a:NLSmCherry (Sasaki et al., 2013), ruling out these hypotheses. Another alternative hypothesis regarding the identity of the ectopic bone is that it is a precocious interopercle, a bone that forms anterior to the opercle later in development (Eames et al., 2013). However, when we examined the interopercle in relation to the operculohyoid ligament in wild types we observed that this bone initiates near the interhyal cartilage, far from the opercle, whereas the ectopic bone in mef2ca^b1086 mutants develops adjacent to the opercle. Moreover, we found that the interopercle develops surrounding the operculohyoid ligament (Fig. S2) unlike the ectopic bone in mutants, which is gained at the expense of the operculohyoid ligament. These data are consistent with the conclusion that the ectopic bone in mef2ca^b1086 mutants, in part, originates from fate-switched ligament cells and is not a precocious interopercle.

A substantial fraction of the variation observed in mef2ca^b1086 mutants is due to stochastic developmental noise (DeLaurier et al., 2014; see also below). However, we observed that when mef2ca^b1086 mutants develop the mild phenotype (wild-type-like opercles) this phenotype was more often bilateral than expected from chance alone (4/50 observed versus 0.72/50 expected by chance alone; P<0.0001, χ² analysis). As developmental noise (Polak, 2003) should not often yield such bilateral pairing, we hypothesized that the genetic background also influenced whether the phenotype would be mild or severe. To test this hypothesis we designed a genetic selection experiment based on the classic selective breeding strategy of family selection (progeny testing) (Falconer et al., 1996; Lush, 1935; Lynch and Walsh, 1998; Prentice, 1935). In our selection experiment, we scored mef2ca^b1086 heterozygous sibling pairs for the penetrance of fate switching in their offspring. Whereas we found that the offspring penetrance was consistent among repeat crosses of the same parental pair (Fig. S3), we discovered a great deal of variation among sibling parental pairs in the first generation, ranging from 2% to 87% penetrance (Fig. 5A, low 1 and high 1). Using this scoring system, we generated two separate selection lines: fish that produced offspring with high penetrance of the ectopic bone phenotype were bred (‘upward’ selection, which should increase the frequency of the phenotype), whereas in the other line breeding pairs that gave low penetrance were selected (‘downward’ selection). We found that we could, in a few generations, shift the distribution of offspring penetrance scores both upward and downward as our selection experiment progressed. This experiment clearly demonstrates heritability of fate-switching penetrance. To examine further the nature of heritability, we hybridized pairs of animals from high or low penetrance lines and examined the F1 hybrid generation (Fig. 5B). We found that the F1 progeny had a mean phenotype between that of the parental lines, although the distribution of penetrance in the two F1 families differed slightly and both were skewed towards higher penetrance.

What is the selectable difference between high and low penetrance strains? Because the function of the mef2ca gene is clearly important for the variable fate choice between ligament versus bone identity, we hypothesized that the expression of mef2ca might also be involved in the variable decision to become a ligament or bone cell. To test this hypothesis, we assayed mef2ca gene expression in mef2ca^b1086 mutant animals from high and low penetrance strains by in situ hybridization. We discovered that mef2ca^b1086 transcript expression appeared to be increased in mutants from the high penetrance strain compared with mutants from the low penetrance strain (Fig. 6A), an observation confirmed by quantitative transcript analyses (Fig. S4).

The severe phenotypes in animals with high quantities of the mutant transcript and mild phenotypes in individuals with low amounts of transcript motivate the hypothesis that strong expression of the mef2ca^b1086 transcript contributes to the severe phenotype in the high penetrance strain. To test this hypothesis directly, we injected mef2ca^b1086 mRNA into mutants from the low penetrance strain (mef2ca^b1086-mild) and found a significant increase in fate-switching penetrance in injected animals compared with uninjected controls (Fig. 6B, P<0.01). Conversely, we knocked down mef2ca^b1086 translation in the high penetrance strain (mef2ca^b1086-severe) with morpholino injection (Miller et al., 2007) and found a significant decrease in fate-switching penetrance in injected animals compared with uninjected controls (Fig. 6B, P<0.01). These findings support the hypothesis that variable expression of mef2ca^b1086 underlies the variable mutant phenotype. Because increased expression of the mutant allele is associated with the more severe phenotype we suggest that the mef2ca^b1086 mutation is antimorphic, as we consider further next.

That high expression of the mef2ca^b1086 allele promotes fate switching suggests that ligament-to-bone switching is a special feature of this allele. To test directly the hypothesis that high penetrance of fate switching is specific to the mef2ca^b1086 allele, we analyzed another mef2ca mutant allele (mef2ca^b631) in which the initiating methionine is destroyed (Miller et al., 2007). We discovered by pair-wise crossing individual animals from an unselected background that mutants homozygous for the mef2ca^b631 allele had overall very low penetrance of fate switching (Fig. 6C), comparable to our mef2ca^b1086 ‘mild’ selected line (Fig. 5A, downward selection strain), and that most opercles resembled wild type.

Although both the mef2ca^b631 and mef2ca^b1086 alleles are completely recessive to the wild-type allele, increased expression of the fate-switching mef2ca^b1086 allele increases the penetrance of the phenotype, characteristic of an antimorph. Because antimorphs are characteristically dominant, we wondered if there is ever a genetic context in which the fate-switching mef2ca^b1086 allele can act as a dominant. We tested the hypothesis that mef2ca^b1086 is dominant to mef2ca^b631 by generating heteroallelic combinations. In mutants heteroallelic for mef2ca^b1086 from the mild selected line (mef2ca^b1086-mild) and mef2ca^b631, the fate-switching penetrance was low, similar to either the mef2ca^b1086-mild or mef2ca^b631 homoallelic mutants. By contrast, in mutants heteroallelic for mef2ca^b1086 from the high penetrance strain (mef2ca^b1086-severe) and mef2ca^b631, the penetrance was significantly increased relative to mef2ca^b1086-mild or mef2ca^b631 homoallelic mutants (Fig. 6B, P<0.01). These results resemble those obtained by crossing the mef2ca^b1086-severe line to the mef2ca^b1086-mild line (Fig. 5B, P0) and suggest that a single copy of the mef2ca^b1086 allele from the severe strain, and any accompanying cis elements, are sufficient to promote fate switching in the mutant context. We interpret these results to indicate that high expression of the mef2ca^b1086 transcript contributes to the severe phenotype beyond a simple loss of function, functioning as an antimorph, as we consider further in the Discussion.

How can we understand heritable, but apparently noisy, variation in transcription of mef2ca? Gene expression can be influenced by nearby transposons (Rebollo et al., 2012). Furthermore, variable epigenetic silencing of transposons near endogenous genes mechanistically cause variable phenotypes by impacting local gene expression (Morgan et al., 1999; Rakyan et al., 2003). Therefore, we searched for a transposable element in the zebrafish genome near the mef2ca locus that could be mechanistically related to the variable mef2ca^b1086 mutant expression and the resulting phenotypic variation. We identified a candidate Class I long terminal repeat-containing retrotransposon annotated in the zebrafish genome (Zv9 and GRCz10) at the 5′ end of the mef2ca locus, juxtaposed to the annotated mef2ca transcriptional start site. Because this transposon, and its epigenetic status, is in a perfect position to impact mef2ca expression, we hypothesized that epigenetic silencing at this repetitive element differs between mutants from the high and low penetrance strains. To test this hypothesis, we directly assayed the DNA methylation status of this transposon by sodium bisulfite sequencing in mef2ca^b1086 mutants from the high and low penetrance strains. We found that this transposon was subject to CpG DNA methylation, and discovered that the overall amount of transposon DNA methylation was significantly reduced in the high penetrance strain compared with the low penetrance strain (P<0.01, χ² and Fisher's exact test) (Fig. 6D,E). Because the variable methylation status of a transposon juxtaposed to the mef2ca locus is associated with high versus low transcription, we propose that a fundamental basis of phenotypic variation in this system is noisy, epigenetic transposon silencing.

Developmental transitions, like cell fate choices, are times of high developmental noise because progenitor cells must destabilize their transcriptional program so that the new program can be adopted. These noisy transition states are probably modulated by the epigenome (Pujadas and Feinberg, 2012). In accordance with this, the fate-switch phenotype in mef2ca^b1086 mutants is especially susceptible to the underlying epigenetic variability. This is a different interpretation from that proposed earlier (DeLaurier et al., 2014), when we speculated that the process of morphogenesis, occurring at a later step in embryonic development, was particularly sensitive to developmental instability. It is intriguing to consider that variable ligament-to-bone fate switching in the hyoid arch might be a deeply conserved phenomenon, as craniofacial ligaments in the hyoid arch of mammals can variably mineralize, exemplified in humans by Eagle's syndrome (Piagkou et al., 2009).

We cannot conclude that mef2ca is a general ligament versus bone fate switch, as we do not observe ligaments outside pharyngeal arch 2 fate switching to bone. Rather, our data are consistent with the previously proposed model that mef2ca is involved in upstream specification of particular ectomesenchymal lineages, at least in part through its role in dorsal-ventral patterning within the pharyngeal arches, well before the stage when specific differentiated fates develop (DeLaurier et al., 2014; Kimmel et al., 2003; Miller et al., 2007). Here, we report the new discovery that a major aspect of mispatterning in mef2ca^b1086 mutants is a fate switch from ligament to bone. The newly identified ligament in question, which we named the operculohyoid, links the opercle bone to the ceratohyal cartilage, similar to the function performed in most adult teleosts by the interopercular bone as a component of the opercular four-bar lever system (Anker, 1974; Muller, 1996). This ligament has thus far escaped detection in any fish, probably owing to the absence of zebrafish ligament markers, until recently (Chen and Galloway, 2014).

Rapid phenotype enrichment in our selective breeding experiment suggests that relatively few heritable factors modify the fate switch. The pattern of inheritance we observed appears to fit that of a threshold character in which a continuous, normally unobservable variable controlled by a combination of genes, environment and developmental noise, underlies a phenotype that is easily observed: ligament or bone. In support of the threshold model, we found that family selection (progeny testing) efficiently shifts the proportion of mef2ca^b1086 mutants that display fate switching, a selection method known to be effective for threshold characters (Curnow, 1984). Further support for the threshold model comes from the analysis of our F1 hybrid animals. We find that different parental strains may appear similar in their outward phenotype (high or low penetrance) but the proportion of affected individuals can vary markedly among F1 families, consistent with classical selective breeding experiments testing threshold models (Wright, 1934). Another well-known example of threshold character selection is Waddington's genetic assimilation experiment. He revealed cryptic variation, which he could select upon by inducing environmental stress to move the threshold (Waddington, 1953). In our study, we revealed selectable cryptic variation by disrupting mef2ca function. Thus, a quantitative variable probably underlying fate switching in our system is the variation in expression of mef2ca^b1086 due to variable epigenetic silencing of the transposon at the mef2ca locus.

Our threshold model predicts that we have selected for genetic changes that affect the probability of a particular epigenetic and transcriptional state. That is, the epigenetic variation in our system is modified by genetic variation. However, our experiments cannot rule out a transgenerational epigenetic inheritance model in which we applied selection on the epialleles themselves. The best-documented example of transgenerational epigenetic inheritance is the inheritance of coat color in the agouti viable yellow (A^vy) mouse. In A^vy mice, a retrotransposon inserted upstream of the agouti locus drives ectopic agouti expression. The transposon is sensitive to silencing by DNA methylation, and the variable DNA methylation pattern results in variable agouti expression and coat color phenotypes that are heritable transgenerationally among isogenic mice (Morgan et al., 1999). The inheritance within an A^vy isogenic mouse background demonstrates transgenerational epigenetic (as opposed to genetic) inheritance. Future work recapitulating the mef2ca^b1086 lesion in clonal zebrafish strains (Streisinger et al., 1981) will allow testing for epigenetic inheritance.

Our data are consistent with the model that high expression of the mef2ca^b1086 mutant transcript, when wild-type function is absent, interferes with development beyond a simple loss of function. The mef2ca^b1086 transcript is predicted to encode a C-terminally truncated protein product consisting of essentially just the N-terminal MADS domain that mediates DNA binding and dimerization (Black and Cripps, 2010; Shore and Sharrocks, 1995). Thus, one possible mechanism for the high penetrance ligament-to-bone transformation phenotype is that the high levels of mutant transcript that we showed to be present in the high penetrance strain results in high levels of truncated protein, which has deleterious activity. We find support for this mechanism in our studies of the mef2ca^b631 allele, in which the initiating methionine is destroyed. One potential model is that the mef2ca^b631 allele is a null, because no translated protein is predicted. However, we cannot rule out the possibility that downstream internal translational start sites are functional in mef2ca^b631. In fact, another methionine is found just 11 amino acids downstream of the mef2ca^b631 mutated methionine, potentially resulting in hypomorphic function. Moreover, the idea that mef2ca^b631 functions as a hypomorph is supported by previous studies with another allele, mef2ca^tn213, and mef2ca morpholino injection, which both produce significantly higher penetrance for cartilage phenotypes than mef2ca^b631 (Miller et al., 2007). Nevertheless, the high penetrance of the severe, ligament-to-bone fate-switch phenotype appears to be specific to mef2ca^b1086 probably owing to deleterious activity of this particular allele. We find further support for the model that C-terminally truncated forms of mef2ca are detrimental in studies with mammalian cell culture, where expression of C-terminally truncated MEF2 protein fragments interferes with wild-type MEF2C activity (Leupin et al., 2007; Maiti et al., 2008; Molkentin et al., 1996; Ornatsky et al., 1997; Rao et al., 1998). We note that these dominant-negative constructs, encoding both the MADS domain and the adjacent MEF domain, are not precisely identical to mef2ca^b1086, which is predicted to encode just the MADS domain.

Given this antimorph model for the high penetrance of the mef2ca^b1086 allele, it is not clear why the allele is not inherited as a dominant. Careful examination revealed no phenotypic difference between heterozygotes and homozygous wild types. One possible explanation is that null or hypomorphic alleles can be buffered by a compensatory network (Rossi et al., 2015) or redundant loci, whereas mef2ca^b1086 cannot. In this scenario, mef2ca^b1086 interfering activity is not strong enough to overcome the normal function of wild-type mef2ca in the heterozygous condition; it is able only to disrupt the less efficient compensatory or redundant network. Precedent for such an unusual allele comes from work in Arabidopsis (Sijacic et al., 2011) in which the authors describe recessive antimorph alleles that can overcome redundant loci but not wild-type alleles. Based on our experimental evidence, we propose that mef2ca^b1086 is a recessive antimorph. With antibodies directed to the N terminus of zebrafish Mef2ca (not currently available) future studies might test directly hypotheses related to the molecular nature of inhibitory activities of the predicted truncated protein.

Recent models propose that the epigenome mediates canalization, or the ability of development to produce a consistent phenotype (Pujadas and Feinberg, 2012). Local repetitive elements and their variable silencing influence the expression of genes (Rebollo et al., 2012), and we envisage a general mechanism whereby epigenetic silencing of repetitive elements functions to canalize development by attenuating transcriptional chatter mediated by these transposons. This canalization is especially important in Mendelian mutants to attenuate expression of potentially disruptive mutant alleles such as mef2ca^b1086. Therefore, Mendelian mutants can be overlaid with complex (epi)genetic traits that contribute to mutant phenotypic variation. Understanding the phenomena that underlie mutant phenotypic variation is not only important for understanding variation in mutant laboratory animals but also human diseases in which Mendelian disorders present with variable penetrance in the clinic (Badano and Katsanis, 2002; Dipple and McCabe, 2000a; Nadeau, 2001).

All fish lines were maintained and staged according to established protocols (Kimmel et al., 1995; Westerfield, 1993). All of our work with zebrafish has been approved by the University of Oregon Institutional Animal Care and Use Committee (IACUC). Assurance number for animal research: A-3009-01. The mef2ca mutant alleles mef2ca^b1086 and mef2ca^b631 and their genotyping protocols have been previously described (Miller et al., 2007). KASP genotyping (LGC) and a StepOnePlus Real-Time PCR System (Applied Biosystems) was also used for mef2ca^b1086. See supplementary Materials and Methods for these and all other primer sequences. The following transgenic lines have been previously described: Tg(sp7:EGFP)b1212 (DeLaurier et al., 2010), sox9a:EGFP^zc81Tg (Eames et al., 2013), Tg(alx4a:DsRed2)pd52 (Nachtrab et al., 2013). To generate Tg(fli1a:NLSmCherry)b1252, the fli1a ectomesenchyme enhancer (p5E-fli1a-F-hsp70I) (Das and Crump, 2012) was recombined upstream of NLS-mCherry using the Tol2 kit (Kwan et al., 2007).

Alcian Blue and Alizarin Red stains of fixed animals and vital staining with Alizarin Red were performed as described previously (Kimmel et al., 2010; Walker and Kimmel, 2007). Whole-mount Thrombospondin 4b immunostaining was modified from Subramanian and Schilling (2014); details can be found in supplementary Materials and Methods. Whole-mount in situ hybridization was carried out as previously described (Talbot et al., 2010) with either fluorescence or immunohistochemical detection. The sp7 probe has been previously described (DeLaurier et al., 2010). Details describing mef2ca and scxa probe generation are in supplementary Materials and Methods.

Full sibling mef2ca^b1086 heterozygous adults were identified by PCR and breeding pairs were isolated and housed separately. Offspring from pairs were stained for cartilage and bone at 6 days post-fertilization (dpf) and mef2ca^b1086 mutant offspring were scored for fate-switch penetrance. Sibling pairs producing offspring with high or low penetrance were selectively bred and phenotypically wild-type animals from these crosses were raised to adulthood for the next generation and next round of selection.

DNA was isolated and bisulfite converted with the EZ DNA Methylation Direct Kit (Zymo). Bisulfite-converted transposon DNA was amplified by PCR with primers designed with MethPrimer (Li and Dahiya, 2002). PCR products were cloned using the TOPO TA cloning kit and sequenced with M13R. Analysis of bisulfite data was performed using QUMA (Kumaki et al., 2008). P-values were calculated using χ² and Fisher's exact tests.

Immunohistochemical in situ-stained animals were imaged on a Zeiss Axiophot 2. Static confocal images were captured using either a Zeiss LSM5 Pascal confocal or a Leica SD6000 spinning disk confocal. Images were assembled in ImageJ, metamorph and Photoshop with any adjustments applied to all panels. For time-lapse movies, animals were imaged as reported (Huycke et al., 2012) and skeletonized as described in supplementary Materials and Methods and Fig. S5.

RNA was extracted from a pool of whole 24 hpf embryos from the high penetrance strain using TRIzol (Thermo Fisher Scientific). cDNA was retrotranscribed using the SuperScript III Kit (Thermo Fisher Scientific). Using the mef2ca BC070007.1 sequence as a reference, we amplified cDNA containing the T7 promoter followed by the mef2ca native Kozak sequence, the mef2ca open reading frame and partial 3′UTR. Amplified fragments were cloned using pCR-Blunt II-TOPO, and mef2ca^b1086 clones were identified by Sanger sequencing. Using Phusion polymerase (NEB), a mef2ca^b1086 cDNA template was amplified from the cloned plasmid and RNA for injection was synthesized using the T7 mMessage Machine transcription kit (Thermo Fisher Scientific). RNA was purified by lithium chloride precipitation. A final RNA dose of 300 pg was injected into one-cell-stage embryos from the mef2ca^b1086-mild low penetrance strain. A final dose of 2-3 ng mef2ca mRNA translation-blocking morpholino (Miller et al., 2007) was injected into one- to four-cell-stage embryos from the mef2ca^b1086-severe high penetrance strain. Injection volumes were estimated using a stage micrometer. mef2ca transcripts were quantified by reverse transcription qPCR; see supplementary Materials and Methods for details.

We thank Thomas Desvignes for help with qPCR, Aniket Gore for help with bisulfite sequencing, and Mckenna Fairey for help scoring skeletal phenotypes. We are grateful for discussions with our colleagues especially John Postlethwait, Bill Cresko, Monte Westerfield, Anne Ferguson-Smith and Mike Miller.