Genetic and Environmental Contributions to Facial Morphological Variation: A 3D Population-Based Twin Study

Introduction Facial phenotype is influenced by genes and environment; however, little is known about their relative contributions to normal facial morphology. The aim of this study was to assess the relative genetic and environmental contributions to facial morphological variation using a three-dimensional (3D) population-based approach and the classical twin study design. Materials and Methods 3D facial images of 1380 female twins from the TwinsUK Registry database were used. All faces were landmarked, by manually placing 37 landmark points, and Procrustes registered. Three groups of traits were extracted and analysed: 19 principal components (uPC) and 23 principal components (sPC), derived from the unscaled and scaled landmark configurations respectively, and 1275 linear distances measured between 51 landmarks (37 manually identified and 14 automatically calculated). The intraclass correlation coefficients, rMZ and rDZ, broad-sense heritability (h2), common (c2) and unique (e2) environment contributions were calculated for all traits for the monozygotic (MZ) and dizygotic (DZ) twins. Results Heritability of 13 uPC and 17 sPC reached statistical significance, with h2 ranging from 38.8% to 78.5% in the former and 30.5% to 84.8% in the latter group. Also, 1222 distances showed evidence of genetic control. Common environment contributed to one PC in both groups and 53 linear distances (4.3%). Unique environment contributed to 17 uPC and 20 sPC and 1245 distances. Conclusions Genetic factors can explain more than 70% of the phenotypic facial variation in facial size, nose (width, prominence and height), lips prominence and inter-ocular distance. A few traits have shown potential dominant genetic influence: the prominence and height of the nose, the lower lip prominence in relation to the chin and upper lip philtrum length. Environmental contribution to facial variation seems to be the greatest for the mandibular ramus height and horizontal facial asymmetry.

the lower lip prominence in relation to the chin and upper lip philtrum length. Environmental contribution to facial variation seems to be the greatest for the mandibular ramus height and horizontal facial asymmetry.

Introduction
Facial morphology attracts interest from a variety of scientific disciplines dealing with craniofacial evolution (anthropology) [1], forensic facial reconstruction and identification (forensic science) [2], facial recognition (computer science) [3], predictions of facial ageing, growth and development (clinical orthodontics) [4], and the perception of facial features in social interactions (sociology and psychology) [5]. The common goal of all these disciplines has been a greater understanding of the factors that shape facial morphology in various populations. Human Genome and Phenome projects revolutionized the way scientists have looked at this topic. For a long time, dysmorphic facial characteristics have been the focus of clinical and genetic studies, but the same could not have been said about normal facial variation until recently [6].
Facial morphology is influenced by genetic and environmental factors and their complex interactions that present the greatest challenge in modern biology [7,8]. Traditionally, exploration of genetic variance in humans has been limited mostly to the use of additive effects estimated using pedigree data [9]. The role of genetics in complex traits has been quantified as heritability, e.g. the proportion of the total phenotypic variance explained by additive genetic variance. However, genetic variance of complex traits also includes non-additive variances (dominance and epistasis) [9,10]. With the development of dense panels of single nucleotide polymorphisms (SNPs), the exploration of genetic variation of complex traits is moving to the dissection of genetic variation at individual loci [11,12].
Various environmental influences are also responsible for specific craniofacial shape, including hormones, nutrition, diseases, trauma, surgery, dentofacial orthopaedics, lifestyle factors (e.g. smoking, alcohol, exercise) and oral functions (mastication, respiration, swallowing) [13]. Environmental factors can influence epigenetic mechanisms (acquired and heritable changes in gene function that occur without a change in the DNA sequence) [14]. This can be accomplished through DNA modification (i.e., methylation), histone modification and posttranscriptional silencing by RNA interference. Thus, transcriptional activity of specific genes is controlled at specific points in time in specific organs [15].
Heritability of craniofacial characteristics has been investigated in family (siblings or parents-offspring design) and twin studies (monozygotic, MZ vs. dizygotic twins, DZ) mainly using two-dimensional techniques (anthropometry, photography and lateral skull radiography) [16][17][18][19][20][21][22][23]. The results of these studies are limited mainly due to: 1) a significant loss of information when faces are studied in two dimensions only; 2) risks associated with radiation hazard, 3) issues arising in determining zygosity (it was not always confirmed formally by genetic testing), 4) relatively small samples and a loss of statistical power. Recent three-dimensional (3D) techniques (such as laser scanning, photogrammetry and magnetic resonance imaging) have enabled non-invasive data collection and a more comprehensive description of facial features [24][25][26][27].
Over the last decade, a few studies have investigated 3D facial characteristics of twins [28][29][30][31]. Using facial colour maps of twins (depicting the deviation between faces), Naini and Moss [29] found that the triangular area encompassing the orbital rims, intercanthal area, and the nose is genetically driven in British twins; MZ twins showed similarities in the shape of the eyebrows, bridge of the nose and infraorbital ridges, whereas other parts of the face such as cheeks, chin, and lips showed significant variation [31]. Similar findings were reported in the preliminary study performed on American twins [32]. These 3D studies dealt with small, convenient samples that were not representative of the respective populations. In addition, there was no robust estimation of heritability. Therefore, our current knowledge on heritability of craniofacial traits mainly stems from two-dimensional studies with some conflicting results on heritability estimates of horizontal and vertical skeletal parameters.
The aim of the present study was to assess genetic and environmental contributions to 3D facial soft tissue morphological variation within the framework of the classical twin study design [33]. The objectives were to: 1) evaluate additive and non-additive genetic effects and 2) evaluate the influence of common and unique environment on facial morphological traits of British adult female twins.

