Abstract Numerous lines of evidence point to a genetic basis for facial morphology in humans, yet little is known about how specific genetic variants relate to the phenotypic expression of many common facial features. We conducted genome-wide association meta-analyses of 20 quantitative facial measurements derived from the 3D surface images of 3118 healthy individuals of European ancestry belonging to two US cohorts. Analyses were performed on just under one million genotyped SNPs (Illumina OmniExpress+Exome v1.2 array) imputed to the 1000 Genomes reference panel (Phase 3). We observed genome-wide significant associations (p < 5 x 10−8) for cranial base width at 14q21.1 and 20q12, intercanthal width at 1p13.3 and Xq13.2, nasal width at 20p11.22, nasal ala length at 14q11.2, and upper facial depth at 11q22.1. Several genes in the associated regions are known to play roles in craniofacial development or in syndromes affecting the face: MAFB, PAX9, MIPOL1, ALX3, HDAC8, and PAX1. We also tested genotype-phenotype associations reported in two previous genome-wide studies and found evidence of replication for nasal ala length and SNPs in CACNA2D3 and PRDM16. These results provide further evidence that common variants in regions harboring genes of known craniofacial function contribute to normal variation in human facial features. Improved understanding of the genes associated with facial morphology in healthy individuals can provide insights into the pathways and mechanisms controlling normal and abnormal facial morphogenesis.

Author Summary There is a great deal of evidence that genes influence facial appearance. This is perhaps most apparent when we look at our own families, since we are more likely to share facial features in common with our close relatives than with unrelated individuals. Nevertheless, little is known about how variation in specific regions of the genome relates to the kinds of distinguishing facial characteristics that give us our unique identities, e.g., the size and shape of our nose or how far apart our eyes are spaced. In this paper, we investigate this question by examining the association between genetic variants across the whole genome and a set of measurements designed to capture key aspects of facial form. We found evidence of genetic associations involving measures of eye, nose, and facial breadth. In several cases, implicated regions contained genes known to play roles in embryonic face formation or in syndromes in which the face is affected. Our ability to connect specific genetic variants to ubiquitous facial traits can inform our understanding of normal and abnormal craniofacial development, provide potential predictive models of evolutionary changes in human facial features, and improve our ability to create forensic facial reconstructions from DNA.

