1. Genetics and Genomics

Relating multivariate shapes to genescapes using phenotype-biological process associations for craniofacial shape

  1. Jose D Aponte
  2. David C Katz
  3. Daniela M Roth
  4. Marta Vidal-García
  5. Wei Liu
  6. Fernando Andrade
  7. Charles C Roseman
  8. Steven A Murray
  9. James Cheverud
  10. Daniel Graf
  11. Ralph S Marcucio  Is a corresponding author
  12. Benedikt Hallgrímsson  Is a corresponding author
  1. Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Canada
  2. School of Dentistry, Faculty of Medicine and Dentistry, University of Alberta, Canada
  3. Department of Biology, Loyola University Chicago, United States
  4. The Jackson Laboratory, United States
  5. Department of Medical Genetics, Faculty of Medicine and Dentistry, University of Alberta, Canada
  6. Department of Orthopaedic Surgery, School of Medicine, University of California, San Francisco, United States
  7. Department of Animal Biology, University of Illinois Urbana Champaign, United States
https://doi.org/10.7554/eLife.68623
Share this article
Cite this article
  1. Jose D Aponte
  2. David C Katz
  3. Daniela M Roth
  4. Marta Vidal-García
  5. Wei Liu
  6. Fernando Andrade
  7. Charles C Roseman
  8. Steven A Murray
  9. James Cheverud
  10. Daniel Graf
  11. Ralph S Marcucio
  12. Benedikt Hallgrímsson
(2021)
Relating multivariate shapes to genescapes using phenotype-biological process associations for craniofacial shape
eLife 10:e68623.
https://doi.org/10.7554/eLife.68623
  2. Download BibTeX
  3. Download .RIS

Abstract

Realistic mappings of genes to morphology are inherently multivariate on both sides of the equation. The importance of coordinated gene effects on morphological phenotypes is clear from the intertwining of gene actions in signaling pathways, gene regulatory networks, and developmental processes underlying the development of shape and size. Yet, current approaches tend to focus on identifying and localizing the effects of individual genes and rarely leverage the information content of high-dimensional phenotypes. Here, we explicitly model the joint effects of biologically coherent collections of genes on a multivariate trait – craniofacial shape – in a sample of n = 1145 mice from the Diversity Outbred (DO) experimental line. We use biological process Gene Ontology (GO) annotations to select skeletal and facial development gene sets and solve for the axis of shape variation that maximally covaries with gene set marker variation. We use our process-centered, multivariate genotype-phenotype (process MGP) approach to determine the overall contributions to craniofacial variation of genes involved in relevant processes and how variation in different processes corresponds to multivariate axes of shape variation. Further, we compare the directions of effect in phenotype space of mutations to the primary axis of shape variation associated with broader pathways within which they are thought to function. Finally, we leverage the relationship between mutational and pathway-level effects to predict phenotypic effects beyond craniofacial shape in specific mutants. We also introduce an online application that provides users the means to customize their own process-centered craniofacial shape analyses in the DO. The process-centered approach is generally applicable to any continuously varying phenotype and thus has wide-reaching implications for complex trait genetics.

Editor's evaluation

This paper offers a new take on multivariate genotype-phenotype mapping that identifies the joint phenotypic effect of genes involved in known biological processes that impact craniofacial variation. More specifically, the work expands on the traditional idea of candidate gene investigations into candidate biological process investigations, grouping multiple genes into a single analysis. In doing so, the authors show the joint effects of three strong candidate processes, chondrocyte differentiation, determination of left/right symmetry, and palate development on multidimensional craniofacial shape in the heterogenous Diversity Outbred mouse population.

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

Introduction

Variation in human craniofacial shape is moderately to highly heritable (~30–70%; Cole et al., 2017; Tsagkrasoulis et al., 2017), and resemblances among close relatives as well as twins underscore the strong relationship between shared genetics and shared phenotype (Johannsdottir et al., 2005; Nakata, 2014). Despite many studies in humans and in mice (Claes et al., 2018; Cole et al., 2016; Shaffer et al., 2016), however, we know very little about the genetic basis for variation in craniofacial shape. This is likely due to genetic complexity (Katz et al., 2019; Richtsmeier and Flaherty, 2013; Visscher, 2008; Wood et al., 2014; Wray et al., 2013). Like many aspects of morphological variation, craniofacial shape is extraordinarily polygenic. Genes with major mechanistic roles in facial development such as Fgf8 often contribute little to observed phenotypic variation (Green et al., 2017) while genetic influences without obvious connections to craniofacial development emerge as significant contributors (Kenney-Hunt et al., 2008; Klingenberg and Leamy, 2001; Maga et al., 2015; Pallares et al., 2015; Pallares et al., 2014). The effects of genetic variants on phenotype often depend on genetic background (Mackay and Moore, 2014; Percival et al., 2017), and many mutations have variably penetrant effects even when background is controlled (Kawauchi et al., 2009; Rendel, 1967). These issues likely arise because genetic influences act through multiple layers of interacting developmental processes to influence phenotypic traits, resulting in complex patterns of epistasis and variance heterogeneity (Hallgrimsson et al., 2019; Hallgrimsson et al., 2014; Kawauchi et al., 2009; Wagner and Zhang, 2011; Gasch et al., 2016). Solutions that go beyond studies of single gene effects are required to overcome these significant challenges in complex trait genetics. Here, we implement an enhanced form of the more general candidate gene approach to evaluate the conjoint effects of multiple genes on a complex trait – craniofacial shape.

There are two basic approaches to mapping genetic effects on to phenotypic variation. A candidate gene approach measures genotypic values with known physiological and biochemical relationships to the phenotypes of interest (Cheverud and Routman, 1996). In contrast, a random marker or genome-wide approach seeks to associate any potential genetic variant with variation in the trait of interest. There are advantages and disadvantages to these two approaches. The candidate gene approach is blind to the unknown – phenotypic variation is often associated with genes not expected to be important. On the other hand, a candidate gene approach allows direct measurement of genotypic values and produces results that are interpretable in terms of trait physiology or development. A genome-wide or random marker approach can produce unexpected insight by revealing novel gene-phenotype associations. However, this comes at a great cost in power (Visscher et al., 2017). For highly polygenic traits, this approach often produces a ‘tip of the iceberg’ effect in which studies reveal a small and often incoherent subset of the genes that actually determine variation in the trait of interest (Broman, 2009).

Several strategies have been developed that partially overcome these tradeoffs. One solution is the use of polygenic risk scores. Polygenic risk scores assess the overall genetic influence on a trait without regard to the genome-wide significance of individual SNP effects (Dudbridge, 2013; Wray et al., 2007). Approaches such as meta-analyses of genome-wide association studies (GWAS) or studies based on extreme phenotypes Morozova et al., 2015 have expanded gene lists for a variety of complex traits. However, lengthy lists of genes or overall genomic risk for specific phenotypes do not necessarily constitute tractable genetic explanations for phenotypic variation. When thousands of genes are required to explain heritable variation in stature, for instance, it is not clear what such lists tell you beyond the obvious fact that stature is heritable and polygenic (Yang et al., 2010; Wood et al., 2014). This tension between hypothesis-driven and hypothesis-free approaches and their attendant tradeoffs between statistical power and interpretability is, arguably, a major issue within complex trait genetics. To resolve this conceptual conflict, approaches are needed that integrate quantitative genetics with biological insights regarding the cellular and developmental processes through which genes influence phenotypic variation.

Existing approaches to complex trait genetics also tend to treat phenotypic traits as singular and one-dimensional. Even for morphological variation, most studies reduce shape variation to linear distances, principal components (PCs), regression scores or measures of size that are then mapped as individual traits (Xiong et al., 2019; Shaffer et al., 2016; Cole et al., 2016). This approach disregards the information content of multivariate phenotypic variation. While univariate traits only vary along one dimension, high-dimensional traits such as craniofacial shape can vary in direction as well as magnitude within a multidimensional shape space. To identify the distinctive axes of gene effects on a multivariate trait, one must model such multiple multivariate relationships directly.

Building on Mitteroecker et al.’s (2016) multivariate genotype-phenotype (MGP) method, we extend the candidate gene framework to evaluate the combined contributions of genes to variation in high-dimensional phenotypic traits such as craniofacial shape. Grouping genes by ontological information such as membership in developmental pathways or other relevant biological hypotheses, our process-centered multivariate approach, herein referred to as process MGP, brings traditional GWAS together with a simplified model of the hierarchical genotype-phenotype (GP) map. Gene Ontology (GO) terms are broadly grouped into three categories:cellular components, molecular functions, and biological processes. Our work focuses on biological process gene annotations because they group known relationships between several genes that contribute to a developmental function. The process MGP approach aims to leverage this knowledge by modeling the joint effects of these genes on craniofacial shape variation.

Understanding the genetic determinants of craniofacial variation, as with most complex traits, represents a many-to-many GP map problem (Lewontin, 1974; Houle et al., 2010). Both phenotypic and genotypic measurements have complex within-set covariance structures. On the genetic side, the covariance structure is represented by pathway/biochemical interactions, as well as chromosomal structure like linkage, chromatin, and 3D chromosomal organization. For shape-related phenotypes, the covariance matrix is structured by the chosen set of landmarks and their resulting coordinates. The functional relationship from genotype to phenotype is then described by a between-set covariance (Klingenberg and Leamy, 2001; Mitteroecker et al., 2016). To dissect these relationships, we use a regularized partial least squares (PLS) (Lorenzo et al., 2019) approach to estimate a low-dimensional mapping from the alleles in our sample to variation in adult mouse craniofacial shape. While PLS is well suited for analysis of covariation between two sets of measurements, regularization is essential for mitigating overfitting when there are many alleles simultaneously modeled. We focus on how allelic variation in processes relevant to craniofacial development maps to craniofacial shape variation. We ask the following four questions:

  1. How much shape variation is communally accounted for by genes contributing to a process, for example, chondrocyte differentiation?

  2. How similar are the effects of different processes on shape? For instance, do cell proliferation genes affect face shape in a similar way to genes in the bone morphogenetic protein pathway?

  3. How similar are mutant model effects and process effects? For example, do chondrocyte mutant effects align with the effects of natural variants in chondrocyte differentiation genes?

  4. Can one use the similarity of a mutational effect to MGP process effects predict unobserved phenotypes associated with that mutation?

Together, these questions demonstrate the ability of the process MGP approach to add meaningful understanding of the complex relationships between genotype and phenotype by quantifying higher-level regularities between complex phenotypic and genomic data. We also demonstrate its potential as a resource for the study of mutational effects on complex traits such as craniofacial shape.

Results

Process MGP mapping

For process MGP analyses, we used the mouse genome informatics database (Bult et al., 2018Kanehisa et al., 2017) to identify genes annotated to a given process. Each annotation term has an associated GO ID. For example, ‘chondrocyte differentiation’ has GO ID GO:000206 (Figure 1, box 1). We cross-reference the GO ID with the Ensemble genome database (GRCm38.p6) to find the name, chromosome, and base pair start/end position for each gene (Figure 1, box 2) annotated to the process. For genes with multiple splice variants, we select the longest transcript. For each gene, we compare marker base pair positions and select the closest upstream and downstream markers to the center of each gene. The eight-state genotype probability is then calculated as the average founder allele probabilities between the two selected markers. (Figure 1, box 3). After marker selection, we fit a regularized PLS model using the founder allele probabilities (eight variables/marker) and full landmark data set (Figure 1, box 4). Regularization penalizes the coefficients such that increasing regularization strength causes more coefficients to have a value of zero. We chose a regularization parameter using 10-fold cross-validation. For each of the example process MGP analyses shown, we chose the regularization strength that best represented the tradeoff between minimizing model error and maximizing interpretability of marker effects and the similarity of phenotypic effects with mouse mutant models. The full cross-validation results are shown in Figure 2, Figure 4, and Figure 5—figure supplement 1.

Figure 1
Download asset Open asset
Process multivariate genotype-phenotype (MGP) schematic.

Once a process is selected, we cross-reference the known gene locations using Ensembl with the locations of the genotyped markers in the Diversity Outbred (DO) sample. The founder probabilities of the nearest upstream and downstream markers are averaged for each gene. The compiled founder probabilities and landmark coordinates are then used in a regularized partial least squares (PLS) model to estimate the axis of greatest covariance between the marker data and craniofacial variation.

Figure 2 with 2 supplements see all
Download asset Open asset
Process multivariate genotype-phenotype (MGP) for chondrocyte differentiation with a regularization parameter of 0.075.

(A) PLS1 genetic loadings are shown for each gene in the model sorted from largest to smallest effects. Individual founder allele effect sizes are colored within each bar. The gene in red text corresponds to the mutant used for comparison of phenotypic effects. (B) The estimated chondrocyte differentiation MGP phenotype is shown with a heatmap. Warm colors represent areas of relative expansion, light green represents areas of little shape effect, and cool colors represent areas with relative contraction. (C) Chondrocyte differentiation MGP effects shown in black vectors multiplied 4× are compared to a Bmpr1b (Alk6) homozygous mutant and are shown with red vectors. The vector correlation between chondrocyte differentiation MGP and Bmpr1b is shown below the phenotypic effects.

We demonstrate process MGP mapping with three examples. The first estimates the primary axis of skull shape covariation with genes annotated to ‘chondrocyte differentiation’ (Figure 2). Differentiation of chondrocytes is one of several key developmental processes involved in endochondral ossification. Endochondral bones form the majority of the cranial base through a cartilage model of bone formation (Percival and Richtsmeier, 2013). There are 38 genes annotated to chondrocyte differentiation in the Ensembl database (Yates et al., 2020). In the figure, genetic effects are shown as zero-centered bars that span the range of estimated allele effects across the eight DO founders; individual founder allele effects – eight per marker – are color-coded within those bars (Figure 2A). We chose a regularization parameter of 0.075 for this analysis (Figure 2—figure supplement 1). Among chondrocyte differentiation genes, Nov/Ccn3, Bmpr1b (Alk6), and Nfib are most implicated in the major axis of pathway covariation with craniofacial shape. The phenotypic effects at each landmark primarily relate to anteroposterior positioning of the zygomatic arches and dorsoventral jugal position (Figure 2B). The chondrocyte differentiation GP map explains 2.15% of the total variance in craniofacial shape. Compared to 10,000 random permutations of the model, chondrocyte differentiation explains substantially more craniofacial variation (Figure 2—figure supplement 2).

Figure 2C compares the direction of the chondrocyte differentiation MGP axis – magnified 4× – to the axis of shape variation of a relevant mutant phenotype. We chose homozygous Bmpr1b mutants for this comparison for three reasons. The first is because Bmpr1b in synergy with other bone morphogenic protein pathway receptors regulates chondrocyte proliferation and differentiation in embryonic cartilage condensations (Yoon et al., 2005, ). The second reason we chose Bmpr1b mutant comparisons is because the marker selected for Bmp1rb in the genomic analysis contains one of the strongest allelic effects associated with the morphological effect. Bmpr1b shows stronger loading/association than Bmpr1a (Figure 2A). While Bmpr1a has well-established roles in craniofacial development (Liu et al., 2005; Liu et al., 2018), the role for Bmpr1b on its own is less clear. Bmpr1b mutants show shorter long bones at birth and overexpression of a dominant negative Bmpr1b using a type 1 collagen promoter showed delayed ossification of the frontal, parietal, and occipital bones (Yoon et al., 2005; Zhao et al., 2002). The overall phenotypic directions of Bmpr1b mutant variation and chondrocyte differentiation variation are moderately correlated at r = 0.371 (t = –5.06, df = 160, p<0.0001), but the direction at landmarks with large effects in mutant and MGP is clearly coincident. Over the landmarks we measured, the chondrocyte differentiation effect is less global than the Bmpr1b effect, likely due to the difference in severity of the mutant phenotype.

