Genome-Wide Association Study Reveals Multiple Loci Influencing Normal Human Facial Morphology.

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.

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 () 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.

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.

Results

In total, we observed seven associations in five traits that exceeded the conventional threshold for genome-wide significance (p < 5 x 10−8, ; Figs ). One of the associations also exceeded our study-wide significance threshold of p < 5 x 10−9, calculated based on 10 independent traits (see Methods for details). Due to the large number of traits, we will limit our presentation of results to genome-wide significant signals. The entire list of meta-analysis associations with p-values < 5 x 10−7 is available in . Manhattan plots showing the meta-analysis results, as well as the results for each sample, are available in supplemental figures .

LocusZoom plots showing genome-wide significant associations observed in the meta-analysis for cranial base width (Fig 1A).

(A) chromosome 14 and (B) chromosome 20. LocusZoom plots show the association (left y-axis; log10-transformed p-values) with facial traits. Genotyped SNPs are depicted by stars and imputed SNPs are depicted by circles. Shading of the points represent the linkage disequilibrium (r2, based on the 1000 Genomes Project Europeans) between each SNP and the top SNP, indicated by purple shading. The blue overlay shows the recombination rate (right y-axis). Positions of genes are shown below the plot.

LocusZoom plots showing genome-wide significant associations observed in the meta-analysis for intercanthal width (Fig 1H).

(A) chromosome 1 and (B) chromosome X. LocusZoom plots show the association (left y-axis; log10-transformed p-values) with facial traits. Genotyped SNPs are depicted by stars and imputed SNPs are depicted by circles. Shading of the points represent the linkage disequilibrium (r2, based on the 1000 Genomes Project Europeans) between each SNP and the top SNP, indicated by purple shading. The blue overlay shows the recombination rate (right y-axis). Positions of genes are shown below the plot.

LocusZoom plots showing genome-wide significant associations observed in the meta-analysis for nasal width (Fig 1K) and nasal ala length (Fig 1N).

(A) nasal width on chromosome 20 and (B) nasal ala length on chromosome 14. LocusZoom plots show the association (left y-axis; log10-transformed p-values) with facial traits. Genotyped SNPs are depicted by stars and imputed SNPs are depicted by circles. Shading of the points represent the linkage disequilibrium (r2, based on the 1000 Genomes Project Europeans) between each SNP and the top SNP, indicated by purple shading. The blue overlay shows the recombination rate (right y-axis). Positions of genes are shown below the plot.

LocusZoom plot showing genome-wide significant association observed in the meta-analysis for upper facial depth (Fig 1B) on chromosome 11.

LocusZoom plots show the association (left y-axis; log10-transformed p-values) with facial traits. Genotyped SNPs are depicted by stars and imputed SNPs are depicted by circles. Shading of the points represent the linkage disequilibrium (r2, based on the 1000 Genomes Project Europeans) between each SNP and the top SNP, indicated by purple shading. The blue overlay shows the recombination rate (right y-axis). Positions of genes are shown below the plot.

We observed two significant associations for cranial base width: one at 14q21.1 (top SNP rs79272428, p = 1.01 x 10−8, ) and the other at 20q12 (top SNP rs6129564, p = 1.65 x 10−9, ). Notably, the chromosome 20 association exceeded our strict threshold for study-wide statistical significance. For intercanthal width, we observed two significant associations: one at 1p13.3 (top SNP rs619686, p = 4.70 x 10−8, ) and the other at Xq13.2 (top SNP rs11093404, 4.16 x 10−8, ). There were also significant associations with nasal width at 20p11.22 (rs2424399, p = 2.62 x 10−8, ) and nasal ala length at 14q11.2 (top SNP rs8007643, p = 3.36 x 10−8, ). We observed a second independent peak on chromosome 20 for nasal width located 371kb upstream of the main peak. The second peak remained (top SNP rs80186620, p = 5.32x10-6, ) after conditional association analysis adjusting for the effects of rs2424399 on nasal width. Finally we observed a significant association with upper facial depth at 11q22.1 (top SNP rs12786942, p = 4.59 x 10−8, ). For all of the above associations, the results were driven primarily by the larger Pittsburgh dataset. The cranial base width (14q21.1), intercanthal width (1p13.3) and upper facial depth associations were at least nominally significant (p < 0.05) in both datasets. Sample-specific association results for all SNPs with p-values less than 5 x 10−7 are listed in for the Pittsburgh sample and for the Denver sample.