Materials and methods Sample
The sample comes from the UK Adult Twin Registry (Twins UK, www.twinsuk.ac.uk), a population-based study of twins with ongoing longitudinal data collection [34]. At the time the study was conducted, the Registry database contained 3D facial images of 1521 twins aged 16.5 to 87.3 years (98.6% females), obtained using a non-invasive stereophotogrammetry device (3dMDface system, 3dMD Inc., Atlanta, GA, USA) proven to be valid and reliable [35][36][37]. Facial images of 1380 female participants aged between 23.6 and 86.5 years (mean 58.8, SD 10.5) were used in this study. There were 263 pairs of MZ twins, 341 pairs of DZ twins and 172 unpaired twins (75 MZ and 97 DZ): 31 individuals whose facial image was not taken and 141 individuals who were excluded for one of the following reasons: 1) zygosity was initially assessed by a 'peas in a pod' questionnaire but not confirmed via subsequent genotyping or genome-wide association studies (n = 19); 2) male participants (n = 15); 3) image quality was not good enough to identify 37 facial landmarks (n = 99); 6) facial expression was not neutral (n = 3); 7) head posture was not natural (n = 3); 8) mouth open and/or eyes closed (n = 2). Ethical approval for the study was obtained from the Health Research Authority UK (more specifically from its National Research Ethics Service Committee in London, Ref: EC04/015). Written consent was obtained from all participants.