The similarity of the chondrocyte differentiation effect with the Bmpr1b mutant and the high loading Bmpr1b allele in the DO genome suggest that Bmpr1b mutants may produce chondrocyte differentiation defects in the developing neurocranium. In response to the process MGP results, quantified cell size and distribution in the intersphenoid synchondroses (ISS) of several mutant and control Bmpr1b mice (Figure 3—figure supplement 1). We found that homozygotes show overall larger cell sizes as well as a differing distribution of cell sizes throughout the width of the ISS (Figure 3A–C; χ2 = 21.23, df = 3, p<0.0001). The presence of larger cell sizes in the homozygote Bmpr1b mutants suggests that the synchondroses possess more hypertrophic chondrocytes. Additionally, Bmpr1b homozygous mutant mice show premature fusion of the coronal suture (Figure 3D). Both features have not been reported in the literature.

Figure 3 with 1 supplement see all
Download asset Open asset
Chondrocyte defects in Bmpr1b mutants.

(A, B) Quantification of cell size in the sections of the intersphenoid synchondrosis shows an increase in relative cell size as well as a change in the distribution of cell sizes throughout the width of the synchondrosis. (C) Sections of intersphenoid synchondroses highlighting the midline and extremes of the synchondroses. (D) Premature fusion of the coronal suture is visible in Bmpr1b homozygous mutants.

The second example quantifies cranial shape covariation with the 81 genes annotated to ‘determination of left/right symmetry.’ We used a regularization parameter of 0.04 (Figure 4—figure supplement 1). There are several high loading alleles that contribute to the determination of left/right symmetry MGP phenotype. In particular, an Fgf10 allele inherited from the Castaneus founder background was among the most important (Figure 4A). FGF10 is a key ligand in early development, directing proliferation as well as differentiation for many craniofacial components, including the palate, teeth, and bones (Hilliard et al., 2005; Prochazkova et al., 2018; Watson and Francavilla, 2018). The phenotype associated with left/right symmetry alleles is predominately related to a larger neurocranium volume relative to the outgrowth of the face (Figure 4B). We also visualized the asymmetry in the phenotypic response, which shows subtle asymmetry, particularly in the position of the anterior zygomatic landmark (Figure 4D). Several Fgf ligands including Fgf10 are integral in the asymmetric distribution of organs (Hecksher-Sørensen et al., 2004). Fgf8 has also previously been shown to produce asymmetric craniofacial phenotypes in zebrafish, but craniofacial asymmetry has not previously been observed in Fgf10 mutants (Albertson and Yelick, 2005). Left/right symmetry loci explain 3.4% of the total variance in craniofacial shape, which exceeds the variance explained by 10,000 randomly permuted L/R symmetry MGP analyses (Figure 4—figure supplement 2). We compared the estimated L/R symmetry MGP effect with the direction of an Fgf10 homozygous mutant because of the relative importance of the allelic effect (Figure 4C). The vector correlation between the Fgf10 mutant and the estimated left/right symmetry effect is 0.672 (t = 12.29, df = 160, p<0.0001). The importance of other Fgf ligands in craniofacial symmetry, as well as the high-loading Fgf10 allele for left/right symmetry MGP along with the similar genomic and mutant phenotypes, suggests that Fgf10 mutants could show directional asymmetry in the cranium. To test this hypothesis, we measured a sample of 8 Fgf10 adult mutant crania for object symmetry and detected significant directional asymmetry (Figure 4D; F = 4.91, df = 52, p<0.0001).

Figure 4 with 2 supplements see all
Download asset Open asset
Process multivariate genotype-phenotype (MGP) for determination of left/right symmetry with a regularization parameter of 0.04.

(A) PLS1 genetic loadings are shown for each gene in the model sorted from largest to smallest effects. Individual founder allele effect sizes are colored within each bar. The gene in red text corresponds to the mutant used for comparison of phenotypic effects. (B) The estimated left/right symmetry MGP phenotype is shown with a heatmap. Warm colors represent areas of relative expansion, light green represents areas of little shape effect, and cool colors represent areas with relative contraction. (C) Estimated left/right symmetry MGP phenotype is shown with black vectors multiplied 4×. An Fgf10 homozygous mutant is shown with red vectors for comparison. The vector correlation between left/right symmetry MGP and the Fgf10 mutant is shown below the phenotypic effects. (D) Visualizations of asymmetry in the L/R MGP response and the Fgf10 homozygous mutant. Asymmetry vectors are magnified 4×.

The final example estimates the shape covariation attributed to the 73 genes annotated to ‘palate development.’ Formation and fusion of the palatal shelves are crucial for proper orofacial development and heavily influence overall facial shape (Greene and Pisano, 2010). We used a regularization penalty of 0.05 because it best balances the vector correlation to the mutant comparison and the reduction of prediction error (Figure 5—figure supplement 1). Several genes contribute strongly to the palate development MGP effect including Ephb2, Gli3, and Lrp6 (Figure 5A). The estimated phenotype shows corresponding variation in palate length as well as strong effects in the majority of the cranial base landmarks (Figure 5B). Palate development MGP loci explain 2.4% of the total variance in cranial shape, which is greater than the variance explained by 10,000 randomly permuted palate development MGP models (Figure 5—figure supplement 2). We compared the palate development phenotype to a heterozygous Ankrd11, neural crest-specific knockout mouse. The Ankrd11 locus is associated with KBG syndrome in humans, which presents with generally delayed bone mineralization as well as craniofacial characteristics including palate abnormalities (Low et al., 2016). While the vector correlation between the palate development MGP effect and the Ankrd11 mutant over the complete set of cranial landmarks is moderate at r = 0.339 (Figure 5C; t = 4.31, df = 160, p=0.0001), the vector correlation for palate landmarks is substantially higher at r = 0.536.

Figure 5 with 2 supplements see all
Download asset Open asset
Process multivariate genotype-phenotype (MGP) for palate development.

(A) PLS1 genetic loadings are shown for each gene in the model sorted from largest to smallest effects. Individual founder allele effect sizes are colored within each bar. The gene in red text corresponds to the mutant used for comparison of phenotypic effects. (B) The estimated palate development MGP phenotype is shown with a heatmap. Warm colors represent areas of relative expansion, light green represents areas of little shape effect, and cool colors represent areas with relative contraction. (C) Estimated palate development MGP phenotype is shown with black vectors multiplied 4×. An Ankrd11 mutant mean is shown with red vectors for comparison. The vector correlation between palate development MGP and the Ankrd11 mutant is shown below the phenotypic effects.

In each case above, we have shown how association of gene sets and phenotypic variation can produce highly informative results that can guide subsequent hypothesis testing. For a given biological process, we identified several genes that load strongly on the primary axis of MGP covariation for which mutant samples were available to us, as well. Future investigations could also use this information about genes with high loadings to generate new mutants for analysis of associated developmental processes. For each example, we focus only on the first PLS axis, so distinct joint gene effect combinations may contribute to novel phenotypic directions in lower PLS axes.

Joint versus single-loci effects

While the process MGP approach focuses on the joint effect of markers on craniofacial shape, it is important to measure the extent that joint effects matter for craniofacial shape. Unlike alternative models such as canonical correlation analysis (CCA), the PLS model used for process MGP does not allow for statistical tests of individual marker effects. However, it is possible to measure the similarity of phenotypic effects after successively removing the most heavily loading markers from the full model. Figure 6 shows the change in variance explained (A) as well as the change in the direction of phenotypic response (B) as markers are increasingly removed from the model. For each example analysis, we remove markers in order from most heavily loaded to least heavily loaded. We found that process MGP analyses with few loci of very large effect, like Ccn3/Nov for chondrocyte differentiation MGP, are very sensitive to the removal of the most highly loaded genes. The vector correlation between the full chondrocyte differentiation MGP model and the model with the Ccn3/Nov marker removed is 0.346.

Figure 6 with 1 supplement see all
Download asset Open asset
Gene drop tests.

For each of the example analyses, we show the effect of removing the most heavily loaded markers from the process multivariate genotype-phenotype (MGP) analysis on the (A) variance explained by the model and (B) vector correlation with the full model. The variance explained as well as vector correlation is relatively stable for both L/R symmetry and palate development MGP models, suggesting that the effect is driven by the coordination of many markers. In contrast, chondrocyte differentiation MGP shows large differences, particularly in the direction of the phenotypic effect as the most heavily loaded markers are removed from the analysis.

For process MGP analyses with a more uniform distribution of marker effects, we found that the phenotypic effect is much more reliant on a multitude of marker alleles. For instance, L/R symmetry MGP with the 10 most heavily loaded markers removed still produced a vector correlation of 0.95 with the full model. The majority of process MGP analyses demonstrate a similar importance to several alleles, highlighting the main strength of process-level analyses over individual marker tests (Figure 6—figure supplement 1). In the following sections, we will examine how process MGP phenotypes relate to each other, as well as the phenotypic directions of several mutant mouse models.

Pairwise comparison of craniofacial development processes

We chose 15 processes integral to craniofacial development and compared the pairwise similarity of effect on craniofacial shape using a heatmap based on clustering of the correlation matrix (R Development Core Team, 2017 ). Processes with similar effects on craniofacial shape will be highly correlated, while processes that affect distinct aspects of craniofacial variation will be uncorrelated to each other. The clustering algorithm resulted in two main blocks of strongly correlated effects (Figure 7A). The largest block of highly correlated phenotypic effects includes neural crest cell migration, epithelial to mesenchymal transition, forebrain development, as well as some of the most general developmental processes like cell proliferation, bone development, apoptosis, A/P pattern specification, and FGFR signaling. In addition, there is a general BMP block, with Bmp signaling, dorsoventral pattern formation, endochondral ossification, and positive regulation of skeletal muscle tissue growth. Interestingly, phenotypic variation associated with cranial suture morphogenesis, neural tube patterning, and intramembranous ossification is largely uncorrelated with the other craniofacial developmental processes included here.

Figure 7 with 1 supplement see all
Download asset Open asset
Pairwise multivariate genotype-phenotype (MGP) vector correlations.

(A) Pairwise correlations of phenotypic effects for 15 process MGP analyses. Scale on the right denotes color correspondences to vector correlation, where yellows are high correlations, greens are moderate, and blues are low. (B) Pairwise process MGP vector correlations as a function of the number of shared genes between the processes. Processes that share less than 10 genes can produce very similar and very disparate phenotypic effects. Processes with substantial numbers of shared genes will tend to show highly correlated responses as they increasingly use similar marker sets.

To assess the stability of the clustering result, we estimated the vector correlation between the cluster distances – also known as the cophenetic distance – and the original correlation matrix (Sneath and Sokal, 1973). A high vector correlation suggests reliable clustering, whereas a low correlation suggests a random clustering result. The correlation between the cophenetic distance matrix and the correlation matrix is 0.648 (t = 8.64, df = 103, p=7.6 × 10–14), suggesting a moderate, though significant, structure in the similarity of effects amongst this set of MGP processes.

The similarity in process MGP effects suggests that processes may coordinate in a limited set of potential directions of phenotypic variation. One reason that we could observe this pattern that is not because of common axes of GP variation is that key genes show up repeatedly within processes and largely drive these patterns of phenotypic variation. Figure 7B shows over 45,000 pairwise process MGP vector correlations as a function of the number of shared genes between the two randomly chosen annotations. While the similarity of phenotypic effects generally increases as the number of shared genes increases, GO processes that share genes are not necessarily strongly correlated. For GO processes that share between 0 and 10 genes, the observed correlations in phenotypic response spanned from no correlation to almost entirely concordant. GO processes that share more than 10 genes show generally higher vector correlations, with the lowest vector correlation we observed at 0.35.

Process effects in the mutant morphospace

To assess the extent to which craniofacial shape variation associated with developmental processes aligns with variation from mutants of major effect, we projected seven process effects onto the first two PCs of a dataset containing the DO sample, and samples from 30 mutant genotypes (Figure 8A). Each black label represents the mean shape score of the listed mutant genotype. The shaded ellipse with an orange border displays the 95% data ellipse of PCs 1 and 2 of DO cranial shape variation. The DO mean shape is contrasted by the mutant variation along PC1. The first PC describes vault size relative to the length of the face. The phenotype shown along the X axis of Figure 8A depicts the maximum positive PC1 shape, while the heatmap drawn on the crania represents the local deformations towards the minimum PC1 shape. The positive direction of PC2 describes coordinated variation that includes a relatively wider vault, narrower zygomatic, and shorter premaxilla (Figure 8A, Y axis margin).

Figure 8
Download asset Open asset
Comparisons of multivariate genotype-phenotype (MGP) and mouse mutant directions.

(A) Seven MGP phenotypes projected onto a principal component analysis (PCA) of the Diversity Outbred (DO) and a sample of 30 mutant mouse genotypes. Mutant means are labeled in black. The directions of MGP effects are shown with orange vectors from the DO mean to the associated process MGP. The range of DO variation on principal components (PCs) 1 and 2 is shown with the shaded ellipse with an orange border. (B) A heatmap of vector correlations between 30 mutant effects and 30 process MGP effects. The scale on the right denotes color correspondences to vector correlation, where yellows are high correlations, greens are moderate, and blues are low.

Process effects – highlighted with orange vectors originating at the DO mean shape – are necessarily of smaller magnitude than the total variation in the DO sample. Therefore, to better compare the direction of process effects the vector magnitudes were magnified 4×. Several process effects align in distinct directions of mutant effects, such as bmp signaling pathway and endochondral ossification in the direction of Shh, Nipbl, and Ift88 mutants. Neurotransmitter transport and Wnt signaling pathway is similar in direction to Kcna1Mceph and B9d1 mutant effects. Execution phase of apoptosis and intracellular transport both show similar effects to a cluster of Bmp mutants. Figure 8A focuses on two PCs, which allows for the contextualization of how process MGP analyses and mutants vary in similar directions and allows us to visualize what those phenotypes look like. This combination of context and visualization is only possible in limited axes and cannot account for differences or similarities in the full multivariate shape space.

To show the similarity of process MGP directions with mutants in the full shape space, we present a heatmap of 30 process MGP effects to 30 mouse mutant models in Figure 8B. The heatmap shows the correlation in direction with yellow/green denoting higher correlation and teal/blue denoting lower correlation. The bottom right of the heatmap (highlighted by a white border) shows a block of mutants for which there are strong process correlations. These are among the most extreme phenotypes along PC1 (Figure 8A) and include mutants for Nosip, Bmp2, Grm1, Bmp2; Bmp7 transheterozygote, Bmp7, Ghrhr, Fgf10, and Papps2. The processes most strongly correlated to these mutants are histone methylation, dendrite morphogenesis, chromosome segmentation, vasodilation, and fibroblast growth factor binding.

