Skip to main content
Download PDF

Developmental constraint through negative pleiotropy in the zygomatic arch

EvoDevo volume 9, Article number: 3 (2018) Cite this article

Abstract

Background

Previous analysis suggested that the relative contribution of individual bones to regional skull lengths differ between inbred mouse strains. If the negative correlation of adjacent bone lengths is associated with genetic variation in a heterogeneous population, it would be an example of negative pleiotropy, which occurs when a genetic factor leads to opposite effects in two phenotypes. Confirming negative pleiotropy and determining its basis may reveal important information about the maintenance of overall skull integration and developmental constraint on skull morphology.

Results

We identified negative correlations between the lengths of the frontal and parietal bones in the midline cranial vault as well as the zygomatic bone and zygomatic process of the maxilla, which contribute to the zygomatic arch. Through gene association mapping of a large heterogeneous population of Diversity Outbred (DO) mice, we identified a quantitative trait locus on chromosome 17 driving the antagonistic contribution of these two zygomatic arch bones to total zygomatic arch length. Candidate genes in this region were identified and real-time PCR of the maxillary processes of DO founder strain embryos indicated differences in the RNA expression levels for two of the candidate genes, Camkmt and Six2.

Conclusions

A genomic region underlying negative pleiotropy of two zygomatic arch bones was identified, which provides a mechanism for antagonism in component bone lengths while constraining overall zygomatic arch length. This type of mechanism may have led to variation in the contribution of individual bones to the zygomatic arch noted across mammals. Given that similar genetic and developmental mechanisms may underlie negative correlations in other parts of the skull, these results provide an important step toward understanding the developmental basis of evolutionary variation and constraint in skull morphology.

Background

The skull is a complex structure that supports and protects tissues critical for survival, including the brain, sense organs, and masticatory apparatus. While a wide range of skull morphology has evolved across mammalian species, fundamental patterns of skull bone integration are generally conserved [1,2,3,4]. This shared morphological pattern reflects conserved tissue origins [5], ossification patterns [6,7,8], and strong selective pressure for an adequately integrated and functional head [2, 9]. The same factors that reinforce the development of an integrated head can serve as developmental constraints on the directions that evolution can take [10, 11]. Understanding the genetic basis for this integration and developmental constraint is critical for illuminating how genes and development can influence processes of evolution. Here, we investigate the genetic basis for negative correlations between adjacent skull bones in a large genetically heterogeneous population of mice. These negative correlations may serve as developmental constraints that support overall integration of the head while also allowing for significant variation in the relative size of individual bones within subregions. Identifying genetic factors underlying this type of developmental constraint is important for understanding the basis of morphological integration and will illuminate critical connections between developmental and evolutionary processes.

Random pairs of linear distances across the skull are expected to display positive correlations driven by the overall growth of the integrated head. Negative correlations between raw linear distance measurements are rare and can be evidence of scale-independent negative pleiotropy, which exists when genetic variation leads to an opposite phenotypic effect on two traits [10, 12]. Within the craniofacial skeleton, a pair of negatively correlated and adjacent bone lengths provides a relatively simple system within which to search for genes underlying negative pleiotropy. Although not as commonly discussed, negative pleiotropy is fundamentally different than antagonistic pleiotropy, which occurs when a genetic factor contributes a positive fitness effect for at least one trait and a negative fitness effect for at least one other trait [13,14,15]. The concept of negative pleiotropy does not require either trait to be associated with a fitness effect. However, the existence of negative pleiotropy itself might have either positive or negative fitness effects, depending on the evolutionary context.

A previous analysis of adult skull variation indicated that linear dimensions of the cranial vault, zygomatic arch, and cranial base vary strongly among the eight inbred mouse founder strains of the Collaborative Cross (CC). Furthermore, that analysis indicated that pairs of bones within each of these regions display negative correlations in relative size [16]. Here, we explicitly tested whether pairs of adjacent linear distances within cranial vault, zygomatic arch, and cranial base are negatively correlated. In cases where they were, we performed genome-wide mapping in a large heterogeneous population of Diversity Outbred (DO) mice to identify genomic regions driving negative pleiotropy. The DO mice are derived from the CC founder strains. Each mouse carries a high degree of heterozygosity and a unique combination of alleles [17, 18], which is ideal for high-resolution genetic mapping. Although genome-wide mapping has been performed on DO mice for other phenotypes such as blood measurements and body composition [17,18,19], this is the first study that addresses craniofacial morphology.

Previous quantitative trait locus (QTL) analyses of mouse craniofacial variation have illustrated the ubiquity of pleiotropy in the genetics of craniofacial form [20,21,22] and identified candidate genes associated with variation in multivariate measures of craniofacial shape in a variety of mouse populations [23,24,25,26,27], including another outbred mouse sample [28]. The DO mice may prove more valuable than previously analyzed mouse crosses or backcrosses, because they represent a strongly genetically heterogeneous population that has been outbred for more than 10 generations, leading to relatively high genomic mapping resolution. In addition, haplotype variation at a given SNP can be tied back to the complete genomic sequence of eight diverse inbred founder strains.

Within our mapping results, we expect that a locus driving negative pleiotropy will have a significant haplotype effect on the linear distances of each adjacent bone. Furthermore, we expect that variation in founder strain haplotypes at this locus will be associated with opposite phenotypic effects for each bone (i.e., one bone gets longer when the other gets shorter). Third, if this locus represents a developmental constraint on the morphology of the combined length of both bones, we expect that local haplotype variation will have no effect on overall regional length. Identifying a genomic region associated with negative pleiotropy and candidate genes within it is a critical first step in understanding an important genetic and developmental basis for both developmental constraint and variation in the relative contribution of bones to regional skull morphology. Although our analyses are regionally specific, identifying the general mechanisms underlying negative pleiotropy of adjacent skeletal elements is an important step toward understanding the developmental basis for evolutionary change in the skull.

Methods

Adult and embryonic sample

