Neanderthals and humans interbred in Europe for 5,400 years

Automatic landmarking identifies new loci associated with face morphology and implicates Neanderthal introgression in human nasal shape

Abstract
We report a genome-wide association study of facial features in >6000 Latin Americans based on automatic landmarking of 2D portraits and testing for association with inter-landmark distances. We detected significant associations (P-value <5 × 10−8) at 42 genome regions, nine of which have been previously reported. In follow-up analyses, 26 of the 33 novel regions replicate in East Asians, Europeans, or Africans, and one mouse homologous region influences craniofacial morphology in mice. The novel region in 1q32.3 shows introgression from Neanderthals and we find that the introgressed tract increases nasal height (consistent with the differentiation between Neanderthals and modern humans). Novel regions include candidate genes and genome regulatory elements previously implicated in craniofacial development, and show preferential transcription in cranial neural crest cells. The automated approach used here should simplify the collection of large study samples from across the world, facilitating a cosmopolitan characterization of the genetics of facial features.

Introduction
Genome-wide association studies (GWAS) of human facial features are contributing importantly to elucidating the genetic basis of variation in facial features in the general population1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17. The genomic regions identified often overlap developmental genes, and have been shown to be enriched in regulatory elements active during craniofacial development4,17. These studies were initially performed in individuals of European descent7,9,10,11. Face GWASs are being gradually extended to the characterization populations of non-European ancestry1,2,3,5,6,16. The increasing characterization of non-Europeans is helping to draw a fuller picture of the genetic architecture of facial variation in humans and further our understanding of the evolution of human facial features.

GWASs of facial variation have used a range of phenotyping approaches, from qualitative assessment of morphological features on 2D photographs1, to measurements based on manual landmarking of 2D photographs2, to semi-automatic analyses of 3D facial images4,15. These approaches vary greatly in cost, informativity, and ease of application. Although 3D images fully represent facial morphology, acquisition of such data requires specialized equipment, complicating their widespread application across the world. Although less informative than 3D imaging, standard 2D photographs have the potential to facilitate the collection of large, diverse study samples. However, manual landmarking of 2D photographs is a slow, labor-intensive task. This has fostered an interest in the application of fully automatic landmarking approaches. However, so far these have enjoyed limited success18,19,20. Most studies based on 2D photographs have therefore been based on entirely manual1 or, at times, semi-automatic landmarking9 (i.e. combining automatic landmarking with manual editing).

Here we report a GWAS of facial features derived from a fully automatic landmarking of 2D frontal photographs from Latin Americans of mixed European, Native American and African ancestry. The association signals detected overlap with previous GWAS findings. In addition, we identify 33 novel signals. For most of the novel signals identified, we find evidence of statistical replication in European, East Asian, or African GWAS data, and one mouse homologous region influences craniofacial morphology in mice. One of the novel regions identified includes a tract of introgression from Neanderthals, which we associate with an increase in nasal height, consistent with the morphological differentiation between Neanderthals and modern humans.

www.nature.com/articles/s42003-023-04838-7

Results
Study sample and phenotyping
The 6486 individuals examined here are part of the CANDELA cohort, collected in five Latin American countries21. This cohort has been previously studied in GWASs of various physical appearance traits1,2,22,23,24. This includes two previous facial morphology GWASs based on 2D photographs: one mainly based on categorical (i.e. morphoscopic) phenotyping, and one based on manual landmarking of lateral (profile) photographs1,2,24. Individuals included in these studies were genotyped on Illumina’s OmniExpress chip (including >700,000 SNPs) and characterized for a set of standard covariates (age, sex, BMI, and genetic ancestry estimated from the chip data)1,2,24.

