Article Figures and data Abstract Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract The human face represents a combined set of highly heritable phenotypes, but knowledge on its genetic architecture remains limited, despite the relevance for various fields. A series of genome-wide association studies on 78 facial shape phenotypes quantified from 3-dimensional facial images of 10,115 Europeans identified 24 genetic loci reaching study-wide suggestive association (p < 5 × 10−8), among which 17 were previously unreported. A follow-up multi-ethnic study in additional 7917 individuals confirmed 10 loci including six unreported ones (padjusted < 2.1 × 10−3). A global map of derived polygenic face scores assembled facial features in major continental groups consistent with anthropological knowledge. Analyses of epigenomic datasets from cranial neural crest cells revealed abundant cis-regulatory activities at the face-associated genetic loci. Luciferase reporter assays in neural crest progenitor cells highlighted enhancer activities of several face-associated DNA variants. These results substantially advance our understanding of the genetic basis underlying human facial variation and provide candidates for future in-vivo functional studies. Introduction The human face represents a multi-dimensional set of correlated, mostly symmetric, complex phenotypes with high heritability. This is illustrated by a large degree of facial similarity between monozygotic twins and, albeit less so, between other relatives (Djordjevic et al., 2016), stable facial features within and differences between major human populations, and the enormous diversity among unrelated persons almost at the level of human individualization. Understanding the genetic basis of human facial variation has important implications for several disciplines of fundamental and applied sciences, including human genetics, developmental biology, evolutionary biology, medical genetics, and forensics. However, current knowledge on the specific genes influencing human facial appearance and the underlying molecular mechanisms forming facial morphology is limited. Over recent years, nine separate genome-wide association studies (GWASs) have each highlighted several but largely non-overlapping genetic loci associated with facial shape phenotypes (Adhikari et al., 2016; Cha et al., 2018; Claes et al., 2018; Cole et al., 2016; Lee et al., 2017; Liu et al., 2012; Paternoster et al., 2012; Pickrell et al., 2016; Shaffer et al., 2016). The overlapping loci from independent facial GWASs, thus currently representing the most established genetic findings, include DNA variants in or close to 10 genes that is CACNA2D3 (Paternoster et al., 2012; Pickrell et al., 2016), DCHS2 (Adhikari et al., 2016; Claes et al., 2018), EPHB3 (Claes et al., 2018; Pickrell et al., 2016), HOXD cluster (Claes et al., 2018; Pickrell et al., 2016), PAX1 (Adhikari et al., 2016; Shaffer et al., 2016), PAX3 (Adhikari et al., 2016; Claes et al., 2018; Liu et al., 2012; Paternoster et al., 2012; Pickrell et al., 2016), PKDCC (Claes et al., 2018; Pickrell et al., 2016), SOX9 (Cha et al., 2018; Claes et al., 2018; Pickrell et al., 2016), SUPT3H (Adhikari et al., 2016; Claes et al., 2018; Pickrell et al., 2016), and TBX15 (Claes et al., 2018; Pickrell et al., 2016), which were identified mostly in Europeans. Indicated by the small effect size of the previously identified face-associated DNA variants, together with the high heritability of facial phenotypes, the true number of genes that influence polygenetic facial traits is likely to be much larger than currently known. Despite the previous work, our understanding of the complex genetic architecture of human facial morphology thus remains largely incomplete. This emphasizes on the need for well-powered GWASs that can be achieved through large collaborative efforts, which has been recently demonstrated for appearance traits involving pigmentation variation (Hysi et al., 2018; Visconti et al., 2018), but is currently missing for the human face. A most recent study (Claes et al., 2018) suggested that variants at genetic loci implicated in human facial shape are involved in the epigenetic regulation of human neural crest cells, suggesting that at least some of the functional variants may reside within regulatory elements such as enhancers. Although this is in line with evidence for other appearance phenotypes such as human pigmentation traits, where experimental studies have proven enhancer effects of pigmentation-associated DNA variants (Visser et al., 2012; Visser et al., 2014; Visser et al., 2015), experimental evidence for understanding the functional basis of DNA variants for which facial phenotype associations have been previously established is missing as of yet. Led by the International Visible Trait Genetics (VisiGen) Consortium and together with its study partners, the current study represents a collaborative effort to identify novel genetic variants involved in human facial variation in the largest set of multi-ethnic samples available thus far. Moreover, it demonstrates the functional consequences of identified face-associated DNA variants based on in-silico work and in-vitro experimental evidence. Results Facial phenotypes The current study included seven cohorts, totaling 18,032 individuals. Demographic characteristics and phenotyping details for all cohorts are provided in Materials and methods and Supplementary file 1 - Table 1. A variety of different phenotyping methods are available for 3D facial image data, where some approaches (Claes et al., 2018) require access to the underlying raw image datasets, which often cannot be shared between different research groups. In order to maximize sample size via a collaborative study, we focused on linear distance measures of facial features that can be accurately derived from facial image datasets. In our cohorts with 3D facial surface images, we focused on distances between a common set of thirteen well-defined facial landmarks (Figure 1). In the Rotterdam Study (RS) and the TwinsUK study (TwinsUK), we used the same ensemble method for automated 3D face landmarking that was specifically designed for large cohort studies with the flexibility of landmark choice and high robustness in dealing with various image qualities (de Jong et al., 2018; de Jong et al., 2016). Facial landmarking in Avon Longitudinal Study of Parents and their Children (ALSPAC) and Pittsburgh 3D Facial Norms study (PITT) was done manually. Generalized Procrustes Analysis (GPA) was used to remove affine variations due to shifting, rotation and scaling (Figure 1—figure supplement 1). A total of 78 Euclidean distances involving all combinatorial pairs of the 13 facial landmarks were considered as facial phenotypes. Figure 1 with 5 supplements see all Download asset Open asset Manhattan plot from meta-analysis discovery GWAS (N = 10,115 Europeans) for 78 Euclidean distances between 13 facial landmarks. Results from meta-analysis of four GWASs in Europeans (N = 10,115), which were separately conducted in RS, TwinsUK, ALSPAC, and PITT for 78 facial shape phenotypes that are displayed in a signal 'composite' Manhattan plot at the bottom of the figure. All loci consisting of SNPs reaching study-wide suggestive significance (p<5 × 10−8, red line) were nominated according to the nearby candidate genes in the associated loci, and the top 10 loci were highlighted in colors. The black line indicates study-wide significance after multiple trait correction (p<1.2×10−9). Novel and previously established face-associated genetic loci were differentiated using colored and gray text, respectively. Annotation of the 13 facial landmarks used to derive the 78 facial phenotypes and the associated facial phenotypes are illustrated in the upper part of the figure. Detailed phenotype analyses were conducted in RS. In RS, almost all (72 out of 78, 92%,) of the facial phenotypes followed the normal distribution as formally investigated by applying the Shapiro-Wilk normality test and obtaining non-significant results (Supplementary file 1 - Table 2). For only six phenotypes, this test revealed significant results; however, their histogram looked largely normal. A highly significant sex effect (min p=2.1 × 10−151) and aging effect (min p=8.0 × 10−44) were noted, together explained up to 21% of the variance per facial phenotype (e.g. between alares and Ls; Figure 1—figure supplement 2A–2C, Supplementary file 1 - Table 3). Intra-organ facial phenotypes (e.g. within nose or within eyes) showed on average higher correlations than inter-organ ones (e.g. between nose and eyes), and symmetric facial phenotypes showed higher correlations than non-symmetric ones. Genetic correlations (Zhou and Stephens, 2014) estimated using genome-wide data were similar to correlations obtained directly from phenotype data, with on average higher absolute values seen based on genotypes (Figure 1—figure supplement 3, Supplementary file 1 - Tables 4 and 5). Unsupervised hierarchical clustering analysis identified four distinct clusters of facial shape phenotypes, for example many distances related to pronasale were clustered together (Figure 1—figure supplement 2D). Twin heritability of all facial phenotypes was estimated in TwinsUK as h2 mean(sd)=0.52 (0.15), max = 0.86 and in Queensland Institute of Medical Research study (QIMR) as h2 mean(sd)=0.55 (0.19), max = 0.91 (Supplementary file 1 - Table 6). The most heritable features were linked to the central face (i.e., the outer- and inter- orbital distances, the nose breadth, and the distance between the subnasion and the inner eye corners; Figure 1—figure supplement 2E and F). These twin-derived heritability estimates were generally higher than those obtained from genome-wide SNPs in the same cohorts (Supplementary file 1 - Table 6). These results were overall consistent with expectations and suggest that the phenotype correlation within the human face may be explained in part by sets of shared genetic components. GWASs and replications We conducted a series of GWASs and meta-analyses to test the genetic association of 7,029,494 autosomal SNPs with 78 face phenotypes in RS, TwinsUK, ALSPAC, and PITT, totaling 10,115 subjects of European descent, used as the discovery dataset. Inflation factors were in an acceptable range in all meta-analyses (average λ = 1.015, sd = 0.006; Figure 1—figure supplement 4), and were thus not further considered. A matrix spectral decomposition analysis estimated the effective number of independent facial traits as 43, which corresponds to a Bonferroni corrected study-wide significant threshold of 1.2 × 10−9. A power analysis showed that under this threshold, our discovery data had over 90% power to detect an allelic effect explaining 0.54% variance for individual phenotypes (Figure 1—figure supplement 5A). In the discovery GWAS meta-analysis, a total of 221 SNPs from six distinct loci showed study-wide significant association (p<1.2 × 10−9) with facial phenotypes, among which two have not been previously reported: 4q28.1 INTU (top SNP rs1250495 p=3.03 × 10−10), and 12q24.21 TBX3 (rs1863716 p=1.47 × 10−10). The other four have been described in previous GWASs of facial variation, including 2q36.1 PAX3, 4q31.3 SFRP2, 17q24.3 CASC17, and 20p11.22 PAX1 (Figure 1 and Supplementary file 1 - Table 7) (Liu et al., 2012; Paternoster et al., 2012; Adhikari et al., 2016; Shaffer et al., 2016; Claes et al., 2018). Considering the traditional genome-wide significant threshold of 5 × 10−8 as our study-wide suggestive significance threshold for replication studies, the power analysis showed that under this threshold, our discovery data had over 90% power to detect an allelic effect explaining 0.43% variance for individual phenotypes (Figure 1—figure supplement 5B). In addition to the study-wide significant ones, 273 SNPs from 18 distinct loci passed the 5 × 10−8 threshold (Supplementary file 1 - Table 7), of which 15 were novel and 3 (1p12 TBX15, 9q22.31 ROR2 and 17q24.3 SOX9) have been reported in previous studies (Cha et al., 2018; Claes et al., 2018; Pickrell et al., 2016). Thus, in total, we identified 24 face-associated genetic loci, of which 17 were not reported in previous face GWASs, that we followed-up via a multi-ethnic replication analysis (see below). In general, no significant heterogeneity was observed for the 24 top SNPs of the study-wide significant and suggestive associated loci, that is only one SNP (rs3403289 at 2q36.1 PAX3) showed nominally significant (Cochran's Q test p=0.04, Table 1) heterogeneity between the discovery cohorts, but was not significant after multiple testing correction. Table 1 SNPs significantly associated (p<5×10−8) with facial shape phenotypes from European discovery GWAS meta-analysis (RS, TwinksUK, ALSPAC, and PITT), and their multi-ethnic replication (UYG, CANDELA and QIMR). Discovery meta-analysisReplication(N = 10,115)(N = 7,917)RegionSNPNearest GeneEAOATraitBetaPQCom. PNovel face-associated loci1p36.22rs143353512CASZ1AGPrn-AlL−0.296.44 × 10−90.730.00011p36.13rs200243292ARHGEF19ITEnR-ChL−0.111.46 × 10−80.540.06401p31.2rs77142479RPE65CAEnL-Sn−0.087.65 × 10−90.900.01212p12rs10202675LRRTM4TCEnR-AlR−0.264.43 × 10−80.400.02262q31.1rs2884836KIAA1715TCExR-ChR−0.093.04 × 10−80.520.10963q12.1rs113663609CMSS1AGEnR-ChR−0.123.46 × 10−80.960.07824q28.1rs12504954INTUAGPrn-AlL−0.093.03 × 10−100.520.00156p22.3rs2225718RNF144BCTEnR-AlR−0.092.97 × 10−80.600.00716p21.2rs7738892KIF6TCEnR-N0.111.34 × 10−80.090.00188p23.2rs1700048CSMD1CAExR-ChL−0.223.94 × 10−80.730.11318q21.3rs9642796DCAF4L2AGAlL-Ls−0.114.52 × 10−80.070.309010q22.1rs201719697SUPV3L1ADExL-AlL−0.283.42 × 10−80.66NA12q24.21rs1863716TBX3ACPrn-AlR−0.091.47 × 10−100.532.45 × 10−513q14.3rs7325564LINC00371TCN-Prn0.091.34 × 10−80.160.351914q32.2rs1989285C14orf64GCN-Prn−0.104.54 × 10−80.820.001916q12.1rs16949899SALL1TGExR-ChL−0.133.35 × 10−80.760.891216q12.2rs7404301RPGRIP1LGAAlL-Ls−0.093.49 × 10−90.432.87 × 10−7Previously reported face-associated loci1p12rs1229119TBX15TCExR-ChR0.082.25 × 10−80.260.12192q36.1rs34032897PAX3GAEnR-N0.141.96 × 10−170.040.00204q31.3rs6535972SFRP2CGEnL-AlL0.123.51 × 10−120.125.89 × 10−129q22.31rs2230578ROR2CTEnL-Sn0.099.02 × 10−90.661.09 × 10−617q24.3rs8077906CASC17AGPrn-EnL−0.093.00 × 10−110.240.050517q24.3rs35473710SOX9GAAlR-ChL0.193.95 × 10−80.370.263420p11.22rs4813454PAX1TCAlL-AlR−0.102.32 × 10−110.580.0002 Except rs2230578 (3'-UTR of ROR2) and rs7738892 (a synonymous variant of KIF6), all SNPs are intronic or intergenic. Gene symbols in bold indicate successful replication after Bonferroni correction of 24 loci (p<2.1×10−3). Meta P in bold indicates study-wide significance in the discovery analysis after multiple trait correction (p<1.2×10−9). Replication P in bold means significant after Bonferroni correction for multiple tests (p<2.1×10−3). Q: p-value from the Cochran's Q test for testing heterogeneity between discovery cohorts. EA and OA, the effect allele and the another allele. Com. P: p-values from a combined test of dependent tests. NA, not available. The association signals at the 24 genetic loci involved multiple facial phenotypes clustered to the central region of the face above the upper lip (Figure 1), and closely resembled the twins-based face heritability map (Figure 1—figure supplement 2E and F). All seven previously reported loci showed largely consistent allele effects on the same or similar sets of facial phenotypes as described in previous GWASs (Supplementary file 1 - Table 8). Four (2q36.1 PAX3, 9q22.31 ROR2, 17q24.3 CASC17, and 20p11.22 PAX1) out of the 24 face-associated loci have been noted in the GWAS catalog for face morphology or related phenotypes that is nose size, chin dimples, monobrow thickness, and male pattern baldness (Supplementary file 1 - Table 9). Five out of the 24 loci showed genome-wide significant association (p<5e-8) with non-facial phenotypes and/or diseases in the GeneATLAS database (Canela-Xandri et al., 2018) (Supplementary file 1 - Table 10), suggesting pleiotropic effects. Six out of the 24 loci showed a high degree of diversity among the major continental groups in terms of both allele frequency differences (>0.2) and Fst values (empirical p<0.05), including 1p31.2 RPE65, 2q36.1 PAX3, 4q28.1 INTU, 16q12.1 SALL1, 16q12.2 RPGRIP1L and 17q24.3 CASC17 (Supplementary file 1 - Tables 11 and 12). Three out of the 24 face-associated loci showed signals of positive selection in terms of the integrated haplotype scores (iHS, empirical p<0.01), including RPE65, RPGRIP1L and CASC17 (Supplementary file 1 - Table 12). None of the 24 loci overlapped with the previously reported signals of singleton density score (SDS, empirical p<0.01) in Europeans (Field et al., 2016), which may indicate that the observed selection signals are unlikely the consequences from the recent selective pressures detectable via SDS analysis. We selected the 24 top-associated SNPs (Table 1), one SNP per for each of the identified 24 genetic loci passing the study-wide suggestive significance threshold (p<5×10−8), for replication studies in a total of 7917 multi-ethnic subjects from three additional cohorts that is the Consortium for the Analysis of the Diversity and Evolution of Latin America (CANDELA) study involving Latin Americans known to be of admixed Native American and European ancestry, the Xinjiang Uyghur (UYG) study from China known to be of admixed East Asian and European ancestry, and the Queensland Institute of Medical Research (QIMR) study from Australia involving Europeans. We used a 'combined test of dependent tests' (Kost and McDermott, 2002) to account for the incompatibility among facial traits from the three replication cohorts, which tests for H0: no association for all facial phenotypes in all replication cohorts vs. H1: associated with at least one facial phenotypes in at least one replication cohorts. After Bonferroni correction of multiple testing of 24 SNPs, significant evidence of replication (p<2.1×10−3) was obtained for 10 SNPs, among which six were novel loci: 1p36.22 CASZ1, 4q28.1 INTU, 6p21.2 KIF6, 12q24.21 TBX3, 14q32.2 C14orf64, and 16q12.2 RPGRIP1L (Figure 2), and four have been described in previous face GWASs: 2q36.1 PAX3, 4q31.3 SFRP2, 9q22.31 ROR2, and 20p11.22 PAX1 (Adhikari et al., 2016; Claes et al., 2018; Paternoster et al., 2012; Pickrell et al., 2016; Shaffer et al., 2016). Figure 2 with 24 supplements see all Download asset Open asset Six novel genetic loci associated with facial shape phenotypes and successfully replicated. The figure is composed by six parts, (A) 1p36.22 CASZ1 region, (B) 4q28.1 INTU region, (C) 6p21.2 KIF6 region, (D) 12q24.21 TBX3 region, (E) 14q32.2 C14orf64 region, and (F) 16q12.2 RPGRIP1L region. Each part includes two figures, Face map (left) denoting all of the top-SNP-associated phenotypes weighted by the statistical significance (-log10P). The dashed line means that the P value was nominally significant (p<0.01) but did not reach study-wide suggestive significance (p>5e-8). LocusZoom (right) shows regional association plots for the top-associated facial phenotypes with candidate genes aligned below according to the chromosomal positions (GRCh37.p13) followed by linkage disequilibrium (LD) patterns (r2) of EUR in the corresponding regions. Similar figures for all 24 study-wide suggestively significant loci are provided in Figure 2—figure supplements 1–24. The allelic effects of the replicated SNPs were highly consistent across all participating cohorts, while the individual cohorts were underpowered (Figure 2—figure supplements 1–6), further emphasizing on the merit of our collaborative study. The associated phenotypes of the replicated SNPs demonstrated a clear pattern of facial symmetry that is the same association was observed on both sides of the face, while such a symmetric pattern was less clear for most phenotypes involved in the non-replicated loci (Figure 2—figure supplements 7–17). Several novel replicated loci showed effects on the similar sets of facial shape phenotypes; for example, 1p36.22 CASZ1 and 4q28.1 INTU (Figure 2A and C) both affected the distances between pronasion and alares. Moreover, several novel replicated loci affected phenotypes in multiple facial regions for example 4q32.1 C14orf64 and 12q24.1 TBX3 (Figure 2D and E) both affected many distances among eyes, alares and mouth. It should be noted however, that the discovery dataset is of European ancestry while the replication dataset is multi-ethnic and the facial phenotypes could not be obtained in the same way as in the discovery set; hence, the replication results were rather conservative, that is the possibility of false negatives in this replication study cannot be excluded. This was evident from the results of the replication attempts of the seven known loci (1p12 TBX15, 2q36.1 PAX3, 4q31.3 SFRP2, 9q22.31 ROR2, 17q24.3 CASC17, 17q24.3 SOX9, 20p11.22 PAX1) (Figure 2—figure supplements 18–24). Out of these 7, only four were replicated (PAX3, SFRP2, ROR2 and PAX1, Table 1). For this reason, we considered all of the 494 SNPs from the 24 loci passing the study-wide suggestive threshold in the subsequent polygenic score analyses and in the functional annotation analyses. Integration with previous literature knowledge In our discovery GWAS dataset, we also investigated 122 face-associated SNPs reported in previously published facial GWASs (Adhikari et al., 2016; Cha et al., 2018; Claes et al., 2018; Cole et al., 2016; Lee et al., 2017; Liu et al., 2012; Paternoster et al., 2012; Pickrell et al., 2016; Shaffer et al., 2016). Out of those, a total of 44 (36%) SNPs showed significant association on either the same set of facial phenotypes or at least one facial phenotypes in a combined test, after Bonferroni correction for the number of SNPs (adjusted p<0.05) (Supplementary file 1 - Table 8). Notably, none of the 5 SNPs reported in a previous African GWAS (Cole et al., 2016) was replicated in our European dataset. This is unlikely explained by the allele frequency differences between Sub-Saharan African (AFR) and Europeans (EUR) (absolute frequency differences between AFR and EUR were less than 0.1 in the 1000 Genomes project) but more likely due to different phenotypes, that is in the previous study these loci were primarily associated with facial size, which was not tested on our study. We also investigated the potential links between facial phenotypes and 60 SNPs in 38 loci previously implicated in non-syndromic cleft lip and palate (NSCL/P) phenotype, which is linked with normal variation of facial morphology as suggested previously (Beaty et al., 2010; Beaty et al., 2011; Birnbaum et al., 2009; Grant et al., 2009; Leslie et al., 2017; Ludwig et al., 2017; Ludwig et al., 2012; Mangold et al., 2010; Sun et al., 2015; Yu et al., 2017). Among these, 17 SNPs in 13 loci showed significant face-associations in the discovery cohorts after Bonferroni correction for the number of SNPs and phenotypes, and many (11 SNPs in eight loci) involving cleft related facial landmarks as expected, for example subnasale, labliale superius, and labiale inferius (Supplementary file 1 - Table 13). Multivariable model fitting and polygenic face scores In the RS cohort, a multiple-regression conditioning on the effects of the 24 lead SNPs identified 31 SNPs showing significant (adjusted p<0.05) independent effects on sex- and age-adjusted face phenotypes (Supplementary file 1 - Table 14). These 31 SNPs individually explain less than 1%, and together explain up to 4.62%, of the variance for individual phenotypes (Figure 3A and B). The top-explained phenotypes clustered to the central region of the face, that is the distances between nasion, subnasale, pronasale, left endocanthion and right endocanthion. The variance explained was potentially overestimated in RS cohort since the SNP selection and model building were conducted using the same dataset. However, potential overfitting is unexpected to have a drastic effect, given that some of the selected SNPs showed weaker effects in RS than in other cohorts used (Supplementary file 1 - Table 7). We then used the fitted models to derive a set of polygenetic face scores for EUR, AFR and East Asian (EAS) subjects (total N = 1,668) from the 1000-Genomes Project. The polygenic face scores based on the 31 SNPs were largely reversely distributed between EUR and AFR, with EAS situated in between (Figure 3C). EUR showed the longest distances among the vertical phenotypes, while AFR had the widest of the horizontal phenotypes. EAS was consistently in-between with the most differential phenotype scores observed for nose wing breadth (AFR > EAS > EUR) and the nose length (EUR > EAS > AFR). An exception was for mouth phenotypes, where AFR had the largest values in both mouth width and lip thickness (Figure 3D). The nose-related score distribution, for example nose breadth and length, thus largely assembled the facial features in major continental groups, consistent with the hypotheses regarding adaptive evolution of human facial morphology (Noback et al., 2011; Weiner, 1954), for example nose width and nose length have been found to correlate with temperature and humidity. It should be noted, however, that the possibility of potential bias cannot be excluded due to the fact that the discovery GWAS was conducted in only Europeans. Notably, the polygenic scores of the NSCL/P-associated SNPs explained up to 1.2% of the age- and sex- adjusted phenotypic variance for phenotypes mostly between nose and mouth, for example right alare – right cheilion and pronasale – right cheilion (Figure 3—figure supplement 1). Figure 3 with 1 supplement see all Download asset Open asset Polygenic scores of 78 facial phenotypes. (A) Facial phenotype variance is explained by a multivariable model consisting of 31 face-associated SNPs with independent effects. (B) Unsupervised clustering reveals the details on the explained phenotype variance, where the SNPs were attributed to the candidate genes in the associated regions. All phenotypes with the total explained variance > 2% are shown. SNPs were colored according to their raw p-values and those passing FDR < 0.05 were indicated in squares. (C) Mean standardized polygenic scores calculated by applying the multivariable model trained in the RS discovery sample to Sub-Saharan Africans (AFR), East Asians (EAS) and Europeans (EUR) in the 1000 Genomes Project data. (D) Examples of the distributions of the face scores in samples of different ancestry origin/country. Gene ontology of face-associated genetic loci A gene ontology enrichment analysis for the 24 face-associated genetic loci (Table 1) highlighted a total of 67 biological process terms and seven molecular function terms significantly enriched with different genes group (FDR < 0.01 or q value < 0.01, Supplementary file 1 - Table 15), with the top significant terms being 'embryonic digit morphogenesis' and 'embryonic appendage morphogenesis' (FDR < 1 × 10−8, Figure 4—figure supplement 1). Biological process terms were categorized in four broader categories, 'morphogenesis', 'differentiation', 'development' and 'regulation'. These results underlined the important role of embryonic developmental and regulatory processes in shaping human facial morphology. Furthermore, cis-eQTL effects of the 494 face-associated SNPs at the 24 loci (Supplementary file 1 - Table 7) were investigated using the GTEx data. Significant multi-tissue cis-eQTL effects (adjusted p<0.05) were found around three genes, ARHGEF19, RPE65, and TBX15 (Supplementary file 1 - Table 16). No significant tissue-specific eQTL effects were observed, likely due to the lack of cell types indexed within the GTEx database that are directly involved in craniofacial morphology. Expression of face-associated genetic loci in embryonic cranial neural crest cells Embryonic cranial neural crest cells (CNCCs) arise during weeks 3–6 of human gestation from the dorsal part of the neural tube ectoderm and migrate into the branchial arches. They later form the embryonic face, consequently establishing the central plan of facial morphology and determining species-specific and individual facial variation (Bronner and LeDouarin, 2012; Cordero et al., 2011). Recently, CNCCs have been successfully derived in-vitro from human embryonic stem cells or induced pluripotent stem cells (Prescott et al., 2015). Taking advantage of the available CNCC RNA-seq data (Prescott et al., 2015) and other public RNA-seq datasets compiled from GTEx (Lonsdale et al., 2013) and ENCODE (ENCODE Project Consortium, 2012), we compared 24 genes, which were nearest to the top-associated SNPs at the 24 face-associated genetic loci identified in our discovery GWAS meta-analysis, with the background genes over the genome for their expression in CNCCs and 49 other cell types. In 89% of 9 embryonic stem cells, 65% of 20 different tissue cells, and 15% of 20 primary cells these 24 genes showed significantly higher expression in all cell types (p<0.05) after multiple testing correction compared to the background genes (Figure 4A), and the most significant difference was observed in CNCCs (p=3.8 × 10−6). Furthermore, these 24 genes also showed preferential expression in CNCCs compared to their expression in other cell types, that is in 65% of 20 different tissue cells, in 60% of 20 primary cells, and in 22% of 9 embryonic stem cells (Figure 4A). Specifically, 16 (66.7%) of these 24 genes showed preferential expression in CNCCs compared to 49 other cell types. The top-ranked three preferentially expressed genes w