Genetic mechanisms underlying East Asian and European Facial differentiation

Facial morphology, the most conspicuous feature of human appearance, is highly heritable. Previous studies on the genetic basis of facial morphology were mainly 3 performed in European populations. Applying a proven data-driven phenotyping 4 and multivariate genome-wide scanning protocol to the largest collection of 3D 5 facial images of an East Asian population to date, we identified 244 leading 6 variants associated with normal-range facial variation, of which 130 are novel. A 7 newly proposed polygenic shape analysis indicates that the effects of the variants 8 on East Asian facial shape can be generalized into the European population. Based on this analysis, we further identified 13 variants mainly related to differences 10 between European and East Asian facial shape. Natural selection analyses suggest 11 that the difference in European and East Asian nose shape is caused by a 12 directional selection, mainly due to a local adaptation in Europeans. Our results 13 expand the knowledge of human facial genetics and illustrates for the first time 14 the underlying genetic basis for facial differences across populations.

Here we performed a genome-wide association study (GWAS) based on the largest 30 collection of 3D facial images in Han Chinese populations to date. Using a proven data-31 driven phenotyping approach, we successfully identified hundreds of associated 32 variants 1,9 . In addition, we identified specific variants distinguishing between European 33 and East Asian facial appearance. We further provided evidence that those population-34 based facial differences, especially for nose shape, are under selection. A schematic 35 overview of our study design can be found in Extended Data Figure 1. 36 37 GWAS on hierarchical facial phenotypes discovered 244 signals 38 To study facial variation from a global to local scale, we used 3D facial surface scans 39 from a large-scale East-Asian population (Methods, Supplementary Table 1) in a 40 discovery (n = 6,968) and replication cohort (n = 2,706), and subsequently combined them in a meta-analysis. A semi-supervised phenotyping procedure defined 63 42 hierarchically arranged facial segments using the discovery cohort (Methods). Next, we 43 performed a canonical correlation analysis (CCA) based GWAS on each facial 44 segment's group of principal components (PC), which defined the linear combination 45 of PCs that are most correlated with each variant. In the replication cohort, the 46 associations were tested based on projections onto the CCA directions followed by 47 Pearson's correlation. Subsequently, we identified 50 independent tests using parallel 48 analysis and permutation test (Methods) 1,20,21 . Thus, besides conventional genome-wide 49 threshold (P = 5×10 -8 ), we set the stricter study-wide threshold to P = 1×10 -9 (P = 5×10 -50 8 /50) after Bonferroni correction. In the discovery datasets, we identified 156 genome-51 wide significant variants (P < 5×10 -8 ) after condition analysis and peak selection 52 (Methods). In which, 81 variants are study-wide significant (P < 1×10 -9 ). Majority of 53 variants (120 out of 156 genome-wide variants, 76 out of 81 study-wide variants) are 54 replicated at nominal threshold (P < 0.05). To increase statistical power, we performed 55 a meta-analysis using Stouffer's method to combine the P-values obtained from the 56 discovery and replication cohort 1,22 . As a result, we identified 244 genome-wide 57 threshold independent variants associated with normal range facial variation ( Fig. 1a; 58 Supplementary Table 2). In which, 151 variants are study-wide significant. According 59 to the anatomical structure, we classified the 244 genome-wide significant variants in 60 ten facial regions (Methods), including forehead, glabella, eye, tempora, zygoma, nose, 61 maxillary, upper mouth, lower mouth, and mandible ( Fig. 1b; Supplementary Table 3). 62 The nose region has most of the variants (107) out of the ten regions.    Shape (PPS) is the mean polygenic shape for a given population (Methods). Using data 176 of EUR (n = 404, including TSI (107), GBR (91), IBS (107) and CEU (99)) and EAS 177 (n = 208, including CHB (103), and CHS (105)) individuals from 1000GP, we 178 calculated the PPS of the two populations for the whole face and ten anatomical facial 179 regions 23 . Thus, the differentiated accumulated genetic effects on facial shape between 180 the two populations then becomes the difference between two polygenic population 181 shapes; − . To visualize and compare this effect to the true population 182 mean facial shapes for each facial region, we constructed EUR and EAS derived shapes 183 by adding and subtracting, respectively, ( − )/2 to and from the 184 overall average facial shape (i.e., a population neutral average face, which was 185 constructed as the average of the EUR and EAS population mean shapes (Methods)). 186 Firstly, we visualized the EUR and EAS derived faces. We used 3D facial scans of EAS 187 and EUR individuals to calculate each population's average face and therefore 188 generated the EUR and EAS mean faces. Compared with the mean face of EUR, we 189 found that EAS had more protrusion in the cheek; more concavity in the forehead, 190 glabella, nose, and mandible ( Fig. 3a-i), which are consistent with a previous study 19 . 191 Interestingly, when we amplify the differentiated accumulated genetic effects five times, 192 the PPS derived faces look remarkably similar to EAS and EUR's actual mean face (Fig.  193 3a-ii). The EUR and EAS derived whole faces show similar facial variations to the 194 ground truth, especially in the glabella and nose region ( Fig. 3a-ii). 195 To test the generalization of the association results from our study to Europeans, we 196 We also performed the same analyses locally in ten anatomical facial regions. The EAS 231 and EUR derived regional facial shapes using leading variants were significantly more 232 similar to the population mean shapes than derived shapes using random variants in the 233 upper mouth, nose, maxillary, glabella, eye, tempora and zygoma in all three 234 measurements of similarity, cosine similarity, Euclidean distance, and EAS-FS (seg 17,  Table 6). As a result, 13 variants that 263 passed filtering are regarded as increasing to EAS-FS (Table1; Supplementary Fig. 4). 264 The 13 variants have a higher FST than the other variants (t.test P < 1.0×10 -16 ) 33 , 265 indicating that these 13 variants have significant allele frequency differences between 266 EUR and EAS. We also conducted an enrichment analysis with Population Branch 267 Statistics (PBS), which identifies alleles that have experienced great changes in 268 frequency in one population (EAS) relative to two reference populations (EUR and 269 YRI) 37 . The mean PBS values of these 13 variants are significantly higher in EAS (P < 270 1.0×10 -16 ) but not in EUR (P = 0.188), which suggests that these variants may be under 271 selection in East Asian. Thus, those variants potentially have some contribution to the 272 morphological differences between Europeans and East Asians. Most of the EAS-273 specific variants might be standing genetic variations as the alternative allele frequency 274 was relatively high, given the evolutionary time, as shown in Table 1. Intriguingly, six 275 variants had FST above 0.5 between EUR and EAS, including rs3827760 in the EDAR 276 gene with the FST of 0.95, rs12632544 close to MRPS22 with the FST of 0.60. 277 Furthermore, three of them affected the glabella facial segment, and two of them 278 affected the nasal region. This result also suggests that local adaption might play a role 279 in forming facial variations between East Asians and Europeans. Among the 13 variants, 280 six were previously reported to be associated with facial shape variation. Well-known 281 facial genes such as EDAR, TBX15 and MRPS22 were associated with craniofacial 282 morphology in many studies 1-4,6,9,13,15,16,38-42 . Our study shows a signal in an intron of 283 TBX15 (rs10923710) contributing to maxillary and tempora shape (seg 19 and 26) in 284 EAS, consistent with the observation that this locus has multiple spatial effects on the 285 face 1 . The variant rs12632544 is an intergenic variant near MRPS22 on chromosome 286 3q23. It was in LD (r 2 = 0.932) with rs12633011, which was reported associated with 287 morphology of eyes in a previous East Asian study 3 . MRPS22 has been reported to be 288 associated with human earlobe size and a mouse skeleton phenotype 16,43,44 . A reanalysis 289 of a GWAS study on cranioskeletal variation in outbred mice revealed that variants in 290 the region of chromosome 9, overlapping with Mrps22, are significantly associated 291 with craniofacial variation (FDR < 5%) ( Supplementary Fig.5). The variants in this 292 region are associated with the protrusion of the maxillary bone region, and shrinkage 293 of the eye and malar bone region. In our study, rs12632544 contributes to EAS-FS in 294 the glabella, eye and tempora (Extended Data Fig. 4b), which is similar to what was 295 seen in the outbred mice. We also identified seven novel variants contributing to EAS-296 FS, six out of which are close to novel genes. Some of these have been reported in the 297 context of craniofacial dysmorphology-for instance, LHX8 and DIS3L2. The variant 298 rs6669519, which contributes to the shape of glabella, is an intergenic variant near 299 LHX8 on chromosome 1p31.1. LHX8, one of the members of the LIM homeobox family 300 of proteins, is a candidate gene for cleft palate 45 . It was also reported to be associated 301 with forebrain neuron development and differentiation 46 . Rs12473319, associated with 302 EAS-FS in glabella (seg 24), is an intronic variant of the DIS3 Like 3'-5' 303 Exoribonuclease 2 (DIS3L2) gene (Extended Data Fig. 4c). DIS3L2 was reported to be 304 associated with a skeleton phenotype in a mouse genome study 44,47 . Besides, DIS3L2 is 305 a candidate gene for the Perlman syndrome, which is presented with craniofacial 306 abnormalities. Moreover, the frequency of the derived C allele is higher in EAS (47.8%) 307 than EUR (2.2%). The estimated allele age of this variant is about 7,875 years old 308 (Table1), when the East Asians and Europeans were already two separated and diverse 309 populations 48 . This suggests that the mutation of this variant happened in East Asians. There are four independent association signals with nose morphology in the SOX9 locus 314 in our EAS GWAS (Extended Data Fig. 7) that are of interest. Although SOX9 is a well-315 known gene contributing to the variation of the nose shape, rs8068343 is a novel 316 independent variant that affects nose shape differences between populations 9,10,13 . In 317 contrast, the other three variants may affect nose shape within populations. Compared 318 with the previously found variants near SOX9, the reference (T) allele of this rs8068343 319 has a higher frequency in EUR (0.96) than in EAS (0.46). Moreover, this variant has a 320 significantly higher FST (0.528) between EUR and EAS than other variants, and the iHS 321 score of this variant is also significantly higher in CHB (2.55). These results indicate 322 this variant may be associated with EUR-EAS nose differentiation.  Fig 4d). The results further prove that nose shape may be under local selection 359 in Europeans rather than East Asians. 360 The analyses above showed that the genetic variants associated with nose shape are 361 more differentiated than expected by random drift. We further analyzed the direction 362 of genetic differentiation. Based on a study of He et al, we firstly estimated and tested 363 differences using the selection coefficients for nose-increasing variants between EUR 364 and EAS 50 . In an enrichment analysis, the EUR population shows higher selection 365 coefficients for nose-increasing variants than the EAS population (P = 4.88×10 -2 ) 366 (Supplementary Table 8, Fig. 4e). Moreover, by comparing the mean genetic prediction 367 of EAS and EUR's facial variation to the expected difference under random drift 368 (Methods), the nose and glabella morphology in EUR is more protruding than EAS and 369 the divergence of the nose is greater than expected under the neutral model (Fig. 4f) 32 . 370 Furthermore, by comparing the PPS using the leading variants with the expected PPS 371 under random drift in the EUR and EAS population, we obtained the direction (and 372 significance) of natural selection on facial morphology in each population 32 . Similar to 373 the population differences, the nose, glabella, and zygoma are under significant natural 374 selection in the European population (Fig. 4g). However, in the East Asian population, 375 the effects of natural selection are not significant (Fig. 4h). These results suggest that 376 facial morphology in European populations undergoes local adaptation, producing a 377 more protruded nose, glabella, and flatter zygoma. 378 Based on the above results, we speculate that the differences between EAS and EUR's 379 facial features are probably due to the adaptive selection that occurred in the European 380 population, which makes Europeans have protruded and narrow noses, significantly 381 different from those of East Asians.