Here we used the Face++ cloud service platform (https://www.faceplusplus.com) to automatically place 106 landmarks on frontal 2D photographs (i.e. portraits) from the CANDELA individuals (Supplementary Fig. 1). Previously, 16 of these landmarks had been placed manually on a small subset of these individuals1 and we used these data to evaluate the robustness of the Face++ landmarking. We also compared Face++ with Dlib, a popular landmarking tool25,26 (Supplementary Table 1). We calculated Interclass Correlation Coefficients (ICCs) and median Euclidean distances between landmarks placed manually, by Face++, or by Dlib. According to both metrics, the landmarks placed by Face++ were very close to the manual landmarks, and the performance of Face++ was superior to Dlib for certain landmarks (Supplementary Table 1).

After Procrustes superposition, we calculated inter-landmark distances (ILDs) between 34 Face++ landmarks (mostly corresponding to well-defined anatomical landmarks, Fig. 1a and Supplementary Table 2)1,8,10,17,27,28,29. Accounting for face symmetry, we obtained 301 distances. Some of the landmarks retained are on the eyebrow edges, making distances based on them sensitive to eyebrow size (Fig. 1b). The distances obtained show considerable variation and are approximately normally distributed (Supplementary Fig. 2). Many distances show a significant correlation with three head angles estimated by Face++ (pitch, roll, and yaw angle), reflecting the effect of head pose (Supplementary Table 3, Supplementary Fig. 3). Consequently, we excluded 76 individuals with extreme head angle values, and included these angles as covariates in the genetic association tests.

Fig. 1: Overview of the facial features GWAS performed here.



a Dots indicate the location of the 34 facial landmarks used for calculation of 301 inter-landmark distances (ILDs, Supplementary Table 2 provides additional information on these landmarks). b Lines represent the 148 ILDs for which we detected significant association with at least one genomic region in the CANDELA data. Black lines refer to ILDs that resulted in replication of previously reported associations (Table 1). Red lines represent ILDs that revealed novel associations (with a darker red highlighting the ILDs associated with five genomic regions discussed in the text; Table 2). c Combined Manhattan plot illustrating all significant GWAS hits (-log(P) > 7.3, red line). Black labels indicate candidate genes from previous GWASs that were replicated here (Table 1). Red labels highlight the main candidate genes at the five novel regions discussed in the text. For visibility, the y-axis has been truncated at a -log(P) of 13.

Follow-up of novel face regions in the mouse
To evaluate the potential effect in the mouse of the face regions newly identified here, we reanalyzed published genome-wide SNP data from outbred mice characterized for craniofacial shape variation34. Of the 33 novel regions identified here, 30 could be successfully mapped onto the mouse genome (Supplementary Table 5). Of these, a region on mouse chromosome 5q (homologous to human 22q12.1) showed significant association for SNPs over a ~1.5 Mb segment, with the index SNP in this region (rs32069343, P-value: 2 × 10−34), impacting on multiple aspects of mouse skull and mandible shape (Fig. 2, Supplementary Movie). In the CANDELA GWAS, SNPs in 22q12.1 are associated with ILD D437 between landmarks 3 and 31 (Fig. 1), a distance sensitive to the height of the lower face (smallest P-value of 1.8 × 10−8 for rs9608473, Fig. 2). In previous studies, SNPs in 22q12.1 have been strongly associated with height35, and suggestively associated with facial features36,37 and cleft lip/palate38.

Fig. 2: Regional association plots in mouse and human.




a At the top are shown P-values for SNPs in the mouse chromosome 5 region homologous to human 22q12.1 . Dot colors reflect LD with the index SNP (rs32069343). Underneath are shown genes annotated in the region. b Effect of rs32069343 on mouse skull shape (expansion/contraction relative to the mean shape is shown in blue/brown). To facilitate visualization, the effect has been magnified ×20 (see also Supplementary Movie). c LocusZoom plot for 22q12.1: the top panel shows SNP P-values (dot color reflecting LD with the index SNP rs9608473). The bottom panel displays genes annotated in the region. d, e Facial landmarks placed by Face++ (dark dots indicate the landmarks retained in ILD calculation). Lines indicate ILDs associated with SNPs in 22q12.1. In d line color reflects association P-value; in e line color reflects the direction and magnitude of the association. The ILD with strongest association P-value (D437) is labelled.