There are a set of mutant phenotypes that have generally low correlations to the set of processes chosen. These mutants include Fgf3, Shh, Nipbl, Disp, Pten, Hhat, and Alk2; Alk3 transheterozygote. Interestingly, this group of mutants varies more along PC2 than PC1 (Figure 8A). Notably, regulation of intracellular protein transport and regulation of cell death are strongly uncorrelated with the majority of mutant directions.

Real-time process GP mapping

Finally, we provide an online tool to visualize process effects and make comparisons to mutant effects in real time. This application is found at https://genopheno.ucalgary.ca/MGP/ and can be used for analyses similar to those described in this paper. When the user selects GO terms, the program searches for genotype markers adjacent to each gene listed and uses the selected markers to fit a regularized PLS model. The result is an estimate of the many-to-many relationship between the selected markers and cranial shape variation. The visual outputs include barplots depicting the relative allele effect sizes for each gene in the process and a 3D plot of the corresponding axis of shape variation. Users can compare the effects of different processes and also compare process effects to mutant effects from a provided database of 30 mutant genotypes.

To illustrate how to use this application, we have provided the graphical user interface used to select the parameters (Figure 9). As an example, in the ‘Process text’ entry field, supply a starting term; we chose ‘brain.’ The GO database is then filtered, returning a user-selectable subset of biological process ontology annotation terms in the ‘Process filter’ field. We chose ‘forebrain morphogenesis,’ which has 11 associated genes. We chose to magnify the process phenotype vectors 4× and compare the effect to a heterozygous Ift88 mutant. Ift88 is a core component of the primary cilia, which are responsible for promoting developmental signals involved in many facets of facial development (Tian et al., 2017). Further, the plots that are generated are interactive. For example, marker loadings can be highlighted and subset by genes of interest (Sievert, 2019). There is further information about using this online tool in the ‘About this app’ tab.

Figure 9 with 1 supplement see all
Download asset Open asset
Example screenshot of web version of process analysis.

Analyses include a barplot of the relative effect sizes of each selected marker and the associated phenotype shown with black vectors at each landmark. If a mutant comparison is selected, the vector correlation is provided and the mutant phenotype is shown with red vectors. Selecting ‘send me the results’ generates an HTML report with an interactive 3D model.

Discussion

A key goal in genomics is to create tractable genetic explanations for phenotypic variation. In this study, we used the process MGP approach to model the joint effects of genomic markers on multivariate craniofacial shape. This approach allows us to address the joint contributions of multiple genes that share ontological characteristic such as pathway membership on craniofacial shape as a multivariate trait. Specifically, we chose markers adjacent to genes annotated under a developmental process of interest. We showed three process MGP analyses in depth, each with distinct phenotypic effects. Each of these comparisons highlighted the integrated structure of phenotypic variation in mouse craniofacial shape. We found that while there are processes with distinct and localized effects, genetic effects generally converge on a limited set of directions in phenotype space. Further, these process effects often correspond with the directions of major mutations known to affect these same processes.

Many recent studies have addressed the genetics of craniofacial shape in humans and mice (reviewed in Roosenboom et al., 2016; Weinberg et al., 2018; White et al., 2021). While these studies are yielding a growing list of genes, suggesting that facial shape is highly polygenic, they have left the vast majority of heritable variation unexplained. Existing studies have either used univariate measures of facial shape such as linear measurements or univariate summaries of multivariate shape (e.g., Procrustes distances or PC scores). In addition, most genomic studies of craniofacial shape quantify the effects of each genomic marker independently, with notable exceptions focusing on epistatic effects (e.g., Varón-González et al., 2019). Our approach shares common features with some predecessor GP mapping strategies in which candidate genes/SNPs are selected a priori because of common involvement in a pathway (or other mechanistic cluster) (Claes et al., 2014; Liu et al., 2012; Wang et al., 2010; Wang et al., 2007). Wang and colleagues selected SNPs based on proximity to genes of interest and effect size to jointly model the pathway-level effects on Parkinson disease data. Their approach is similar to gene set enrichment analysis, weighing overrepresentation of statistical effects related to case-control group membership. In contrast, our approach focuses on estimating a multivariate set of continuous craniofacial responses. Importantly, our approach jointly identifies GP axes that maximally covary. This differs significantly from approaches that determine phenotypes for analysis a priori or based on a predetermined method of data reduction such as the principal component analysis (PCA). Our implementation also differs from similar methods like CCA that was used to associate single-locus effects with a multivariate phenotype (Claes et al., 2018). In comparison with the process MGP approach, CCA has the advantage of allowing for a parametric hypothesis test, whereas PLS analyses are limited to permutation-based hypothesis testing. A distinguishing feature of the process MGP approach is the ability to penalize the model with regularization. This is ideal for models with many simultaneous genetic effects in order to mitigate the effects of overfitting. Regularization is not unique to PLS as applications of ridge penalties to CCA have been used for genomic analyses (Waaijenborg and Zwinderman, 2009; Le Floch et al., 2012).

A key finding of our application of the MGP method to craniofacial shape is that multivariate phenotypic variation aligns nonrandomly to genetic markers associated with pathways or developmental processes. Process MGP effects that are generally not driven by single loci of large effect are possible, like with the chondrocyte MGP analysis (Figures 2—5, Figure 6—figure supplement 1). These covarying effects represent the joint genetic effects of multiple contributors to phenotypic variance. While these patterns of MGP covariation may include genetic variants that do not actually affect the phenotype, many other high-loading alleles will be contributors that we lack statistical power to detect under a typical univariate approach (Pitchers et al., 2019; Varón-González et al., 2019). Here, the overall pattern of GP covariance is the level of genetic explanation for phenotypic variation. When such patterns involve genes that are ontologically linked in meaningful ways, they can provide novel insights into the coordination of genetic effects on phenotypic variation and bolster existing hypotheses from developmental studies.

Another valuable asset that arises from the process MGP approach is the ability to generate testable hypotheses or predictions from MGP observations. The chondrocyte differentiation MGP analysis suggested differentiation defects in the Bmpr1b mutant that could contribute to craniofacial variation. We followed up the MGP analysis with histological analysis of Bmpr1b mutants and showed premature suture fusion as well as atypical distribution of hypertrophic chondrocytes in the ISS. Similarly, the analysis of left/right symmetry genes suggested that Fgf10 alleles can contribute to directional asymmetry. A follow-up morphometric analysis of symmetry showed that Fgf10 mutants do display significant craniofacial asymmetry (Figure 4D). Process MGP can also be used to test existing hypotheses about GP relationships. The relative importance of the Ankrd11 locus in the palate development analysis and the similarity between the genomic and mutant phenotype further validate the role of Ankrd11 in palate development. These examples illustrate the additional insights that a process MGP analysis of a mutational effect can provide. Given that such comparisons can be run quickly with our web application, this creates a tool with the potential for hypothesis generation and initial screening for hypotheses about process-level effects on craniofacial variation in mice.

The explicit modeling of multivariate relationships between phenotypes and genotypes also allows a focus on pleiotropy. Developmental studies in mice demonstrate widespread craniofacial morphological effects from localized developmental perturbations (Martínez-Abadías et al., 2012; Stelzer et al., 2007; Young et al., 2010) Perturbations to specific processes in development generally produce effects on multiple aspects of phenotype due to knock-on effects at later stages or to interactions at the level of tissues or anatomical structures (Hallgrímsson et al., 2007; Hallgrímsson et al., 2009). A change in cartilage growth in basicranial synchondroses produces a global change in craniofacial form, for example (Parsons et al., 2015). Remarkably, enhancers with highly specific temporospatial effects on gene expression also produce global rather than localized changes in craniofacial shape (Attanasio et al., 2013). Given that pleiotropy is likely ubiquitous (Churchill et al., 2012; Wagner et al., 2007), explicitly multivariate approaches to understanding GP maps are clearly needed.

This convergence of genetic effects on axes of covariation is reflected in our finding that mutations to major developmental genes produce effects that tend to align with the directions of effect associated with the corresponding broader pathways or ontological groups (Hallgrímsson et al., 2009). These results suggest that perturbations that are developmentally similar tend to move the phenotype in the same direction in multivariate space (Figure 8). Even so, both mutational and higher-level pathway/process effects tend to converge on a few directions of variation, suggesting that multiple pathways and processes lead to common developmental outcomes (Houle and Fierst, 2013; Uller et al., 2018; Hallgrímsson and Lieberman, 2008; Gonzalez et al., 2013). This conclusion is further supported by our finding that the genetic axes of covariance for individual processes/pathways can align with multiple directions of mutational effect. For example, the process MGP phenotypes clustered in the bottom right of Figure 8B are all highly correlated with a set of BMP and growth hormone-related mutants.

In some cases, mutants and MGP map directions do not correspond. There are several ways this can occur. The first is that the DO population may simply lack alleles as deleterious as found in mutant lines. A small effect allele in the DO may not align with the direction of a mutant almost completely lacking expression of the target gene. Further, there are many examples where a mutation may have different and sometimes even opposite effects depending on genetic background (Mackay, 2013; Percival et al., 2017). Mutations of major effect may also differ in direction from variants in related genes that have smaller phenotypic effects due to underlying nonlinearities in development (Green et al., 2017). Investigating how variants in genes that are functionally related vary in phenotypic effect is an important avenue of inquiry that is revealed by analyses such as those we have performed here. Additionally, relationships between process and mutant effects may stimulate hypotheses about previously unknown or unvalidated interactions between loci or pathways.

A second potential reason that MGP effects may not correspond to major mutation effects is the use of only one PLS axis for each process analysis. With only one axis, we only show the phenotypic direction with greatest covariance with genetic marker variation. If there are multiple large marker effects that do not covary, the weaker marker effect will be masked in the analysis. For instance, there may be a PLS axis for ‘chondrocyte differentiation’ that corresponds more strongly with the Bmp2 mutant phenotype. This phenomenon may be particularly prominent for pathways with substantially different mutant effects, like FGF (Figure 8A). While we did not delve into the directions outside of the first PLS axis, we have facilitated the selection of lower axes in the web application for users to explore and compare with mutants of interest.

Finally, our analysis shares the limitation of all approaches based on gene annotation data. Incomplete annotation may lead to faulty or incomplete groupings of genes when defining pathway/process hypotheses. Gene annotation is a huge undertaking, and there is substantial variation in the completeness of different process annotations. Many process annotations are manually assigned using inference from the literature, while most are a combination of automated efforts based on transcript similarity and human curation (Mudge and Harrow, 2015; Finger et al., 2017). Related to this, we assign gene annotation data to genetic markers based on the closest protein-coding region. While this is a reasonable proxy, there will be regulatory sites that affect genes other than the one immediately adjacent and this is a potential source of uncertainty in our analysis (Forrest et al., 2014; Yue et al., 2014).

In addition, this approach does not currently model the temporal and spatial aspects of gene function throughout development. As a result, alleles of high importance in an MGP analysis do not necessarily produce craniofacial variation through the selected process. A strong allelic effect like Ccn3 can load heavily in several processes, like ‘chondrocyte differentiation,’ ‘fibroblast migration,’ and ‘negative regulation of inflammatory response.’ We do not know the mechanism through which any individual allele contributes to variation from an MGP analysis alone. Genomic data with more fine grain measurements of variation in expression and utility of individual loci may be better suited to teasing out the potential mechanisms that alleles produce variation.

The MGP method represents a deliberate decision to trade higher-level insight from GP association data at the expense of statistical certainty about the significance of individual gene effects. The current implementation of the method also does not allow for quantification of individual epistatic effects. Epistasis occurs when the genotypic trait value for a locus is altered by the genotype of a different locus. Such effects generate nonlinear GP maps, but when considered genome-wide, contribute mainly to additive variance (Cheverud and Routman, 1996; Hill, 2017). The MGP method is additive in that it models only the linear effects of genes. However, since it captures the covariances among genotypic effects, much of this ‘additive’ variation is likely epistatic in origin.

Complex traits present a massive challenge in genomics because so many are turning out to be enormously polygenic. To generate tractable explanations of the genetic basis for such traits, methods are needed that extract higher-level representation of GP relationships than those that emerge from single-locus-focused approaches. Here, we present a process-driven framework for deriving such higher-level genetic explanations for phenotypic variation. Our approach leverages the biological tendency for developmental processes to produce covariation among aspects of a multivariate phenotypic trait (Kawauchi et al., 2009; Wagner et al., 2007). The underlying assumption in this approach is that there are latent variables within high-dimensional GP data that correspond to developmental architecture. We believe that analyses aimed at defining and characterizing such latent variables represent a level of genetic explanation for phenotypic variation that is complementary to genetic analyses designed to establish the significance of single-locus effects. Pursuing such questions will help bridge the gap between emerging mechanistic accounts of morphogenesis and our growing understanding of the genetics of morphological variation.

Materials and methods

Mice

We use a sample (n = 1145) of DO mice (Jackson Laboratory, Bar Harbor, ME) to map GP relationships for craniofacial shape (Churchill et al., 2012; 2004). The DO is a multiparental outcross population derived from the eight founding lines of the Collaborative Cross (CC). Each animal’s genome is a unique mosaic of the genetic diversity found in the CC – more than 45 million segregating SNPs (Collaborative Cross Consortium, 2012). Random outcrossing over many DO generations maintains this diversity and, with recombination, increases mapping resolution . Discussions of recommended sample sizes in univariate DO studies can be found in Churchill et al., 2012. Both studies recommend a sample size greater than 800 mice for small univariate effect sizes (1–5% variance explained). Further, there are inherent power advantages to our approach because multivariate responses represent maximized differences in phenotype given a set of genotypic measurements. In contrast, univariate approaches such as analyses of individual PCs can only detect effects along those predefined axes that may not have clear biological significance.

Our DO sample was sourced from three separate laboratories and seven DO generations. 386 are from the Jackson Laboratory (JAX), 287 from the University of North Carolina (UNC), and 472 come from the Scripps Research Institute. Figure 10—figure supplement 1 shows the distribution of the sample by lab source and generation of breeding. Imaging of mice at the University of Calgary was performed under IACUC protocol AC13-0268. Ankrd11 and Bmpr1b mutant mice were bred at the University of Alberta by the Graf lab under Animal Use and Care Committee protocol AUP1149 in accordance with guidelines of the Canadian Council of Animal Care.

Genotyping

Request a detailed protocol

Genotyping was performed by Neogen (Lincoln, NE). Ear clippings were used to extract DNA for all samples. Mice from generations 9, 10, and 15 were genotyped using the MegaMUGA genotyping array (77,808 markers); mice from generations 19, 21, 23, and 27 were genotyped using the larger GigaMUGA array (143,259 markers) (Morgan et al., 2015). To pool the genotype data from these two SNP arrays with differing numbers of markers, we imputed markers between the two genotyping arrays using the ‘calc_genoprob’ function in the qtl2 package (Broman et al., 2019). The function uses a hidden Markov model to estimate genotype probabilities and missing genotype data (Gatti et al., 2014). After imputation, the merged genetic dataset consists of 123,309 SNPs that vary among CC founders. Each animal’s genetic record is a 123,309 * 8 matrix of estimated diplotype contributions of each CC founder to each marker.