Discussion 383
In summary, as the first large-scale East Asian facial GWAS using a data-driven global-384 to-local phenotyping, our study greatly expanded the knowledge of craniofacial 385 genetics outside the traditionally investigated European populations. Our study also provided many insights into the genetics basis of the facial shape 414 difference between Europeans and East Asians. In addition to identifying 13 primary 415 variants contributing to the European-Asian facial difference, we provided a method to 416 investigate the genetic factors associated with inter-population phenotypic variations. 417 These 13 variants all have a positive and larger than a subtle effect on the EUR-EAS 418 facial difference, shaping the faces of East Asians to be more EAS-FS. Again, 419 corresponding with the PPS results, due to the innate limitation of GWAS, our study 420 overlooks rare or fixed variants in the East Asian population, which are also associated 421 with facial variation, and which also generate more EAS-FS. If the same methods are 422 used in a European population, additional variants affecting EAS-FS can be discovered. 423 Moreover, for those rare or fixed variants with opposite alleles between EUR and EAS, 424 a GWAS on a single population could never identify these, and an admixture population 425 instead is needed. 426 Due to the large number of significant variants associated with the nose region, our 427 natural selection analysis further supported the hypothesis made in the study of Zaidi 428 et al that human nose shape has evolved in response to selection pressures 51 . Again, 429 independent PBS studies using our signals and signals from a European study indicated 430 that the nose shape difference between Europeans and East Asians is mainly due to 431 natural selection in Europeans rather than in East Asians 1 . 432 In conclusion, this study presents the largest East Asian population GWAS on 3D facial 433 shape, using a data-driven global-to-local phenotyping. Our study revealed a large 434 number of novel variants associated with normal range facial shape variation.

