Introduction

The face is a highly visible trait system that has a pervasive role in human life, including social interaction, medical diagnosis1 and potentially forensic issues.2 A genetic basis of facial morphology is evident from heritability studies3 and from anecdotal evidence about familial resemblance (most notably in monozygotic twins) that has a role in social interaction as well as clinical genetics.4 Links between pathological and normal variation are plausible, because many syndromes show a clear connection between genetic alteration and typical facial gestalt,4 implying that genes involved in affected individuals might also contribute to normal facial variation. However, such a link has not been established as yet.

With a birth prevalence of around 1.7 in 1000 lifebirths, non-syndromic cleft lip with or without cleft palate (NSCL/P) is one of the most common congenital disorders worldwide,5 and is defined by a pathological facial trait. Risk to NSCL/P has long been conjectured to be connected to normal facial variation; earlier observations have established deviating facial measurements in relatives of NSCL/P patients as compared with control populations.6, 7 A recent meta-analysis highlighted nose width and facial width as most prominently deviating traits.8 Deviations were found in width measurements of the face, which is in line with the mechanism of NSCL/P formation in embryonic development being characterized by failure of palatal fusion, or narrowing of lateral facial structures.9, 10 These findings suggest that relatives with phenotypic facial deviations share alleles with NSCL/P patients, implying variable or modified penetrance of identical alleles spanning the phenotypic spectrum of normal facial variation to NSCL/P.

NSCL/P has a multifactorial etiology encompassing both genetic and environmental components. Candidate-gene studies have reported a number of genes associated with NSCL/P. However, to date, only the IRF6 gene has shown a convincing degree of consistency across studies with respect to association with NSCL/P.11 Recent genome-wide association studies (GWASs) have identified five chromosomal regions showing replicated associations with NSCL/P as confirmed genetic risk factors.12, 13, 14

In the present study, we empirically tested the hypothesis that genetic variants associated with pathological NSCL/P phenotype are also involved in determining normal facial morphology by investigating the six replicated genetic risk factors for NSCL/P, together with three additional suggestive loci from recent GWASs12, 13, 14 for any association with variation of facial and nasal width in two samples (Essen, Germany and Rotterdam, The Netherlands) from the general population. Additionally, we assessed the predictive power of these genetic markers on normal facial morphology, which is relevant for practical applications of established knowledge, for example, in future forensic studies, where estimating facial appearance of unknown persons from DNA left behind is becoming important for finding these persons.2

Material and methods

Samples and phenotype extraction, Essen, Germany

We have recruited a random sample of 529 individuals of German European descent aged between 18 and 35 years in the Essen area. Two-dimensional digital portraits of all persons were made using a standardized procedure described elsewhere.1 The sample consisted of 106 sibships and 346 singletons, and the female to male ratio was 2.3:1. All individuals provided written informed consent and the local ethics committee approved the study.

Graphs were superimposed on images by the bunch graph method.1 These graphs were used to standardize width measurements relative to facial image area.15 Graphs were readjusted manually before standardization. The nose width is defined as the distance between the left and right alar landmarks, and for the facial width phenotype, the distance between the left and right zygoma landmarks was used, in the following called bizygomatic distance (Figure 1a).

Figure 1
figure 1

The bizygomatic distance and the nose width as defined in the two-dimensional photos (a) and the three-dimensional MRI (b).

Full size image

Samples and phenotype extraction, Rotterdam, The Netherlands

Furthermore, 2497 subjects of Dutch European descent aged ≥45 were selected from the population-based Rotterdam Study.16 These subjects were scanned on a 1.5 T General Electric MRI unit (GE Healthcare, Milwaukee, WI, USA), using an imaging protocol including a three-dimensional T1-weighted sequence with an interpolated voxel size of 0.49 × 0.49 × 0.80 mm3. More details on image acquisition can be found elsewhere.17 The Rotterdam Study has been approved by the Medical Ethics Committee of the Erasmus University Medical Center, Rotterdam, The Netherlands, and all participants provided written informed consent. The nose width and bizygomatic distance phenotypes were extracted from these magnetic resonance images (MRI) with an automated technique, which can transfer predefined landmarks from a limited set of annotated images to an unmarked image.