Data collection
During the taking of 3D facial images, the participants were instructed to do the following: sit up straight in front of the device and look at the distance; to keep the head in the natural head posture; to maintain neutral facial expression with mouth closed, lips at rest (if competent) and eyes open; and to be still. In addition, they were asked to take off their glasses and tie up their hair. Image acquisition took approximately 3 minutes per person (1.5 ms shot followed by automated transfer, processing and triangulation of the data). The quality of the images was checked using the proprietary software 3dMD Patient. The images were saved as TSB files (a 3dMD Patient format) and converted to OBJ format for further analysis with Rapidform 2006 (INUS Technology Inc., Seoul, South Korea).
Processing, normalising and landmarking of facial images was performed using in-house developed subroutines for Rapidform 2006 [38]. First, the background objects (if any were captured), pieces of clothes, hair and most of the neck were removed using a freehand cutting tool. The area of the neck approximately 2 cm under the mandibular body and ramus was left to enable correct identification of two landmarks: 'gnathion' and 'gonion'. Secondly, the faces were normalized by fitting them into the same reference frame. The origin of this reference system (with the coordinates x = 0, y = 0, z = 0) was the point halfway between the inner canthi of the eyes. Coronal plane (xy) was determined by the cylinder that fitted all data points of structure made from the original face and its mirror reflection. Sagittal plane (yz) was identified as the symmetry plane of that structure. Transverse plane (xz) was perpendicular to the previous two planes and connected the inner canthi of the eyes [38]. The normalization procedure was performed to facilitate landmarking and construct average faces [27].
Stereophotogrammetry enables capturing of thousands of points across the facial surface that contribute to facial variation [24][25][26][27]. In this study, the analyses were limited to landmark configurations and linear distances based on 37 anthropometric landmarks (14 bilateral and 9 mid-sagittal). The definitions of eight mid-sagittal and 10 bilateral landmarks were identical to a widely-known anthropometric definitions by Farkas [39]. One mid-sagittal and four bilateral landmarks were newly defined by the authors of the study (Fig 1, Table 1). All landmarks were manually identified on all the images by one operator. For every participant, 111 x, y and z coordinates were saved.

Facial parameters
The first step in obtaining facial parameters for the analyses was the generalized Procrustes registration of landmark configurations. It was performed on the sample of 1380 participants (263 pairs of MZ twins, 341 pairs of DZ twins and 172 unpaired twins) using the R software and programming language [40]. Since both facial shape and size variations were considered important for the study, two registrations were done: without scaling (taking into consideration both facial size and shape) and with scaling (preserving only information about facial shape) [41].
The second step was to import the Procrustes registered coordinates into SPSS (v.20.0, SPSS Inc, Chicago, Illinois, USA) and perform two Principal Component Analyses: the first one was done on the unscaled dataset and the results (principal components) obtained using the covariance matrix and the second one was done on the scaled dataset and the results based on the correlation matrix [42][43][44]. The principal components were further validated by replicating the results on two subsamples consisting of 690 individuals.
The third step was to calculate linear distances between manually identified 37 landmarks. Additional distances were included by automatically calculating the position of 14 mid-points of bilateral landmarks. These 51 landmarks determined 1275 linear distances (51x50/2). However, only 1247 of them show different results. The reason is that the 14 pairs of bilateral landmarks and their respective mid-points produce 42 distances, with only 14 of them being independent. For instance, distances 'alare left'-'alare right', 'alare left'-'mid-alare', and 'alare right'-'mid-alare' are not independent.
The intraclass correlation coefficients, r MZ and r DZ , were calculated for each facial parameter for 1208 twins who had a pair (263 MZ and 341 DZ). These formed the bases for further analyses on genetic and environmental contribution to facial morphology, as explained below.

Glabella (g)
The most prominent midpoint between the eyebrows Nasion (n) Midline point between the nasal root and nasofrontal suture, above the line that connects the two inner canthi The majority of the landmarks were identified according to the definitions given by Farkas [39]. The landmarks marked with an asterisk (*) were defined by the authors of the study.

Error estimates
To estimate the variance due to errors in scanning and landmarking, a sample of 73 faces (5% of the original database sample) was randomly selected for an intra-examiner reliability analysis. The faces were landmarked on two occasions six months apart by placing 37 landmark points. The reliability of a landmark was evaluated as the mean 3D Euclidean distance between its positions determined at two landmarking sessions. For rough estimates, it was assumed that all the landmarks were independent. The total variance of landmark coordinates was evaluated as the sum of variances of all individual coordinates across the sample. After Procrustes registration of either set of 73 landmark configurations, the variances of the 111 individual coordinates were calculated. The total variance (V) for one set was 519.1 mm 2 and 485.9 mm 2 for the other, with the mean equal to 502.5 mm 2 .
According to Lübbers et al. [36], the mean global error of a 3dMD stereophotogrammetry system was 0.2 mm for a mannequin head. For live subjects, this should be at least doubled. The total variance due to scanning errors (V SE ) was calculated as 37 x 0.4 2 = 5.9 mm 2 , which makes up about 1% (5.9/502.5) of the total variance.
The variance due to landmarking errors can be estimated as follows. In the reliability sample, each landmark was placed twice. If x 1 and x 2 are two measurements of the x-coordinate of a landmark, its true value and the error are estimated as (x 1 +x 2 )/2 and |x 1 −x 2 |/2, respectively. By summing up the variances of the individual errors of all coordinates, we get the variance due to landmarking errors (V LE ) equal to 25.1 mm 2 , which makes up about 5% (25.1/502.5) of the total variance.