Materials and methods 565
Sample and recruitment details 566 The samples in this study were collected from three independent cohorts, the National 567 Since we used two different genotyping platforms in the discovery and replication 602 datasets, we chose to impute the two data sets separately. Prior to imputation, samples 603 with a genotyping missing rate > 0.05 were excluded. Haplotypes were estimated from 604 the genotypes using SHAPEIT2 2 . Then, the samples were imputed to the 1000 605 Genomes Project Phase 3 reference panel using IMPUTE2 3 . After imputation, variants 606 with an INFO score < 0.8 or certainty score < 0.9 were excluded. The discovery and 607 replication datasets were then filtered by variant missingness (--geno 0.05), minor allele 608 frequency (--maf 0.02), and Hardy Weinberg equilibrium (P < 1×10 -6 ). After post-609 imputation quality control, 8,018,212 shared variants between discovery dataset and 610 replication dataset were obtained for analysis. 611

3D image registration and quality control 613
After image acquisition, the 3D images were imported into MeshMonk, a 3D 614 registration software, from wavefront.obj format files to perform a spatially dense 615 surface registration process 4 . 616 The nose tip landmark was manually identified as the point of origin in each 3D image 617 to trigger an initial but crude pose alignment as input to MeshMonk. Afterwards, a 618 symmetric anthropometric mask of 7,906 landmarks was non-rigidly mapped onto all 619 3D images. In each dataset (discovery and replication), generalized Procrustes analysis 620 (GPA) was performed on a stack of the mapped 3D images and their reflection so that 621 any difference in position, orientation, and size of both original and reflected shapes 622 was eliminated. The average of an original and its reflected facial shape resulted in 623 symmetrized facial shape. In this study, we were only interested in symmetrized facial 624 shape, although the asymmetry of the human face is of great interest for future work. 625 After GPA and symmetrization, we investigated every mapped image manually and 626 identified outlier images, typically exposed by locally inconsistent triangles on the 627 surface, which are stretched and compressed irregularly in the images. Outlier images 628 can be caused by poor quality or large noise (i.e., isolated pieces, captured position, 629 clothes) of facial images. We removed four images with poor quality and performed 630 the whole registration process on the noisy images again after removing the noise. In segmentation based on the phenotypic correlation between facial landmarks using the discovery dataset 5,6 . To calculate the phenotypic correlations, we first corrected the 638 symmetrized facial shapes for the covariates of age, age squared, sex, BMI, and four 639 SUGIBS components using a partial least-squares regression (PLSR, function 640 plsregress from MATLAB TM 2018a) to obtain the residual facial shapes 7 . We then 641 calculated the pairwise RV coefficients as the phenotypic 3D correlations between the 642 landmarks of the residual facial shapes. In addition to the RV similarity matrix, a scaled 643 Euclidean distances matrix (the closest points scaled to 1; the most far points scaled to 644 0) between landmarks on the anthropometric mask was calculated. This favors the 645 clustering of points closer together in 3D and therefore avoids isolated points lying 646 outside but clustered with specific facial regions, or, i.e., non-coherent facial regions. 647 Finally, we performed a hierarchical spectral clustering on a combined matrix, as 648 0.9×RV similarity matrix + 0.1×distance matrix, up to level five, resulting in a total of 649 63 facial segments ( Supplementary Fig. 1). 650 In each segment, we calculate each individual's shape by adding his/her residual shape 651 from PLSR to the average shape. For every variant, the meta-analysis described above yielded 63 P-values representing 691 63 facial segments. In the conditional analysis and peak selection, we selected the 692 lowest P-value for each variant. For the initial selection, we selected the variants with 693 P-value below the genome-wide threshold (P = 5×10 -8 ) and calculated the pairwise 694 r 2 between these variants. In each chromosome, we grouped the selected variants 695 consecutively in a way that the r 2 between every two neighbor selected variants in the 696 group is larger than 0.05, which resulted in 230 groups. Then for each group, we 697 performed a conditional analysis to identify the independent signal. In detail, we first 698 selected the variant with the lowest P-value as the conditional variant. Then, we 699 performed association tests of the remaining variants on the condition of the conditional 700 variant. The variant with the most significant P-value still lower than the genome-wide 701 threshold was then added into the list of the conditional variants. We repeated these two 702 steps until the remaining variants were not significant. Finally, we obtained 244 lead 703 SNPs from all groups. To determine the study-wide Bonferroni P-value threshold, we 704 calculated the number of independent tests by both the eigenvalues of the correlation 705 matrix of the segments 13 and the permutation analysis scheme used in the study of 706 White et al 6,14 . The numbers of independent tests obtained from the eigenvalues of the 707 correlation matrix and the permutation analysis are 50 and 45.62, respectively. Here, 708 we used the more stringent threshold 5×10 -8 /50 = 1×10 -9 . 709