Scanning and landmarking

Request a detailed protocol

We used micro-computed tomography to acquire 3D scans of the full heads of the mice. Scanning was done at the University of Calgary at 0.035 mm voxel resolution (Scanco vivaCT40). One of us (WL) then acquired 54 3D landmarks (Figure 10) manually on each volume using Analyze 3D. A discussion of the error associated with manual landmarking can be found in Katz et al., 2019. In addition to the DO phenotype data, the mutant mouse data used for comparisons were collected, scanned, and landmarked between the Hallgrimsson and Marcucio labs.

Figure 10 with 1 supplement see all
Download asset Open asset
54 3D landmark configuration.

(A) Sagittal view of representative scan with landmarks shown as red spheres. (B) Dorsal view of landmark configuration. (C) Ventral view of landmark configuration.

Landmark registration and analysis

Request a detailed protocol

We symmetrized landmarks along the midline of the skull using Klingenberg et al.’s method for object symmetry that configures landmark pairs into a common orientation with reflection and subsequently removes variation associated with translation, scale, and rotation, using Generalized Procrustes Analysis (Adams et al., 2013; Klingenberg et al., 2002; Mardia, 2000 et al.; Schlager, 2017). We tested for directional asymmetry using the Procrustes ANOVA approach described in Klingenberg et al., 2002. To focus on shared, within-generation patterns in our multigenerational DO sample without sex effects, we regressed symmetric shape on DO generation and sex and used the residual shapes with the grand mean added as the observations for analysis.

Genetic relatedness

Request a detailed protocol

Adjustment of phenotypes for the influence of genetic relatedness is a common approach in genomic studies to prevent spurious associations. However, it is not necessary in all cases, such as situations with low genetic relatedness and little variation in relatedness. We evaluated whether accounting for genetic relatedness was important for our sample. To do so, we estimated a kinship matrix based on DO genotype correlations (Cheng and Palmer, 2013; Broman et al., 2019). The kinship values in our sample have a mean of 0 and a standard deviation of.047. As a result of these findings, we performed all subsequent analyses on the within-generation symmetric shape data, without an adjustment for relatedness.

Regularized PLS analysis

Request a detailed protocol

MGP methods for explicitly modeling multivariate phenotypes and for overcoming the limitations of simple linear regression are increasingly common in mapping studies. One example of a multivariate genomic approach is found in Claes et al., 2018, where the authors used CCA to quantify individual SNP effects for a multivariate measurement of facial shape. CCA returns a vector of the linear combination of phenotypic effects that maximally correlates to the alleles at a given locus. Mitteroecker et al., 2016 developed a similar multivariate strategy around a singular value decomposition (SVD) of GP covariance matrices (versus decomposition of correlation matrices in CCA). PLS describes a family of approaches that use SVD to decompose cross-covariance matrices (Lee et al., 2011; Mitteroecker et al., 2016; Singh et al., 2016). PLS is increasingly used with large genetic datasets in order to model how genomic effects extend to multiple traits (Bjørnstad et al., 2004; Mehmood et al., 2011; Tyler et al., 2017). However, its implementation for MGP mapping is, thus far, much more limited.

SVD decomposes the covariance matrix into three matrices:

Y=UDV′

where Y is the mean-centered covariance matrix, U denotes the left singular vectors, a set of vectors of unit length describing the relative weighting of each variable on each axis, and D denotes the variance along each axis. V denotes the set of right singular vectors. For a full (square, symmetric) covariance matrix, U and V are identical, and the decomposition is equivalent to PCA. For a non-symmetric matrix of covariances, that is, one describing covariance between two distinct blocks of traits, each successive column of U and V provides a pair of singular vectors describing the best least-squares approximation of covariance between the two blocks, in order of greatest covariance explained to least.

PLS is most often used to find low-rank linear combinations that maximize covariance between two sets of features. Here, we use the data-driven regularized PLS model implemented in the mddsPLS package to find paired axes that maximize covariance between allelic and shape variation (Lorenzo et al., 2019). The model uses a lasso penalty to minimize the coefficients (loadings) towards zero to prevent overfitting (James et al., 2013). Overfitting can occur when many genotypic markers are included in the model, particularly when markers are colinear. The genotype block is composed of the full set of DO founder probabilities for each selected marker. Thus, an analysis of 20 markers would estimate 160 genotype coefficients. The phenotype block consists of the full set of 54 3D landmarks (162 phenotype coefficients). In all biological process analyses undertaken herein, we used a regularization parameter of 0.06 and report only the first paired axes of the PLS model, that is, the genotype and phenotype axes that explain the most covariance.

We generate graphical displays of process results using the R packages ggplot2 (Wickham, 2016) and Morpho (Schlager, 2017). An example script to reproduce the analyses is provided at https://github.com/J0vid/MGP_shiny/tree/main/analyses (copy archived at swh:1:rev:61fb597c48ced306dd588e289e69c3e3d8f9ce15, Aponte, 2021).

Statistical results and comparisons

Request a detailed protocol

We estimate the magnitude and direction of MGP process effects using R2 and vector correlations, respectively. R2 is calculated as the ratio of trace of the predicted model covariance to the trace of the phenotypic covariance matrix. We contextualize the MGP process R2 by comparing it to the R2 value of 10,000 randomly drawn marker sets of the same size. For instance, a process annotated with 40 genes would be compared to 10000 40-gene MGP analyses with random markers selected in each iteration. Random marker selection for permutation is constrained to follow similar patterns of linkage disequilibrium to the observed marker set of interest. The null expectation in this scenario is that gene annotation does not provide better information about coordinated marker effects than a randomly selected set of markers.

Vector correlations between process MGP effects are calculated by taking the Pearson product-moment correlation of the two sets of process PLS1 phenotypic loadings. Vector correlations between process effects and mutant effects are calculated by taking the correlation between the process PLS1 phenotypic loadings and mutant MANOVA coefficients. MANOVA was used to compare the mutant group phenotype against the DO sample specified as the reference group. The coefficients of MANOVA describe the relative weights of each landmark coordinate difference between the DO mean shape and the mutant mean shape.

Chondrocyte morphometrics

Request a detailed protocol

Chondrocyte morphometrics were performed using a novel technique developed by the Marcucio Laboratory. Images of the ISS were stained with H&E, SafO, or picrosirius red and were captured and imported into ImageJ (2–6 sections from at least four mice/genotype/synchondrosis; Rueden et al., 2017). Landmarks were placed in a defined order (left, right, top, bottom) of visible chondrocytes in the synchondrosis using ImageJ’s multi-tool. Data points were then exported as XY coordinates and imported into Microsoft Excel for calculation of major and minor axes relative to overall width of synchondrosis. Area of individual cells was determined from height and width values based on the assumption that each cell is roughly ellipsoidal. An example of major and minor axis measurements and ellipsoidal area measurements on a slide is provided in Figure 3—figure supplement 1.

We compared differences in the distribution of cell sizes along normalized synchondroses between Bmpr1b mutants and controls with a mixed effects model approach. We used ellipsoidal area of cell size (in microns) as our dependent variable. For fixed effects, we modeled the normalized synchondrosis position (first and second order), where a value of 0 represents the relative midline of the synchondrosis and values of –1 and 1 represent the most distant cells in that synchondrosis. We also modeled genotype as a fixed effect as well as a genotype by cell position interaction (both first- and second--order interactions). For each individual within each genotype, we measured multiple histological sections. These repeated and nested measurements of cell size in multiple sections for each individual were modeled as random effects. To test for cell size differences between genotypes, we used a likelihood ratio test to compare the full model to a reduced model with the fixed effect of genotype and all genotype interactions removed.

Visualization tools

Request a detailed protocol

We introduce an interactive web application that allows the user to select processes of interest with a graphical user interface; see the resulting craniofacial effect at https://genopheno.ucalgary.ca/MGP/. The web apps were written using the shiny package in R (Chang et al., 2018; R Development Core Team, 2017). The application dynamically filters the MGI GO database based on the initial user input. Queries will only list GO terms with exact matches. For example, ‘chond’ will return GO terms that incorporate either ‘chondrocyte’ and ‘mitochondria’.

Multiple queries can be selected. An analysis of ‘chondrocyte differentiation’ and ‘chondrocyte hypertrophy’ will select the joint gene set of both processes. Processes with different names can be jointly queried with the pipe operator ‘|,’ which is interpreted as an OR (union) operator. For example, to generate the list of GO terms associated with either apoptosis or WNT, we used the ‘apoptosis|WNT’ query and selected the processes ‘Wnt signaling pathway’ and ‘execution phase of apoptosis’ to perform the analysis on the joint set of associated genes (Figure 9—figure supplement 1).

Several other parameters can be specified by the user including the type of plot to be generated for the genetic loadings, the amount of magnification applied to the phenotype effect vectors, the regularization parameter, and the option to overlay a mutant phenotype for comparison. The comparative database currently includes craniofacial shape contrast data (wild-type vs. mutant) for 30 mutant genotypes. If a mutant comparison is selected, the full set of DO specimens are registered with the mutants added (with size removed). We then provide the vector correlation between the process effect and the mutant effect (see Figure 9). The database also includes PC1 of the DO sample for comparison.

The app enables users to save results. A save request will generate and download an HTML report of the analysis that includes several versions of the genetic effect plot and an interactive 3D model of the estimated phenotypic effect. If a mutant comparison is selected, it will also appear in the report.

The application tracks recent searches by the user for their reference. A heatmap of process vector correlations of the PLS phenotype loadings is also available under the ‘recent searches’ tab. The user can select between a heatmap of the processes in their search history or a random assortment of process correlations from past anonymous user searches.

Finally, we provide programmatic access to our model for both process MGP analyses as well as custom gene lists over the web through an application programming interface (API). Queries can be formatted using curl commands as well as request URLs and return results in JavaScript object notation (JSON) format. Documentation for the available functions and their parameters, as well as examples for queries, can be found at https://genopheno.ucalgary.ca/api/__docs__/. The API was written using the plumber package (Schloerke and Allen, 2021) in R, with code available at https://github.com/J0vid/MGP_shiny/tree/main/MGP_API.

Data availability