The method operates in three steps. In the first step, the left/right zygoma and alare were marked in the images of 18 subjects by a trained observer (we shall refer to these images as the atlases). To assess intra-observer error, the measurements were repeated approximately 6 months later. In the MRI, the zygoma is defined as the most lateral point of the zygomatic arch of the skull when the head is aligned to the Frankfurt horizontal plane. The alare is defined as the most lateral point of the nose. In the second step, each unmarked image (which we shall call the target) was non-rigidly registered to all atlas images. Image registration is the process of determining a coordinate transformation, which relates any position in the atlas image to the anatomically corresponding position in the target image.18 The registrations were performed with the Elastix software (http://elastix.isi.uu.nl/)19 by first estimating an affine transformation, followed by a non-rigid transformation. Finally, the zygoma and alare positions in the 18 atlas images were transformed to the target image coordinate system, resulting in 18 estimates of each landmark's position. These estimates were then combined by finding the geometric median. The nose width and bizygomatic distance were computed based on these median alare and zygoma landmarks, respectively. Visual inspection revealed 935 scans (37%) in which the alare was not visible. The associations involving nose width were therefore based on 1562 subjects. Missing the nose in numerous individual MRI scans is explained by the scanning procedure targeting the brain for other study purposes. Although bias due to informative missingness cannot be entirely excluded, it is unlikely as measurements were standardized by head size and a strong correlation between factors leading to missingness (ie, field of measurement) and nose width would be a prerequisite. An example of a typical facial image as reconstructed from the MRI data is shown in Figure 1b, with nasal width and bizygomatic distance indicated.

The accuracy of width estimation was evaluated by leave-one-out experiments on the 18 atlas images. In each atlas image, the four landmark positions were estimated using the other 17 atlases. Resulting width estimates were then compared with the distance between the manually positioned landmarks. For the bizygomatic distance, a mean error of 0.9±2.2 mm (relative: 0.7±1.8 %) was measured, and for the nose width −0.2±0.9 mm (−0.7±2.4%). As a reference, we also estimated the mean intra-observer errors of the same distances. These were 0.3±0.5 mm (0.2±0.4%) for the bizygomatic distance, and 0.7±0.4 mm (1.9±1.0%) for the nose width.

In order to normalize the two measures for differences in global head size, we obtained the global volume change of the affine transformation, which relates each target image to an average atlas constructed from the 18 atlas images.17 As we are only interested in the volume change, it was not necessary to derive this transformation explicitly. Instead, we first determined the separate volume changes of the 18 affine transformations from the target to each atlas by computing the determinant of the transformation matrix. We then computed the average volume change by calculating the geometric mean and inverted it.

SNP ascertainment and genotyping

The initial single nucleotide polymorphism (SNP) ascertainment was based on the most comprehensive study of the IRF6 locus11 and the GWASs performed in Central Europeans.12, 13 For the IRF6 locus (1q32), the functional variant rs642961 was chosen; for the three replicated loci at 8q24, 10q25 and 17q22, the most strongly associated SNPs (rs987525, rs7078160, rs227731) were chosen; and for the 8q24 and the 17q22 loci, two additional variants showing independent effects were selected (rs16903544, rs17760296).12, 13 To cover the five most associated loci from both studies,13 two more SNPs (rs9574565, rs1258763) from 13q31 and 15q13 were included. During the preparation of this manuscript, a GWAS in samples of European and Asian ancestry reported two additional genome-wide significant loci at 1p21 and 20q12.14 These loci are represented by rs560426 and rs13041247, which were also included in our study; however, genotypes for these two markers were only available in the Rotterdam sample.

The Sequenom MALDI-TOF mass-spectrometer (MassArray system, Sequenom Inc., San Diego, CA, USA) was used for genotyping in the Essen sample, and the data were analyzed using the Spectrodesigner Software package (Sequenom, San Diego, CA, USA). Primer sequences are available upon request. Three distinct clusters were analyzed with the Typer Analyzer 4.0.1 software (Sequenom). For the Rotterdam samples, SNP genotypes were taken from a previously established genome-wide SNP dataset in the Rotterdam Study, using the Human 610 Quad Arrays (Illumina Inc., San Diego, CA, USA).20 All SNPs in both studies were in Hardy–Weinberg-equilibrium (Supplementary Table S1).

Statistical analysis

The two horizontal distances as depicted in Figure 1 were analyzed for genetic association. Variances of these phenotypes were standardized to 1 before further analysis. We used a standard linear model to evaluate genetic association. An additive genetic model by encoding genotypes as 0, 1 and 2 (counting minor alleles), as well as a genotypic model with dummy variables for genotypes were analyzed. Normality of phenotypes was confirmed by visual inspection of the residual distribution. In the Essen sample, we used a generalized estimation equation approach to account for familial dependencies. Sibs were treated as exchangeable within families. P-values in the Essen sample are based on 2 × 104 permutations, because some genotype frequencies were less than 10%. Permutations were performed within families and between families of identical structure. We also stratified analyses by sex to detect possible sex-specific effects. P-values of both studies were combined by a random effects meta-analysis.

The prediction model was constructed using a penalized regression model with L1-penalty (LASSO).19 The tuning parameter λ was optimized using a 40-fold cross validation based on mean-squared error minimization. For families, the mean-squared error was averaged within families and then across families. The regression model included all SNPs in the study, all of which were either coded for the additive or genotype model. Except for the discussion section, all reported P-values are nominal. All analyses were conducted using software package R version 2.10.1 (http://www.r-project.org/).20 The generalized estimation equations were fitted using package geepack. Penalized regression was conducted using package penalized, and meta-analysis made use of package meta.

Results

Association results

Results for an additive model of association with the traits nose width and bizygomatic distance for the 11 SNPs tested are given in Table 1. The most prominent result was an association between rs1258763 at locus 15q13 with nose width in the Essen sample (P=6 × 10−4). Although association at this marker was not statistically significant in the Rotterdam sample (P=0.2), effects pointed in the same direction (−0.24 and −0.13 for Essen and Rotterdam, respectively) and the combined P-value was statistically significant (P=0.0005). Consistent with this observation, the G-allele at rs1258763 was over-represented among controls as compared with NSCL/P patients in the previous GWAS.13 Suggestive evidence for nose width determination was observed for SNPs rs17760296 and rs227731 (17q22) in the Essen sample (P<0.1), but effect size directions were different in the Rotterdam study data.

Table 1 Association tests for width measurements of the face (phenotype) with an additive model
Full size table

For the bizygomatic distance phenotype, only the Rotterdam sample showed statistically significant associations that were strongest for rs987525 (8q24, P=0.0171), which was supported by a same direction, albeit small effect in the Essen sample; however, the combined analysis was not significant (P=0.26). The minor risk allele A is associated with a narrowing of the bizygomatic distance. Interestingly, this SNP showed the largest effect size in one of the previous GWASs.12 A second significance for bizygomatic distance was seen for rs7590268 in the Rotterdam sample (2p21, P=0.035), where the minor allele leads to reduction in bizygomatic distance. Furthermore, rs17760296 at 17q22 showed borderline significance for bizygomatic distance in the Rotterdam sample (P=0.0504). Although it provided suggestive evidence for association in the Essen sample (P=0.0883), the effect directions were not agreeing. Discrepancies between studies and effect sizes may have their origin in aspects of the phenotype collection as discussed below.

The nose width effect of rs1258763 (15q13) was best explained by a recessive model of the major allele, whereas rs987525 in CCDC26 seemed to be best explained by an additive model (Supplementary Table S2). All analyses included age and sex as covariates. As sex-specific effects were suggested from previous findings,7 we also conducted a stratified analysis according to sex. The rs1258763 (15q13) association appeared stronger in males as compared with females (Table 2) in both cohorts. This finding concurs with measurements in unaffected relatives of NSCL/P patients, showing a larger deviation in nose width in males compared with females.7 For rs987525 (8q24), females exhibited a stronger effect than males in the Rotterdam sample, whereas the male effect in Essen was close to zero. Full information on effect sizes and corresponding P-values for the additive and genotype models are given in Supplementary Tables S3 and S4.

Table 2 Sex specific effects for the additive model
Full size table

Prediction results

All SNPs were entered into a penalized regression analysis for, both, additive and genotype models (Table 3). These SNPs explained 2.04 and 0.57% of the age- and sex-adjusted variance of nose width in the Essen cohort under the additive and genotypic models, respectively. Bizygomatic distance phenotype had 0.28% of the variance explained in the Rotterdam sample, under the additive model. All other models showed very little explanatory power. Coefficients of the prediction models are given in Supplementary Table S5. The nose width phenotype was mainly explained by SNP rs1258763 (15q13) and, albeit to some lesser extent, by SNPs rs7590268 (2p21) and rs17760296 (17q22). The bizygomatic distance prediction model did not show as strong a concentration of explanatory power on a few SNPs. SNPs with the largest coefficients included rs17760296 (17q22), rs987525 (8q24), rs7590268 (2p21) and rs7078160 (10q25) with absolute effect sizes ranging between 0.23 and 0.32. All SNPs had non-zero coefficients for bizygomatic distance phenotype.

Table 3 Explained variance for prediction models of the phenotypes nose width and bizygomatic distance in the Rotterdam Study and Essen samples
Full size table

Discussion

In this study, we investigated a hypothesized link between risk alleles for NSCL/P and variation in normal facial morphology in the general population. Such a connection is obvious for traits where disease definition depends on a threshold in quantitative traits like blood pressure.21 For NSCL/P, this connection has been documented through measurements of facial geometry in unaffected relatives of NSCL/P patients.7, 8, 22 Population-based samples that would allow to verify this hypothesis were identified and approached with the idea of a study regarding a restricted number of SNPs. Our decision on which facial phenotypes to study was based on two more recent studies. A meta-analysis showed nasal cavity and maximum head width as the most distinguishing features.8 Whereas nasal cavity width was found increased, the most prominent feature in a later study was nasal width reduction in males.7 A recent study found the naso-labial region to be ‘pinched’ in males.22 These findings and the fact that they are both measurable in two-dimensional facial pictures and three-dimensional head MRI data led us to focus on nose width and bizygomatic distance. However, the two phenotype data sets were not completely comparable to each other based on differences in the initial data sets, but also in ways of phenotype extractions. The bizygomatic distance phenotype was well defined for the head MRI data in the Rotterdam Study sample; however, for two-dimensional facial pictures in the Essen sample, the bizygomatic distance was more indirectly defined by neighboring landmarks of the face.15 Differences in phenotype definition are one possible explanation for the differing findings in these two cohorts. The most significant finding, SNP rs1258763 (15q13), showed the same effect size in both studies, corroborating the expectation for the nose width phenotype to be measured very similarly in both cohorts. Another influential factor differing between cohorts was age. Rotterdam Study participants were substantially older on average at time of MRI, compared with the Essen subjects at time of photography (Figure 2). To the authors’ knowledge, there are no longitudinal studies investigating age effects on facial traits and there were no other studies establishing a quantitative relationship (perhaps with the exception of ear lobe crest association with coronary heart disease).23 Anecdotal knowledge, however, suggests a strong influence of age on facial traits in general, although it is not clear to what degree this would apply to the two distance measures used here. Age may have a potentially confounding effect on the genetic association of facial characteristics, which may explain some of the discrepant findings between the Essen and Rotterdam Study data. Alternatively, these discrepant findings may reflect false positives due to multiple testing. Figure 2 shows that the age distribution is substantially different between the study population samples. Supplementary Table S6 summarizes the significance levels of covariates for both cohorts over all statistical tests performed. On average, sex was not a significant covariate in the Essen cohort, but was significant in the Rotterdam Study data (P=10−3). The effect of age was marginally significant in Essen (P=0.048), whereas in the Rotterdam Study sample, age was a highly significant covariate (P=10−253). In the prediction model, age and sex explained ∼10% of phenotypic variation in the Rotterdam Study, compared with ∼2% in the Essen sample. Most likely, the facial phenotypic trajectory with age acts as a strong confounder. However, potentially, markers could also influence facial traits in an age specific manner. To conclude, age matching seems to be desirable in studies on facial traits.

Figure 2
figure 2

Age distribution of probands in the Essen and Rotterdam Study sample.

Full size image

In the initial GWASs,12, 13 all six SNPs showed allele frequencies between 0.2 and 0.5 for risk alleles. These findings point to a low-penetrance disease mechanism in which effects on normal variation could follow a different mechanism. NSCL/P patients show broader mid-facial measurements in general,24 whereas unaffected relatives tend to show reduced width for some measurements, for example, nose width.7, 8, 22 Our findings for the nose width phenotype support a widening effect of the NSCL/P risk allele, and therefore, contradict a simple disease model where the allele action in the normal population could be equated with that in unaffected relatives. However, the sex specificity of findings in unaffected relatives (major male effect) is confirmed. Taking into account that unaffected relatives constitute a selected subgroup with a potentially specific subset of modifiers/environmental exposure, our findings support a complex genotype/phenotype relationship for the two facial traits investigated here.

In this study, we focused on testing a specific biological hypothesis and have therefore reported nominal P-values. Nevertheless, multiplicity of testing had to be accounted for. Nine to eleven SNPs were evaluated in these two population samples using two phenotypes. Hence, taking a conservative stance and correcting for 44 tests in individual samples and 22 tests for combined results, the nose width finding for SNP rs1258763 (15q31) was significant after a strict Bonferroni correction, with a P-value of 0.026 in the Essen sample and 0.01 in the combined analysis. Considering nose width as the primary phenotype resulted in P-values of 0.013 (Essen) and 0.005 (combined). However, associations for SNPs rs7590268 (2p21), rs987525 (8q24) and rs17760296 (17q22) with the bizygomatic distance phenotype in the Rotterdam data were no longer significant after Bonferroni correction. Although no SNP was statistically significant in both studies simultaneously and thereby replicated, we believe that our results, at least for the nose width phenotype, warrant the interpretation of a genuine association.

For both cephalometric and anthropometric measures, high heritability has been reported,3 suggesting a strong genetic determination of facial traits. We therefore explored the explanatory power of our genetic markers to predict the nose width and the bizygomatic distance phenotypes from DNA variants. Genetic effects were strongly superimposed by age and sex effects, especially in the Rotterdam sample. Nevertheless, some variation could be explained by the genetic markers, namely 2% for the nose width phenotype in the Essen sample and 0.57% for the bizygomatic distance phenotype in the Rotterdam Study sample. Although these numbers appear low, they were driven by few signals (eg, 17q22, 15q13, 2p21, nose phenotype Essen, additive model). When compared with body height, for which a recent study found ca. 3% of the variance explained by ca. 20 associated SNPs,25 these numbers seem comparable. We speculate that the human face is determined by a large number of genes whose effects are spread over multiple subtraits, making it potentially easier to dissect genetic contributions to facial morphology as compared with body height.

In conclusion, we have demonstrated that association with one marker could explain ca. 2% of nose width variation, and a tentative association between bizygomatic distance and other markers could account for about 0.5% of variation. Finally, our study represents the first approach to understanding genetic control of facial morphology, demonstrating that predicting facial distance traits from genetic markers is not nearly as straightforward as it is for human eye20 and hair color,26 and that further genetic research will be needed to identify predictive genetic markers, which could achieve the accuracy needed for practical applications such as future forensics.