Anatomical facial regions 711
Following the data-driven facial segmentation, we manually selected 10 non-712 overlapping segments, which well fit with facial anatomy, composing the whole face. 713 We named these segments as the following anatomical regions: forehead, glabella, eye, 714 tempora, zygoma, nose, maxillary, upper mouth, lower mouth, and mandible. In the ten 715 regions, forehead and mandible are located in the third layer of the 63-segmentation 716 pattern, the remaining eight regions are located in the fifth layer. In further analysis, we 717 accumulated the other segments variants into one of the 10 regions (Supplementary  718   Table 3). 719 720

Gene annotation 721
We used three gene-mapping criteria implemented in Functional Mapping and 722 Annotation (FUMA) to identify the candidate gene for each lead variant 15 . First, we 723 tried to map each variant to genes based on physical distance (within a 10,000 base pair 724 window) from the known protein-coding genes in the human reference assembly. 725 Second, we also included the genes which have a significant eQTL association with the 726 leading variants. We extracted the cis-eQTLs analysis, which uses a linear regression 727 evaluated association between leading variants and expression levels of nearby genes 728 (1 Mb distance to the leading variant), using 10 tissue types from the GTEx v8 729 database 16-18 . We used a false discovery rate (FDR) of 0.05 to define significant eQTL 730 associations. Finally, we also identified candidate genes for each leading variant if there 731 is a 3D DNA-DNA interaction between the leading variant region and the gene region, 732 or in other words, chromatin interaction mapping. To further prioritize candidate genes, 733 we limited interaction-mapped genes to those who interact with a predicted enhancer 734 region identified in any of the 111 tissues or cell types from the Roadmap Epigenomics upstream to 500 bp downstream of the TSS and also predicted by the ROADMAP to 737 be a promoter region) 19 . We expected that the resulting candidate genes are more likely 738 to have a plausible biological function. We used an FDR of 1×10 -6 to define significant 739 interactions. 740 To further narrow down the candidate genes, we investigated whether any gene in the 741 window was previously associated with craniofacial development or morphology 742 through normal-range facial association studies, genetic disorders with facial 743 dysmorphology as a symptom, or animal models. We calculated genome-wide natural selection signatures based on XP-EHH using 823 rehh2 25 . The genome-wide XP-EHH z-scores were standardized through normalization 824 within each derived allele frequency bin (bin widths = 0.01). We estimated two-tailed 825 P-values of the variant according to the normalized z-scores. We calculated the FST and 826 estimated from the 1000GP Phase 3. EUR (n = 404, EUR: TSI (107), GBR (91), IBS 828 (107), and CEU (99)), EAS (n = 208, EAS: CHB (103) and CHS (105)) and YRI (n = 829 103) individuals were used. Using these frequencies, we estimated pairwise FST 830 between groups for the PBS test. The estimator is calculated as follows: 831