All diversity outcross microCT scan and QTL data have been deposited with Facebase (https://doi.org/10.25550/1-731C). Scripts are available at https://github.com/j0vid/MGP_shiny in the analyses folder, (copy archived at https://archive.softwareheritage.org/swh:1:rev:61fb597c48ced306dd588e289e69c3e3d8f9ce15) and the associated online tool is available at https://genopheno.ucalgary.ca/MGP/.

The following previously published data sets were used

References

  46. Book
    1. Lewontin R
    (1974)
    The Genetic Basis of Evolutionary Change
    Columbia University Press.
  71. Book
    1. Rendel JM
    (1967)
    Canalization and Gene Control
    London: Logos Press.
    1. Schlager S
    (2017)
    Statistical Shape and Deformation Analysis
    217, Chapte9, Statistical Shape and Deformation Analysis, Amsterdam, Netherlands, Elsevier.
  80. Book
    1. Sneath PHA
    2. Sokal RR
    (1973)
    Numerical Taxonomy: The Principles and Practice of Numerical Classification
    Freeman.
    1. Wood AR
    2. Esko T
    3. Yang J
    4. Vedantam S
    5. Pers TH
    6. Gustafsson S
    7. Chu AY
    8. Estrada K
    9. Luan J
    10. Kutalik Z
    11. Amin N
    12. Buchkovich ML
    13. Croteau-Chonka DC
    14. Day FR
    15. Duan Y
    16. Fall T
    17. Fehrmann R
    18. Ferreira T
    19. Jackson AU
    20. Karjalainen J
    21. Lo KS
    22. Locke AE
    23. Mägi R
    24. Mihailov E
    25. Porcu E
    26. Randall JC
    27. Scherag A
    28. Vinkhuyzen AAE
    29. Westra HJ
    30. Winkler TW
    31. Workalemahu T
    32. Zhao JH
    33. Absher D
    34. Albrecht E
    35. Anderson D
    36. Baron J
    37. Beekman M
    38. Demirkan A
    39. Ehret GB
    40. Feenstra B
    41. Feitosa MF
    42. Fischer K
    43. Fraser RM
    44. Goel A
    45. Gong J
    46. Justice AE
    47. Kanoni S
    48. Kleber ME
    49. Kristiansson K
    50. Lim U
    51. Lotay V
    52. Lui JC
    53. Mangino M
    54. Leach IM
    55. Medina-Gomez C
    56. Nalls MA
    57. Nyholt DR
    58. Palmer CD
    59. Pasko D
    60. Pechlivanis S
    61. Prokopenko I
    62. Ried JS
    63. Ripke S
    64. Shungin D
    65. Stancáková A
    66. Strawbridge RJ
    67. Sung YJ
    68. Tanaka T
    69. Teumer A
    70. Trompet S
    71. van der Laan SW
    72. van Setten J
    73. Van Vliet-Ostaptchouk JV
    74. Wang Z
    75. Yengo L
    76. Zhang W
    77. Afzal U
    78. Ärnlöv J
    79. Arscott GM
    80. Bandinelli S
    81. Barrett A
    82. Bellis C
    83. Bennett AJ
    84. Berne C
    85. Blüher M
    86. Bolton JL
    87. Böttcher Y
    88. Boyd HA
    89. Bruinenberg M
    90. Buckley BM
    91. Buyske S
    92. Caspersen IH
    93. Chines PS
    94. Clarke R
    95. Claudi-Boehm S
    96. Cooper M
    97. Daw EW
    98. De Jong PA
    99. Deelen J
    100. Delgado G
    101. Denny JC
    102. Dhonukshe-Rutten R
    103. Dimitriou M
    104. Doney ASF
    105. Dörr M
    106. Eklund N
    107. Eury E
    108. Folkersen L
    109. Garcia ME
    110. Geller F
    111. Giedraitis V
    112. Go AS
    113. Grallert H
    114. Grammer TB
    115. Gräßler J
    116. Grönberg H
    117. de Groot L
    118. Groves CJ
    119. Haessler J
    120. Hall P
    121. Haller T
    122. Hallmans G
    123. Hannemann A
    124. Hartman CA
    125. Hassinen M
    126. Hayward C
    127. Heard-Costa NL
    128. Helmer Q
    129. Hemani G
    130. Henders AK
    131. Hillege HL
    132. Hlatky MA
    133. Hoffmann W
    134. Hoffmann P
    135. Holmen O
    136. Houwing-Duistermaat JJ
    137. Illig T
    138. Isaacs A
    139. James AL
    140. Jeff J
    141. Johansen B
    142. Johansson Å
    143. Jolley J
    144. Juliusdottir T
    145. Junttila J
    146. Kho AN
    147. Kinnunen L
    148. Klopp N
    149. Kocher T
    150. Kratzer W
    151. Lichtner P
    152. Lind L
    153. Lindström J
    154. Lobbens S
    155. Lorentzon M
    156. Lu Y
    157. Lyssenko V
    158. Magnusson PKE
    159. Mahajan A
    160. Maillard M
    161. McArdle WL
    162. McKenzie CA
    163. McLachlan S
    164. McLaren PJ
    165. Menni C
    166. Merger S
    167. Milani L
    168. Moayyeri A
    169. Monda KL
    170. Morken MA
    171. Müller G
    172. Müller-Nurasyid M
    173. Musk AW
    174. Narisu N
    175. Nauck M
    176. Nolte IM
    177. Nöthen MM
    178. Oozageer L
    179. Pilz S
    180. Rayner NW
    181. Renstrom F
    182. Robertson NR
    183. Rose LM
    184. Roussel R
    185. Sanna S
    186. Scharnagl H
    187. Scholtens S
    188. Schumacher FR
    189. Schunkert H
    190. Scott RA
    191. Sehmi J
    192. Seufferlein T
    193. Shi J
    194. Silventoinen K
    195. Smit JH
    196. Smith AV
    197. Smolonska J
    198. Stanton AV
    199. Stirrups K
    200. Stott DJ
    201. Stringham HM
    202. Sundström J
    203. Swertz MA
    204. Syvänen AC
    205. Tayo BO
    206. Thorleifsson G
    207. Tyrer JP
    208. van Dijk S
    209. van Schoor NM
    210. van der Velde N
    211. van Heemst D
    212. van Oort FVA
    213. Vermeulen SH
    214. Verweij N
    215. Vonk JM
    216. Waite LL
    217. Waldenberger M
    218. Wennauer R
    219. Wilkens LR
    220. Willenborg C
    221. Wilsgaard T
    222. Wojczynski MK
    223. Wong A
    224. Wright AF
    225. Zhang Q
    226. Arveiler D
    227. Bakker SJL
    228. Beilby J
    229. Bergman RN
    230. Bergmann S
    231. Biffar R
    232. Blangero J
    233. Boomsma DI
    234. Bornstein SR
    235. Bovet P
    236. Brambilla P
    237. Brown MJ
    238. Campbell H
    239. Caulfield MJ
    240. Chakravarti A
    241. Collins R
    242. Collins FS
    243. Crawford DC
    244. Cupples LA
    245. Danesh J
    246. de Faire U
    247. den Ruijter HM
    248. Erbel R
    249. Erdmann J
    250. Eriksson JG
    251. Farrall M
    252. Ferrannini E
    253. Ferrières J
    254. Ford I
    255. Forouhi NG
    256. Forrester T
    257. Gansevoort RT
    258. Gejman PV
    259. Gieger C
    260. Golay A
    261. Gottesman O
    262. Gudnason V
    263. Gyllensten U
    264. Haas DW
    265. Hall AS
    266. Harris TB
    267. Hattersley AT
    268. Heath AC
    269. Hengstenberg C
    270. Hicks AA
    271. Hindorff LA
    272. Hingorani AD
    273. Hofman A
    274. Hovingh GK
    275. Humphries SE
    276. Hunt SC
    277. Hypponen E
    278. Jacobs KB
    279. Jarvelin MR
    280. Jousilahti P
    281. Jula AM
    282. Kaprio J
    283. Kastelein JJP
    284. Kayser M
    285. Kee F
    286. Keinanen-Kiukaanniemi SM
    287. Kiemeney LA
    288. Kooner JS
    289. Kooperberg C
    290. Koskinen S
    291. Kovacs P
    292. Kraja AT
    293. Kumari M
    294. Kuusisto J
    295. Lakka TA
    296. Langenberg C
    297. Le Marchand L
    298. Lehtimäki T
    299. Lupoli S
    300. Madden PAF
    301. Männistö S
    302. Manunta P
    303. Marette A
    304. Matise TC
    305. McKnight B
    306. Meitinger T
    307. Moll FL
    308. Montgomery GW
    309. Morris AD
    310. Morris AP
    311. Murray JC
    312. Nelis M
    313. Ohlsson C
    314. Oldehinkel AJ
    315. Ong KK
    316. Ouwehand WH
    317. Pasterkamp G
    318. Peters A
    319. Pramstaller PP
    320. Price JF
    321. Qi L
    322. Raitakari OT
    323. Rankinen T
    324. Rao DC
    325. Rice TK
    326. Ritchie M
    327. Rudan I
    328. Salomaa V
    329. Samani NJ
    330. Saramies J
    331. Sarzynski MA
    332. Schwarz PEH
    333. Sebert S
    334. Sever P
    335. Shuldiner AR
    336. Sinisalo J
    337. Steinthorsdottir V
    338. Stolk RP
    339. Tardif JC
    340. Tönjes A
    341. Tremblay A
    342. Tremoli E
    343. Virtamo J
    344. Vohl MC
    345. Amouyel P
    346. Asselbergs FW
    347. Assimes TL
    348. Bochud M
    349. Boehm BO
    350. Boerwinkle E
    351. Bottinger EP
    352. Bouchard C
    353. Cauchi S
    354. Chambers JC
    355. Chanock SJ
    356. Cooper RS
    357. de Bakker PIW
    358. Dedoussis G
    359. Ferrucci L
    360. Franks PW
    361. Froguel P
    362. Groop LC
    363. Haiman CA
    364. Hamsten A
    365. Hayes MG
    366. Hui J
    367. Hunter DJ
    368. Hveem K
    369. Jukema JW
    370. Kaplan RC
    371. Kivimaki M
    372. Kuh D
    373. Laakso M
    374. Liu Y
    375. Martin NG
    376. März W
    377. Melbye M
    378. Moebus S
    379. Munroe PB
    380. Njølstad I
    381. Oostra BA
    382. Palmer CNA
    383. Pedersen NL
    384. Perola M
    385. Pérusse L
    386. Peters U
    387. Powell JE
    388. Power C
    389. Quertermous T
    390. Rauramaa R
    391. Reinmaa E
    392. Ridker PM
    393. Rivadeneira F
    394. Rotter JI
    395. Saaristo TE
    396. Saleheen D
    397. Schlessinger D
    398. Slagboom PE
    399. Snieder H
    400. Spector TD
    401. Strauch K
    402. Stumvoll M
    403. Tuomilehto J
    404. Uusitupa M
    405. van der Harst P
    406. Völzke H
    407. Walker M
    408. Wareham NJ
    409. Watkins H
    410. Wichmann HE
    411. Wilson JF
    412. Zanen P
    413. Deloukas P
    414. Heid IM
    415. Lindgren CM
    416. Mohlke KL
    417. Speliotes EK
    418. Thorsteinsdottir U
    419. Barroso I
    420. Fox CS
    421. North KE
    422. Strachan DP
    423. Beckmann JS
    424. Berndt SI
    425. Boehnke M
    426. Borecki IB
    427. McCarthy MI
    428. Metspalu A
    429. Stefansson K
    430. Uitterlinden AG
    431. van Duijn CM
    432. Franke L
    433. Willer CJ
    434. Price AL
    435. Lettre G
    436. Loos RJF
    437. Weedon MN
    438. Ingelsson E
    439. O’Connell JR
    440. Abecasis GR
    441. Chasman DI
    442. Goddard ME
    443. Visscher PM
    444. Hirschhorn JN
    445. Frayling TM
    (2014) Defining the role of common variation in the genomic and biological architecture of adult human height
    Nature Genetics 46:1173–1186.
    https://doi.org/10.1038/ng.3097
    1. Yue F
    2. Cheng Y
    3. Breschi A
    4. Vierstra J
    5. Wu W
    6. Ryba T
    7. Sandstrom R
    8. Ma Z
    9. Davis C
    10. Pope BD
    11. Shen Y
    12. Pervouchine DD
    13. Djebali S
    14. Thurman RE
    15. Kaul R
    16. Rynes E
    17. Kirilusha A
    18. Marinov GK
    19. Williams BA
    20. Trout D
    21. Amrhein H
    22. Fisher-Aylor K
    23. Antoshechkin I
    24. DeSalvo G
    25. See L-H
    26. Fastuca M
    27. Drenkow J
    28. Zaleski C
    29. Dobin A
    30. Prieto P
    31. Lagarde J
    32. Bussotti G
    33. Tanzer A
    34. Denas O
    35. Li K
    36. Bender MA
    37. Zhang M
    38. Byron R
    39. Groudine MT
    40. McCleary D
    41. Pham L
    42. Ye Z
    43. Kuan S
    44. Edsall L
    45. Wu Y-C
    46. Rasmussen MD
    47. Bansal MS
    48. Kellis M
    49. Keller CA
    50. Morrissey CS
    51. Mishra T
    52. Jain D
    53. Dogan N
    54. Harris RS
    55. Cayting P
    56. Kawli T
    57. Boyle AP
    58. Euskirchen G
    59. Kundaje A
    60. Lin S
    61. Lin Y
    62. Jansen C
    63. Malladi VS
    64. Cline MS
    65. Erickson DT
    66. Kirkup VM
    67. Learned K
    68. Sloan CA
    69. Rosenbloom KR
    70. Lacerda de Sousa B
    71. Beal K
    72. Pignatelli M
    73. Flicek P
    74. Lian J
    75. Kahveci T
    76. Lee D
    77. James Kent W
    78. Ramalho Santos M
    79. Herrero J
    80. Notredame C
    81. Johnson A
    82. Vong S
    83. Lee K
    84. Bates D
    85. Neri F
    86. Diegel M
    87. Canfield T
    88. Sabo PJ
    89. Wilken MS
    90. Reh TA
    91. Giste E
    92. Shafer A
    93. Kutyavin T
    94. Haugen E
    95. Dunn D
    96. Reynolds AP
    97. Neph S
    98. Humbert R
    99. Scott Hansen R
    100. De Bruijn M
    101. Selleri L
    102. Rudensky A
    103. Josefowicz S
    104. Samstein R
    105. Eichler EE
    106. Orkin SH
    107. Levasseur D
    108. Papayannopoulou T
    109. Chang K-H
    110. Skoultchi A
    111. Gosh S
    112. Disteche C
    113. Treuting P
    114. Wang Y
    115. Weiss MJ
    116. Blobel GA
    117. Cao X
    118. Zhong S
    119. Wang T
    120. Good PJ
    121. Lowdon RF
    122. Adams LB
    123. Zhou X-Q
    124. Pazin MJ
    125. Feingold EA
    126. Wold B
    127. Taylor J
    128. Mortazavi A
    129. Weissman SM
    130. Stamatoyannopoulos JA
    131. Snyder MP
    132. Guigo R
    133. Gingeras TR
    134. Gilbert DM
    135. Hardison RC
    136. Beer MA
    137. Ren B
    138. The Mouse ENCODE Consortium
    (2014) A comparative encyclopedia of DNA elements in the mouse genome
    Nature 515:355–364.
    https://doi.org/10.1038/nature13992

Decision letter

  1. Cheryl Ackert-Bicknell
    Reviewing Editor; University of Colorado, United States
  2. George H Perry
    Senior Editor; Pennsylvania State University, United States
  3. Laura Saba
    Reviewer
  4. Peter Claes
    Reviewer; Medical Imaging Research Center, Belgium

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Shapes and genescapes: Mapping multivariate phenotype-biological process associations for craniofacial shape" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and George Perry as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Laura Saba (Reviewer #1); Peter Claes (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions

All three reviewers emphatically stated they found the discussion and conclusions of this paper to be overstated, and not completely supported by the data. That said, they unanimously felt that a substantial revision of this work could deal with this problem and hence a revision opportunity was suggested over an outright rejection. Overall, the requested revisions do not require new data, but do require some reanalysis of the data. These issues are presented in greater detail in the point by point revisions listed below. In addition, when the reviewers tried to use the web app on the browsers Chrome and Firefox, just when the results appeared, they received an error that it disconnected from the server. Please check the web app stability and where necessary provide guidance on browser use for both Mac and PCs.

1. The use of the word 'process' throughout was not initially clear. It would help if it were clearly defined early on.

2. There is some inconsistency between the results, Figure 2A, and the discussion about how Nov/Ccn3 is annotated.

3. For the likelihood ratio test for differences in the distribution of cell sizes between Bmpr1b mutants and controls, it is important that all the interaction effects that included genotype were eliminated as well.

4. The matric used for judging the stability of the cluster results is often referred to as the cophenetic correlation coefficient. This common terminology may help to clarify the first sentence of the paragraph that starts on line 287. It isn't clear why a t-statistics was used to determine if the correlation coefficient was significantly different from zero. Perhaps a reference is needed.

5. In the 'Process effects in mutant morphospace', the authors point out similarities in the direction of the process vector to some of the mutants. It would be impressive if those process/mutant pairs annotated and if those genes had high weights in the MGP genetic vector. Similar comment for the gene/process pairs in the heatmap.

6. The conclusions starting on line 462 and ending on line 466 were not tested rigorously.

7. More detail about how the mutant MANOVA coefficients were calculated is needed.

8. It is not clear how the regularization parameter was chosen.

9. The paper is presented as the development of a new approach; however, it is not written in that way. The methods section "Regularized PLS analysis" is not accessible for a general public. If the reader is not already familiar with multivariate PLS, they will have no chance of understanding what this paper is all about. The results don't mention the actual genotype-phenotype approach, a description of the process-based MGP mapping should be done first thing in the results. Accordingly, the first figure should be current figure 10. And hopefully the legend will include more details.

10. The approach provides the phenotypic effect of genetic variation in already known pathways but it does not result in new genotype-phenotype associations; this is acknowledged in the text. However, the manuscript suggest that the results generate testable hypothesis which is just not supported by the provided data. The predictions of phenotypic effects are done on genes that are previously known to be involved in pathways that result in such phenotypes. For example, the Bmpr1b gene is annotated in the 'Chondrocyte Differentiation' biological process and genetic variation in such gene is therefore used for the genotype-phenotype analysis, but then, the authors suggest that given the results of the mapping, Bmpr1b should show chondrocyte development phenotypes, which is a circular argument. This example is also true for the other two predictions described in the manuscript. A way of showing that the approach results in testable hypothesis will be to run the genotype-phenotype mapping on pathways that are not previously known to affect craniofacial shape, and then test whether the high loading genes indeed affect the phenotype in the expected way. Until that is shown, the claim that this approach results in new testable hypothesis should be taking with caution and this paper must be revised to reflect this reality.

11. The approach is not new. Multiple other methods perform versions of multivariate genotype-phenotype mapping (this is discussed in the manuscript), but more importantly, the authors state to be using a method previously published by Mitteroecker et al. 2016 with the twist of restricting the analysis to known biological processes. It is not clear in the manuscript how much of their approach is actually new and how much is Mitteroecker's applied to a subset of markers. To highlight some lines where it sounds as if the authors did develop the method are found on lines 390-392 and lines 423-425, but rereading of Mitteroecker et al. 2016 suggests this is not the case. This is a crucial point that is not explained nor discussed at all.

12. One of the challenges with multivariate analyses of this type is how to measure success of the model. In this case, the authors compared their genotype-phenotype results to phenotype results from genetically manipulated mice. As written, it is not clear where the information from the mutant mice came from. Although not explicitly stated, it appears that the mutant data is assigned to a process simply based on the inclusion of the mutated gene in the process being tested. There was no consideration of the genetic loadings of that particular gene in the original MPG model. If the mutated gene is in the process but has very weak genetic loadings, the assumption that the mutant phenotype should be similar to the phenotype associated with the process genetic loading is suspect. Theoretically, the MPG predicted phenotype in the mutated rodent would be more accurate if the mutated gene had a stronger genetic loading. Furthermore, if the mutant data was derived from the literature than it likely includes some publication bias. All genes included likely have some effect on cranial morphology or their results likely would not have been published.

13. Within the manuscript, there is an emphasis on the concordant direction of association between the process MGP axis and the axis of shape variation of a relevant mutant phenotype. However, some ambiguity remains about the relevance of the direction of association between the direction of the process MGP axis and the axis of shape variation of a relevant mutant phenotype. The loadings per gene are presented as negative/positive values in all figures, but it is never explained in the methods what does the sign mean. From the description in the methods, the direction of the correlation between the process PLS1 phenotypic loadings and the mutant MANOVA coefficients cannot be predicted unless you knew that the mutation of the gene was similar to one or more of the 8 founder haplotypes. This would require functional knowledge about the markers/haplotype included in the genetic loadings. Explain this in the methods and figure legends to help the reader actually understand what the PLS is doing.

14. There was also some confusion surrounding the rationale and methodology for estimating the null distribution of the R2 values. More detail is needed on how random marker selection for permutation was constrained to follow similar patterns of linkage disequilibrium to the observed marker set. Also, it is unclear what the rationale for randomly selecting markers rather than randomly selecting genes to generate the null distribution.

15. In general, the section 'Comparison of processes to principal component' lacks statistical rigor and conclusions are supported by the data. It is not relevant to compare all 1000 processes with the principal components without first filtering based on R2 values from the MPG analysis. The overall value of the analysis in this section is unclear.

16. In the discussion, the authors do not make a compelling case for how these types of analyses could be used to generate hypotheses for further studies.

17. As stated in the introduction this is a study of complex relationships between genotypes and phenotypes, and it is important to keep in mind that the technique used here as many others in similar endeavors, remain linear in nature, and therefore limited in exposing potentially non-linear relationships.

18. While the three example processes are interesting and easy to understand or follow in terms of, how this is of interest. I do struggle with the interpretation of the follow-up analyses, including correlating effects of processes, comparing of processes to principal component directions and process effects in mutant morpho-space. It is unclear if all these analyses make sense, and more importantly, the dimensionality into which these explorations are being made is rather low. First, any process effect is modelled by the first PLS-component only. Therefore, is the process and its effect on craniofacial shape modelled completely? Probably not. Second, the morpho-space in the last two mentioned analyses here above, is restricted to the first two principal components only. Therefore, is craniofacial variability modelled completely? Probably not. As a result, Figure 7 A, especially, represents only a sub-dimensional view of process effects versus mutant phenotypes. Discrepancies observed, e.g., none of the Bmp mutants align with the bmp signaling pathway, can be due to the limited linear dimensionality, or because of other underlying genuine biological explanations. However, it is hard to tell.

19. It is suggested by one of the reviewer that a better overlap between PLS and CCA and the use of both in genotype-phenotype associations is welcome. E.g., the discussion simply states, "Our approach also differs from methods that associate single locus effects with a multivariate phenotype (Claes et al., 2018)" without any explanation or insight how and what is different. In fact, the difference is not that great, besides the difference in underlying paradigm (GWAS versus candidate selection) and looking across multiple genes instead, which is the forte of this work. The incorporation of gene-based or haplotype-based MV GWAS into the discussion or introduction is also encouraged.

20. The reviewers need to make a stronger case about the benefit of groups of genes, because of the regularization and the selection of a 1-dimensional latent variable. I.e., It is of interest, for at least the three examples given, to run a gene-by-gene based analysis alongside, and to verify the benefit of the process-based analysis.

21. It is suggest that the authors investigate the effect of the user-defined regularization on their results. Additional cross-validation based results to see how good the first PLS component generalizes is needed. It is recognized that collecting a completely independent replication cohort is outside the scope or typical research resources. However, a more elaborated investigation using cross-validations might help.

22. Somewhat related to the point above, the authors need to provide more information on the dimensionality reductions performed, e.g., please present CV results using multiple PLS components, and provide insights on how many components are expected to be relevant in describing a process and its effect on craniofacial variation. The same applies for the 2-PC morphospace, of interest is to provide insight into choosing only 2PCs and no more?

23. The reviewers found the assessment cluster stability to be circular. Please include some references illustrating this approach or to investigate alternative measures for clustering techniques in machine learning, e.g., the silhouette score.

24. The Comparison of processes to principal component directions, ends with listing processes correlated to PCs. It was hard finding out if these groups listed made sense to be grouped as such. Additional and more explicit information for the reader in those directions is welcome here and in other places where such listing of processes (e.g., the clustering) is done.

25. How are the authors estimating % total phenotypic variance explained by the PLS genotype vector (lines 193, 226)? Is this the R2 they talk about in the methods? If so, please use the same terminology in results and methods sections.

26. The authors show that the %var explained in the first three examples (chondrocytes, palate, and symmetry) is higher than randomly chosen set of genes. But maybe, an additional biologically important point to test is whether pathways that are known to be involved in craniofacial development explain significantly higher phenotypic variance than pathways that are expected not to contribute to such phenotypes. This, besides being a testable hypothesis, will shed light into our understanding of the GO terms annotated to affect craniofacial development, and probably highlight pathways previously not associated with such processes that could be followed up in future work.

27. It is highlighted as a main finding that the phenotypic effect of different biological processes correlate with each other and with single mutants, however it is not discussed how much of that correlation is driven by different processes sharing genes. Could you elaborate on that aspect of the analysis? How many genes are shared and whether that correlates with the correlation of vector directions? I wouldn't be sure about how solid the main finding is unless it is clearly shown that it is not coming from shared genetic effects.

28. Supp Figure 4 shows several covariates for the data, however the methods only described how 'generation' was accounted for. What about sex? Mouse source? I don't have a way of knowing whether you are using them as covariates in your model, or whether you can actually add covariates to the model, because the model is not described in the methods.

29. The Results section is dry and doesn't provide any conclusion for the different sections. This could be improved if the conclusions of the analyses that are currently described in the Discussion section be moved to the relevant Results sections. This will also benefit the discussion that is currently very long and most of it is an actual description of results.

30. The information necessary to understand the figures should be in the corresponding figure legend not in the main text (e.g. l 327-329, l 331-333, l346-349). And, all figures showing results of the MGP analysis will benefit from ranking the genes according to variance of the loadings instead of alphabetical order – so it becomes easier to asses which genes are more important.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Relating Multivariate Shapes to Genescapes Using Phenotype-Biological Process Associations for Craniofacial Shape" for further consideration by eLife. Your revised article has been reviewed by 3 peer reviewers and the evaluation has been overseen by George Perry as the Senior Editor, and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below. The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. It is recognized that this is the second request for revisions on this paper.

Essential revisions:

A lengthy discussion was had among the reviewers about clarifying what is meant by new hypothesis. It was agreed that this issue must be addressed, but not does not necessarily need new data or experiments. Two of the three points below directly center on this issue.

1) In line 327-329 (of the tracked changes version) the authors write" "The majority of process MGP analyses demonstrate a similar importance to many alleles, highlighting the main strength process-level analyses over individual marker tests". This is an interesting claim, but what data supports this statement exactly? It would be great if this is supported by some density plot with for example # of genes ranking high (whatever threshold you choose) for each process. The authors have tested ~30 or so processes for some of the analyses. Please clarify if this claim coming from there?