Citation: Shaffer JR, Orlova E, Lee MK, Leslie EJ, Raffensperger ZD, Heike CL, et al. (2016) Genome-Wide Association Study Reveals Multiple Loci Influencing Normal Human Facial Morphology. PLoS Genet 12(8): e1006149. https://doi.org/10.1371/journal.pgen.1006149 Editor: Gregory S. Barsh, Stanford University School of Medicine, UNITED STATES Received: March 5, 2016; Accepted: June 8, 2016; Published: August 25, 2016 Copyright: © 2016 Shaffer et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All of the phenotypic measures and genotypic markers used here are available to the research community through the dbGaP controlled access repository (http://www.ncbi.nlm.nih.gov/gap) at accession number: phs000949.v1.p1. The raw source data for the phenotypes – the 3D facial surface models – are available for the 3D Facial Norms dataset through the FaceBase Consortium (www.facebase.org). Finally, searchable results datasets of the p-values from the studies reported here are available through the FaceBase Human Genomics Analysis Interface (facebase.sdmgenetics.pitt.edu). Funding: The National Institute for Dental and Craniofacial Research (http://www.nidcr.nih.gov/) provided funding through the following grants: U01-DE020078 (SMW, MLM); U01-DE020057 (JCM, MLM); R01-DE016148 (MLM, SMW); U01-DE020054 (RAS); U01-DE024425 (MLM); K99-DE02560 (EJL). Funding for genotyping was provided by the National Human Genome Research Institute (https://www.genome.gov/): X01-HG007821 (MLM, SMW, EF). Funding for initial genomic data cleaning by the University of Washington (CAL and CCL) was provided by contract #HHSN268201200008I from the National Institute for Dental and Craniofacial Research (http://www.nidcr.nih.gov/) awarded to the Center for Inherited Disease Research (CIDR). Additional funding for was provided by the National Institute of Justice (http://www.nij.gov/Pages/welcome.aspx): 2013-DN-BX-K005 (RAS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Numerous lines of converging evidence indicate that variation in facial morphology has a strong genetic basis. These include the results of human heritability studies using twin and parent-offspring designs [1–5], Mendelian craniofacial syndromes [6], transgenic animal models with distinctive craniofacial phenotypes [7–9], and studies mapping QTLs for craniofacial shape in several mammalian models [10–14]. However, we still have little understanding of how genetic variation relates to the diversity of normal facial traits commonly observed in humans. Understanding the genetic basis for normal facial variation has important implications for human health. Many genetic syndromes with dysmorphic facies are characterized by relatively subtle morphological changes, often involving quantitative traits with continuous distributions [6]. The range of variation for any given facial trait often displays substantial overlap between affected and healthy individuals. Thus, understanding the genetic factors that contribute to normal facial trait variation may provide valuable insights into the causes of craniofacial dysmorphology, including common craniofacial birth defects such as orofacial clefts [15,16]. To date, only a few studies have explicitly tested for associations between aspects of normal human facial morphology and common genetic variants. Among these, two genome-wide association (GWA) studies have been carried out on healthy individuals of European ancestry using 3D facial imaging and a combination of traditional and more advanced morphometric methods to derive phenotypes [17,18]. Between these two studies, a handful of intriguing genome-wide significant signals were reported, although they were largely non-overlapping. Notably, both studies reported an association between PAX3 variants and anatomical changes in interorbital region, an intriguing finding given that mutations in PAX3 cause Waardenburg Syndrome type 1 which is characterized by hypertelorism among other morphological abnormalities. Both studies also reported significant associations with measures of nasal projection in their discovery cohorts, although different genomic regions were implicated. In addition, several more focused candidate gene studies of loci implicated in craniofacial syndromes or in developmental pathways involved in craniofacial development have connected one or more craniofacial dimensions or aspects of shape with a small number of common genetic variants [19–28]. At least three candidate gene studies [20,25,28] have reported modest associations between common variants in FGFR1 and normal variation in craniofacial morphology, but in each case a different constellation of traits was involved. It is notable that none of the genes from these studies, including FGFR1, were identified in the two previous GWA studies of facial morphology. Thus, while prior studies have detected a handful of biologically plausible genes associated with variation in craniofacial features, it is clear that these efforts are just scratching the surface and the potential for additional discovery is great. In the current study, we performed GWA analyses on a set of 20 craniofacial measurements commonly used in clinical assessment (Fig 1) derived from 3D surface images in two well-characterized samples of unrelated White individuals of European ancestry from the USA: a sample comprised of 2447 individuals collected through the University of Pittsburgh (i.e., the Pittsburgh sample) and an independent sample of 671 individuals collected under the direction of the University of Colorado (i.e., the Denver sample). All participants were genotyped using the same SNP array (Illumina OmniExpress+Exome v1.2), which included just under one million SNPs, and were imputed to the 1000 Genomes reference panel (Phase 3). We conducted association tests in each sample separately and combined the results using a meta-analytic approach. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Set of 20 linear distance measurements used in the current study. (A) Cranial base width, (B) Upper facial depth*, (C) Middle facial depth*, (D) Lower facial depth*, (E) Morphological facial height, (F) Upper facial height, (G) Lower facial height, (H) intercanthal width, (I) Outercanthal width, (J) Palpebral fissure length*, (K) Nasal width, (L) Subnasal width, (M) Nasal Protrusion, (N) Nasal ala length*, (O) Nasal height, (P) Nasal Bridge Length, (Q) Labial fissure length, (R) Philtrum length, (S) Upper lip height, and (T) Lower lip height. Measurements with an asterisk (*) are bilateral, but only the left side is shown in the figure. https://doi.org/10.1371/journal.pgen.1006149.g001