In an attempt to replicate the main findings from the prior two GWA studies in Europeans, we tested previously implicated SNPs against traits from our Pittsburgh dataset that capture similar aspects of morphology. This was not possible for every prior genotype-phenotype association given differences in the measurements available. With that limitation in mind, the Pittsburgh dataset was chosen for comparison because it was the larger of our two datasets and the phenotyping protocol was most similar to prior GWA studies. For Paternoster et al. [17], we attempted to test three of their four genome-wide significant associations, two of which involved nasal ala length (). In our data, nasal ala length showed a nominally significant association (p = 0.018) with rs1982862, an intronic variant in CACNA2D3. Conversely, we found no evidence of association between this measure and rs11738462, an intronic variant in C5orf64. The previously observed association between PAX3 and the position of nasion relative to the orbits could not be tested directly. However, we found no evidence of association between the implicated SNP rs7559271 and intercanthal width, which captures aspects of interorbital septum morphology. As a further exploratory analysis we also looked at the association between rs7559271 and several vertical or projective measurements involving nasion, but no significant associations were found for any of these traits. For Liu et al. [18], we attempted to test each of their six previously reported genome-wide significant associations (). We observed a strong association (p = 1.70 x 10−5) between nasal ala length and rs4648379, an intronic variant in PRDM16. We also observed an association between rs6555969, a SNP near C5orf50 and upper facial depth (p = 0.005), which is a reasonable approximation of the zygion-nasion distance reported by Liu et al. [18]. To test the association between interorbital distance and rs17447439, an intronic variant in TP63, we used measures of intercanthal and outercanthal width; however, we did not observe an association with either measure. Finally, Liu et al. [18] reported associations between SNPs in PAX3, C5orf50 and COL17A1 and the position of nasion relative to the orbits. We tested these three SNPs in our dataset against intercanthal width, a trait involving roughly similar anatomical components. Notably, we found associations between rs974448 (PAX3, p = 0.002) and rs6555969 (C5orf50, p = 0.049) and intercanthal width.

Discussion

Based on meta-analysis, we observed seven associated loci for five facial traits: cranial base width (), intercanthal width (), nasal width (), nasal ala length (), and upper facial depth (). 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 . 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 (). 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 and listed in . 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 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 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 .

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 . 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).

Associations with p-values < 5 x 10−7 for all 20 traits from genome-wide meta-analysis.

(XLSX)

Click here for additional data file.

Associations with p-values < 5 x 10−7 for all 20 traits from genome-wide analysis of the Pittsburgh sample.

(XLSX)

Click here for additional data file.

Associations with p-values < 5 x 10−7 for all 20 traits from genome-wide analysis of the Denver sample.

(XLSX)

Click here for additional data file.

Descriptive statistics of the two study samples.

(DOCX)

Click here for additional data file.

List of linear distance measurements.

(DOCX)

Click here for additional data file.

Number of SNPs omitted and retained for each quality filter.

(DOCX)

Click here for additional data file.

Final number of SNPs used in our analysis.

(DOCX)

Click here for additional data file.

Manhattan plots for cranial base width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for upper facial depth.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for middle facial depth.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for lower facial depth.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for morphological facial height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for upper facial height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for lower facial height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for intercanthal width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for outercanthal width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for palpebral fissure width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for nasal width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for subnasal width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for nasal protrusion.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for nasal ala length.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for nasal height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for nasal bridge length.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for labial fissure width.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for philtrum length.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for upper lip height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

Manhattan plots for lower lip height.

(A) meta-analysis results, (B) Pittsburgh sample results, and (C) Denver sample results. Lines for p-value thresholds set at 5 x 10−8 for genome-wide significance and 5 x 10−7 for suggestive significance.

(PDF)

Click here for additional data file.

LocusZoom plots comparing the (A) discovery analysis, and (B) conditional analysis for the observed genetic association of nasal width near PAX1 at 20p11.22. Genetic association (left y-axis; log10-transformed p-values) is shown for genotyped SNPs depicted as stars and imputed SNPs depicted as circles. Shading of the points represent the linkage disequilibrium (r2) between each SNP and the rs2424399 (the top SNP from the discovery analysis), indicated by purple shading. The blue overlay shows the recombination rate (right y-axis). Positions of genes are shown below the plot. In the discovery analysis, a possible second peak in low-LD with the rs2424399 was observed around chromosomal position 22.0 Mb. After conditioning on rs2424399, variants at position 22.0 Mb showed some independent evidence of association, although not meeting genome-wide or suggestive thresholds for significance.

(PDF)

Click here for additional data file.

3D facial surface model showing the location of the 24 standard landmarks used to generate the linear distances.

Landmarks shown in frontal view (A) are n = nasion; prn = pronasale; sn = subnasale; ls = labiale superius; sto = stomion; li = labiale inferius; sl = sublabiale; gn = gnathion; en = endocanthion; ex = exocanthion; al = alare; sbal = subalare; cph = crista philtra; ch = chelion (for bilateral points only right side labeled). Landmarks shown in the lateral view (B) are ac = alar curvature point and t = tragion (only left landmark shown for these two bilateral points).