Statistical analyses
The intraclass correlation coefficients with the corresponding p-values and 95% confidence intervals for all facial parameters (principal components and linear distances) of 1208 twins (263 MZ and 341 DZ pairs) were evaluated using a non-parametric approach implemented in R programs [47]. A bootstrapping technique and a resampling technique with 100,000 random permutations each were used for the calculations [48]. The level of statistical significance was set at 0.05.

Results
The intra-examiner reliability for 31 out of 37 landmarks was less than 2 mm (Table 2). Six remaining landmarks had reliability between 3.0 and 6.4 mm. The sample size was adequate and all the landmarks were used for the subsequent analyses.
The findings of the principal component analysis for the total sample were replicated in both subsamples (S1 Table). Nineteen uPC (accounting for 87% of the total variance in facial form) were found for the unscaled dataset and 23 sPC (accounting for 83% of the total variance in facial shape) were found for the scaled dataset (Table 3). Although each principal component was contributed by all 111 coordinates, there were only relatively few that showed significant contribution (S1 Table). The anatomical interpretation of the above principal components (Table 4) was based on only those coordinates that contributed the most (factor loadings over 0.5 in magnitude). The effect of a PC on the face may be illustrated using a sequence of average faces [27] corresponding to different PC scores (Fig 2).
Genetic and environmental contributions to facial principal components are presented in Table 5. Broad-sense heritability (h 2 ) was statistically significant for majority of principal components: 13 out of 19 uPC (ranging from 38.8% to 78.5%) and 17 out of 23 sPC (ranging from 30.5% to 84.8%). Heritability of 50% and above was found for 12 uPC and 10 sPC. Morphology of the nose (width, height and prominence), prominence of the lips, inter-ocular distance and facial size had the highest heritability (all above 70%). Common environment significantly contributed to the variation in uPC11 and sPC12 (vertical position of the zygoma). Potential dominant genetic effect was found for sPC1 (prominence of the upper lip relative to the chin). The contribution of unique environment was statistically significant for 17 out of 19 uPC (ranging from 23% to 83%) and 20 out of 23 sPC (ranging from 24.6% to 87.7%). The highest unique environment contribution was found for uPC13 (horizontal asymmetry of the nasal root, lips and chin), uPC19 (representing sagittal position of zygoma and mandibular angle) and uPC10 (horizontal asymmetry of the nose, philtrum and lips). High unique environment contributions of over 50% was also found for uPC18, uPC16 and uPC11. These principal components were associated with the vertical position of the zygoma and the mandibular angle (could be interpreted as mandibular ramus height) and vertical position of the landmark 'nasion' (i.e., nasal root). In the scaled dataset, these corresponded to sPC8, sPC11, sPC18 and sPC20.
The genetic and environmental contributions to facial linear distances are presented in S2 Table. Thirty most heritable distances are presented in Table 6. Statistically significant genetic contribution was found for 1222 out of 1247 investigated independent distances. Broad-sense heritability of 70% and above was found for 476 distances (38.2%) and a heritability of 90% and above for 41 (3.3%) distances.     Common environment was found to have a significant positive contribution to 53 out of 1247 distances (4.3%). These distances are between: 1) bilateral landmarks 'zygion' and different landmarks located on the tip of the nose, lips, mandibular angle (landmark 'gonion') and chin and 2) eye and lip landmarks and middle points between bilateral 'zygion' and 'gonion'.
Possible dominant genetic effect was found for 11 distances (0.9%). These were: three distances between right eye landmarks and the tip of the nose ('prn'); one distance between the bilateral 'exocanthion' midpoint and the tip of the nose ('prn'); two distances between the tip of the nose and two nose base landmarks ('sn', 'phnr'); two distances between the most prominent point of the chin ('pg') and two lower lip landmarks ('li', 'mmll') and three distances between upper lip philtrum points ('mphl', 'mphlr' and 'dpc') and nose base points ('sbalr', 'msbal'). These distances represent the prominence and height of the nose, prominence of the lower lip in relation to the chin and length of the upper lip philtrum. The contribution of unique environment was significant in 1245 out of 1247 distances.  r MZ and r DZ , intra-class correlation coefficient in MZ and DZ twins respectively; h 2 , broad-sense heritability; c 2 , relative common environment contribution; e 2 , relative unique environment contribution. S2