2) The reviews all agreed that there was some level of circularity in the examples shown by the authors. The reviewers strongly felt that the response to comment #10 was not convincing, and the main text was not sufficiently modified to clarify this point. The reviewers agreed that if they remained confused on this point, readers of this manuscript likely would be too. The following is a second attempt at explaining the reviewers point again: all three mutants tested (Bmpr1, Fgf10, and Ankrd11) have previous published data showing their effects on craniofacial shape, sometimes even directly related to the biological process the authors are testing. So, given a) they are functionally annotated to be part of that process, b) there are previous publications -possibly the annotation is derived from those studies, although not necessarily- showing their phenotype effect on craniofacial shape, it does not sound as the authors are generating new hypothesis with their method. This would sound more convincing for the readers if the authors explained better for each example what is their rational for calling each of them a new hypothesis generated by the method. Explicitly list what is previously known or not known about those genes and their effect on craniofacial shape (they reference previous studies, but not in the context of framing how this is different from what they are doing), and be emphatically clear about what is the new "thing" they are trying to show. The reviewers are not saying it is not great to see those genes highly ranked within the process MPG vector, but they are just saying that strong claims of this being a new hypothesis are being made when working with genes for which an effect on craniofacial shape has been shown before. This exact point was also raised in query #16 but there was no new text in the Discussion addressing this point.

3. Related to the point above, there were several comments in by the reviewer previously on how, with the current data, the authors could approach more directly the 'testable new hypotheses' that this method generates, but they were not really implemented. This of course does not diminish the value of the paper, but the opposite would have increased its ability to prove that new hypotheses are indeed generated and proven true (not related to genes obviously related to craniofacial shape). For example, point #26 suggests to-test effects on processes that are not related/or expected to be related to craniofacial shape. The authors reply that is difficult to delineate such processes. The reviewers wondered about considering or at least mentioning 'gonad /sperm / egg development'? Social / mating behavior? Circadian rhythm? Etc.

https://doi.org/10.7554/eLife.68623.sa1

Author response

Essential revisions

All three reviewers emphatically stated they found the discussion and conclusions of this paper to be overstated, and not completely supported by the data. That said, they unanimously felt that a substantial revision of this work could deal with this problem and hence a revision opportunity was suggested over an outright rejection. Overall, the requested revisions do not require new data, but do require some reanalysis of the data. These issues are presented in greater detail in the point by point revisions listed below. In addition, when the reviewers tried to use the web app on the browsers Chrome and Firefox, just when the results appeared, they received an error that it disconnected from the server. Please check the web app stability and where necessary provide guidance on browser use for both Mac and PCs.

1. The use of the word 'process' throughout was not initially clear. It would help if it were clearly defined early on.

The reviewer is correct. We now realize that this word occurs many times in the paper and is not used consistently. In the majority of contexts, we have used the term to essentially mean “gene annotation” and that has now been corrected this in the revision. We now define process, in the context of gene annotation and that is now corrected in the paper. Where we use “process” we are now clear about the context and whether we are referring to genes linked by membership in a pathway or genes with known associations to a subcellular, cellular or tissue level biological activity.

2. There is some inconsistency between the results, Figure 2A, and the discussion about how Nov/Ccn3 is annotated.

This has been corrected.

3. For the likelihood ratio test for differences in the distribution of cell sizes between Bmpr1b mutants and controls, it is important that all the interaction effects that included genotype were eliminated as well.

The reduced model for the likelihood ratio test did not include interaction terms for genotype and distribution. This has been clarified in the manuscript.

4. The matric used for judging the stability of the cluster results is often referred to as the cophenetic correlation coefficient. This common terminology may help to clarify the first sentence of the paragraph that starts on line 287. It isn't clear why a t-statistics was used to determine if the correlation coefficient was significantly different from zero. Perhaps a reference is needed.

We thank the reviewer for this suggestion and have revised the text to use this term. The reason for the t-test is to determine whether the correlation coefficient is real – that is different from the range one might expect under the null hypothesis of no correlation.

5. In the 'Process effects in mutant morphospace', the authors point out similarities in the direction of the process vector to some of the mutants. It would be impressive if those process/mutant pairs annotated and if those genes had high weights in the MGP genetic vector. Similar comment for the gene/process pairs in the heatmap.

Our results show that one can get directions of change in morphospace that resemble those of mutations of major effect via many different combinations of SNPs. This implies that these directions are determined by developmental processes that can be influenced in multiple ways to produce the same direction of effect. This is, in fact, the central point of the paper. We have clarified the text to make this easier to understand.

6. The conclusions starting on line 462 and ending on line 466 were not tested rigorously.

We agree with the reviewers that testing the significance of process-associated directions in morphospace is currently theoretically intractable. Although we believe there is value in the analysis we presented, we agree that the patterns cannot be tested rigorously. For this reason, we have removed these analyses from the paper. Instead, we have followed the reviewers suggestions to focus more on testing the significance of single vs joint SNP effects on directions of variation. We believe that this change strengthens rather than detracts from the paper.

7. More detail about how the mutant MANOVA coefficients were calculated is needed.

We have provided more detail to clarify this in the statistical results and comparisons section in the methods.

8. It is not clear how the regularization parameter was chosen.

As the reviewer is probably intimating, the choice of regularization parameter represents a tradeoff between the risk of overfitting and the risk of not seeing real patterns in the data. In this case, the choice of parameter was driven by our goal to minimize model error without overly penalizing the phenotypic response in comparison to mutant phenotypes. We have clarified this in the results describing the MGP method and have added a supplemental figure showing a range of regularization parameters for all of the example analyses. The regularization parameter can also be changed in the app for users to experiment with settings and determine the effect on the results.

9. The paper is presented as the development of a new approach; however, it is not written in that way. The methods section "Regularized PLS analysis" is not accessible for a general public. If the reader is not already familiar with multivariate PLS, they will have no chance of understanding what this paper is all about. The results don't mention the actual genotype-phenotype approach, a description of the process-based MGP mapping should be done first thing in the results. Accordingly, the first figure should be current figure 10. And hopefully the legend will include more details.

This is an excellent point. We have moved a general description of the method to the introduction and have moved Figure 10 to Figure 1. Figure 10 is actually placed where it is because of the eLife format which has methods coming after results. This is easily fixed and we agree with the reviewer on this.

10. The approach provides the phenotypic effect of genetic variation in already known pathways but it does not result in new genotype-phenotype associations; this is acknowledged in the text. However, the manuscript suggest that the results generate testable hypothesis which is just not supported by the provided data. The predictions of phenotypic effects are done on genes that are previously known to be involved in pathways that result in such phenotypes. For example, the Bmpr1b gene is annotated in the 'Chondrocyte Differentiation' biological process and genetic variation in such gene is therefore used for the genotype-phenotype analysis, but then, the authors suggest that given the results of the mapping, Bmpr1b should show chondrocyte development phenotypes, which is a circular argument. This example is also true for the other two predictions described in the manuscript. A way of showing that the approach results in testable hypothesis will be to run the genotype-phenotype mapping on pathways that are not previously known to affect craniofacial shape, and then test whether the high loading genes indeed affect the phenotype in the expected way. Until that is shown, the claim that this approach results in new testable hypothesis should be taking with caution and this paper must be revised to reflect this reality.