(PDF)

Click here for additional data file.

Plot showing population stratification across the first two principal components of ancestry (EV1 and EV2).

The proportion of total genetic variation explained by each principal component of ancestry is indicated on the axis.

(PDF)

Click here for additional data file.

Table 1

Genome-wide significant meta-analysis results for five traits.

				Pittsburgh Sample				Denver Sample				Meta-Analysis
Trait	SNP	Locus	Minor allele	MAF	N	Beta (se)a	p	MAF	N	Beta (se)a	p	p
Cranial base width (Fig 1A)	rs17106852	14q21.1 (38038468)	G	0.106	2368	-1.104 (0.205)	7.91 x 10−8	0.097	671	-0.858 (0.406)	3.51 x 10−2	1.01 x 10−8
	rs6129564	20q12 (38904203)	A	0.116	2368	-1.210 (0.193)	4.07 x 10−10	0.100	671	-0.449 (0.412)	2.77 x 10−1	1.65 x 10−9
Intercanthal width (Fig 1H)	rs619686	1p13.3 (110218761)	G	0.056	2426	-0.763 (0.163)	3.12 x 10−6	0.052	671	-0.536 (0.186)	4.12 x 10−3	4.70 x 10−8
	rs11093404	Xq13.2 (72289467)	A	0.243	2426	0.427 (0.075)	1.31 x 10−8	0.248	671	0.073 (0.075)	3.32 x 10−1	4.16 x 10−8
Nasal width (Fig 1K)	rs2424399	20p11.22 (21632545)	C	0.235	2429	0.377 (0.070)	9.53 x 10−8	0.300	671	0.177 (0.098)	7.04 x 10−2	2.62 x 10−8
Nasal ala length (Fig 1N)	rs8007643	14q11.2 (21365801)	T	0.067	2426	1.064 (0.186)	1.19 x 10−8	0.069	671	0.221 (0.216)	3.07 x 10−1	3.36 x 10−8
Upper facial depth (Fig 1B)	rs12786942	11q22.1 (101394765)	T	0.119	2368	1.429 (0.317)	6.95 x 10−6	0.105	670	1.674 (0.523)	1.43 x 10−3	4.59 x 10−8

a Beta values determined based on minor allele

Table 2

Testing of previously identified genome-wide significant SNPs in European samples.

	Published GWAS Results						Current Results
Study	Locus	Candidate gene	SNP	Minor allele (MAF)	Beta (p)	Associated phenotype(s)	Closest phenotype(s) in our dataset	Minor allele (MAF)	Beta (p)
Paternoster et al. [17]	12q21.3	TMTC2	rs10862567	T (0.31)	0.181 (4.4x10-8)	Position of the right endocanthion point	N/A	T (0.32)	N/A
	2q36	PAX3	rs7559271	G (0.38)	0.169 (2.2x10-10)	Nasion to mid-endocanthion distance	Intercanthal width (Fig 1H)	G (0.37)	0.066 (0.392)
	3p14.3	CACNA2D3	rs1982862	A (0.12)	-0.257 (1.8x10-8)	Pronasale to alare distance	Nasal ala length (Fig 1N)	A (0.16)	-0.297 (0.018)
	5q12	C5orf64	rs11738462	A (0.17)	-0.204 (1.8x10-8)	Pronasale to alare distance	Nasal ala length (Fig 1N)	A (0.18)	0.028 (0.818)
Liu et al. [18]	1p36.3	PRDM16	rs4648379	T (0.28)	-0.260 (1.1x10-8)	Pronasale to alare distance	Nasal ala length (Fig 1N)	T (0.28)	-0.447 (1.7x10-5)
	2q36	PAX3	rs974448	G (0.17)	0.290 (1.6x10-8)	Nasion to orbit distance a	Intercanthal width (Fig 1H)	G (0.17)	0.171 (0.002)
	3q28	TP63	rs17447439	G (0.04)	-0.910 (4.4x10-8)	Distance between the orbits a	Intercanthal width (Fig 1H)	G (0.04)	0.043 (0.758)
							Outercanthal width (Fig 1I)	G (0.04)	0.261 (0.273)
	5q35.1	C5orf50	rs6555969	T (0.33)	0.410 (1.2x10-9)	Nasion to zygion distance	Upper facial depth (Fig 1B)	T (0.33)	0.615 (0.005)
					0.260 (2.3x10-9)	Nasion to orbit distance a	Intercanthal width (Fig 1H)	T (0.33)	0.159 (0.049)
	10q25.1	COL17A1	rs805722	T (0.19)	0.290 (4.0x10-8)	Nasion to orbit distance a	Intercanthal width (Fig 1H)	T (0.18)	0.047 (0.628)

a orbit landmark measured from MRI at the approximate location of the pupil