Discussion
Phenotypic variation in humans is produced through a complex interplay between genotype and environment, but the characterization of phenotype lags behind the characterization of the genotype [49]. In facial research, there is no uniform approach to the analysis of facial phenotype due to: 1) variable facial data acquisition (two-dimensional or three-dimensional techniques); 2) lack of a standardized way to quantify spurious expressions (relying on subjective opinion of the examiner); 3) variable selection and/or definition of phenotypic traits; 4) the amount of information analysed (sparse or spatially-dense points across the facial surface that can be identified manually or automatically) and 5) statistical analyses (univariate/multivariate).
In this study, facial phenotype was characterized by principal components and linear distances based on 37 anthropometric landmarks manually identified on the 3D facial images. Most of the landmarks could be reliably identified (within 2 mm), except for soft tissue 'zygion' (difficult to describe in anatomical terms and traditionally identified by trial measurement), 'gonion' (normally identified by palpation), 'pogonion' and 'gnathion' (these were difficult to identify in vertical plane if a chin had a flat surface) [39]. However, as the sample was quite large (1380 twin faces), it was decided to include these landmarks in the subsequent analyses.
We found that soft tissue facial traits in adult British female twins have moderate to high heritability, which is in general agreement with previous family and twin studies. In Koreans [6], the heritability values ranged from 0.25 (lower facial height, sn-gn) to 0.61 (intercanthal width; en-en) for linear measurements derived from digital photographs of family members; the authors also identified three factors with heritability estimates from 0.45 (total face height, upper face height and nose height) to 0.55 (lower face height and width of the mandible, mouth and nose). High correlations between parents and offspring and siblings were also found in Indian families [50] for mandibular position, chin prominence, nasal prominence, nasal width, lip length at philtrum, lip prominence and facial height. However, heritability was not calculated. In Belgian families [51], heritability estimates for soft tissue facial parameters ranged from 0.46 (nose height) to 0.72 (external biocular breadth). Cephalometric analysis of soft tissue parameters in Turkish siblings [21] revealed high heritability (0.72 to 1.0) for soft tissue chin thickness, soft tissue facial angle, Merrifield angle (formed by the Frankfort plane and profile line joining the chin and the more prominent lip, usually the upper) and Holdaway angle (formed by a line tangent to the chin and upper lip with the cephalometric NB line). Differences in estimates can be explained by different methodology: photos [6,50,51] or lateral cephalograms [21]; measurement errors; small samples and different ethnicities. In addition, the components of genetic variance for a given trait vary from population to population [51].
The first 3D twin study focusing on soft tissues [28] was performed on British twins (10 pairs of MZ and 8 pairs of same-sex DZ twins aged 9 to 17 years) using stereophotogrammetry. Significant differences between faces of MZ and DZ twins were found for intercanthal width, right eye width, nose width, nose height, mouth width, and upper lip height. This coincides with our results, despite the obvious differences in the sample size and age. In [29], faces of British twins (10 pairs of MZ and 10 pairs of same-sex DZ twins) were also analysed using laser scanning. It was apparently the first time that facial surfaces were suggested for analysis apart from analyzing linear distances. Significant genetic determination was revealed for midfacial parameters, especially left eye width, intercanthal width, nose height, and nose width. However, the study by Naini and Moss [29] did not show any significant differences in mouth width and upper lip height, which contradicts the study done by Burke [28] as well as our findings, which show a strong genetic contribution above 60% for these traits. The differences can be explained by the small sample used in [29], wide age range of participants (6-42 years) and their mixed ethnicity, as well as no formal heritability calculation.
In [38], a well-known British population study (ALSPAC) was the source of 37 twin pairs, whose faces were captured using laser scanning. Configurations of 21 facial landmarks as well as facial surfaces were compared between 19 MZ and 18 DZ twin pairs aged 15 years. Procrustes analysis did not reveal any significant difference in facial landmark configurations. On the other hand, average female MZ and DZ twin faces coincided in the eyes, supraorbital and infraorbital ridges, philtrum and lower part of the cheeks. In the absences of heritability estimates, the findings of that study indirectly show that central facial structures are the most heritable ones.
A preliminary study performed on American twins [32], 10 MZ and 11 same-sex DZ twins 5-12 years of age, found that only three out of nine extracted principal components showed statistically significant genetic contribution. These were related to the horizontal separation between the eyes, the length, breadth and projection of the nose, and the height and projection of the upper lip. Heritability estimates approached 1.0 and the authors explained this over-estimation by the small sample size, small number of landmarks and a very crude calculation of heritability. In our study, a significantly greater number of landmarks were used and therefore more principal components were extracted. Most of these components showed an evidence of statistically significant genetic contribution.
Various landmark-based traits (e.g.,distances, angles, ratios and principal components) as well as surface-based traits (e.g., geodesic distances and curvatures) can be used to seek genes responsible for normal facial morphological variation. Only a few genome-wide association studies have been conducted so far, which have revealed a relatively small number of associations between certain facial traits and single-nucleotide polymorphisms (SNPs). The findings are in agreement with our results on highly heritable traits, especially those associated with the morphology of the nose and lips. The SNPs found in the PAX3 gene [52][53][54] are associated with the nasal root morphology in Europeans and Latin Americans. In addition, multiple intronic SNPs in the PRDM16 gene are associated with nose width and nose height [53]. A SNP close to the C5orf50 gene is associated with the position of the landmark 'nasion' [53]. An intronic SNP in TP63 gene is associated with the inter-ocular distance and a missense SNP in COL17A1 is associated with the distance between the eyes and 'nasion' [53]. Significant associations were recently found for three more nose-related traits: columella inclination (4q31), nose bridge breadth (6p21) and nose wing breadth (7p13 and 20p11) [54]. The rs642961 SNP in the IRF6 gene (a known risk factor of non-syndromic cleft lip/palate) was found to strongly predict normal lip shape variation in Han Chinese females but not males [55].
The finding that facial asymmetry is not genetically driven complies with the results of two studies, [38] and [56]. In the first one, the amount of three-dimensional asymmetry was calculated after superimposing (registering) the original face with its mirror reflection and measuring the average distances between the two facial surfaces. There was no statistically significant difference in the amount of facial asymmetry between MZ and DZ twins [38]. In the other study, the relationship between facial asymmetry (evaluated from nine mid-facial landmarks) and genetic variation at 102 SNP loci (recently associated with facial shape variation) was investigated. The authors failed to identify any SNP relating to either fluctuating or total asymmetry [56].
Our finding on the mandibular ramus height is in agreement with a recent cephalometric study conducted on 141 adult Lithuanian twin pairs with completed mandibular growth and DNA confirmed zygosity [57]. The authors estimated the significance of additive (A) and nonadditive (D) genetic factors as well as shared (C) and unique environment (E) using a maximum likelihood genetic structural equation. Their results indicate that the shape and sagittal position of the mandible is under stronger genetic control than is its size and vertical relationship to cranial base. For linear measurements, such as mandibular body length, ramus width and ramus height, the best-fitting model was found to be ACE, indicating low genetic determination.
The present study has certain advantages and limitations, which will be discussed below. One of the merits is a large sample size, which enabled us to use additional landmarks (in comparison to previous studies on facial morphology) and contributed to the validity of heritability estimates (without compromising statistical power). In addition, the sample was homogenous in terms of ethnicity and came from a well-designed population-based study. Finally, zygosity was confirmed by genetic testing, hence avoiding the misclassification of twins.
On the other hand, the limitations of the study were due to: 1) limitations of the classical twin design, 2) complexity of facial morphology and 3) inclusion of only female individuals in the sample. The observed facial morphological variation is a combination of: 1) genetic contribution (encompassing additive and non-additive genetic effects as well as gene interactions), 2) environmental contribution (consisting of common and unique environment), 3) gene-environment interaction and 4) measurement errors (due to scanning and landmarking).
The trait correlation between twins is the result of their genetic similarity and sharing common environment. However, the classical twin design does not allow for the determination of any gene-environment interactions. In addition, the scanning and landmarking errors (usually included in the unique environment component) may have some effect (generally negligible) on the heritability values of some traits, especially those that just reach statistical significance (p-values less than but close to 0.05). Genetic contributions calculated here are likely to be slightly overestimated because the model disregards gene-environment interactions. Furthermore, the imperfections of the classical twin model can lead to heritability values (h 2 ) over 1, which are demonstrated by the following results: the heritability (h 2 ) of sPC1 was evaluated as 1.015 (95%CI: 0.755 to 1.279), h 2 of the linear distance 'prn-men' was 1.002 (95%CI: 0.798 to 1.217) and h 2 of the 'prn-mex' distance was 1.048 (95%CI: 0.832 to 1.271). However, these do not diminish the importance of our findings. Instead of focusing on the actual value of an estimate, it is more important to reveal which facial traits demonstrate the most compelling evidence of heritability and use that knowledge for future genome-wide association studies of normal facial morphology.
The second limitation is related to the complexity of studying facial morphology (as explained in the introductory paragraph of the discussion). The third limitation is due to inclusion of only females in the sample. This reflects the prevalence of females in the TwinsUK register [34]. In the data available to us, there were 15 males, which were excluded intentionally as they were too few for a separate statistical analysis; the inclusion of male individuals in a predominantly female sample could have affected the findings. This issue needs to be addressed in future studies. It would be interesting to look at a sufficiently large male twin sample and compare the results with those for female twins. In addition, the sample we dealt with was ethnic specific and further research is needed on other populations, since environmental effects and gene alleles frequencies may differ between populations [51,52].

Conclusions
The study provides the estimates of genetic and environmental contributions to three groups of landmark-based facial traits in adult female twins. Based on the analysis of principal components, statistically significant genetic influence on the facial form was found to range from 38.8% to 78.5%, whereas that on the facial shape accounts for 30.5% to 84.8% of the total phenotypic variance. Genetic factors can explain more than 70% of the phenotypic variance in 7 principal components related to facial form, 5 principal components related to facial shape and 474 linear distances. These facial parameters represent: facial size (height), nose (width, prominence and height), lips prominence and inter-ocular distance. A few traits show potential dominant genetic influence, namely the prominence and height of the nose, the prominence of the lower lip in relation to the chin and length of the upper lip philtrum. The highly heritable traits are likely candidates for genome-wide association studies. Environmental contribution to facial variation is the greatest in the mandibular ramus height and horizontal facial asymmetry. This heritability study may inform future genetic studies which facial traits should be focused on.
Supporting Information S1 Table. The results of the Principal Component Analysis (PCA) for the total sample and two sub-samples (unscaled and scaled datasets). Writingreview & editing: JD AIZ SR.