We don’t agree with the logic of this comment. The point here, as we understand it, as that there is circularity or ascertainment bias introduced by the fact that gene annotations are based on known phenotypic effects. The reviewer would be correct in this if the phenotypic outcome in our model was the same thing as the annotation. In other words, the observation that the annotation, “craniofacial shape” reveals a differs in craniofacial shape would be uninteresting and circular. However, for the vast majority of process annotations, the relationship between craniofacial shape and the process represented by that annotation is not known. For the example given, there is no known general association between “chondrocyte proliferation” and craniofacial shape. Most of those annotations would have been derived from studies of cartilage or bone development. It is doubtful that any derive primarily from studies of mouse craniofacial morphology. Thus, we don’t see how the circularity that is stated by the reviewer would actually arise in our data.

11. The approach is not new. Multiple other methods perform versions of multivariate genotype-phenotype mapping (this is discussed in the manuscript), but more importantly, the authors state to be using a method previously published by Mitteroecker et al. 2016 with the twist of restricting the analysis to known biological processes. It is not clear in the manuscript how much of their approach is actually new and how much is Mitteroecker's applied to a subset of markers. To highlight some lines where it sounds as if the authors did develop the method are found on lines 390-392 and lines 423-425, but rereading of Mitteroecker et al. 2016 suggests this is not the case. This is a crucial point that is not explained nor discussed at all.

Mitteroecker explains a framework for many-to-many genomic models. We have taken this framework and applied it to a gene subset strategy. Furthermore, we have introduced regularization to mitigate the overfitting effects commonly observed with models with many parameters. In addition, we have applied this method on a dataset with much higher genotyping SNP resolution and have provided an app that allows users to pursue their own analyses quickly. The combination of these factors are novel, despite previous use of PLS in the genomics literature. We have edited the language in the manuscript to more carefully detail what contributions of this work are novel.

12. One of the challenges with multivariate analyses of this type is how to measure success of the model. In this case, the authors compared their genotype-phenotype results to phenotype results from genetically manipulated mice. As written, it is not clear where the information from the mutant mice came from. Although not explicitly stated, it appears that the mutant data is assigned to a process simply based on the inclusion of the mutated gene in the process being tested. There was no consideration of the genetic loadings of that particular gene in the original MPG model. If the mutated gene is in the process but has very weak genetic loadings, the assumption that the mutant phenotype should be similar to the phenotype associated with the process genetic loading is suspect. Theoretically, the MPG predicted phenotype in the mutated rodent would be more accurate if the mutated gene had a stronger genetic loading. Furthermore, if the mutant data was derived from the literature than it likely includes some publication bias. All genes included likely have some effect on cranial morphology or their results likely would not have been published.

Aside from annotation data, the data presented in this paper are not derived from the literature. The mutants that are presented in this paper and that are available for comparison in the app represent data collection by the Hallgrimsson and Marcucio labs. All of the data were obtained from scans in the Hallgrimsson lab and none of these data were obtained from the literature. We have clarified the origin of the data in the methods section.

13. Within the manuscript, there is an emphasis on the concordant direction of association between the process MGP axis and the axis of shape variation of a relevant mutant phenotype. However, some ambiguity remains about the relevance of the direction of association between the direction of the process MGP axis and the axis of shape variation of a relevant mutant phenotype. The loadings per gene are presented as negative/positive values in all figures, but it is never explained in the methods what does the sign mean. From the description in the methods, the direction of the correlation between the process PLS1 phenotypic loadings and the mutant MANOVA coefficients cannot be predicted unless you knew that the mutation of the gene was similar to one or more of the 8 founder haplotypes. This would require functional knowledge about the markers/haplotype included in the genetic loadings. Explain this in the methods and figure legends to help the reader actually understand what the PLS is doing.

The loadings represent the relative strength of allelic effects in the Diversity Outbred genome. There is no clear interpretation of individual effects of an ordination approach like PLS. However, we have now added an analysis of the effect of the removal of highly-loaded genes in the analysis. This result is now presented in figure 6 of the manuscript.

Except for the rare case of novel mutations that arose during the breeding of the DO mice, all of the SNP variants represent genomic sequences of the original 8 founder strains. It is not clear to us what is meant by similarity to the 8 founder strains. The SNP-phenotype associations are what they are, as is the case in any other genetic study. However, it is in the mutant to MGP comparisons that knowledge of the functional effect is relevant. In those cases, we are testing the similarity between genotype-phenotype associations as determined in the DO data to known function to phenotype associations. Those latter associations are based on studies of those mutants and in one case, we have actually predicted a functional effect from a genotype-phenotype association in the DO data. Beyond this, a SNP by SNP functional analysis of the genomic diversity of the DO mice is obviously intractable.

14. There was also some confusion surrounding the rationale and methodology for estimating the null distribution of the R2 values. More detail is needed on how random marker selection for permutation was constrained to follow similar patterns of linkage disequilibrium to the observed marker set. Also, it is unclear what the rationale for randomly selecting markers rather than randomly selecting genes to generate the null distribution.

We have changed the methodology for estimating the null distribution of R2 values to permute the phenotypes instead. This avoids the disadvantages of random marker selection and attempting to match the original linkage disequilibrium pattern of the observed marker set. These changes have been reflected in the methods as well as the figures.

15. In general, the section 'Comparison of processes to principal component' lacks statistical rigor and conclusions are supported by the data. It is not relevant to compare all 1000 processes with the principal components without first filtering based on R2 values from the MPG analysis. The overall value of the analysis in this section is unclear.

We agree and this section has been removed. See response to point 6.

16. In the discussion, the authors do not make a compelling case for how these types of analyses could be used to generate hypotheses for further studies.

This is an excellent point. We have included a paragraph in the discussion highlighting how MGP analyses can generate hypotheses for further studies.

17. As stated in the introduction this is a study of complex relationships between genotypes and phenotypes, and it is important to keep in mind that the technique used here as many others in similar endeavors, remain linear in nature, and therefore limited in exposing potentially non-linear relationships.

We agree with the reviewers that linearity is a drawback to this approach. The discussion has been expanded to highlight this caveat more clearly. This same criticism could be levied at the majority of genotype-phenotype mapping approaches.

18. While the three example processes are interesting and easy to understand or follow in terms of, how this is of interest. I do struggle with the interpretation of the follow-up analyses, including correlating effects of processes, comparing of processes to principal component directions and process effects in mutant morphospace. It is unclear if all these analyses make sense, and more importantly, the dimensionality into which these explorations are being made is rather low. First, any process effect is modelled by the first PLS-component only. Therefore, is the process and its effect on craniofacial shape modelled completely? Probably not. Second, the morphospace in the last two mentioned analyses here above, is restricted to the first two principal components only. Therefore, is craniofacial variability modelled completely? Probably not. As a result, Figure 7 A, especially, represents only a sub-dimensional view of process effects versus mutant phenotypes. Discrepancies observed, e.g., none of the Bmp mutants align with the bmp signaling pathway, can be due to the limited linear dimensionality, or because of other underlying genuine biological explanations. However, it is hard to tell.

Although the analyses presented in this work may be held to reduced dimensions, they are ultimately driven by the g-p relationships in the data. This is in contrast to a priori reduction of the phenotype dimensionality followed by genomic analysis. We agree that this caveat needs to be made clear to the reader and have added detail to that end, particularly focusing on figure 8 (previously Figure 7).

19. It is suggested by one of the reviewer that a better overlap between PLS and CCA and the use of both in genotype-phenotype associations is welcome. E.g., the discussion simply states, "Our approach also differs from methods that associate single locus effects with a multivariate phenotype (Claes et al., 2018)" without any explanation or insight how and what is different. In fact, the difference is not that great, besides the difference in underlying paradigm (GWAS versus candidate selection) and looking across multiple genes instead, which is the forte of this work. The incorporation of gene-based or haplotype-based MV GWAS into the discussion or introduction is also encouraged.

We thank the reviewers for this suggestion. We have expanded on the similarity between CCA and PLS in the methods as well as the Discussion section. Additionally, we included review of existing gene-based GWAS attempts in the introduction.

20. The reviewers need to make a stronger case about the benefit of groups of genes, because of the regularization and the selection of a 1-dimensional latent variable. I.e., It is of interest, for at least the three examples given, to run a gene-by-gene based analysis alongside, and to verify the benefit of the process-based analysis.

It should be noted that no direct comparison between a many-to-many approach and a gene-by-gene is possible. That is because each gene-by-gene analysis solves for a unique phenotypic effect in multivariate space that must be interpreted individually. From this perspective alone, a many-to-many approach is advantageous because it returns a single set of phenotypic and genotypic responses that can be interpreted in unison. In addition, the many-to-many approach explicitly models the coordination of genetic effects at the process level. One point of entry into such a comparison is the gene drop analysis that we have implemented in the revision. See response to point 6.

21. It is suggest that the authors investigate the effect of the user-defined regularization on their results. Additional cross-validation based results to see how good the first PLS component generalizes is needed. It is recognized that collecting a completely independent replication cohort is outside the scope or typical research resources. However, a more elaborated investigation using cross-validations might help.

We’ve done cross validation for regularization and have included that result in supplemental figure 1. We thank the reviewer for this suggestion.

22. Somewhat related to the point above, the authors need to provide more information on the dimensionality reductions performed, e.g., please present CV results using multiple PLS components, and provide insights on how many components are expected to be relevant in describing a process and its effect on craniofacial variation. The same applies for the 2-PC morphospace, of interest is to provide insight into choosing only 2PCs and no more?

The 2PC figure is meant to visually introduce the concept of the similarity of directions of effects. This has the advantage of interpreting the entire context for those first 2PCs. For example, figure 7A shows “wnt signaling pathway” MGP has a similar direction of effects to many mutants, like Fgf10, and the reader knows the corresponding phenotypic effect using the heatmap figure along the x-axis. In contrast, figure 7B shows compares the full vector correlations for several processes to the same set of mutants, with the caveat that we cannot provide visualizations of what those phenotypes look like. We have clarified the intent of figure 7 and its components in the text.

23. The reviewers found the assessment cluster stability to be circular. Please include some references illustrating this approach or to investigate alternative measures for clustering techniques in machine learning, e.g., the silhouette score.

The reviewers may be over-interpreting the significance of the clusters. We point of the analysis in Figure 7 is that there is a higher-level structure to the vector correlations among processes and that individual processes don’t tend to point in independent directions. We are not interpreting the specific higher level structure here in terms of the number of clusters or their biological underpinnings. As the reviewers intimate, that would be an entirely new level of analysis in which is beyond the scope of the present paper. However, we did pursue the methodological suggestion of the reviewer and calculate a silhouette score. This actually returns a structure of five clusters which is very similar to what we show in Figure 7. Given that the biological nature of these clusters is not a focus of the current paper, however, we have chosen not to include this result in the paper but provide it as Author response image 1. In terms of the reviewer’s suggestion to add references supporting the use of hierarchical clustering, we now provide an appropriate supporting reference (Sneath and Sokal, 1973) and note that this is a widely used method for clustering multivariate data.

Author response image 1
Download asset Open asset

24. The Comparison of processes to principal component directions, ends with listing processes correlated to PCs. It was hard finding out if these groups listed made sense to be grouped as such. Additional and more explicit information for the reader in those directions is welcome here and in other places where such listing of processes (e.g., the clustering) is done.

This section has been removed. See point 6.

25. How are the authors estimating % total phenotypic variance explained by the PLS genotype vector (lines 193, 226)? Is this the R2 they talk about in the methods? If so, please use the same terminology in results and methods sections.

The R2 is calculated in the same way as described in the methods. This has been clarified in the text.

26. The authors show that the %var explained in the first three examples (chondrocytes, palate, and symmetry) is higher than randomly chosen set of genes. But maybe, an additional biologically important point to test is whether pathways that are known to be involved in craniofacial development explain significantly higher phenotypic variance than pathways that are expected not to contribute to such phenotypes. This, besides being a testable hypothesis, will shed light into our understanding of the GO terms annotated to affect craniofacial development, and probably highlight pathways previously not associated with such processes that could be followed up in future work.

We thank the reviewers for this suggestion. It is certainly of interest to see if known craniofacial processes affect craniofacial variation in comparison to processes with no obvious relationship to craniofacial development. However, it is difficult to delineate between processes that are known to be involved in craniofacial development vs processes that are not known to be involved in craniofacial development. For example, cell cycle process related gene sets may not directly influence craniofacial shape, but have indirect effect through phenomena like allometry. Virtually any mutation will have differential effects across cell types, and to the extent that the relevant cell types are present in the developing skull and face, they will contribute to craniofacial variation.

27. It is highlighted as a main finding that the phenotypic effect of different biological processes correlate with each other and with single mutants, however it is not discussed how much of that correlation is driven by different processes sharing genes. Could you elaborate on that aspect of the analysis? How many genes are shared and whether that correlates with the correlation of vector directions? I wouldn’t be sure about how solid the main finding is unless it is clearly shown that it is not coming from shared genetic effects.

We thank the reviewers for this suggestion. We have added an analysis of the relationship between the number of genes 2 process MGP analyses share and the similarity of the resultant phenotype. This is presented in figure 7.

28. Supp Figure 4 shows several covariates for the data, however the methods only described how ‘generation’ was accounted for. What about sex? Mouse source? I don’t have a way of knowing whether you are using them as covariates in your model, or whether you can actually add covariates to the model, because the model is not described in the methods.

We thank the reviewers for pointing out these missing details. They have been clarified in the methods section.

29. The Results section is dry and doesn’t provide any conclusion for the different sections. This could be improved if the conclusions of the analyses that are currently described in the Discussion section be moved to the relevant Results sections. This will also benefit the discussion that is currently very long and most of it is an actual description of results.

These are great suggestions and we have moved the summarized descriptions of the analyses to the relevant subsections of the results.

30. The information necessary to understand the figures should be in the corresponding figure legend not in the main text (e.g. l 327-329, l 331-333, l346-349). And, all figures showing results of the MGP analysis will benefit from ranking the genes according to variance of the loadings instead of alphabetical order – so it becomes easier to assess which genes are more important.

We have expanded on figure legends to provide more detail. We have also revised the figures to sort the loadings by variance.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Essential revisions:

A lengthy discussion was had among the reviewers about clarifying what is meant by new hypothesis. It was agreed that this issue must be addressed, but not does not necessarily need new data or experiments. Two of the three points below directly center on this issue.

1) In line 327-329 (of the tracked changes version) the authors write" "The majority of process MGP analyses demonstrate a similar importance to many alleles, highlighting the main strength process-level analyses over individual marker tests". This is an interesting claim, but what data supports this statement exactly? It would be great if this is supported by some density plot with for example # of genes ranking high (whatever threshold you choose) for each process. The authors have tested ~30 or so processes for some of the analyses. Please clarify if this claim coming from there?

We have included a supplemental figure to figure 6 to better support this claim.