Discussion Based on meta-analysis, we observed seven associated loci for five facial traits: cranial base width (Fig 1A), intercanthal width (Fig 1H), nasal width (Fig 1K), nasal ala length (Fig 1N), and upper facial depth (Fig 1B). The most significant of these, meeting the strict study-wide threshold for significance (i.e., p < 5 x 10−9), was the association of cranial base width at 20q12 410kb downstream of MAFB, a transcription factor previously implicated in orofacial clefts [29] and facial characteristics in cleft families [30]. However, the MAFB SNP associated with clefting was 250kb away and not in LD with the SNP observed here. In addition to orofacial clefting, mutations in MAFB cause multicentric carpotarsal osteolysis syndrome, which includes mild facial anomalies. These phenotypes are consistent with the developmental role of MAFB in regulating the migration of cranial neural crest cells during the patterning of skeletomuscular features of the head [31]. Altogether, these lines of evidence suggest a possible role for MAFB in normal facial variation. Another association for cranial base width was observed at 14q21.1 in the vicinity of PAX9, SLC25A2, MIPOL1, and FOXA1. The homeodomain protein-coding PAX9 is important for craniofacial development in mice [32,33] and dental development in humans [34]. Using in situ hybridization, Peters et al. [32] showed that Pax9 is expressed throughout the developing cranial base in mice at E13.5. Biological evidence for other genes in this region also suggests possible roles in facial variation including MIPOL1, which has been observed to be affected by chromosomal aberrations in patients with craniofacial phenotypes, including holoprosencephaly [35]. Because holoprosencephaly involves alteration in the horizontal spacing of facial structures, variants in genes associated with this condition may also influence measures of craniofacial width in healthy subjects. Taken together, these previous observations point to the plausibility of genetic variants in this region influencing normal facial variation. Two genetic associations were observed for intercanthal width. One of these was at 1p13.3 within the gene GSTM2, which codes an enzyme involved in detoxification of compounds. Among the genes within 250kb of the peak signal are two potentially relevant candidate genes, GNAI3 and ALX3. Mutations in GNAI3, which encodes a G protein subunit involved in pharyngeal arch patterning, cause auriculocondylar syndrome, a rare craniofacial disorder [36,37], although hyper- or hypotelorism have not specifically been described. ALX3 is a homeobox gene essential for head and face development. Mutations in ALX3 result in frontonasal dysplasia 1 [38] in humans and nasal clefts in mice [39]. Ocular hypertelorism is a prominent feature of frontonasal dysplasia and is believed to result from disruptions of the Hedgehog signaling pathway [40,41]. The second association with intercanthal width was observed for a broad 900kb LD block on Xq13.2. The peak of the diffuse association signal is over PABP1C1L2A, which encodes an uncharacterized poly-A binding protein. However, at the edge of the LD block, roughly 500kb centromeric to the peak signal, is HDAC8, which encodes a histone deacetylase involved in epigenetic gene silencing during craniofacial development [42]. Mutations in HDAC8 cause Cornelia de Lange syndrome [43,44], a developmental disorder characterized by facial dysmorphology including hypertelorism. A mutation in HDAC8 has also been described in a family with an X-linked intellectual disability syndrome with distinctive facial features, which included hypertelorism [45]. A number of other genetic associations with facial traits were observed at loci harboring genes with relevant biological roles. For example, an association with nasal width was observed at 20p11.22 near the PAX1 gene. Mutations in PAX1 cause otofaciocervical syndrome [46], characterized by facial dysmorphology, including specific nasal features such as a sunken nasal root and excessive narrowing. PAX1 plays a role in chondrocyte differentiation [47], which may explain its association with nasal width, a measure of the distance between the left and right cartilaginous nasal alae. Nevertheless, a study of Pax1 expression in mice showed expression in the pharyngeal arches at E11.5, but not in the developing olfactory placodes [48], so it is unclear how this gene influences nasal development. An association with nasal ala length was observed at 14q11.2 in a region containing an RNase gene cluster plus at least 25 other genes (within about 400kb of the association peak). Among the many genes in region are ZNF219, which encodes a transcriptional partner of Sox9 essential for chondrogenesis in mice [49], and CHD8, mutations in which are associated with autism spectrum disorder in conjunction with macrocephaly and distinct facial features including a broad nose [50]. A similar story pertains to the association between SNPs on 11q22.1 and upper facial depth. The peak signal occurs within TRPC6, which encodes a cation channel subunit mutated in hereditary renal disease [51]. TRPC6 has no known connection to craniofacial development, but other genes in the region have reported craniofacial expression, including YAP1 [52]. In aggregate, we observed a number of genetic associations near genes with biologically plausible roles in facial variation. A common theme was that associated loci harbored genes involved in syndromes with craniofacial phenotypes. This result fits with a long-standing hypothesis about the relationship between Mendelian syndromes and complex traits in which common variants near genes causing Mendelian syndromes are involved in related common, complex diseases and traits, including normal phenotypic variation [53]. That being said, for any of the observed associations, it is not clear which variant might be functional, though we hypothesize that the functional variants underlying the statistical signal will be regulatory. Moreover, it is not clear which genes they regulate. Thus, references to implicated genes should always be treated with appropriate caution. While none of our genome-wide or suggestive (p < 5 x 10−7) signals included SNPs implicated in either of the previous two European-focused GWA studies [17,18], we nevertheless found evidence of association when we tested the top SNPs from these studies against comparable phenotypes from our data. Strongest among these was nasal ala length, a lateral projective measure of the nose extending from alar cartilage to the nasal tip, previously associated with 1p36.32 (rs4648379, PRDM16) [18] or 3p14.3 (rs1982862, CACNA2D3 [17]. We found at least nominal associations with both of these regions in our data, with one (rs4648379, PRDM16) showing evidence at p = 1.70 x 10−5. Both prior GWA studies reported an association between SNPs at 2q35 (PAX3) and morphology of the interorbital septum. We tested these SNPs and found an association between rs974448 and intercanthal width (p = 0.002), lending some additional support to the claim that common variants in PAX3 might influence aspects of normal facial morphology. Our ability to test previously reported genetic associations was limited by a lack of directly comparable phenotypes, which is related to differences in data collection methods and the type and number of measurements available. In addition, the prior two European GWA studies each used imaging modalities different from the kind used here. Similar factors may also explain some of the discrepancies in association results observed between our two study cohorts. Although care was taken to generate the same set of distance measures in both cohorts, the different 3D cameras and landmarking protocols used could result in different patterns of association. Despite these differences, the measurements from both cohorts were found to be very similar in their overall distributions. Alternative phenotypes, such as multivariate measures of facial shape, can also be used in these types of studies. However, prior attempts to extract shape variation from landmark and surface data in similarly sized samples (e.g., using Procrustes–based methods) have not yielded statistically significant associations [17,18]. One reason for this may be that the effect of any single gene is diluted because the resulting phenotypes represent such a complicated mix of local and global shape features. The problem is highly complex and there is presently little consensus on the most prudent approach to complex facial phenotyping [54]. Fortunately, several promising approaches are on the horizon, such as the BRIM methods being developed by Claes et al. [28]. It is likely that samples an order of magnitude larger than anything available at the moment will be required before we can begin to exploit the richness contained in human 3D facial datasets. Despite these limitations, we have found evidence of genetic association between chromosomal regions containing genes with important roles in facial development and quantitative traits that characterize key features of the normal human craniofacial complex. In addition to improving our general knowledge of the factors that underlie the diversity of facial forms we see in humans, these genotype-phenotype associations may help us better understand the wide range of phenotypic expression and severity seen in some rare genetic syndromes. Common variants in a number of different genes or regulatory elements may contribute to the expression of dysmorphic phenotypes present in these conditions. Moreover, such associations in healthy individuals may aid the search for clues to the etiology of much more common craniofacial anomalies. For example, three of the traits reported here (cranial base width, nasal width and intercanthal width) have all been previously implicated as potential phenotypic risk factors for orofacial clefting, the most common craniofacial birth defect in humans [55]. Weinberg et al. [15,56] have shown that the unaffected, but genetically at-risk, relatives of cleft-affected individuals exhibit increased breadth through the middle and upper face. The identification of the genes that influence these traits may help us identify important risk-modifiers for clefting [16]. Testing the SNPs implicated here for associations in our cleft families is a future goal of our research group.

Materials and Methods Ethics statement Institutional ethics (IRB) approval was obtained at each recruitment site and all subjects gave their written informed consent prior to participation (University of Pittsburgh Institutional Review Board #PRO09060553 and #RB0405013; UT Health Committee for the Protection of Human Subjects #HSC-DB-09-0508; Seattle Children’s Institutional Review Board #12107; University of Iowa Human Subjects Office/Institutional Review Board #200912764 and #200710721; Colorado Multiple Institutional Review Board #09–0731; UCSF Human Research Protection Program Committee on Human Research #10–00565). Study samples Our study included two independent samples, each comprised of unrelated self-described White individuals of European ancestry from the United States. The Pittsburgh sample included 2447 unrelated individuals ranging in age from three to 49 years. The majority of these participants (n = 2272) were recruited at research centers in Pittsburgh, Seattle, Houston and Iowa City as part of the FaceBase Consortium’s 3D Facial Norms Dataset, described in detail by Weinberg et al. [57]. The remaining subjects were recruited as healthy controls for a separate study at Pittsburgh on orofacial cleft genetics. The Denver sample included 671 unrelated individuals ranging in age from three to 12 years. These participants were recruited from Denver and San Francisco as part of a separate FaceBase Consortium study of facial shape genetics in multiple ethnic populations [58]. The basic demographic features of these samples are provided in S4 Table. In both samples, subjects were excluded if they had a personal history of facial trauma, a personal history of facial reconstructive or plastic surgery, a personal history of orthognathic/jaw surgery or jaw advancement, a personal history of any facial prosthetics or implants, a personal history of any palsy, stroke or neurologic condition affecting the face, a personal or family history of any facial anomaly or birth defect, and/or a personal or family history of any syndrome or congenital condition known to affect the head or face. Phenotyping 3D facial surfaces were captured using digital stereophotogrammetry, a standard imaging method resulting in high-density, geometrically accurate point clouds representing the surface contours of the human body [59]. Facial surfaces in the Pittsburgh sample were collected with 3dMD imaging systems (3dMD, Atlanta, GA). Facial surfaces in the Denver sample were imaged using the Creaform Gemini camera system (Quebec, Canada). A common set of 24 standard facial soft-tissue landmarks [60] was collected on each 3D facial surface and the xyz coordinate locations recorded (S22 Fig). Landmarks were collected manually in the Pittsburgh sample as described in Weinberg et al. [57]. An automated landmark collection method was used in the Denver sample. From these landmarks, we calculated 20 linear distances that correspond to craniofacial measurements commonly used in clinical assessment [61]. These measurements are shown in Fig 1 and listed in S5 Table. For bilateral measurements, values were summed across the left and right sides in order to minimize the number of traits tested. Trait values were inspected for outliers by computing age- and sex-specific z-scores. Genotyping, quality checks, imputation, and population structure For each study sample, genotyping was performed by the Center for Inherited Disease Research (CIDR). Saliva samples were used to genotype 3,186 participants for 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. In addition, 70 duplicates and 72 HapMap control samples were genotyped for quality assurance purposes. Data cleaning was performed by the University of Washington Genetics Coordinating Center (UWGCC) using standard analysis pipelines implemented in the R Environment for Statistical Computing, as previously described [62]. These analyses include interrogating samples for genetic sex, chromosomal anomalies, relatedness among participants, missing call rate, and batch effects, and interrogating SNPs for missing call rate, discordance between duplicate samples, Mendelian errors (as measured in HapMap control parent-offspring trios), Hardy-Weinberg equilibrium, and differences in allele frequency and heterozygosity between sexes (for autosomal and pseudo-autosomal SNPs). Supplemental S6 Table shows the number of SNPs omitted and retained for each quality filter. Imputation was performed to capture information on unobserved SNPs as well as sporadically missing genotypes among genotyped SNPs, using haplotypes from the 1000 Genomes Project [63] Phase 3 reference panel of 2,504 samples from 26 worldwide populations. First, pre-phasing was performed in SHAPEIT2 [64], and then imputation of 34,985,077 variants was performed in IMPUTE2 [65,66]. Masked variant analysis–that is, imputation of genotyped SNPs as though they were unobserved in order to assess imputation quality–showed high concordance between imputed and observed genotypes (0.998 for SNPs with MAF < 0.05 and 0.982 for SNPs with MAF ≥ 0.05) indicating high quality imputation. Imputed SNPs were included in analyses if the minimum genotype probability for a given variant was greater than 50%. Principal component analysis using 96,700 autosomal SNPs pruned from the total panel based on call rate (> 95%), MAF (> 0.05), and LD (pairwise r2 < 0.1 in a sliding window of 10 Mb), was used to assess population structure. Supplemental S23 Fig depicts the observed genetic structure of the population across the first two principal components of ancestry (i.e., eigenvectors from the PCA). Linear regression was used to test the association between each PC, as the dependent variable, and each SNP in the genome. These analyses confirmed that none of the first 20 principal components were due to local variation in specific genomic regions. Statistical approach Prior to genetic analysis, each of the 20 linear distance measures was adjusted for the effects of sex, age, age2, height, weight, and facial size (calculated as the geometric mean of the linear distance measures) using linear regression in order to generate 20 adjusted phenotypes (i.e., residuals). The inclusion of age and age2 as covariates was done in an effort to adjust for both linear and non-linear aspects of age on the phenotypes. After model fitting different sets of covariates, including more complicated spline functions, we settled on a combination of age and age2 as the most reasonable approach based on akaike information criterion values calculated across age-adjustment models. Linear models were then used to test genetic association between each phenotype and each SNP, under the additive genetic model, while simultaneously adjusting for the first four principal components of ancestry. For SNPs on the X chromosome we coded hemizygous males as 0/2 so they are on the same scale as 0/1/2 females. Analyses were performed separately for the Pittsburgh and Denver cohorts, and combined via inverse variance-weighted meta-analysis. To appropriately model SNP effects, we required that the minor allele be present in at least 30 participants, corresponding to a MAF threshold of 0.6% in the Pittsburgh cohort and 2% in the Denver cohort. SNPs meeting the minor allele count criterion in both Pittsburgh and Denver cohorts were included in the meta-analysis. The final number of genotyped SNPs available for analysis after minor allele filtering was 659,955 for the Pittsburgh sample, 638,772 for the Denver sample, and 637,391 for the meta-analysis. The number of imputed and total (genotyped plus imputed) SNPs is available in S7 Table. Given the large number of tests, we used the conventional threshold of p < 5 x 10−8 (i.e., Bonferroni correction for 1 million tests) for genome-wide statistical significance. Because we expect many of our traits to be correlated, we used the eigenvalue method described by Li and Ji [67] to determine that the effective number of independent traits was 10. Thus, we set the threshold for study-wide statistical significance at p < 5 x 10−9 (i.e. p < 5 x 10−8 divided by 10). Because these thresholds are very conservative, we also reported “suggestive” evidence of association of p < 5 x 10−7 in S1–S3 Tables. Phenotypes were generated using the R Environment for Statistical Computing, and genetic association was performed using PLINK [68]. Availability of data All of the phenotypic measures and genotypic markers used here are available to the research community through the dbGaP controlled access repository (http://www.ncbi.nlm.nih.gov/gap) at accession number: phs000949.v1.p1. The raw source data for the phenotypes–the 3D facial surface models–are available for the 3D Facial Norms dataset through the FaceBase Consortium (www.facebase.org). Finally, searchable results datasets of the p-values from the studies reported here are available through the FaceBase Human Genomics Analysis Interface (facebase.sdmgenetics.pitt.edu).

Acknowledgments The authors thank all the individuals who participated in this study as well as the following individuals for their help with recruitment and data management: Pooja Gandhi, Carla Sanchez, Beth Emanuele, Rebecca DeSensi, Laura Stueckle, Linda Peters, Erik Stuhaug, Eden Palmer, Trylla Tuttle, Lidia Leonard, Maria Elena Serna, Rosa N. Martinez, Syed Hashmi, Jennifer Rigdon, Samantha Stachowiak and Paulene Holland.

Author Contributions Conceived and designed the experiments: SMW JRS EF RAS MLM. Analyzed the data: JRS EO MKL EJL EF MLM SMW. Contributed reagents/materials/analysis tools: ZDR CLH MLC JTH CHK NLN LMM GLW JCM CAL CCL JC TF OK WM BH. Wrote the paper: SMW JRS EO MKL EJL EF CAL CCL SS BH RAS MLM.