Associations Between Genetic Data and Quantitative Assessment of Normal Facial Asymmetry

Human facial asymmetry is due to a complex interaction of genetic and environmental factors. To identify genetic influences on facial asymmetry, we developed a method for automated scoring that summarizes local morphology features and their spatial distribution. A genome-wide association study using asymmetry scores from two local symmetry features was conducted and significant genetic associations were identified for one asymmetry feature, including genes thought to play a role in craniofacial disorders and development: NFATC1, SOX5, NBAS, and TCF7L1. These results provide evidence that normal variation in facial asymmetry may be impacted by common genetic variants and further motivate the development of automated summaries of complex phenotypes.


INTRODUCTION
The ability to make connections between genetic and phenotypic variation, hinges on phenotypic descriptions that are sufficiently detailed to capture the traits of interest. Biomedical imaging creates very high dimensional datasets that can be analyzed and used to extract phenotype descriptions. Traditional phenotyping from images consists of 2D and 3D measurements of landmarks manually placed on the image. Landmark data is typically sparse and is likely insufficient to capture the complexity necessary for an association with genetic data. A recent study, testing the relationship of facial asymmetry, estimated from nine mid-facial landmarks, with genetic variation at 102 single nucleotide polymorphism (SNP) loci, recently associated with facial shape variation, was unable to identify any SNP relating to asymmetry (Windhager et al., 2014). Methods for automatically phenotyping images and incorporating complex shape information, will be key in understanding the genetic basis of morphology. New approaches such as the BRIM method, developed by Claes et al., have shown the promise of summarizing morphological differences in novel ways to identify genes affecting normal morphology (Claes et al., 2014). The aim of this study is to use automated phenotyping to produce a score of facial asymmetry that incorporates local morphological measurements and their spatial distribution to investigate the genetic basis of facial asymmetry.
Previous analysis of symmetry in 3D facial images has included manual landmarks (Devlin et al., 2007;Stauber et al., 2008), automated measurements (Mercan et al., , 2018, plane of symmetry calculation (Linden et al., 2017), and dense surface registration of a 3D image with a mirrored version (Yu et al., 2009;Demant et al., 2010;Darvann et al., 2011;Djordjevic et al., 2012). Surface registration-based methods show particular promise due to their independence from the plane of symmetry and ability to provide dense shape information across the surface of the face.
Recent applications of surface-registration based methods have been validated by comparison to traditional landmark methods and used quantified asymmetry in individuals using the average transform magnitude or root mean squared error from predefined regions (Claes et al., 2011;Kornreich et al., 2016;Öwall et al., 2016;Verhoeven et al., 2016) and principle modes of variation (Lanche et al., 2007).
In previous work, our group have developed voxel-based deformable morphology analysis methods capable of quantifying facial development in embryos and postnatal animals, from 3D imaging modalities with high precision. Using compact feature representation of image differences, facilitates the comparisons between individuals and across groups (Rolfe et al., 2011(Rolfe et al., , 2013(Rolfe et al., , 2014. In this work we introduce a surface-registration based method to quantify bilateral symmetry in individuals and a metric to summarize how an individual's facial asymmetry and its spatial distribution compares to asymmetry in a healthy, control population.
In this study, we preform GWA analysis on two facial asymmetry scores using a sample of 3186 healthy subjects. Highly significant genetic associations were identified for one of our scores, including genes known to play a role in craniofacial disorders: NFATC1, SOX5, NBAS or likely to play a role in craniofacial development: TCF7L1.

Data
The datasets used in this work were previously collected as part of the FaceBase Consortium's 3D Facial Norms Dataset, described in detail by Weinberg et al. (2016). This study was one of the purposes, under informed consent, and IRB approval was obtained for their use in this work. The dataset consisted of 3D photographic facial surface scans and genetic data from 3186 3D facial meshes from healthy subjects of European Caucasian ancestry, between the ages of 3 and 40 years old. Error screening and quality control measures were followed to reduce variability, due to factors such as facial expression and poor image quality. Subjects were screened for many confounding environmental factors, including: (1) a personal history of facial trauma; (2) a personal history of facial reconstructive or plastic surgery; (3) a personal history of orthognathic/jaw surgery or jaw advancement; (4) a personal history of any facial prosthetics or implants; (5) a personal history of any palsy, stroke, or neurologic condition affecting the face; (6) a personal or family history of any facial anomaly or birth defect; and/or (7) a personal or family history of any syndrome or congenital condition known to affect the head and/or face (Weinberg et al., 2016). To demonstrate that the age range in this dataset did not disrupt the results, we also ran a GWAS excluding pre-pubertal individuals (under 14). The genes identified as significant on the whole dataset still met our threshold for significance on the restricted dataset. These results are reported in Figure S1 and Table S1.
All image data used in this project were acquired using the 3dMD imaging systems (3dMD, Atlanta, GA). These commercial stereo-photography systems incorporate multiple camera viewpoints to provide a 3D mesh of the human face, at no risk to the subject, with the high level of anatomical integrity required for medical research. Several recent studies have assessed the amount of noise or variability that may be present in 3D meshes acquired using the 3dMD system, compared to alternative methods such as direct anthropometry and digital photogrammetry (Dindaroglu et al., 2016), or highaccuracy industrial "line-laser" scanning (Zhao et al., 2017). The findings from these studies suggest that the number of errors likely to be present in a 3dMD dataset is similar to, or an improvement over more traditional methods. The facial surface scans were stored as 3D meshes that were not aligned and may contain extraneous objects such as hair and clothing. Prior to analysis, images were preprocessed to remove noise, cropped to extract the facial region and aligned using custom software developed by our research group (Wu et al., 2014).
A standard set of 24 facial landmarks was collected for each 3D facial mesh. In this study, a subset of 18 landmarks was selected to minimize the number of subjects excluded, due to missing the landmark points. A diagram of the landmarks used for analysis is shown in Figure S2. Details on the procedures used to identify the landmarks on 3D facial surfaces can be found on the "Technical Notes" section of the 3DFN website (https://www. facebase.org/facial_norms/notes).
The genotype data consists of 964,193 SNPs on the Illumina (San Diego, CA) OmniExpress+Exome v1.2 array plus 4,322 custom SNPs chosen in regions of interest based on previous studies of the genetics of facial variation.

Facial Asymmetry Score
Most attempts to summarize image characteristics rely on global features that describe the image as a whole, or local features calculated point-wise across the image. Previous work evaluating asymmetry in facial images has tended toward a local, point-wise approach (Claes et al., 2011;Kornreich et al., 2016;Öwall et al., 2016). While these features have been shown to be effective, we propose a method to produce a richer phenotype description by scoring an individual's relationship to a model of normal asymmetry using both global and local differences. In this work, the global assessment of facial asymmetry is restricted to the region below the eyes (defined by the right and left endocathion landmarks). This restriction limits noise caused by eyelashes, eyebrows, and the hairline and is consistent with landmark-based analysis for this data set, as landmarks were not collected in the forehead region (Weinberg et al., 2016).
In our score assignment system, two local asymmetry metrics were defined to produce independent scores of asymmetry. The local metrics were assessed at each point on the surface of each image. For each asymmetry metric, a statistical model of asymmetry was calculated and each image was scored by its distance from the average model using our novel similarity measure to combine global distribution information with local point-wise correspondences, to produce a summary of the local and global differences. The block diagram of the asymmetry score assignment system is shown in Figure 1.

Local Asymmetry Metrics
The 18 manually-placed landmarks shown in Figure S2 were used to align each subject mesh in a common orientation. After alignment, a base mesh was chosen and corresponding points in each source mesh were found for each point in the base mesh, using the dense point correspondence method developed by Hutton et al. (2003). The locations of corresponding points for each point in the base mesh were averaged over the group to generate the average mesh. Each subject mesh was mirrored across the mid-line and the original and mirrored image were densely mapped to the average mesh. For each point on the average image, an asymmetry flow vector was defined by the difference in position between the corresponding points on the mirror and original images, representing the transformation due to asymmetry. This is illustrated in Figure 2.
We defined two properties of local morphology calculated at each point on an individual facial mesh to capture independent aspects of facial asymmetry.
1. Angle of surface orientation: angle between the normal vectors at corresponding points on and the mirror image. This value quantifies the asymmetry in surface orientation at each point on the image. 2. Angle of deformation: angle between asymmetry flow vector and surface normal on the original image. This value quantifies the direction of the transformation between an image and its mirrored copy at each corresponding point.
These local asymmetry features are illustrated in Figure 3. These angle-based features capture one aspect of asymmetry and are independent of the magnitude of asymmetry.
The magnitude of the deformation can also be used as a local feature of asymmetry using our method. It is defined at each FIGURE 2 | Example of corresponding point mapped from an average image (C) to subject mesh (A) and mirrored copy (B). The asymmetry flow vector is defined between corresponding points on the subject mesh and its mirrored copy. point on an individual facial mesh as the length of the 3D vector between that point and the corresponding point on the mirror image. Results from this approach are included in Figure S3 and Table S5.

Average Model of Normal Asymmetry
Some asymmetry is expected in normal human facial features and the type and amount expected varies with location on the face. For example, asymmetry in the corners of the lips and eyes is more common than asymmetry in the nasal tip. To take into account these spatial differences, each asymmetry score was based on the distance between an individual and an average model of normal asymmetry rather than the absolute asymmetry of the face.
The asymmetry heat maps were used to create an average model of normal asymmetry for each feature. For every point on the average mesh, the average and standard deviation of each feature distribution over all corresponding points in the dataset were calculated. The average and standard deviation heat maps for the angle of surface orientation feature are shown in Figure 4.

Distance From Average Model of Normal Asymmetry
To assess the similarity between two feature heat maps the following questions must be addressed: 1. What feature values are present in the image? 2. Where are regions of similar feature values approximately located?
To simultaneously address these two questions, we developed a similarity metric that combines information about the global feature distribution and point-wise differences. Histograms of image features provide a robust description of global image data that has proven to be powerful in detecting similarity. However, the use of histogram representations of features presents two primary drawbacks: the loss of spatial distribution information and the loss of information due to quantization. To address this, histograms can be augmented by the inclusion of additional spatial information and other local properties (Birchfield and Rangarajan, 2005;Lyons, 2009;Prabhu and Kumar, 2014;Zeng et al., 2015). In previous work, our group developed a method to simultaneously assess similarities in feature values and their regional distribution based on spatial histograms (Rolfe et al., 2014).
Intuitively, the spatial histogram, or spatiogram is an image histogram where the distribution of values is spatially weighted by the similarity of the spatial positions of the values in each bin. Typically, this is done by modeling the spatial location of the contents of each histogram bin with a single or mixture of Gaussian distributions. In this application the known point correspondences between images in the data set calculated in section 2.2.1 are leveraged to provide a more precise score of spatial matching between histogram bins. The spatial information is incorporated as the set of coherent feature regions in a histogram bin. For an image I, the histogram of I is defined as: where n b is the number of pixels with values assigned to the b-th bin and B is the total number of bins. Our spatially augmented histogram is defined as: where n b is the number of points with values assigned to the bin b, and R b is the set of m coherent regions r b1 , ...r bm where r bi is a vector of j point indexes < x 1 , ...x j >. Coherence of regions is determined by computing connected components. A connected region r bi is the set of mesh points such that for any point x, x ′ ∈ r bi , there is a path in r bi to from x to x ′ . A threshold for coherence can be set for a connected region with greater than τ mesh points. For this study regions with τ < 20 (corresponding to less than 0.1% of the image) were classified as incoherent. An example of a feature heat map and coherent regions extracted from the histogram is shown in Figure 5. In Figure 5A, feature values from a feature heat map are grouped into histogram bins. Figure 5B shows the original feature heat map and the extraction of coherent image regions assigned to Bin 4 of the histogram in Figure 5A.
The distance metric between two augmented histograms is typically based on the Bhattacharyya distance between histograms, weighted by the spatial similarity of the contents of bin b as in Birchfield and Rangarajan (2005). The difference between spatial histograms h and h ′ is expressed as: The spatial weighting term m b expresses the similarity of the m spatial regions in bin b. Previously, the Mahalanobis distance, or number of standard deviations between the means of the Gaussian distributions in each bin, has been used to weight the spatial similarity. In this work, we utilized the spatial weighting term to incorporate the point-wise similarity between corresponding points from two feature heat maps. This modification addressed both the need for spatial information and the loss of information due to histogram quantization. We defined the spatial weighting term as the mean feature error between histogram regions, normalized by the standard deviation Frontiers in Genetics | www.frontiersin.org at that point, calculated from the average model of asymmetry. This is expressed as: where w bi is the weight of the ith coherent region in bin b, A(x j ) and A ′ (x j ) are the feature values from the two feature maps at the corresponding point j and σ¯j is the standard deviation at point j. This spatial weighting term represents average error between feature maps, measured in standard deviations, for each coherent region. To achieve a symmetric distance measure, the total distance between h and h ′ was defined as: This distance ρ(h, h ′ ) was applied to assess the similarity of each individual feature heat map to the average feature heat map. This provided a hybrid local-plus-global summary of the abnormality of asymmetry of an individual and was assigned as our score of asymmetry. Average feature heat maps for subjects with the lowest (lower 10 percent of the data set) and highest (upper 10 percent of the data set) asymmetry scores for the angle of surface orientation feature are shown in Figure 6. In the average heat map from the high asymmetry score group in Figure 6B, regions with high values contributed the most to the score in individuals with high levels of asymmetry.

Genetic Association Analyses
Whole genome association with each phenotype score was done using PLINK (Purcell et al., 2007). SNPs with the minor allele present in less than 5 subjects were removed, resulting in 747,780 remaining SNPs. The first four principle components of the genetic data were used as covariates to adjust for the effects of ancestry. A linear model was used to test genetic association between our phenotype scores and each SNP, controlling for the effects of age and gender. The Benjamini-Hochberg procedure was used to adjust the original p-values globally over both phenotype scores in order to control the false discovery rate FIGURE 6 | Group average heat maps from subjects with low asymmetry (angle of surface orientation score in lower 10 percent of data set) (A) and subjects with high asymmetry (angle of surface orientation score in upper 10 percent of data set) (B). Regions with low asymmetry are blue and regions with high asymmetry are red.
(FDR) (Benjamini and Hochberg, 1995). Genome-wide Complex Trait Analysis (GCTA) was used to estimate the proportion of variance in each phenotype score explained by all GWAS SNPs, i.e., heritability (Yang et al., 2011). Each phenotype score was tested for associations with age and sex. The Pearson correlation coefficient was used to test for an association with age. An association with sex was tested using the Kendall rank correlation coefficient (tau). The Kendall test does not rely on the assumption of normally distributed data and so is more appropriate for dichotomous data such as sex. The correlations found between the asymmetry scores, age and sex were weak, though the correlations had high levels of significance in terms of the p-values, as reported in the Tables S2, S3. We speculate that this effect is likely due to the large sample size (i.e., statistical power), which made it possible to detect the significant associations.

Angle of Surface Orientation Score
The top 10 SNPs significantly associated with the angle of surface orientation scores (p-value < 5 × 10 −8 ) are listed in Table 1, and the Manhattan plot is shown in Figure 7. Of these SNPs with highly significant associations, three are located on genes with known links to craniofacial abnormality and asymmetry (NFATC1, SOX5, and NBAS) and one (SNX6) is on a gene with a potential link. NFATC1 encodes a transcription factor that plays a role in mandibular development and the Wnt signaling pathway, which is instrumental to facial morphogenesis (Winslow et al., 2006;Brugmann et al., 2007;Doraczynska-Kowalik et al., 2017). Mutations in NFATC1 are linked to Cherubism, a disorder characterized by abnormal bone tissue in the lower part of the face and a characteristic facial phenotype (Kadlub et al., 2016). A recent GWA study of morphological measurements has also suggested a possible link between this gene and measurements of the mouth (Lee et al., 2017). SOX5 encodes a transcription factor involved in the regulation of embryonic development that is thought to play a role in chondrogenesis. SOX5 is linked to Lamb-Shaffer Syndrome, which can cause an abnormal craniofacial phenotype including a facial asymmetry, depressed and/or broad nasal bridge, and bulbous nasal tip (Lamb et al., 2012). Mutations in NBAS are associated with Pelger-Huet Anomaly, which has a phenotype including facial asymmetry, long face, and straight nose (Segarra et al., 2015). It is also linked to Feingold Syndrome 1, which can result in craniofacial dismorphology including asymmetry, triangular shaped face, and flat nasal tip (Chen et al., 2012). SNX6, a member sorting nexin family, has not been definitively linked to craniofacial disorders, however multiple studies have suggested it as a candidate gene for holoprosencephaly, the most common developmental field defect in patterning of the human prosencephalon and associated craniofacial structures (Kamnasaran et al., 2005;Segawa et al., 2007). Also of interest is TCF7L1, which encodes for a transcription factor that mediates the Wnt signaling pathway and has been found to have high expression in the developing murine palate (Potter and Potter, 2015). Genes with known or potential association with craniofacial abnormality and asymmetry are boldfaced.
The angle of surface orientation phenotype scores were assessed for heritability using GCTA and were found to have a proportion of variance consistent with a substantial heritability. Detailed results are reported in Table S4.

Angle of Deformation Score
The angle of deformation phenotype scores showed less significance than the facial asymmetry scores. The top 10 SNPs associated with the angle of deformation scores are reported in Table 2, and the Manhattan plot is shown in Figure 8. While many of the p-values for the SNPs associated with this phenotype score are not considered significant, it is possible that the multiple testing correction might have been overly conservative when significant linkage equilibrium was present. As there are a number of genes which have known or potential links to facial development or morphology, we reported the genes associated with these SNPs of interest, though the associations are weak.
AMBRA1 encodes a protein that regulates different steps of the autophagic process and is an important regulator of embryonic development. Its mutation or inactivation in mice was shown to result in embryonic malformations (Fimia et al., 2007). Rare deletions in NRXN3 was linked to autism spectrum disorder (Vaags et al., 2012). While there is as yet no consensus on facial phenotypes associated with autism spectrum conditions (ASC), there is evidence to suggest that there are morphologically distinct subgroups within ASC that correspond with different cognitive and behavioral symptomatology (Boutrus et al., 2017). Two SNPs of interest are located on the gene FANCC, which encodes a DNA repair protein with a role in the maintenance of normal chromosome stability. FANCC is implicated in Gorlin syndrome that has a phenotype including broad nasal root, cleft lip, and cleft palate (Reichert et al., 2015). FANCC is also linked to Fanconi anemia that has a phenotype including craniosynostosis, microencephaly and small eyes (de Winter et al., 2000). FTO is a protein coding gene associated with growth retardation, developmental delay, and facial dysmorphism (Boissel et al., 2009;Daoud et al., 2015). The associated phenotype includes skull asymmetry, coarse facial features, abnormal positioning of the maxilla or mandible, prominent alveolar ridge, and cleft palate. The retinoid acid receptor-responsive gene RARRES1 contains two SNPs of interest. This gene is thought likely to play a role in embryonic morphogenesis (Oldridge et al., 2013).
The angle of deformation phenotype scores were assessed for heritiability using GCTA and were found to have a proportion of variance suggesting minimal heritiability and a p-value suggesting low significance. This is a possible explanation for the low levels of significance observed. These results are detailed in Figure S4.

Comparison to Asymmetry Scores Based on the Deformation Vector Magnitude
Angle-based measurements capture one aspect of asymmetry, which may be relevant to specific biological processes. Deformation magnitude, defined as the magnitude of the distance between each point on a facial image and its corresponding point on a mirrored image, is another common choice for meshbased shape analysis. For comparison, we implemented our asymmetry score using the deformation magnitude as the local asymmetry feature. This local property was then used to calculate an overall score of asymmetry following the procedure outlined in the Methods section 2.2. The GWAS results from our magnitude-based asymmetry score are reported in Figure S3 and Table S5.
Since the average value of the deformation magnitude over an image surface is a metric frequently used in other studies, we also implemented an established measure from the literature to compare to our deformation magnitude asymmetry scores (Verhoeven et al., 2016). In this work, local asymmetry is defined as the magnitude of the distance between each point on a facial image and its corresponding point on a mirrored image. The measure of total facial asymmetry was calculated using the average of these distances over the face. This method was selected because it is similar to those used by several other groups and the results were validated on a data set with known ground truth. The GWAS results from this comparable deformable morphology approach are detailed in Figure S4.
Both magnitude-based methods we tested had lower significance and did not identify genes known to result in facial asymmetry. One gene of interest identified by both methods, MYO10, has been linked to craniofacial development in zebrafish. The genes identified by these two magnitude-only methods overlapped, but our magnitude-based asymmetry score showed higher levels of significance.

Comparison to Asymmetry Scores Based on Landmark Measurements
The motivation for developing the mesh-based methods in this work was to provide more complex phenotypes for genetic association than standard landmark-based approaches. Subtle differences in asymmetry that may be scientifically interesting are  The genes with known or potential association with facial development or morphology are boldfaced.
unlikely to be captured by landmark data, which is usually very sparse. While landmark-based methods may identify associations with genes of interest, they may identify different pathways than mesh-based analysis as they do not use data between the landmark points. To compare our method to GWAS using a traditional, landmark-based approach of measuring asymmetry, a score of facial asymmetry was defined using the Procrustes distance between an image and its mirrored copy (Bookstein, 1997). Each image was rigidly aligned with its mirrored copy. A subset of the 12 bilaterally paired landmarks was selected from the original 18 landmarks shown in Figure S2 and the Euclidean distance between right/left landmark pairs was measured. The facial asymmetry score was calculated using the average of these distances. Using this method, no SNPs were found to meet the threshold of genome-wide significance of p = 5 × 10 −8 , as detailed in Figure S5. These results provide additional motivation for the use of mesh-based analysis, in addition to the improvements in precision and reproducibility,

DISCUSSION
Asymmetry is the topic of a large number of studies investigating how genetic and environmental factors influence normal development. It is likely to be influenced by complex and interelated factors, which can be difficult to control for in human studies, presenting significant challenges for analysis. Asymmetry, especially fluctuating asymmetry, has been hypothesized to be closely linked to developmental instability and many studies have interpreted it as a marker of environmental stress during development (Klingenberg and McIntyre, 1998;DeLeon, 2007;Ozener, 2010). However, several recent studies have called these findings into question and have suggested a stronger role of heredity (Quinto-Sánchez et al., 2015. Further studies using genotype and phenotype data are needed to better understand how the developmental processes leading to asymmetry are impacted by environmental factors. While subjects in our study were screened for a number of possible environmental influences on facial asymmetry, as detailed in section 2.1, many other potential confounding factors remain, such as twinning status and smoking behavior that are unknown or could not feasibly be controlled for in this study. Despite these limitations, we have applied a data-driven approach to evaluate methods for quantifying aspects of asymmetry that may be related to biological processes resulting in facial asymmetry. Consistent with other recent findings on the genetic basis of normal facial variation, several of the genes associated with variation in normal asymmetry are involved in syndromes with craniofacial phenotypes. This supports the hypothesis that common variants near the genes related to Mendelian syndromes are implicated in normal phenotypic variation (Shaffer et al., 2016). While we are cautious about interpreting the results from a angle of deformation asymmetry score, due to the weak associations, several of the genes of interest identified are associated with embryonic morphology and development and craniofacial abnormality. The genes identified by the angle of deformation score do not overlap with the genes identified by the angle of asymmetry score. This indicates the possibility that the two aspects of asymmetry quantified may be useful for identifying different biological pathways impacting facial asymmetry.
The heat maps of local asymmetry features provide information about the regions of the face that contribute most to the asymmetry scores. Figure 6B shows the average feature heat map for subjects with the highest asymmetry scores (top 10 percent of the data set). This heat map shows higher levels of asymmetry than the average feature heat map in Figure 4 and also suggests the relative importance of the nasal tip, nasal bridge, upper lip, and chin regions in subjects with high levels of normal asymmetry.
Questions still remain about the ability of complex phenotypes to be accurately associated with genetic data as the genotypephenotype map for facial morphology is likely to be incredibly complex (Hallgrimsson et al., 2014). A single gene can result in local or global shape differences and be intertwined with environmental factors. Despite these challenges, we have demonstrated that our hybrid local-to-global score of abnormal asymmetry was able to find associations with genes known to play a role in craniofacial morphology and asymmetry. While we do not have an assurance that our automated phenotyping method is the optimal strategy to summarize phenotypes for genetic association, the significance of the results motivates its further development. One limitation of this study was our lack of a comparable dataset with which to replicate our findings. If one becomes available in the future, applying these methods to identify an overlapping set of genes would significantly strengthen the findings in this work.
In future work, new local morphological metrics can be investigated using this framework. This method can also be implemented to compare subjects to an average model of a group of interest, rather than a control population, to assess similarity to a known phenotype. Taking a data-driven approach to optimizing phenotypic descriptors, guided by the significance of the genetic associations uncovered, will contribute to both our understanding of the genetic basis of human facial variation and the creation of new metrics for biologically relevant phenotype data.

ETHICS STATEMENT
The study was carried out in accordance with the recommendations of University of Washington IRB #42874 with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The University of Washington IRB approved the protocol.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study were obtained from FaceBase (www.facebase.org), and were generated by projects U01DE020078 and U01DE020054. The FaceBase Data Management Hub (U01DE020057) and the FaceBase Consortium are funded by the National Institute of Dental and Craniofacial Research. The phenotype scores developed for this work will be made available on request.

AUTHOR CONTRIBUTIONS
SR carried out the main efforts on the research including developing the theory behind the angle of surface orientation score and the angle of deformation score as well also carrying out all the experiments, and writing the paper. S-IL acted as a consultant on the analysis of the GWAS results and helped to write the Results and Discussion sections of the paper. LS served as the primary adviser to SR in this work.

FUNDING
Research reported in this publication was supported by the National Institute of Dental and Craniofacial Research (NIDCR) of the National Institutes of Health under award numbers: 5F32DE025519 and U01-DE020050. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.