Our measurements of the Collaborative Cross (CC) founder strains and their F1 crosses derive from craniofacial landmarks collected from 1211 specimens for a previous analysis [16]. The Diversity Outbred (DO) mice are the result of multiple generations of random outcrossing of 175 breeding pairs from partially inbred CC lines [18] at the Jackson Laboratory (JAX; Bar Harbor, ME). Our primary sample consisted of mice raised at the University of North Carolina (UNC), and at JAX. The 287 adult specimens raised at UNC were male and female sibling pairs from outbreeding generation 10 that were raised in previously described conditions [19, 29] under approval and conduced in accordance with the guidelines set forth by the Institutional Animal Care and Use Committee (IACUC) at the University of North Carolina at Chapel Hill. The 277 adult specimens raised at Jackson Labs (JAX IACUC #99066) were females of outbreeding generations 9, 10, and 15.

A sample of 472 DO mice of generations 19, 21, and 23, raised at the Scripps Research Institute (IACUC #08-0150-3), were used to validate relevant QTLs identified by genome-wide association in our primary sample. This validation sample was chosen because it includes mice from later DO generations with greater recombination, which might allow for the refinement of any validated QTL (i.e., shorter confidence intervals). After collection, adult specimens were stored at −20 °C. Receipt and imaging of specimens from other institutions was conducted in accordance with approved IACUC protocol #AC13-0268 at the University of Calgary.

Embryonic specimens from A/WySnJ (AWS), C57BL/6J (C57), and WSB/EiJ (WSB) inbred backgrounds were collected at the University of Calgary at embryonic day (E) 11.5 and processed as recently described [30]. Briefly, embryos were fixed overnight in PaxGene tissue fix solution (Qiagen, PreAnalytics, cat #765312) then stored at −20 °C in the PaxGene tissue stabilizer solution prepared to manufacturer specification (Qiagen, PreAnalytics, cat #765512). Embryonic collection and processing were conducted in accordance with approved IACUC protocol AC13-0267 at the University of Calgary.

Linear distances versus relative linear dimensions

Negative relationships between the length of bones making up the zygomatic arch, sagittal cranial vault, and posterior cranial base were previously identified among the eight CC founder strains from plots of relative linear dimensions that differ strongly between a given founder strain and several of the other founder strains [16]. Unlike standard linear distances, these relative dimensions were calculated from landmark coordinates after they were scaled to remove overall size variation and then transformed to remove the linear component of static allometry. These steps were necessary in the previous analysis of craniofacial shape between genotypes with a wide range of sizes, including small wild-derived strains and the very large New-Zealand Obese strain. However, because most raw linear distances are positively correlated with head size, removing variation that covaries with overall size can create artefactual negative associations between many of the resulting linear dimensions. For instance, two randomly selected linear distances are both strongly positively correlated with head size (centroid size) (Fig. 1a, b). They are also strongly positively correlated with each other (Fig. 1c). After scaling landmark coordinates underlying these measurements by overall skull size during Procrustes superimposition, variation that is correlated with overall scale is removed, regardless of whether that variation is mechanistically or developmentally related to overall size variation. This frequently leads to an artificial negative correlation between the resulting linear dimensions (Fig. 1d). Correlation coefficients from many linear distance pairs illustrates how scaling during Procrustes superimposition changes an asymmetric distribution of correlation coefficients between raw linear distances into a symmetric distribution of correlation coefficients centered on zero (Fig. 1e). Because this scaling procedure magnifies aspects of negative covariation, we completed our current analysis of potential developmental constraints on adjacent bones using raw linear distances to make certain that our evidence for negative relationships between these skull lengths are genuine.

Fig. 1
figure 1

Linear distance correlations and scaling. Plots of two randomly chosen linear distances calculated from our CC founder and F1 sample illustrating the change in correlation direction after scaling measures to overall skull size (centroid size). Plots of raw linear distance a L15–R22 and b R32–R2 versus skull centroid size. Plots of linear distance L15–R22 versus R32-R2 c before and d after Procrustes superimposition-based scaling (based on centroid size). e Correlation coefficients from 10,000 randomly permuted linear distance pairs calculated from raw linear distances (left) and after (right) Procrustes superimposition and scaling

Full size image

Phenotype measurement

Micro-computed tomography (µCT) images of heads were obtained in the 3D Morphometrics Centre at the University of Calgary with a Scanco vivaCT40 scanner (Scanco Medical, Brüttisellen, Switzerland) at 0.035 mm voxel dimensions at 55 kV and 145 µA. Three dimensional coordinates of 54 adult landmarks (8 midline, 46 bilateral), as previously defined [16], were collected by a single observer from minimum threshold defined bone surfaces within Analyze 3D (www.mayo.edu/bir/).

We calculated linear distances associated with zygomatic arch length, sagittal cranial vault length, and midline cranial base length from raw landmark coordinates collected on both CC founder/F1 specimens and DO specimens (Fig. 2). Full zygomatic arch length was calculated as the linear distance between landmarks L(R)3 and L(R)32, with the length of zygomatic process of the maxilla defined between L(R)3 and L(R)24 and zygomatic bone length defined between L(R)24 and L(R)32. Full cranial vault length was defined between M21 and M27, while the length of the frontal and parietal bones were defined between M21–M26 and M26–M27, respectively. Posterior cranial base length was defined between M36 and M38, while the lengths of the basioccipital and sphenoid bones were defined between M36–M37 and M37–M38, respectively. We completed our analysis of the zygomatic arch lengths on the average of the left and right sides to help control for any stochastic landmarking error or stochastic developmental variation between the left and right sides. We were unable to do this for the other measurements, because they are found along the midline of the skull.

Fig. 2
figure 2

Linear distances for analysis. Pairs of linear distances that were tested for a negative correlation within CC founder/F1 and DO samples, plotted between previously defined craniofacial landmarks [16]. a Sagittal cranial vault lengths of the frontal bone (red) and parietal bone (blue) plotted on a superior view, b posterior cranial base lengths of the basioccipital (blue) and sphenoid body (red plotted on an interior view, c zygomatic arch lengths of the zygomatic bone (blue) and zygomatic process of the maxilla (red) on a lateral view

Full size image

The Pearson’s correlation coefficient (r) for each pair of component linear distances (e.g., zygomatic bone length vs. zygomatic process of maxilla length) and between each component linear distance and the overall trait length (e.g., zygomatic bone length vs. total zygomatic arch length) was calculated to identify correlation direction and strength. A t test of whether the r is different from 0 was completed for each pair to determine whether their correlation is significant, after Bonferroni correction to account for multiple testing (α = 0.0028). The coefficient of determination (R2) is interpreted as a measure of how much variance in one length in a pair is explained by the other.

Genotyping and association mapping

Tail biopsies were taken from mice at 6 weeks of age, and DNA was either extracted from the tissue using the QIAGEN DNeasy kit per manufacturer’s instructions or sent to NeoGen GeneSeek for DNA extraction. The primary DO sample was genotyped using the MegaMUGA SNP array [GeenSeek (Neogen), Lincoln, NE] [31], while the validation sample was genotyped using the GigaMUGA array (GeenSeek (Neogen), Lincoln, NE) [32]. We used a subset of 57,977 MegaMUGA SNPs or 120,789 GigaMUGA SNPs that distinguish among the genotypes of the eight CC founder strains and their heterozygous F1 offspring and have a quality tier of 1 or 2 [32]. The probability that each of the eight founder strains contributed to a given SNP maker was calculated for each DO specimen based on array intensity values using the DOQTL package [17] within R [33].

Association mapping was performed using DOQTL [17] for both components and the overall trait length of linear distance pairs that displayed a significant negative relationship within our primary DO mouse sample. First, we completed genome scans for the primary sample using an additive haplotype model for regression of the specimen founder genotype dosage on an individual linear distance, with age at sacrifice and sex as covariates. Peaks indicating regions of the genome where haplotype covaries with a given trait were identified as those with LOD scores above a genome-wide significance threshold determined with permutation tests (1000 iterations, indicating a LOD score ~ 7.76). The support intervals under these significant peaks were identified as continuous regions including 2 LOD scores below the LOD score value of that significant peak. The results of these additive haplotype regressions also indicate which founder haplotypes are associated with positive or negative increases in linear distance length. We completed association mapping using an additive SNP-based regression model [17] across the genomic intervals of interest that were identified in the previous step. A permutation test of the SNP-based regression model across the genome was used to determine the LOD score significance threshold for these tests. Second, to validate significant peaks on chromosome 17, we completed the same association mapping steps with our validation sample for zygomatic arch distances across chromosome 17 (rather than across the whole genome).

RT-PCR

The E11.5 maxillary process was chosen for RT-PCR because both the maxilla and zygomatic bones are derived from mesenchymal condensations within this process (see “Discussion”), starting at approximately E11.5. Specifically, this age was chosen because it is the approximate time when morphogenesis and differentiation of the facial skeleton begins [34, 35]. Mouse embryos were dissected in ice cold PBS and immediately preserved using the PaxGene tissue system (Qiagen, PreAnalytics cat #765312, 765512). Embryos were fixed overnight in the fix at room temperature with rocking, then transferred to the stabilization buffer and stored at −20 °C. Maxillary processes were subsequently micro-dissected from five A/WySnJ (AWS), five C57BL/6J (C57), and six WSB/EiJ (WSB) embryos that had been stored in PaxGene tissue stabilization buffer and were stored in fresh stabilization buffer at −20 °C until extraction. AWS embryos were used in place of A/J embryos, because they were available at the University of Calgary and because the two are closely related inbred strains.

Following a recently published analysis [30], RNA was extracted using the PaxGene RNA extraction kit (Qiagen, PreAnalytics cat #766134), which includes a DNA removal step. RNA was analyzed using a NanoDrop 1000 (ThermoFisher). While RIN analysis was not performed on these samples, a similar group of samples processed during the same time span (Agilent BioAnalyzer) had RIN scores in the area of 7.9–8.3, which is in line with kit specifications. 500 ng of RNA was converted to cDNA using the Maxima First Strand Kit (ThermoFisher, Cat #K1641) in a 25 µl reaction. Real-time PCR was performed on an Applied Biosystems Quantiflex Studio 6 using standard cycling conditions with the low volume (10 µl) setting. Reactions were performed using the 2× PrimeTime gene expression mastermix from Integrated DNA Technologies with low ROX, PrimeTime assays (Gapdh—Mm.PT.39a.1, Camkmt—Mm.PT.58.7890215, Six3os1—Mm.PT.58.43925739, Six2—Mm.PT.58.22007192), an ABI Taqman assay (18s—Mm0477571_s1), and custom Six3 primers (Probe: 5′-CAAACTTCGCCGATTCTCACCACTGCT-3′, Forward primer: 5′-TCTCTATTCCTCCCACTTCTTGTTG-3′, Reverse primer: 5′-GCCGCTACTCGCCAGAAGTA-3′) [36]. Additional primer sequence details are found in Additional file 1. Normalization was done using the arithmetic average of the deltaCT from Gapdh and 18 s RNA runs. Reference genes were selected based on stability from previous experiments and RNAseq data from this region of the face.

Since C57 mice had an intermediate phenotypic effect, mean C57 RNA expression was used as the baseline upon which to compare the expression of all specimens (measured as fold change). One-way ANOVA tests of fold-change values between genotypes were completed for Camkmt, Six2, Six3, and Six3OS1 using Graphpad Prisim (Version 6) software. If there were differences in expression between genotypes, we looked for similarities between variation in RNA expression and the phenotypic effects of A-strain, C57, and WSB haplotypes on relative zygomatic bone length, which might indicate that variation in expression of these candidate genes is associated with the identified negative correlation in zygomatic bone lengths. This was done by using a post-test for linear trend.

Results

We explicitly tested whether the lengths of adjacent bones within the cranial base, cranial vault, and zygomatic arch were negatively correlated. While we expected that most linear distances in the skull would be positively correlated, a negative correlation is evidence for a developmental constraint in how component bones (e.g., frontal and parietal bones) contribute to a larger overall trait (e.g., cranial vault length). All component bone lengths were significantly positively correlated with corresponding overall trait lengths (Table 1). The lengths of component bones of the zygomatic arch were negatively correlated within the CC Founder/F1 and DO samples, while cranial vault components were negatively correlated within the DO sample. There was no evidence of a negative association between components of the posterior cranial base.

Table 1 Linear distance correlations
Full size table

Association mapping

Given the negative correlations between the length of bones contributing to the zygomatic arch and cranial vault, we completed genome-wide association mapping to look for evidence of a genomic region that might drive these correlations based on the mechanism of negative pleiotropy. Our association mapping analysis revealed two intervals on chromosome 17 that were associated with zygomatic bone length variation (42.04829–45.95447; 85.30648–85.88324 Mb), while the second of these intervals (85.30648–85.88324 Mb) was also associated with zygomatic process of the maxilla length variation (Fig. 3b, c). The phenotypic effects of founder haplotypes under the second peak were in opposite directions for the two components (Fig. 4), meaning that a founder haplotype associated with an increase in zygomatic length was also associated with a decrease in zygomatic process of the maxilla length. This closely matches our expectation for a gene underlying negative pleiotropy between two components of a larger trait. Furthermore, this interval displays a significant LOD score for both components but not for total zygomatic arch length (Fig. 3a), which meets our expectation that a gene underlying a developmental constraint on zygomatic arch morphology will have opposite effects on the relative contribution of the components to zygomatic arch length without effecting overall arch length. Chromosome specific association mapping with our validation sample confirmed the second zygomatic peak on chromosome 17 (84.37429–85.86122 for zygomatic bone; 83.71674–85.88897 for zygomatic process of maxilla), the general phenotypic effects of haplotypes under this peak (Additional file 2), and suggested a few other significant peaks on chromosome 17 related to zygomatic arch length (Fig. 5). In addition, a single significant peak on the X chromosome was noted for total cranial vault length (98.220635–100.409455) in our primary sample, but it did not reach significance in our validation sample. This vault length peak was not further pursued in this study because it did not reach significance for either cranial vault component.

Fig. 3
figure 3

Genome-wide scan results. The results of genome-wide scan using an additive haplotype model to identify genomic regions significantly associated with a full zygomatic arch length (landmark 3–33), b zygomatic process of the maxilla length (landmark 3–24), and c zygomatic bone length (landmark 24–33). The red line is a permutation-based (1000 iterations) significant LOD score of 7.76

Full size image
Fig. 4
figure 4

Founder haplotype effects. CC founder strain specific phenotype coefficients (above) and LOD scores from genome-wide scan (below) for significant association between haplotype and a zygomatic process of the maxilla length and b zygomatic bone length across Chromosome 17. Phenotypic coefficients are the effect of having a certain founder strain genotype at a specific genomic location on a measurement

Full size image
Fig. 5
figure 5

Validation sample chromosome 17 scan results. The results of chromosome 17 specific scan using an additive haplotype model to identify genomic regions significantly associated with a full zygomatic arch length (landmark 3–33), b zygomatic process of the maxilla length (landmark 3–24), and c zygomatic bone length (landmark 24–33). The red line is a permutation-based (1000 iterations) significant LOD score of 5.75

Full size image

Association mapping across the support interval of the second peak associated with zygomatic arch variation on chromosome 17 (85.30648–85.88324 Mb) indicates there are 19 known or predicted genes in this region (Fig. 6). These include three protein-coding genes (Six2, Six3, Camkmt) and one well-studied non-coding RNA (Six3os1). We noted that the WSB and A/J founder haplotypes are associated with opposite phenotypic effects for zygomatic bone and zygomatic process of the maxilla lengths, while the C57 haplotype effect is intermediate (Fig. 4). Therefore, if differences in the expression level of a protein-coding gene were responsible for founder haplotype associated variation in zygomatic bone length, we expected that A-strain and WSB expression levels would be most different, while C57 expression would be intermediate.

Fig. 6
figure 6

Association mapping results. Results of association mapping under the significant peak that is found in the genome-wide scan results for both zygomatic bone length and zygomatic process of maxilla length

Full size image

RT-PCR

To test whether the expression of these four candidate genes met this expectation, we collected maxillary prominences from embryonic day (E) 11.5 embryos of AWS, C57, and WSB mice, which had been stage matched by tail somite number. We then completed RT-PCR on these tissue samples using three replicates for each sample to quantify RNA expression levels for Camkmt, Six2, Six3, and Six3os1. All fold-change values were compared to the C57 mean as a baseline. One-way ANOVA results indicate genotype identity significantly contributes to Camkmt and Six2 RNA expression (Table 2). In both cases, a post-test for linear trends is significant when genotypes are ordered as WSB, C57, then AWS. WSB displays relatively high mean Camkmt RNA expression levels and relatively low Six2 levels. AWS displays relatively high mean Six2 levels and intermediate mean Camkmt levels (Fig. 7). No significant trends are noted for either Six3 or Six3OS1.

Table 2 ANOVA post-test for a linear trend in RNA expression, where the alternate hypothesis is that the level of expression is ordered by genotype with C57 being intermediate to the other two genotypes (AWS and WSB)
Full size table
Fig. 7
figure 7

RT-PCR results. RNA expression levels from E11.5 maxillary processes of three inbred mouse genotypes, relative to C57 expression levels

Full size image

Discussion

Within a sample of DO mice, we confirmed negative correlations between the lengths of bones contributing to the cranial vault and found strong evidence for negative pleiotropy between the length of the zygomatic bone and the zygomatic process of the maxilla. A genomic interval on chromosome 17 (85.3–85.9 Mb) met all our expectations for a genetic basis of negative correlation between these adjacent zygomatic arch bones. Specifically, there were significant and opposite haplotype effects on zygomatic and zygomatic process length (Figs. 3a, b, and 4), but no significant haplotype association with overall zygomatic arch length (Fig. 3c). This example of negative pleiotropy shows how skull shape can be conserved while the individual bone contributions to that shape can vary. The development of an integrated skull that fits together well enough for masticatory, sense organ, breathing, and brain function is necessary for survival [1, 11, 12, 37]. Just as with any skull region, the range of possible zygomatic arch phenotypes that allow for proper skull integration and function is limited. Our results illustrated a mechanism of developmental constraint that supports skull integration while allowing for variation in how specific bones contribute to a fundamentally conserved mammalian skull morphology.

Phenotypic impact of candidate region variation

There is a significant association between haplotype variation under a genomic region on mouse chromosome 17 and variation in the lengths of the zygomatic process of the maxilla and the zygomatic bone. While we are confident that there is a causal factor in this region, zygomatic arch element lengths are likely highly polygenic as is the norm for skull morphology [22, 38]. As with most skull bone lengths, the length of the zygomatic bone and the zygomatic process of the maxilla also correlate positively with skull size. In fact, the amount of variation explained by the correlation with skull size (as measured by R2) is greater than the amount of variation explained by the negative correlation between the two zygomatic components (Table 1). Additionally, the greatest difference between haplotype specific phenotypic effects in our genomic interval is about 0.4 mm (Fig. 4), which is approximately one standard deviation for these bone length measurements across our DO mice. Although the identified negative pleiotropy plays an important role in limiting overall zygomatic arch morphology, system-wide growth factors play a stronger role in determining all zygomatic arch bone lengths.

In addition to the interval displaying negative pleiotropy, another peak on chromosome 17 between 42 and 46 Mb met the genome-wide significance level for association with zygomatic bone length. With larger samples and different measurements of zygomatic arch morphology, other regions of interest would also be identified (as in [39]). Furthermore, although a mouse with the WSB haplotype under our interval of interest had, on average, a 0.3 mm shorter zygomatic bone than other DO mice (Fig. 4), WSB inbred founder mice don’t all have a shorter zygomatic bone than other inbred founder strains. This is because WSB alleles in other regions of the genome also contribute to the total WSB founder strain phenotype. Variation in even small-scale skull morphologies is produced by the combination of numerous factors, some acting globally across an organism and some acting more locally [11, 22]. However, determining the genetic factor on mouse chromosome 17 that drives negative pleiotropy within the zygomatic arch may help to reveal an important basis of developmental constraint and evolutionary change within the skull.

Candidate genes

Three protein-coding genes (Camkmt, Six3, and Six2) and one non-coding RNA with known function (Six3os1) are found under our candidate region of interest (Fig. 6). Our real-time PCR results indirectly support Camkmt and Six2 as candidate genes. Although we cannot definitively rule out other genetic factors under and near this genomic interval as candidates, we speculate that changes in the expression of at least one of these four identified factors might be responsible for the noted negative pleiotropy in zygomatic arch element length.

CAMKMT (calmodulin-lysine N-methyltransferase) is expressed across a wide range of tissues and plays a pivotal role in the methylation of calmodulin, which changes across developmental stages and varies in a tissue specific manner [40]. Deletion of a genomic region including CAMKMT in humans has been associated with micrognathia, dolichocephaly, and cleft palate, although the specific loss of CAMKMT has been associated with intellectual disability and muscle fiber abnormalities instead of these craniofacial phenotypes [41]. Since CAMKMT regulates calmodulin (CaM) function, it is also important to note that calmodulin has been linked to variation in beak length in Darwin’s finches and chicks [42]. CAMKMT is critical to basic physiological function across the body and has tentatively been associated with severe craniofacial birth defects.

SIX2 (sine oculus-related homeobox 2) is known to play a significant role in the skeletal development of pharyngeal arch derivatives. SIX2 is upregulated in neural crest-derived cranial mesenchyme in mice at E8.5 [43] and E9.5 [44], becoming localized to mesenchymal cells of nasal prominences, midline, and developing skull vault, as well as olfactory epithelium by E11.5 [43, 45, 46]. Downregulation of SIX2 can lead to loss or reduction in facial skeletal elements [47], reduced cranial base length, and cleft palate [48]. Later in development, SIX2 loss has been linked to increased rates of cartilage replacement by bone during endochondral ossification of the presphenoid, leading to an abnormal cranial base morphology [49]. SIX2 regulates the formation of bones from the first pharyngeal arch, which includes both the maxilla and zygomatic bones.

SIX3 (sine oculus-related homeobox 3) interacts with BMP, WNT, and NODAL, is critical during eye development [50,51,52], and during anterior neural plate specification [53]. Although SIX3 is not expressed in the maxillary arch or other craniofacial mesenchyme during embryonic development [44], genetic variation in Six3 can lead to craniofacial dysmorphology of maxillofacial bones [54], particularly in the facial midline. Mutations in SIX3 can cause holoprosencephaly in humans, a condition associated with forebrain malformation, intellectual disability, ophthalmological abnormalities, and craniofacial features including cyclopia, nasal dysmorphology, and cleft lip/palate [55]. A combination of forebrain loss and modified ectoderm/mesenchyme interactions may underlie the associated craniofacial dysmorphology (e.g., [56]). SIX3OS1 (Six3 opposite strand) is an antisense non-coding RNA that is independently coexpressed with SIX3 in the forebrain and eye after E8.5 in mice [57]. SIX3OS1 likely acts as a transcriptional scaffold for SIX3 and modulates the ability of Six3 to regulate expression of target genes in retinal cells [58]. SIX3 and SIX3OS1 may indirectly regulate facial bone development.

Although all four candidate genes have been previously associated with craniofacial dysmorphology, Six2 is particularly tantalizing because it is a major player in facial bone ossification and is associated with RNA expression level differences between founder strain maxillary prominences. Assuming that variation in Six2 or one of its cis-regulatory factors is responsible, we speculate about how variation in SIX2 expression might modify development to produce zygomatic arch variation in our DO mice and across mammalian clades.

Developmental mechanisms

All bones that contribute to the mouse zygomatic arch form from neural crest-derived mesenchyme [59, 60] within the first pharyngeal arch [61,62,63]. The maxilla and zygomatic bone derive from neural crest cells that migrate from the posterior mesencephalic region, while much of the squamous temporal neural crest mesenchyme probably originates in the first couple of rhombomeres [5, 64]. Our results indicate that a gene or regulatory element within our candidate interval determines the location of the border between the maxilla and zygomatic within the adult zygomatic arch. Because the zygomatic bone is the last remaining dermally ossified circumorbital bone in mammals [6, 65] and appears to be the only dermatocranial element that does not develop in proximity to chondrocranial elements [66], it may ossify in response to a different signal than the maxilla. Although Dlx genes and related factors are associated with determining upper and lower jaw fates within the first pharyngeal arch [67], what determines the individual bone fate of mesenchymal cells within the maxillary portion of the first pharyngeal arch remains unknown.

Heuzé and colleagues [68] recently suggested that the zygomatic progenitor cells experience developmental cues analogous to Dlx expression patterns that distinguish them from maxillary cell populations. It is possible that one of these cues is the spatiotemporal pattern of Six2 expression. Supporting this perspective, loss of SIX1 results in a partial transformation of the zygomatic process of the maxilla into a mandible, possibly through downstream effects on Dlx gene expression [69]. This partial morphological transformation includes a significant increase in the length and volume of the zygomatic process of the maxilla and a significant reduction of the zygomatic bone or the fusion of zygomatic and maxillary portions of the zygomatic arch together [69]. As two members of the same gene family that are expressed within the developing pharyngeal arches, it is possible that SIX1 and SIX2 regulate development through similar pathways and mechanisms.

If a mutation in SIX2 leads to a change in how regional segmentation genes like the Dlx genes are expressed, it is possible that the variation in the contribution of maxilla and zygomatic to the zygomatic arch (Fig. 8a) may be based on a change in the location of a regulatory gene expression border between their presumptive cell populations (Fig. 8b). Tissue boundary definition is critical in many developmental contexts [70]. One well-documented craniofacial example of a gene regulatory border prevents mesenchymal and osteogenic cells from crossing the presumptive coronal suture between neural crest and mesodermally derived mesenchymal cells [70,71,72] (although, see [73]). A change in the relative location of a gene regulatory boundary may serve to define the final location of the zygomaticomaxillary suture.

Fig. 8
figure 8

Potential developmental mechanisms underlying zygomatic arch variation. Schematic models of developmental mechanisms by which an a adult phenotype of relative skeletal contribution to the zygomatic arch might occur. Options include b external regulatory definition of a presumptive suture, c differences in initial mesenchymal condensation size, and d differences in relative growth or ossification rates between the two bone primordia

Full size image

There are other ways that a mutation in SIX2 might modify developmental processes to lead the measured zygomatic arch length variation. The site of the zygomaticomaxillary suture may not be defined prior to osteogenesis but may simply occur wherever the growing bones meet. Six2 expression has been associated with increased mesenchymal cell proliferation in the developing head and renal system [48, 49, 74]. Recent results indicate that Six2 mRNA and protein levels are highest in palatal tissues during the period of initial palatal shelf outgrowth and suggest that later spatiotemporal expression patterns are responsible for local increases in mesenchymal cell proliferation [48]. It is possible that genetic variation under our candidate region leads to a change in the timing, location, or level of mesenchymal precursor cell populations.

A change in proliferation within either the maxillary or zygomatic mesenchymal condensations may result in size variation of that condensation and the resulting bones (Fig. 8c). Within chicken eyes, it has been shown that the largest and most widely spaced intramembranously ossified scleral ossicles tend to be those with the earliest forming mesenchymal condensation precursors. In addition, if one ossicle fails to form, the adjacent ossicles fill in the extra space [75]. In another relevant example, a Fuz ciliopathy mutation leads to the formation of a single frontal bone pair at the expense of parietal bones in mice, perhaps because an excess proliferation of precursor cells leads to an unusually wide frontal bone mesenchymal condensation [76]. It is possible that a similar change in the proliferation rate within the maxillary or zygomatic bone mesenchymal condensations may cause them to become larger at the expense of the other. While Six2 is a tantalizing candidate gene that may drive the negative pleiotropy in zygomatic arch bone lengths among DO mice, further work is required to confirm this.

Evolutionary implications

The jugal bone (homologous with zygomatic bone) is first noted in the fossil record as one of the circumorbital bones within the dermal skeleton of agnathans. Within the presumed ancestral tetrapod, the jugal was a narrow bone of the inferior orbital margin that articulated with facial bones including the lacrimal, maxilla, and squamosal (reviewed by [65, 68]). The zygomatic bone of the last common ancestor of mammals was likely a linear bone connecting the maxilla and squamosal bones that lacked a postorbital connection to the frontal bone. This view is supported by the existence of similar morphologies in living monotremes, marsupials and many other mammalian clades. Among mammals, zygomatic arches vary in their width, breadth, height, length, degree and direction of curvature, among other characteristics. A complete postorbital border between frontal and zygomatic bones is noted in some clades (e.g., cervids, equids, primates), is sometimes represented by zygomatic and frontal processes that do not touch (e.g., carnivorans, lagomorphs), and completely absent in other mammals (e.g., bats, insectivores, rodents) [77,78,79]. Arch thickness and robusticity may vary in response to selective pressure based on mechanical loading requirements [80,81,82,83], with robust arches found in beaver and rabbits, thinner arches in mice and moles, and a practically complete loss in shrews.

Although a wide range of morphologies occur across mammals, our results specifically relate to variation in the relative contribution of bones to zygomatic arch length between the main bodies of the maxilla and squamosal bones. The zygomatic process of the maxilla, the zygomatic bone, and the zygomatic process of the squamous temporal typically contribute to form the zygomatic arch. Because only a short squamous temporal portion exists in laboratory mice, our analysis focused on the zygomatic and maxillary bone contributions to overall arch length. Mouse strain haplotype variation noted under a genomic interval on mouse chromosome 17 (85.3–85.9 Mb) is responsible for significant variation in the relative contribution of these bones to arch length.

Analogous genetic variation may underlie variation in the relative contribution of zygomatic arch elements to total arch length across mammalian species. This sort of phenotypic variation is common across mammalian taxa. For example, carnivorans including felids and canids typically have a very short zygomatic process of the maxilla, an anteriorly placed zygomatic bone, and a relatively long zygomatic process of the temporal [78, 84] (Fig. 9a). On the other hand, rodents typically have a long zygomatic process of the maxilla, a more posteriorly placed zygomatic bone, and a short zygomatic process of the temporal [78, 85] (Fig. 9b). We propose that the identified pattern of negative pleiotropy contributes to this mammalian variation in relative bone length, but not that it explains all variation in zygomatic arch length or shape.

Fig. 9
figure 9

Mammal zygomatic contributions. a snapshot of mammalian variation in the relative contributions of zygomatic arch bones, including a feral domestic cat (Felis silvestris catus), b a wild caught common rat (Rattus norvegicus), and c a European mole (Talpa europae). Outlines based on images found on DigiMorph (digimorph.org). Colored lines representing the relative length of the zygomatic process of the maxilla (red), the zygomatic bone (blue), and the zygomatic process of the squamous temporal (green) match the colors noted in Figs. 2 and 8

Full size image

Two clades that were commonly thought to lack a zygomatic bone illustrate how an extreme imbalance in two zygomatic arch elements might manifest. First, moles (Talpidae) may represent the logical extreme of the proposed mechanism. Although moles lack a separate zygomatic bone as adults [78, 86] (Fig. 9c), one of the multiple small zygomatic arch ossification centers [77] may represent the zygomatic bone [87] and fuse with the zygomatic process of the maxilla quite early in development. The mole zygomatic arch is a complete arch, with a minor zygomatic bone contribution, and a lack of zygomatic arch sutures. The fossil order of Multiturberculata was long considered unique as mammals with a robust zygomatic arch, but lacking a zygomatic bone. However, careful work by Hopson and colleagues [88] indicated the zygomatic bone is found medial to the maxillary and temporal portions of the arch rather than in between.

If the modified expression of a factor like SIX2 leads to the expansion of one arch bone at the expense of another, a fusion of the arch bones together and/or a displacement of the smaller bone may occur. In fact, SIX1 null mutant mice with the enlarged zygomatic process of the maxilla also display either a smaller displaced zygomatic bone or fusion of zygomatic and maxillary portions of the developing arch [69]. It is temping to speculate that this reduction and fusion of the zygomatic bone might serve as a foundation for total zygomatic bone loss. In moles, where the existence of a zygomaticomaxillary suture is not functionally necessary, a complete loss of the zygomatic bone ossification center would result in the same morphology as long as the growing temporal and maxillary bones expanded further to fill in the gap as occurs among chick scleral ossicles [75]. However, in the case of the Multiturberculata, the continued presence of the displaced zygomatic bone may serve to reinforce the zygomatic arch response to mechanical forces [88]. Regardless of the responsible factor, similar changes in developmental processes [89] may underlie the variation in zygomatic arch bone contributions among our DO mice and variation that has arisen during mammalian evolution.

Conclusions

Investigating how genetic factors constrain the range of possible skull variation is critical for identifying mechanisms of integration and developmental constraint. Out of three cranial regions studied, we identified a genomic region underlying negative pleiotropy within the zygomatic arch. Association mapping and subsequent RT-PCR analysis identified candidate genes that might underlie this pattern. Further study is required to determine how the responsible genetic factor modifies developmental processes to limit overall zygomatic arch length while allowing for variation in the relative length of contributing zygomatic bones. This pattern of negative pleiotropy may have contributed to the evolution of mammalian zygomatic arch diversity. Identifying the particular developmental basis for this negative pleiotropy has implications beyond the zygomatic arch, because changes in similar instances of negative pleiotropy underlie significant evolutionary variation in the relative contribution of adjacent bones to larger morphological features in other regions like the cranial vault and upper jaw. These results provide a significant toe-hold in unraveling an example of negative pleiotropy and developmental constraint, which is an important step toward understanding the developmental basis for evolutionary change in the skull.

References

  1. Porto A, de Oliveira FB, Shirai LT, De Conto V, Marroig G. The evolution of modularity in the mammalian skull I: morphological integration patterns and magnitudes. Evol Biol. 2009;36:118–35.

    Article  Google Scholar 

  2. Jamniczky HA, Hallgrímsson B. A comparison of covariance structure in wild and laboratory muroid crania. Evolution. 2009;63:1540–56.

    Article  PubMed  Google Scholar 

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

    Article  PubMed  Google Scholar 

  4. Steppan SJ. Phylogenetic analysis of phenotypic covariance structure. I. Contrasting results from matrix correlation and common principal component analyses. Evolution. 1997;51:571–86.

    PubMed  Google Scholar 

  5. Gross JB, Hanken J. Review of fate-mapping studies of osteogenic cranial neural crest in vertebrates. Dev Biol. 2008;317:389–400.

    Article  CAS  PubMed  Google Scholar 

  6. Schoch RR. Skull ontogeny: developmental patterns of fishes conserved across major tetrapod clades. Evol Dev. 2006;8:524–36.

    Article  PubMed  Google Scholar 

  7. Sánchez-Villagra MR, Goswami A, Weisbecker V, Mock O, Kuratani S. Conserved relative timing of cranial ossification patterns in early mammalian evolution. Evol Dev. 2008;10:519–30.

    Article  PubMed  Google Scholar 

  8. Koyabu D, Werneburg I, Morimoto N, Zollikofer CP, Forasiepi AM, Endo H, et al. Mammalian skull heterochrony reveals modular evolution and a link between cranial development and brain size. Nat Commun. 2014;5:1–9.

    Article  Google Scholar 

  9. Young NM, Hu D, Lainoff AJ, Smith FJ, Diaz R, Tucker AS, et al. Embryonic bauplans and the developmental origins of facial diversity and constraint. Development. 2014;141:1059–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Cheverud JM. Quantitative genetics and developmental constraints on evolution by selection. J Theor Biol. 1984;110:155–71.

    Article  CAS  PubMed  Google Scholar 

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

    Article  PubMed  PubMed Central  Google Scholar 

  12. Cheverud JM. Developmental integration and the evolution of pleiotropy. Am Zool. 1996;36:44–50.

    Article  Google Scholar 

  13. Williams GC. Pleiotropy, natural selection, and the evolution of senescence. Evolution. 1957;11:398–411.

    Article  Google Scholar 

  14. Williams PD, Day T. Antagonistic pleiotropy, mortality source interactions, and the evolutionary theory of senescence. Evolution. 2003;57:1478–88.

    Article  PubMed  Google Scholar 

  15. Carter AJ, Nguyen AQ. Antagonistic pleiotropy as a widespread mechanism for the maintenance of polymorphic disease alleles. BMC Med Genet. 2011;12:160.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Percival CJ, Liberton DK, Pardo-Manuel de Villena F, Spritz R, Marcucio R, Hallgrímsson B. Genetics of murine craniofacial morphology: diallel analysis of the eight founders of the Collaborative Cross. J Anat. 2016;228:96–112.

    Article  PubMed  Google Scholar 

  17. Gatti DM, Svenson KL, Shabalin A, Wu L-Y, Valdar W, Simecek P, et al. Quantitative trait locus mapping methods for diversity outbred mice. G3: genes| Genomes|. Genetics. 2014;4:1623–33.

    Google Scholar 

  18. Svenson KL, Gatti DM, Valdar W, Welsh CE, Cheng R, Chesler EJ, et al. High-resolution genetic mapping using the mouse diversity outbred population. Genetics. 2012;190:437–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Smallwood TL, Gatti DM, Quizon P, Weinstock GM, Jung K-C, Zhao L, et al. High-resolution genetic mapping in the diversity outbred mouse population identifies Apobec1 as a candidate gene for atherosclerosis. G3: genes| Genomes|. Genetics. 2014;4:2353–63.

    Google Scholar 

  20. Kenney-Hunt JP, Wang B, Norgard EA, Fawcett G, Falk D, Pletscher LS, et al. Pleiotropic patterns of quantitative trait loci for 70 murine skeletal traits. Genetics. 2008;178:2275–88.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Leamy LJ, Routman EJ, Cheverud JM. Quantitative trait loci for early-and late-developing skull characters in mice: a test of the genetic independence model of morphological integration. Am Nat. 1999;153:201–14.

    Google Scholar 

  22. Wolf JB, Leamy LJ, Routman EJ, Cheverud JM. Epistatic pleiotropy and the genetic architecture of covariation within early and late-developing skull trait complexes in mice. Genetics. 2005;171:683–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Burgio G, Baylac M, Heyer E, Montagutelli X. Exploration of the genetic organization of morphological modularity on the mouse mandible using a set of interspecific recombinant congenic strains between C57BL/6 and mice of the Mus spretus species. G3: genes |Genomes|. Genetics. 2012;2:1257–68.

    Google Scholar 

  24. Burgio G, Baylac M, Heyer E, Montagutelli X. Nasal bone shape is under complex epistatic genetic control in mouse interspecific recombinant congenic strains. PLoS One. 2012;7(5):e37721. https://doi.org/10.1371/journal.pone.0037721.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Burgio G, Baylac M, Heyer E, Montagutelli X. Genetic analysis of skull shape variation and morphological integration in the mouse using interspecific recombinant congenic strains between C57BL/6 and mice of the mus spretus species. Evolution. 2009;63:2668–86.

    Article  PubMed  Google Scholar 

  26. Maga AM, Navarro N, Cunningham ML, Cox TC. Quantitative trait loci affecting the 3D skull shape and size in mouse and prioritization of candidate genes in-silico. Front Physiol. 2015;6:92. https://doi.org/10.3389/fphys.2015.00092.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Pallares LF, Harr B, Turner LM, Tautz D. Use of a natural hybrid zone for genomewide association mapping of craniofacial traits in the house mouse. Mol Ecol. 2014;23:5756–70.

    Article  CAS  PubMed  Google Scholar 

  28. Pallares LF, Carbonetto P, Gopalakrishnan S, Parker CC, Ackert-Bicknell CL, Palmer AA, et al. Mapping of craniofacial traits in outbred mice identifies major developmental genes involved in shape determination. PLoS Genet. 2015;11:1–25.

    Article  Google Scholar 

  29. Shorter J, Huang W, Beak JY, Hua K, Pardo-Manuel de Villena F, Pomp D, et al. Quantitative trait mapping in Diversity Outbred mice identifies two genomic regions associated with heart size. Mamm Genome. 2017. https://doi.org/10.1007/s00335-017-9730-7.

    PubMed  Google Scholar 

  30. Green RM, Leach CL, Hoehn N, Marcucio RS, Hallgrímsson B. Quantifying three-dimensional morphology and RNA from individual embryos. Dev Dyn. 2017;246:431–6.

    Article  CAS  PubMed  Google Scholar 

  31. Welsh CE, Miller DR, Manly KF, Wang J, McMillan L, Morahan G, et al. Status and access to the Collaborative Cross population. Mamm Genome. 2012;23:706–12.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Morgan AP, Fu C-P, Kao C-Y, Welsh CE, Didion JP, Yadgary L, et al. The mouse universal genotyping array: from substrains to subspecies. G3: genes |Genomes|. Genetics. 2016;6:263–79.

    Google Scholar 

  33. R Developmental Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2008. http://www.R-project.org.

  34. Kaufman MH. The atlas of mouse development. London: Academic Press; 1992.

    Google Scholar 

  35. Miyake T, Cameron AM, Hall BK. Detailed staging of inbred C57BL/6 mice between Theiler’s [1972] stages 18 and 21 (11–13 days of gestation) based on craniofacial development. J Craniofac Genet Dev Biol. 1996;16:1–31.

    CAS  PubMed  Google Scholar 

  36. Geng X, Acosta S, Lagutin O, Gil HJ, Oliver G. Six3 dosage mediates the pathogenesis of holoprosencephaly. Development. 2016;143:4462–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  Google Scholar 

  38. Mackay TF. The genetic architecture of quantitative traits. Annu Rev Genet. 2001;35:303–39.

    Article  CAS  PubMed  Google Scholar 

  39. Wood AR, Esko T, Yang J, Vedantam S, Pers TH, Gustafsson S, et al. Defining the role of common variation in the genomic and biological architecture of adult human height. Nat Genet. 2014;46:1173–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Magen S, Magnani R, Haziza S, Hershkovitz E, Houtz R, Cambi F, et al. Human calmodulin methyltransferase: expression, activity on calmodulin, and Hsp90 dependence. PLoS ONE. 2012;7:e52425.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Bartholdi D, Asadollahi R, Oneda B, Schmitt-Mechelke T, Tonella P, Baumer A, et al. Further delineation of genotype–phenotype correlation in homozygous 2p21 deletion syndromes: first description of patients without cystinuria. Am J Med Genet Part A. 2013;161:1853–9.

    Article  CAS  Google Scholar 

  42. Abzhanov A, Kuo WP, Hartmann C, Grant BR, Grant PR, Tabin CJ. The calmodulin pathway and evolution of elongated beak morphology in Darwin’s finches. Nature. 2006;442:563–7.

    Article  CAS  PubMed  Google Scholar 

  43. Fogelgren B, Kuroyama MC, McBratney-Owen B, Spence AA, Malahn LE, Anawati MK, et al. Misexpression of Six2 is associated with heritable frontonasal dysplasia and renal hypoplasia in 3H1 Br mice. Dev Dyn. 2008;237:1767–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Brunskill EW, Potter AS, Distasio A, Dexheimer P, Plassard A, Aronow BJ, et al. A gene expression atlas of early craniofacial development. Dev Biol. 2014;391:133–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Ohto H, Takizawa T, Saito T, Kobayashi M, Ikeda K, Kawakami K. Tissue and developmental distribution of Six family gene products. Int J Dev Biol. 2002;42:141–8.

    Google Scholar 

  46. Brodbeck S, Besenbeck B, Englert C. The transcription factor Six2 activates expression of the Gdnf gene as well as its own promoter. Mech Dev. 2004;121:1211–22.

    Article  CAS  PubMed  Google Scholar 

  47. Garcez RC, Le Douarin NM, Creuzet SE. Combinatorial activity of Six1-2-4 genes in cephalic neural crest cells controls craniofacial and brain development. Cell Mol Life Sci. 2014;71:2149–64.

    CAS  PubMed  Google Scholar 

  48. Okello DO, Iyyanar PPR, Kulyk WM, Smith TM, Lozanoff S, Ji S, et al. Six2 plays an intrinsic role in regulating proliferation of mesenchymal cells in the developing palate. Front Physiol. 2017;8:1–11.

    Article  Google Scholar 

  49. He G, Tavella S, Hanley KP, Self M, Oliver G, Grifone R, et al. Inactivation of Six2 in mouse identifies a novel genetic mechanism controlling development and growth of the cranial base. Dev Biol. 2010;344:720–30.

    Article  CAS  PubMed  Google Scholar 

  50. Del Bene F, Tessmar-Raible K, Wittbrodt J. Direct interaction of geminin and Six3 in eye development. Nature. 2004;427:745–9.

    Article  PubMed  Google Scholar 

  51. Lang RA. Pathways regulating lens induction in the mouse. Int J Dev Biol. 2004;48:783–91.

    Article  CAS  PubMed  Google Scholar 

  52. Liu W, Lagutin OV, Mende M, Streit A, Oliver G. Six3 activation of Pax6 expression is essential for mammalian lens induction and specification. EMBO J. 2006;25:5383–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Gestri G, Carl M, Appolloni I, Wilson SW, Barsacchi G, Andreazzoli M. Six3 functions in anterior neural plate specification by promoting cell proliferation and inhibiting Bmp4 expression. Development. 2005;132:2401–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Lagutin OV, Zhu CC, Kobayashi D, Topczewski J, Shimamura K, Puelles L, et al. Six3 repression of Wnt signaling in the anterior neuroectoderm is essential for vertebrate forebrain development. Genes Dev. 2003;17:368–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Lacbawan F, Solomon BD, Roessler E, El-Jaick K, Domené S, Velez JI, et al. Clinical spectrum of Six3-associated mutations in holoprosencephaly: correlation between genotype, phenotype, and function. J Med Genet. 2009. https://doi.org/10.1136/jmg.2008.063818.

    Google Scholar 

  56. Hu D, Marcucio RS. A SHH-responsive signaling center in the forebrain regulates craniofacial morphogenesis via the facial ectoderm. Development. 2009;136:107–16.

    Article  CAS  PubMed  Google Scholar 

  57. Geng X, Lavado A, Lagutin OV, Liu W, Oliver G. Expression of Six3 Opposite Strand (Six3OS) during mouse embryonic development. Gene Expr Patterns. 2007;7:252–7.

    Article  CAS  PubMed  Google Scholar 

  58. Rapicavoli NA, Poth EM, Zhu H, Blackshaw S. The long noncoding RNA c acts in trans to regulate retinal development by modulating Six3 activity. Neural development. 2011;6:1.

    Article  Google Scholar 

  59. Jiang X, Iseki S, Maxson RE, Sucov HM, Morriss-Kay GM. Tissue origins and interactions in the mammalian skull vault. Developmental Biology. 2002;241:106–16.

    Article  CAS  PubMed  Google Scholar 

  60. Noden DM, Trainor PA. Relations and interactions between cranial mesoderm and neural crest populations. J Anat. 2005;207:575–601.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Depew MJ, Tucker AS, Sharpe PT. Craniofacial development. In: Rossant J, Tam PPL, editors. Mouse development, patterning, morphogenesis, and organogenesis. San Diego: Academic Press; 2002. p. 421–98.

    Google Scholar 

  62. Cerny R, Lwigale P, Ericsson R, Meulemans D, Epperlein H-H, Bronner-Fraser M. Developmental origins and evolution of jaws: new interpretation of “maxillary” and “mandibular”. DevBiol. 2004;276:225–36.

    CAS  Google Scholar 

  63. Lee S-H, Bédard O, Buchtová M, Fu K, Richman JM. A new origin for the maxillary jaw. Dev Biol. 2004;276:207–24.

    Article  CAS  PubMed  Google Scholar 

  64. Santagati F, Rijli FM. Cranial neural crest and the building of the vertebrate head. Nat Rev Neurosci. 2003;4:806–18.

    Article  CAS  PubMed  Google Scholar 

  65. Gai Z, Yu X, Zhu M. The evolution of the zygomatic bone from Agnatha to Tetrapoda. Anat Rec. 2017;300:16–29.

    Article  Google Scholar 

  66. Kawasaki K, Richtsmeier JT. Association of the chondrocranium and dermatocranium in early skull formation. In: Percival CJ, Richtsmeier JT, editors. Building bones: bone development and formation in anthropology. 1st ed. Cambridge: Cambridge University Press; 2017. p. 52–78.

    Chapter  Google Scholar 

  67. Jeong J, Li X, McEvilly RJ, Rosenfeld MG, Lufkin T, Rubenstein JL. Dlx genes pattern mammalian jaw primordium by regulating both lower jaw-specific and upper jaw-specific genetic programs. Development. 2008;135:2905–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Heuzé Y, Kawasaki K, Schwarz T, Schoenebeck JJ, Richtsmeier JT. Developmental and evolutionary significance of the zygomatic bone. Anat Rec. 2016;299:1616–30.

    Article  Google Scholar 

  69. Tavares AL, Cox TC, Maxson RM, Ford HL, Clouthier DE. Negative regulation of Endothelin signaling by SIX1 is required for proper maxillary development. Development. 2017;144:2021–31.

    Article  CAS  PubMed  Google Scholar 

  70. Yen H-Y, Ting M-C, Maxson RE. Jagged1 functions downstream of Twist1 in the specification of the coronal suture and the formation of a boundary between osteogenic and non-osteogenic cells. Dev Biol. 2010;347:258–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Merrill AE, Bochukova EG, Brugger SM, Ishii M, Pilz DT, Wall SA, et al. Cell mixing at a neural crest-mesoderm boundary and deficient ephrin-Eph signaling in the pathogenesis of craniosynostosis. Hum Mol Genet. 2006;15:1319–28.

    Article  CAS  PubMed  Google Scholar 

  72. Ting M-C, Wu NL, Roybal PG, Sun J, Liu L, Yen Y, et al. EphA4 as an effector of Twist1 in the guidance of osteogenic precursor cells during calvarial bone growth and in craniosynostosis. Development. 2009;136:855–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Holmes G, Basilico C. Mesodermal expression of Fgfr2S252W is necessary and sufficient to induce craniosynostosis in a mouse model of Apert syndrome. Dev Biol. 2012;368:283–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Self M, Lagutin OV, Bowling B, Hendrix J, Cai Y, Dressler GR, et al. Six2 is required for suppression of nephrogenesis and progenitor renewal in the developing kidney. EMBO J. 2006;25:5214–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Jabalee J, Hillier S, Franz-Odendaal T. An investigation of cellular dynamics during the development of intramembranous bones: the scleral ossicles. J Anat. 2013;223:311–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Tabler J, Rice CP, Liu KJ, Wallingford J. A novel ciliopathic skull defect arising from excess neural crest. Dev Biol. 2016;417:4–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. De Beer GR. The development on the vertebrate skull. Chicago: The University of Chicago Press; 1985.

    Google Scholar 

  78. Elbroch M. Animal skulls: a guide to North American species. Mechanicsburg: Stackpole Books; 2006.

    Google Scholar 

  79. Márquez S, Pagano AS, Schwartz JH, Curtis A, Delman BN, Lawson W, et al. Toward understanding the mammalian zygoma: insights from comparative anatomy, growth and development, and morphometric analysis. Anat Rec. 2017;300:76–151.

    Article  Google Scholar 

  80. Herring SW, Teng S, Huang X, Mucci RJ, Freeman J. Patterns of bone strain in the zygomatic arch. Anat Rec. 1996;246:446–57.

    Article  CAS  PubMed  Google Scholar 

  81. Witzel U, Preuschoft H, Sick H. The role of the zygomatic arch in the statics of the skull and its adaptive shape. Folia Primatol. 2004;75:202–18.

    Article  CAS  PubMed  Google Scholar 

  82. Dumont ER, Piccirillo J, Grosse IR. Finite-element analysis of biting behavior and bone stress in the facial skeletons of bats. Anat Rec. 2005;283:319–30.

    Article  Google Scholar 

  83. Franks EM, Holton NE, Scott JE, McAbee KR, Rink JT, Pax KC, et al. Betwixt and between: intracranial perspective on zygomatic arch plasticity and function in mammals. Anat Rec. 2016;299:1646–60.

    Article  Google Scholar 

  84. Gilbert M. Mammalian osteology. Springfield: Missouri Archaeological Society; 1980.

    Google Scholar 

  85. Popesko P, Rajtová V, Horák J. The colour atlas of anatomy of small laboratory animals. London: Saunders; 1992.

    Google Scholar 

  86. Sánchez-Villagra MR, Horovitz I, Motokawa M. A comprehensive morphological analysis of talpid moles (Mammalia) phylogenetic relationships. Cladistics. 2006;22:59–88.

    Article  Google Scholar 

  87. Goswami A, Prochel J. Ontogenetic Morphology and Allometry of the Cranium in the Common European Mole (Talpa europaea). J Mammal. 2007;88:667–77.

    Article  Google Scholar 

  88. Hopson JA, Kielan-Jaworowska Z, Allin EF. The cryptic jugal of multituberculates. J Vertebr Paleontol. 1989;9:201–9.

    Article  Google Scholar 

  89. Hallgrímsson B, Lieberman DE. Mouse models and the evolutionary developmental biology of the skull. Integr Comp Biol. 2008;48:373–84.

    Article  PubMed  Google Scholar 

Download references

Authors’ contributions

CJP contributed to study design, collected computed tomography images, completed morphometric and genome-wide mapping analysis, and was the primary author of the manuscript. RG collected embryonic tissue, completed RT-PCR analysis, and contributed substantially to manuscript. CCR contributed to study design, the interpretation of statistical results, and critically revised the manuscript. DMG performed haplotype probability calculations necessary for genome-wide mapping and critically revised the manuscript. JMM, KMP, KH, DP, JLM, SAM, and LRD contributed to study design, mouse management, sample collection, genotyping, and critically revised the manuscript. RM and BH contributed to study design and critically revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank Wei Liu for completing some of the computed tomography imaging and all craniofacial landmark collection at the University of Calgary. Other phenotypes of some of our mice were collected using the Animal Metabolism Phenotyping core facility within UNC’s Nutrition and Obesity Research Center funded by NIH DK056350. We thank Liyang Zhao and Kuo-Chen Jung for assistance with UNC mouse experiments. We thank Robyn Humphreys for assistance with embryo collection at the University of Calgary. Charles Farber is gratefully acknowledged for assistance with processing of mouse carcasses.

Competing interests

DMG, JLM, SAM, and LRD are employed by The Jackson Laboratory, which sells DO mice like those that make up the sample analyzed here. All other authors declare that they have no competing interests.

Availability of data and materials

Linear distance data and relevant covariates for the CC/F1 specimens and all DO specimens are included within this article’s additional files (Additional files 3, 4, 5). Genotype probabilities of our primary sample are available publicly at http://churchill-lab.jax.org/website/data-percival-etal-2017, while genotype probabilities of our validation sample are available from KMP on reasonable request. Bulk genotypes for the primary sample of DO mice can be found at http://churchill.jax.org/research/cc/do_data/megamuga/raw.

Ethics approval and consent to participate

All mouse procedures were approved and guided by the appropriate institutional animal ethics committees, including at Jackson Labs (#99066), University of North Carolina at Chapel Hill, Scripps Institute (IACUC #08-0150-3), and University of Calgary (IACUC #AC13-0267, AC13-0268).

Funding

This work was supported in part by the NIH via NIDCR (U01DE020054, R01DE019638, R01DE021708 to BH and RM), NIGMS (R01GM070683 to Gary Churchill), NIEHS (R21ES024485 to KMP), and NIDDK (DK076050, DK087346 to DP). This work was also supported in part by NSERC (238992-12) and as a postgraduate fellowship from Alberta Innovates Health Solutions to CJP. Funding for JLM, SAM, & LRD provided by March of Dimes 414-01. None of these funding bodies played a role in study design, collection, analysis, interpretation of data, or in writing the manuscript.

Publisher’s Note

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

Author information

Authors and Affiliations

  1. Department of Anthropology, Stony Brook University, Stony Brook, NY, USA

    Christopher J. Percival

  2. Alberta Children’s Hospital Institute for Child and Maternal Health, University of Calgary, Calgary, AB, Canada

    Rebecca Green & Benedikt Hallgrímsson

  3. The McCaig Bone and Joint Institute, University of Calgary, Calgary, AB, Canada

    Rebecca Green & Benedikt Hallgrímsson

  4. Department of Cell Biology and Anatomy, University of Calgary, Calgary, AB, Canada

    Rebecca Green & Benedikt Hallgrímsson

  5. Program in Ecology Evolution and Conservation Biology, University of Illinois, Urbana, IL, USA

    Charles C. Roseman

  6. The Jackson Laboratory, Bar Harbor, ME, USA

    Daniel M. Gatti, Judith L. Morgan, Stephen A. Murray & Leah Rae Donahue

  7. Department of Molecular Medicine, The Scripps Research Institute, La Jolla, CA, USA

    Jessica M. Mayeux & K. Michael Pollard

  8. Department of Genetics, University of North Carolina Medical School, Chapel Hill, NC, USA

    Kunjie Hua & Daniel Pomp

  9. The Orthopaedic Trauma Institute, Department of Orthopaedic Surgery, UCSF School of Medicine, San Francisco, CA, USA

    Ralph Marcucio

Authors
  1. Christopher J. Percival
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Rebecca Green
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Charles C. Roseman
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Daniel M. Gatti
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Judith L. Morgan
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Stephen A. Murray
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Leah Rae Donahue
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jessica M. Mayeux
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. K. Michael Pollard
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Kunjie Hua
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Daniel Pomp
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Ralph Marcucio
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Benedikt Hallgrímsson
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding authors

Correspondence to Christopher J. Percival or Benedikt Hallgrímsson.

Additional files

Additional file 1.

A.csv file with sequence details and identification of primers used during RT-PCR analysis.

Additional file 2.

Validation Sample Haplotype Effects CC founder strain specific phenotype coefficients (above) and LOD scores from genome-wide scan (below) for significant association between haplotype and A) zygomatic process of the maxilla length and B) zygomatic bone length across Chromosome 17. Phenotypic coefficients are the effect of having a certain founder strain genotype at a specific genomic location on a measurement.

Additional file 3.

A.csv phenotype file with linear distance measures and covariates for the CC Parental/F1 sample discussed in this analysis.

Additional file 4.

A.csv phenotype file with linear distance measures and covariates for the primary DO sample discussed in this analysis.

Additional file 5.

A.csv phenotype file with linear distance measures and covariates for the validation DO sample discussed in this analysis.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Percival, C.J., Green, R., Roseman, C.C. et al. Developmental constraint through negative pleiotropy in the zygomatic arch. EvoDevo 9, 3 (2018). https://doi.org/10.1186/s13227-018-0092-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13227-018-0092-3

Keywords

EvoDevo

ISSN: 2041-9139