2) The reviews all agreed that there was some level of circularity in the examples shown by the authors. The reviewers strongly felt that the response to comment #10 was not convincing, and the main text was not sufficiently modified to clarify this point. The reviewers agreed that if they remained confused on this point, readers of this manuscript likely would be too. The following is a second attempt at explaining the reviewers point again: all three mutants tested (Bmpr1, Fgf10, and Ankrd11) have previous published data showing their effects on craniofacial shape, sometimes even directly related to the biological process the authors are testing. So, given a) they are functionally annotated to be part of that process, b) there are previous publications – possibly the annotation is derived from those studies, although not necessarily – showing their phenotype effect on craniofacial shape, it does not sound as the authors are generating new hypothesis with their method. This would sound more convincing for the readers if the authors explained better for each example what is their rational for calling each of them a new hypothesis generated by the method. Explicitly list what is previously known or not known about those genes and their effect on craniofacial shape (they reference previous studies, but not in the context of framing how this is different from what they are doing), and be emphatically clear about what is the new "thing" they are trying to show. The reviewers are not saying it is not great to see those genes highly ranked within the process MPG vector, but they are just saying that strong claims of this being a new hypothesis are being made when working with genes for which an effect on craniofacial shape has been shown before. This exact point was also raised in query #16 but there was no new text in the Discussion addressing this point.

We thank the reviewers for this point and agree that the flow from MGP analyses to follow up hypothesis testing could be clearer. The paragraphs pertaining to these analyses have been expanded to more clearly lay out (a) what is known about these genes and their involvement in craniofacial development and (b) how each MGP analysis led to follow up hypothesis testing. This has been addressed in the results, as well as the discussion.

For two of the selected genes (Fgf10 and Bmpr1B) background knowledge might suggest that MGP-informed direction of analysis would add little to their functional understanding. As explained below, much of this knowledge is inferred and more importantly the functional involvement of these genes in the chosen processes was not evident from existing literature. As a matter of fact, it had been overlooked for many years. For the third gene (Ankred11) there were only two reports on craniofacial phenotypes in the literature (one is from Daniel Graf who is an author on this paper). MGP informed us to probe for additional, so far overlooked, phenotypes. The demonstration that MGP is able to lead to the discovery of novel phenotypic involvement of ‘well’-studied genes lends strong support to its general applicability as a discovery platform.

Regarding Bmpr1B, we understand the reviewers’ argument but feel that it is actually based on a misunderstanding. In the case of ‘Bmpr1’, we believe that the reviewer is confusing the gene of focus (Bmpr1B) with Bmpr1A. There are well-described craniofacial phenotypes for Bmpr1A (Saito et al. 2012, Hayano et al. 2015, Pan et al. 2017, Liu et al. 2018, Maruyama et al. 2021). In contrast, Bmpr1b is a ‘minor’ Bmp receptor, its expression levels on a tissue level is about two orders of magnitude lower than Bmpr1a (Pan et al. 2017) and to our knowledge there are no reports showing an independent role for Bmpr1b in craniofacial development. There is one paper that describes delayed cranial bone growth in Bmpr1b mutants (CITE) and the study by Yoon et al. 2005 alludes to a craniofacial phenotype, but this is not described in detail nor is its basis analyzed. Both papers are now discussed and cited in the results. Similarly, here is comparatively little known about its role during chondrogenesis. Initial reports described Bmpr1A and Bmpr1B as having overlapping functions (Yoon et al. 2005). More recently, differential requirements have been reported such as the differential requirement of Bmpr1b for hypertrophic chondrocyte differentiation in vitro (Mang et al. 2020). Given this ‘minor’ role, it was surprising to find a strong association of Bmpr1b in the MGP analysis. It is one thing to note that a mouse has an altered craniofacial shape and another entirely to pinpoint the specific cellular basis for that phenotype. Dr. Graf had been working with this particular mutant for years and had never examined the synchondroses until the MGP analysis produced this result. So, this is a clear case of predicting a phenotype that, upon examination, turns out to be there.

A similar point can be made about the Fgf10 mutant. The prediction we made here is that this mutant has altered craniofacial symmetry. This had not previously been reported and there is no annotation for asymmetric skull associated with this mutant. Again, craniofacial shape effects are diverse and so just the fact that it was known that these mice have a craniofacial phenotype would not predict that they would also have directional asymmetry phenotype. It is known that many mutations are associated with elevated fluctuating (random) asymmetry. Directional asymmetry phenotypes, however, are rare and so to observe this after predicting it from the MGP analysis is a pretty clear example of its utility.

Finally, for the Ankrd11 mutant, we are not predicting a craniofacial shape phenotype per se and we agree with the reviewer that this would have been previously known. Rather, we predicted a palatal shape phenotype in Ankrd11 heterozygous mutants based on the MGP analysis. We do agree that there was some circumstantial evidence that might have pointed in this direction in that human KBG syndrome patients do sometimes (20 percent of cases) have palatal abnormalities. Previous and our own work did not report or notice this phenotype, however. We were motivated to look for it based on the MGP analysis. Again, this confirms the strength of the MGP approach.

Based on above explanations, we are really at a loss to understand how the reviewer comes to the conclusion that our analysis is circular, especially for the first two examples. It is worth pointing that a very large proportion of mutations have an effect on craniofacial shape. Facial shape dysmorphism, for example, is reported in 40% of all described human syndromes. But craniofacial shape effects are very diverse and can have myriad possible underlying causes. A reported annotation of “craniofacial” does not make our approach circular because in the vast majority of those cases, the developmental basis of that effect is unknown and the specific nature of that craniofacial phenotype is not quantified or described. We strongly disagree with the contention of the reviewer that readers would come to the conclusion that the MGP approach is circular as we have presented this work multiple times and discussed it with many colleagues and this issue has never come up.

3. Related to the point above, there were several comments in by the reviewer previously on how, with the current data, the authors could approach more directly the 'testable new hypotheses' that this method generates, but they were not really implemented. This of course does not diminish the value of the paper, but the opposite would have increased its ability to prove that new hypotheses are indeed generated and proven true (not related to genes obviously related to craniofacial shape). For example, point #26 suggests to-test effects on processes that are not related/or expected to be related to craniofacial shape. The authors reply that is difficult to delineate such processes. The reviewers wondered about considering or at least mentioning 'gonad /sperm / egg development'? Social / mating behavior? Circadian rhythm? Etc.

As our genetic data consist of a SNP array, any analysis of non-craniofacial processes will still include marker information for genes involved in craniofacial development to the extent that genes are used and reused myriad processes throughout development. For example, if we take the 8 genes annotated to the GO term “spermatid nucleus differentiation” and reverse the search to see what other processes the gene set is associated with, we get back a list of 144 GO terms. Several of these terms would have plausible associations with craniofacial variation, like “negative regulation of apoptotic process” and “developmental growth”. With a SNP array, we lack the resolution to delineate the potential differences in the expression and use of a gene like pygo2, that is annotated to both “developmental growth” and “spermatid nucleus differentiation”. This has been made clearer in the discussion via the addition of a new paragraph.

https://doi.org/10.7554/eLife.68623.sa2

Article and author information

Author details

  1. Jose D Aponte

    Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Calgary, Canada
    Contribution
    Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1608-8612

  2. David C Katz

    Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Calgary, Canada
    Contribution
    Conceptualization, Investigation, Methodology, Software, Supervision, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared

  3. Daniela M Roth

    School of Dentistry, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, Canada
    Contribution
    Conceptualization, Investigation, Methodology, Supervision, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8156-4681

  4. Marta Vidal-García

    Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Calgary, Canada
    Contribution
    Investigation, Visualization, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7617-7329

  5. Wei Liu

    Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Calgary, Canada
    Contribution
    Investigation, Visualization, Writing – review and editing
    Competing interests
    No competing interests declared

  6. Fernando Andrade

    Department of Biology, Loyola University Chicago, Chicago, United States
    Contribution
    Formal analysis, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  7. Charles C Roseman

    Department of Biology, Loyola University Chicago, Chicago, United States
    Contribution
    Conceptualization, Formal analysis, Methodology, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  8. Steven A Murray

    The Jackson Laboratory, Bar Harbor, United States
    Contribution
    Conceptualization, Methodology, Project administration, Resources, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  9. James Cheverud

    Department of Biology, Loyola University Chicago, Chicago, United States
    Contribution
    Methodology, Project administration, Resources, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  10. Daniel Graf

    1. School of Dentistry, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, Canada
    2. Department of Medical Genetics, Faculty of Medicine and Dentistry, University of Alberta, Edmonton, Canada
    Contribution
    Conceptualization, Investigation, Methodology, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  11. Ralph S Marcucio

    Department of Orthopaedic Surgery, School of Medicine, University of California, San Francisco, San Francisco, United States
    Contribution
    Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – review and editing
    For correspondence
    ralph.marcucio@ucsf.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0537-818X

  12. Benedikt Hallgrímsson

    1. Department of Cell Biology & Anatomy, Alberta Children’s Hospital Research Institute and McCaig Bone and Joint Institute, Cumming School of Medicine, University of Calgary, Calgary, Canada
    2. Department of Animal Biology, University of Illinois Urbana Champaign, Urbana, United States
    Contribution
    Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    bhallgri@ucalgary.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7192-9103

Funding

National Institutes of Health (2R01DE019638)

  • Benedikt Hallgrimsson

Natural Sciences and Engineering Research Council of Canada (238992-17)

  • Benedikt Hallgrimsson

Natural Sciences and Engineering Research Council of Canada (RGPIN-2014-06311)

  • Benedikt Hallgrimsson

Canadian Institutes of Health Research (159920)

  • Benedikt Hallgrimsson

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

Acknowledgements

Grants: NIH-2R01DE019638 to RM, BH, and JC, NSERC 238992-17, CIHR Foundation grant 159920 to BH, CFI grant #36262 to BH and NSERC RGPIN-2014-06311 to DG.

MVG was supported by an Alberta Innovates Postdoctoral Fellowship in Health Innovation.

JDA is supported by an Eyes High fellowship, an Alberta Children’s Hospital Research Institute scholarship, and a MITACS graduate fellowship.

Ethics

The work was performed according to protocols approved and reviewed by animal care committees at the University of Calgary (AC13-0268) and the University of Alberta (AUP1149).

Senior Editor

  1. George H Perry, Pennsylvania State University, United States

Reviewing Editor

  1. Cheryl Ackert-Bicknell, University of Colorado, United States

Reviewers

  1. Laura Saba
  2. Peter Claes, Medical Imaging Research Center, Belgium

Version history

  1. Preprint posted: November 12, 2020 (view preprint)
  2. Received: March 21, 2021
  3. Accepted: November 12, 2021
  4. Accepted Manuscript published: November 15, 2021 (version 1)
  5. Version of Record published: November 30, 2021 (version 2)

Copyright

© 2021, Aponte et al.

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

Metrics

  • 1,399
    Page views
  • 144
    Downloads
  • 5
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jose D Aponte
  2. David C Katz
  3. Daniela M Roth
  4. Marta Vidal-García
  5. Wei Liu
  6. Fernando Andrade
  7. Charles C Roseman
  8. Steven A Murray
  9. James Cheverud
  10. Daniel Graf
  11. Ralph S Marcucio
  12. Benedikt Hallgrímsson
(2021)
Relating multivariate shapes to genescapes using phenotype-biological process associations for craniofacial shape
eLife 10:e68623.
https://doi.org/10.7554/eLife.68623

Categories and tags

Research organism

Further reading

    1. Evolutionary Biology
    2. Genetics and Genomics

    The genomic footprint of social stratification in admixing American populations

    Alex Mas Sandoval, Sara Mathieson, Matteo Fumagalli
    Research Article

    Cultural and socioeconomic differences stratify human societies and shape their genetic structure beyond the sole effect of geography. Despite mating being limited by sociocultural stratification, most demographic models in population genetics often assume random mating. Taking advantage of the correlation between sociocultural stratification and the proportion of genetic ancestry in admixed populations, we sought to infer the former process in the Americas. To this aim, we define a mating model where the individual proportions of the genome inherited from Native American, European and sub-Saharan African ancestral populations constrain the mating probabilities through ancestry-related assortative mating and sex bias parameters. We simulate a wide range of admixture scenarios under this model. Then, we train a deep neural network and retrieve good performance in predicting mating parameters from genomic data. Our results show how population stratification shaped by socially constructed racial and gender hierarchies have constrained the admixture processes in the Americas since the European colonisation and the subsequent Atlantic slave trade.

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics

    Chromosome-specific maturation of the epigenome in the Drosophila male germline

    James T Anderson, Steven Henikoff, Kami Ahmad
    Research Article

    Spermatogenesis in the Drosophila male germline proceeds through a unique transcriptional program controlled both by germline-specific transcription factors and by testis-specific versions of core transcriptional machinery. This program includes the activation of genes on the heterochromatic Y chromosome, and reduced transcription from the X chromosome, but how expression from these sex chromosomes is regulated has not been defined. To resolve this, we profiled active chromatin features in the testes from wildtype and meiotic arrest mutants and integrate this with single-cell gene expression data from the Fly Cell Atlas. These data assign the timing of promoter activation for genes with germline-enriched expression throughout spermatogenesis, and general alterations of promoter regulation in germline cells. By profiling both active RNA polymerase II and histone modifications in isolated spermatocytes, we detail widespread patterns associated with regulation of the sex chromosomes. Our results demonstrate that the X chromosome is not enriched for silencing histone modifications, implying that sex chromosome inactivation does not occur in the Drosophila male germline. Instead, a lack of dosage compensation in spermatocytes accounts for the reduced expression from this chromosome. Finally, profiling uncovers dramatic ubiquitinylation of histone H2A and lysine-16 acetylation of histone H4 across the Y chromosome in spermatocytes that may contribute to the activation of this heterochromatic chromosome.

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics

    SUMOylation of Bonus, the Drosophila homolog of Transcription Intermediary Factor 1, safeguards germline identity by recruiting repressive chromatin complexes to silence tissue-specific genes

    Baira Godneeva, Maria Ninova ... Alexei Aravin
    Research Article

    The conserved family of Transcription Intermediary Factors (TIF1) proteins consists of key transcriptional regulators that control transcription of target genes by modulating chromatin state. Unlike mammals that have four TIF1 members, Drosophila only encodes one member of the family, Bonus. Bonus has been implicated in embryonic development and organogenesis and shown to regulate several signaling pathways, however, its targets and mechanism of action remained poorly understood. We found that knockdown of Bonus in early oogenesis results in severe defects in ovarian development and in ectopic expression of genes that are normally repressed in the germline, demonstrating its essential function in the ovary. Recruitment of Bonus to chromatin leads to silencing associated with accumulation of the repressive H3K9me3 mark. We show that Bonus associates with the histone methyltransferase SetDB1 and the chromatin remodeler NuRD and depletion of either component releases Bonus-induced repression. We further established that Bonus is SUMOylated at a single site at its N-terminus that is conserved among insects and this modification is indispensable for Bonus’s repressive activity. SUMOylation influences Bonus’s subnuclear localization, its association with chromatin and interaction with SetDB1. Finally, we showed that Bonus SUMOylation is mediated by the SUMO E3-ligase Su(var)2–10, revealing that although SUMOylation of TIF1 proteins is conserved between insects and mammals, both the mechanism and specific site of modification is different in the two taxa. Together, our work identified Bonus as a regulator of tissue-specific gene expression and revealed the importance of SUMOylation as a regulator of complex formation in the context of transcriptional repression.