Additionally, three direct cranial measurements: maximum cranial breadth (Euryon - Euryon), maximum cranial length (Gonion - Opisthocranium) and maximum cranial height (Vertex - Gnathion), were collected manually using a digital spreading calliper (Paleo-Tech Concepts, USA). Based on the craniofacial and body height measurements, three craniofacial ratios were calculated: Cephalic index: (eu-eu)/(g-op), Head width - Craniofacial height index: (eu-eu)/(v-gn) and Head - Body height index: (v-gn)/(body height), as summarised in Table 1. The linear and angular facial distances were used to calculate 10 principal components (PCs). Additionally, the details on the eyelid and earlobe morphology (single/double and attached/detached respectively) were recorded. In total, 104 craniofacial traits for each individual were recorded and used for genetic association analyses. All the facial measurements were recorded by a single analyst to reduce potential variability.

Candidate genes and SNPs selection

Two complementary strategies were used to generate a preliminary list of candidate genes and genetic markers. The first focused on searching the literature and web resources for candidate genes involved either in normal craniofacial variation or in craniofacial malformations in humans and model organisms (Supplemental Table S2).

The search for candidate genes focused not only on specifically defined craniofacial disorders, but also on genetic syndromes with various manifestations of craniofacial malformations, such as Down syndrome, Noonan Syndrome, Floating-Harbor Syndrome and others, as detailed in the Supplemental Table S2. The main resources for locating candidate genes in the animal models were Mouse Genome Informatics [30] and AmiGo databases [31] The main resources for identifying candidate genes in the human genome were OMIM [32] and GeneCards websites [33]. The online search was performed using the keywords such as: "craniofacial", "craniofacial mutants", "craniofacial anomalies", "craniofacial dimorphism" and "facial morphology" (as summarized in the Supplemental Appendix S1).

The second approach initially implemented a broad search for high Fst SNPs, such as ancestry informative markers (AIMs), with the rationale that many genes affecting craniofacial traits would express significantly different allele frequencies across populations, following non-neutral evolutionary pressure. AIMs were selected from a variety of published and online resources [34-43].

The relevant genes obtained by both approaches were subsequently analysed for potential involvement in the craniofacial embryogenesis, limb development and bilateral body symmetry. It should be noted however, that the final candidate gene list was not limited to craniofacial genes and included high Fst SNPs in genes with unknown function as well as markers located in intergenic regions, potentially possessing regulatory functions.

The resulting set of SNPs was further screened for high Fst SNPs (≥ 0.45) in three '1000 genomes' populations (CAU, ASW, CHB) using ENGINES browser [44] as well as potentially functional polymorphisms, such as non-synonymous SNPs [45], markers in transcription factor binding sites [46] and splicing sites [47] using various web resources, as detailed in the Supplemental Appendix S1 and reviewed elsewhere [48]. The candidate markers search resulted in identification of 1,319 SNPs, located in 177 genes and intergenic regions, as discussed in the Results section. In addition to the craniofacial candidate markers, 464 pigmentation and ancestry - informative polymorphisms were selected from the literature (Supplemental Appendix S1).

The chromosomal locations of final candidate markers (n = 1,783) were submitted to the custom Ampliseq primer design pipeline (Life Technologies), according to manufacturer recommendations. Following primer design failure for 881 markers, the marker list was redesigned to include alternative SNPs in high linkage disequilibrium with the original markers, resulting in 1,670 candidate polymorphisms. The final custom Ampliseq panel was manufactured as two separate pools of 849 and 847 primer pairs, with each amplicon length between 125 bp and 225 bp, covering 15.78 kb of the human genome. This panel comprised of 902 originally targeted candidate markers and 4,732 alternative tag markers, as a replacement for candidate SNPs that failed primer design pipeline. Subsequent genomic annotation of the total amplicon sequences for SNPs with MAF ≤ 1% resulted in 8,518 markers as detailed in Supplemental Table S1 (Supplemental Appendix S2).

DNA extraction and quantification

DNA was purified from buccal swabs using the Isohelix DDK isolation kit (Cell Projects, Kent, UK) according to the manufacturer instructions. DNA samples were quantified using a Real Time quantitative PCR (q-PCR) method using a Bio-Rad CFX96 (Bio-Rad, Gladesville, Australia). This assay amplified a 63bp region of the OCA locus. The primer sequences were 5'-GCTGCAGGAGTCAGAAGGTT-3' (forward primer) and 5'- CATTTGGCGAGCAGAATCC-3' (reverse primer) at a final concentration of 200 mM. All DNA samples were additionally quantified using Qubit 2.0 fluorimeter (Invitrogen) prior to library construction as per manufacturer recommendations.

SNP genotyping and data analysis

Multiple DNA libraries were constructed from sets of 32 Ion Xpress^™ (Life Technologies) barcoded samples using the Ion AmpliSeq^™ library Kit 2.0 (Life Technologies) in conjunction with two custom primer mixes that were pooled according to manufacturer recommendations. Libraries were quantified using the Ion Library Quantitation kit (Life Technologies) and pooled in equal amounts for emulsion PCR, which was performed using the OneTouch^™ 2 instrument (Life Technologies) according to manufacturer recommendations. All DNA samples (n = 587) were genotyped by massively parallel sequencing on the Personal Genome Machine (PGM) (Life Technologies) using the Sequencing 200 v2 kit and 316 Ion chips (Life Technologies).

Raw sequencing data were collected and processed on the Torrent Suite Server v3.6.2 using default settings. Alignment and variant calling were performed against the human genome reference (hg19) sequence at low stringency settings. Binary alignment map (BAM) files were generated and exported into Ion Reporter^™ (IR) cloud-based software for SNP annotation along with the reference hotspot file. The IR analysis resulted in generation of the individual variant caller files (VCF) with genotype calls for each sample as well as sequencing performance statistics.

To reduce potential bias of the self-reported ancestry, ancestry inferences were obtained by 3,302 markers using STRUCTURE version 2.3.4 with default parameters as per software developer recommendations [49]. SNPs in long-range Linkage Disequilibrium (> 100,000 bp) were excluded from the STRUCTURE run. The ancestry was estimated based on four predefined population clusters: Europeans, East Asians, South Asians and Africans, according to software developer recommendations. Relative allele calls for four predefined HapMap population clusters (CEU, YRI, CHB and JPT) were used as reference populations [50]. The ancestry origin was estimated as a single (unmixed) source where the main ancestry cluster could be affiliated with at least 80% of the total mixed ancestry. The samples with mixed ancestry (> 20% admixture) were assigned to an 'Admixture' cluster.

Association analyses were performed using SNP & Variation Suite v7 (SVS) (Golden Helix, Inc., Bozeman, MT) and replicated using PLINK v1.07 software [51]. Statistical analyses in both software programs were performed using linear regression with quantitative phenotypes, and logistic regressions with binary phenotypes under the assumption of an additive genetic model, while each genotype was numerically encoded as 0, 1 or 2. Population stratification correction, incorporated by EIGENSTRAT program was implemented in the analyses [52,53]. In order to reduce any potential confounding effects, all the craniofacial traits association analyses were performed using sex, BMI and EIGENSTRAT ancestry clusters as covariates. In PLINK, p-values were adjusted using the '-adjust' option. The final reported association results are based on the PLINK statistical analyses with the EIGENSTRAT PCA clusters, BMI and gender as covariates.

Annotation analysis of the significantly associated genes was performed using GeneCards, ENTREZ and UniProtKB web portals [33,54]. The MalaCards web site was used to detect association between the genes and hereditary syndromes [55]. The GeneMania web site was used to identify a functional network among the genes and encoded proteins [56]. Gene ontology web resource was used to find orthologs of human genes in other organisms [31,87]. The MGI database was used to search for the phenotype in relevant craniofacial mouse gene mutants [30]. The dbSNP, 1000 genomes, SNPnexus and Alfred websites were used for SNP annotations [58-61].

The SNP Annotation and Proxy Search (SNAP) web portal was used to find SNPs in linkage disequilibrium (LD) and generate LD plots, based on the CEU population panel from the 1000 genomes data set, within a distance of up to 500 kb and an r² threshold of 0.8 [62].

The Regulome database and potentially functional database (PFS) searches were implemented to annotate SNPs with known and predicted regulatory elements in the intergenic regions of the H. sapiens genome [47,63].

Results

Candidate genes search and sequencing data filtering

The online database search for candidate SNPs revealed a list of 150 genes with 310 annotations. In parallel, a search for high Fst markers, resulted in identification of additional 938 targets, for a total of 1,088 genes and intergenic regions (a detailed list of online resources is summarized in the Supplemental Appendix S1).

Manual examination of the candidates list revealed no apparent link with normal craniofacial development or malformations for 592 of these genes and these were therefore excluded. The remaining 496 regions were further screened for non-synonymous and potentially functional SNPs, as well as SNPs with high population differentiation, which resulted in a shortlist of 269 genes and intergenic regions.

Subsequent analysis of these 269 genes/regions for functional annotation, resulted in 177 candidate genes/intergenic regions, possessing 1,319 genetic polymorphisms potentially involved in various stages of human embryonic development, including: embryonic morphogenesis, sensory organ development, tissue development, pattern specification process, tissue morphogenesis, ear development, neural tube morphogenesis, epithelium development, chordate embryonic development and morphogenesis of an epithelium (Supplemental Appendix S1). Notably, the majority of these markers are located in introns and intergenic regions, including 5' and 3' UTRs. This outcome is quite expected, based on the initial assumption that most of the genetic variants influencing phenotypic traits lie in the noncoding part of the human genome, as demonstrated by numerous genome - wide association studies [64]. Such sequence polymorphisms are mostly rare (MAF < 0.1) and have tiny effect size, while thousands of them required to jointly affect complex traits, such as facial morphology or body height [65-68].

In terms of molecular function, AmiGO showed that craniofacial candidate markers might be involved in a range of regulatory activities including: protein dimerization activity, chromatin binding, regulatory region DNA binding, sequence-specific DNA binding RNA polymerase II transcription factor activity, sequence-specific distal enhancer binding activity, heparin binding, RNA polymerase II core promoter proximal region sequence-specific DNA binding transcription factor activity involved in positive regulation of transcription, BMP receptor binding and transmembrane receptor protein serine/threonine kinase binding (Supplemental Appendix S1). The final candidate marker list annotation revealed that almost 99% of the genomic regions (which may cover multiple markers) are associated with one or two genes with 62% of genomic regions located up to 500 kb downstream of a transcription start site (data not shown). This is again an anticipated outcome, given that non-coding variants are likely to reside within cis-regulatory genomic elements, such as enhancers [69,70].

Subsequent analysis of candidate SNPs for mouse phenotype associations confirmed that orthologous candidate markers were previously detected in mouse models displaying abnormal morphology of the skeleton, head, viscerocranium and facial area, as well as specific malformations of the eye, ear, jaw, palate, limbs, digits and tail (data not shown).

Targeted massively parallel sequencing of the 587 samples resulted in 8,518 genetic markers, with the majority of markers (> 5,000) represented by rare polymorphisms with MAF ≤ 1% (data not shown). This outcome was a result of identification of novel and rare markers in individual DNA samples. Three of the DNA samples did not produce high quality genotypes because of poor DNA quality or unsuccessful library and template preparation.

Data quality control steps were performed by removing markers of low genotype quality (GQ > 10) and sequencing depth (DP > 10X). Further filtering of markers using a 2% MAF cut-off resulted in 3,075 markers (Supplemental Appendix S2). Additional filtering based on the HWE threshold of p-value ≥ 0.01 resulted in 2,332 markers. The mean sequencing depth for significantly associated markers was 58 fold (± 48.9 SD).

Genetic association analysis

The association analyses were performed using a linear regression model, incorporating EIGENSTRAT-generated PCAs as well as sex and BMI as covariates. The use of covariates in the statistical analysis aimes to reduce the risk of introducing confounding effects, which can result in false positive associations. While sexual dimorphism in the craniofacial morphology is well-known [71], BMI is likely affect certain craniofacial traits, since the soft facial tissue may change significantly with weight gain or loss. Notably, this potential confounding factor has to date been disregarded in association studies of normal craniofacial morphology. Age was not considered a significant covariate, given that average age of the subjects in this study was 27 (± 8.9 SD). Nevertheless, the potential effect of age as a cofactor was assessed on three craniofacial traits and found not significant (data not shown).

While the majority of current GWA studies rely on a p-value < 5.00E-08 significance threshold, some publications suggest this threshold may be too stringent, especially for complex traits that are regulated by a large number of small effect alleles [72,73]. In contrast to GWAS, candidate gene studies undertake a more focused genetic strategy, concentrating on a relatively limited number of putative markers. As this study analysed a significantly lower number of SNPs than usual GWA-studies, we could use a higher p-value cut-off, since the smaller sample size means the probability of false positive at extremely low p-values is itself lower. Nevertheless, we decided to keep the traditional GWAS p-value significance threshold (< 5.00E - 08) in order to reduce the possibility of detecting false positive results.

In addition, we subsequently applied a more stringent Bonferroni - corrected threshold in order to minimize the chance of detecting spurious associations. Following the association analysis of 104 craniofacial phenotypes with 2,332 genetic markers, the significance threshold based on the Bonferroni correction with a desired α of 0.05 would be 2.06E - 07 (= 0.05/(2,332^*104)).

However, it should be emphasized that the Bonferroni correction is widely considered over-conservative, especially in the case of complex phenotypic traits with small individual effects of each allele. Considering that our results confirm the previously published findings, we believe the GWAS p-value threshold is conservative enough to avoid or at least significantly reduce potentially spurious associations. Following this rationale, we report all the variants, which met the unadjusted genome-wide association p-value threshold, as suggestive. We believe these findings are useful for the future studies focusing on genetics of normal craniofacial morphology.

The results of the association analyses of the craniofacial traits are summarized in Table 2 and Supplemental Figures S1-S19. In general, following the application of a stringent Bonferroni-corrected GWAS threshold (adjusted p-value < 1.6E - 07), we observed five craniofacial traits being associated with six genomic markers. Specifically, nasal width with rs8035124 (p-value 1.74E - 07, Beta = 1.366, SE = 0.209), cephalic index with rs16830498 (p-value 8.67E - 08, Beta = 3.005, SE = 0.4518), nasal index with rs37369 (p-value 1.43E - 07, Beta = 4.025, SE = 0.6124), transverse nasal prominence angle with rs59037879 (p-value 6.07E - 09, Beta = 4.765, SE = 0.6685) and rs10512572 (p-value 1.57E - 08, Beta = 1.505, SE = 0.2171), and principal component (EV = 1391.99) with rs37369 (p-value 2.85E - 08, Beta = -0.021, SE = 0.003079) and rs390345 (p-value 8.55E - 08, Beta = -0.0184, SE = 0.002768). The polymorphisms: rs16830498, rs59037879 and rs390345 are intronic variants in CACNB4, ZEB1 and FOXN3 respectively; rs37369 is a missense mutation in the AGXT2 gene and rs8035124 and rs10512572 are intergenic variants in 15q12.2 and 17q21.33 chromosomal locations respectively.

However, given the over-conservative nature of the Bonferroni correction and also assuming that polygenic traits are likely to be dominated by numerous alleles with small causal effect we also report suggestive associations reaching the unadjusted 5.00E - 08 p-value threshold. In summary, two linear distances (nasal width and nasal tip protrusion), two angular distances (nasolabial angle and transverse nasal prominence angle), three indices (cephalic index, nasal index and nose-face width index) and one principal component revealed significant associations with 33 SNPs in 28 genes and intergenic regions (Table 2).

These factors can be arbitrarily divided into three main categories based on their cellular function: 1) Genes with known roles in the craniofacial morphogenesis and/or mutated in various hereditary syndromes displaying craniofacial abnormalities; 2) Genes or pseudo-genes without known function in the craniofacial morphology regulation or previously uncharacterized genes; and 3) Non-protein coding genes, such as lncRNA class genes. There are also a number of significant variants that are located in the intergenic regions, with or without proximity to open reading frames (ORFs).

The majority of associated markers (n = 21) are located in 17 protein-coding genes and pseudo-genes such as AGXT2, CACNB4, COL11A1, EGFR, EYA1, EYA2, FAM49A, FOXN3, LHX8, LMNA, MYO5A, PAX3, PCDH15, RTTN, SMAD1, XXYLT1 and ZEB1. Five variants are present in RNA-coding (lncRNA) genes, which include AC073218.1, RP11-494M8.4, RP11-408B11.2, RP11-785H20.1 and TEX41. The rest of the markers (n = 6) are found in the intergenic regions, near the following genes and pseudogenes: BMP4, HAS2-AS1, LOC124685, LOC100131241, MRPS36P3, PTCH1, SLC25A5P2, SV2B and TRNAY16P.

Functional annotation of the significantly - associated markers revealed that one variant represent a missense mutation (rs37369), one SNP is a synonymous transversion (rs2290332), while 20 markers are found in intronic sequences and 11 markers located in intergenic regions (Table 1). As stated above, the noncoding genomic location for most of the markers makes sense as they expected to have mild affect on a normal craniofacial phenotype than exonic polymorphisms. This conclusion however, could be biased as a priori, most candidate markers were selected in the non-coding regions (Supplemental appendix S1). As expected, the majority of significantly associated SNPs (n = 27) are found in the regulatory elements of the genome, such as in transcription factor (TF) binding sites, and represent potentially functional SNPs (pfSNPs). These variants may be involved in regulation of gene expression by "fine tuning" of the normal craniofacial phenotype as part of the enhancer/silencer mechanisms, as has been recently suggested [69].

The nasal area measurements, using either "n", "prn", "sn" or "al" landmarks, produced the majority of the total number of significant associations (6 out of 8). These measurements include nasal width (al-al), nasal tip protrusion (sn-prn), nasolabial angle (prn-sn-ls), transverse nasal prominence angle (t_l-prn-t_r), nasal index (al-al/n-sn), and nose face width index (al-al/zy-zy). The apparent overrepresentation of associations with the nasal area is likely a result of the easier allocation and consequent superior reproducibility of the nasal area landmark measurements on 3D images. It may also be a result of a focused selection of candidate genes from the JAX mice database resource, which concentrated on mutants with nasal abnormalities.

The analysis of the direct cranial measurements and their relative indices revealed significant associations for cephalic index (CI) with 9 SNPs (Table 2). The association analysis of the principal components (PC) representing all the craniofacial measurements, revealed one principal component (explaining 73.3% of all the craniofacial phenotypes) that was associated with 14 genetic markers (Table 2).

In contrast to most other craniofacial association studies that focused on a specific homogeneous population group (mostly Europeans), this study included samples from several population groups, which enabled investigation of the genetic factors influencing normal craniofacial morphology in different ethnicities [74]. Self-reported ancestry however, cannot be considered fully reliable, as demonstrated previously [75,76]. In order to address this issue we assessed the self-reported ancestry using STRUCTURE with 186 SNPs removed due to long-range disequilibrium [49]. Following the rationale that the best ancestry estimates are obtained using a large number of random markers [77], we used all the available markers (after MAF ≤ 2% filtering) in STRUCTURE analysis. The STRUCTURE analysis resulted in clusters of 367 Europeans, 51 East Asians, 43 South Asians and 16 Africans, with 107 samples designated as admixed ancestry (Figure 3). Of the samples tested with STRUCTURE, 459 (89%) were assigned the same ancestry cluster (sole or mixed origin) as the self-reported information. Of the remaining 57 individuals, 39 were estimated as 'admixture' (based on up to 20% admixture threshold) and 18 were assigned a single ancestry, different to the self-reported ancestry (Figure 3).

The risk of detecting false positive results due to population stratification was carefully assessed and further reduced by applying an EIGENSTRAT correction. Specifically, EIGENSTRAT's smartpca.perl script was used to perform PCA-clustering in comparison to reference populations from HapMap reference clusters. The Q-Q plots of the associated traits confirmed the expected distribution of data after applied correction (Supplemental Figures S1-S8).

Due to a limited number of genotyped markers typed, genomic allele imputations were not performed.

Association analyses of the non-craniofacial traits

In our attempt to identify genetic markers influencing normal variation in craniofacial traits, we initially incorporated 464 markers previously associated with human pigmentation and ancestry - informative markers, while many of these overlap. These markers were included to validate the statistical methods used for the craniofacial traits association analysis.

The application of the Hardy-Weinberg equilibrium (HWE) threshold resulted in filtering 25% of the total number of SNPs. These markers included almost all the "key" SNPs, previously associated with pigmentation traits, such as rs12913832, rs1129038, rs8039195 and rs16891982. This is not surprising, since population-specific markers are likely not being in HWE 'a priori'. Another explanation for this observation is potential bias from partially uncorrected heterogeneous ancestry, since the ancestry correction algorithm can only minimize, rather than completely remove spurious associations [52]. In fact, the association analyses of the HWE non-filtered genotyping data with pigmentation traits (eye, skin and hair colour), demonstrated highly significant associations, fully concordant with the previously published sources (Supplemental Table S3). It should be noted however, that these results may not necessarily confirm the craniofacial markers associations.

Craniofacial genes and SNPs annotations

Due to space limitations, the functional annotation of the significantly associated SNPs and their potential role in regulating normal craniofacial phenotype has been provided in a Supplemental Appendix S3.

Discussion

This study focused on the identification of genetic markers in a set of candidate genes potentilly involved in normal craniofacial development. We report significant associations (Bonferroni adjusted p-value threshold of < 1.6E - 07) between 5 craniofacial traits (nasal width, cephalic index, nasal index, transverse nasal prominence angle and principal component) and 6 SNPs (rs8035124, rs16830498, rs37369, rs59037879, rs10512572 and rs390345) located in 6 genes and intergenic regions (15q26.1, 17q24.3, CACNB4, AGXT2, ZEB1 and FOXN3 respectively). We also report all the significant associations that met the less stringent GWAS threshold (unadjusted p-value < 5.00E - 08) between 8 craniofacial traits and 33 polymorphisms in 28 genes and intergenic regions, as Bonferroni correction is generally considered over-conservative, especially when analysing complex traits such as craniofacial morphology, which is likely to be influenced by a large number of alleles with relatively small individual effect, similar to body height [65-67,78-80]. In fact, most of the significant SNPs (31 out of 33) are noncoding variants, located in the regulatory regions such as: promoters, enhancers and splicing sites, as discussed below. These associations are suggestive and would require additional studies for potential confirmation.

The association of the PAX3 gene and the COL11A1 gene with the transverse nasal prominence angle and cephalic index respectively, confirms previous findings [11,22,23,80]. In fact, an intronic SNP rs11164649 that was associated with cephalic index in the current study, was recently associated with normal-range effects in various craniofacial traits and used for their prediction [80], while the other variants in COL11A1 (rs4908280) and in PAX3 (rs2289266) have not been reported previously. The rest of the significantly - associated markers are also novel and include 20 polymorphisms in protein-coding genes and pseudo-genes, such as AGXT2, CACNB4, EGFR, EYA1, EYA2, FAM49A, FOXN3, LHX8, LMNA, MYO5A, PCDH15, RTTN, SMAD1, TEX41, XXYLT1 and ZEB1. Interestigly, at least 8 of these genes (EYA1, EYA2, FAM49A, FOXN3, LHX8, SMAD1, TEX41, ZEB1) represent transcription factors, affecting numerous signaling cascades during craniofacial embryonic development (discussed below and in more detail in the Supplemental Appendix S3). Additional 7 significantly-associated SNPs are found in intergenic regions adjacent to several loci, such as BMP4, LOC124685, LOC100131241 and PTCH1. Some of these genes and chromosomal segments were previously linked to craniofacial embryogenesis, while others represent novel associations.

Six genetic variants were found in lncRNA genes, which have not been previously linked to craniofacial morphogenesis. These findings may suggest an yet unexplored level of epigenetic regulation affecting craniofacial morphology. lncRNAs represent a recently discovered class of regulatory factors, whose expression is thought to be important for regulation of gene expression through direct competition with transcription factors in promoter regions, as well as by indirectly affecting gene expression via interaction with miRNAs, histone modifying complexes and other cellular factors [81,82]. The role of these ncRNA molecules in the craniofacial development clearly requires additional research.

Another interesting example of potential epigenetic regulation is the gene EYA1, which encodes a histone phosphatase protein and transcriptional coactivator for SIX gene family. This factor regulates transcription during organogenesis of various craniofacial domains, such as branchial arches, eye and ear and known being mutated in Otofaciocervical Syndrome 1 and Branchiootorenal Syndrome 1 (discussed in detail in the Supplemental Appendix S3). Based on ENCODE and Regulome analysis, at least two of the significantly associated SNPs in the intronic region of this gene (rs1481800 and rs73684719) are potentially functional variants, located distal to EYA1 transcription start site. These polymorphisms may essentially cis-regulate its expression, as well as trans-regulate expression of other craniofacial genes via chromatin modification mechanisms. The comprehensive role of epigenetic gene regulation in general, and in the craniofacial embryonic development in particular, is poorly understood. Our results support the importance of epigenetic regulation behind normal craniofacial morphology and require further research.

Notably, there is a limited number of recent studies showing existence of thousands of enhancer sequences, potentially active in the developing craniofacial complex in mice and in humans [69,83]. Given the extremely high complexity of the craniofacial morphology, its developmental regulation and consequent phenotypic variety is likely influenced by both cis and trans epistatic interactions between genetic and epigenetic regulatory factors, which require further investigation.

Even though a relatively large number of phenotypes were studied (92 linear and angular measurements and indices), this may still represent an oversimplification of the complexity of the human face. Despite the importance of the association between specific 3D measurements and SNPs demonstrated in this study, the association of facial shapes, represented by the principal components might better represent facial morphology complexity. However, only one of the 10 principle components analysed in this study, showed significant associations at the GWAS threshold level. While the explanation of this observation is unclear, it is consistent with other similar studies [22,23]. The specific anthropometric measurements on the other hand, produced numerous significant associations, identifying many genes and intergenic regions that appear to play important roles in the development of normal human facial appearance.

We also note the relatively small sample size used in this study, which led to a major limitation: lack of replication of these results in a separate set of samples. However, the confirmation of the two previously associated genes (PAX3 and COL11A1) supports the validity of our findings.

Conclusions

The high complexity of facial anatomy, and the composite nature of genetic regulation that affects its development, requires novel comprehensive approaches of capturing facial morphology. A number of such methods has recently revealed additional genes with a few polymorphisms associated with the normal craniofacial variation range in the Caucasian population [82,80,84]. Further studies may involve the use of similar or even more advanced methods to capture the majority of craniofacial variation and benefit from larger and more diverse population samples.

A recent attempt to predict facial appearance was performed using only 24 SNPs [85]. This approach has promise, although it is largely based on reconstruction of a 'facial composite image' via prediction of ancestry, sex, pigmentation and human perception of faces. While absolutely reasonable, this approach highlights the need for additional association studies looking at specific craniofacial traits. Genetic association studies of a large scope of anthropometric measurements and other facial phenotypes are essential to uncover the dense genetic and epigenetic landscape, affecting normal facial apperance and may subsequently being useful in its prediction.

Additional association studies utilizing large sample sizes and incorporating dense SNPs and DNA methylation panels or even whole genome sequencing approaches, in conjunction with morphologically adequate representation of the craniofacial characteristics would be a valuable adjunct to the promising results obtained in this study. These studies would likely require extensive collaboration between research groups and greatly improve our understanding of the dynamic gene expression influencing craniofacial morphology, which in turn enable better prediction of the visual appearance of a person from a DNA sequence.

Acknowledgments

We would like to thank the volunteers who participated in this study without whom we could not have performed this research. We would like to thank Olga Kondrashova who helped with the sample collection. We also thank Technical Support Working Group (TSWG) and Pelerman Holdings Pt Ltd. for their generous support of this project.

Ethics and Consent to Participate

The participants provided their written informed consent to participate in this study, which was approved by the Bond University Ethics committee (RO-510).

Competing Interests

The authors declare that they have no competing interests.

Authors' Contributions

MB designed the study, carried out the molecular genetic studies, carried out the data analysis, participated in the statistical analysis and drafted the manuscript. PB performed the statistical analysis and drafted the manuscript. AvD participated in the design of the study and drafted the manuscript. All authors read and approved the final manuscript.

Availability of Data and Materials

The data supporting the conclusions of this article are included within the article and its additional files.

Funding

The funding for this research was provided by the Technical Support Working Group (Award Number: IS-FB-2946) and Pelerman Holdings Pt Ltd.

References

1. Kohn LAP (1991) The Role of Genetics in Craniofacial Morphology and Growth. Annual Review of Anthropology 20: 261-278.
2. Sperber GH, Sperber SM, Guttmann GD (2010) Craniofacial embryogenetics and development. (2^nd edn). Shelton, CT: People's Medical Pub. House USA.
3. Richtsmeier JT, Cheverud JM (1986) Finite element scaling analysis of human craniofacial growth. J Craniofac Genet Dev Biol 6: 289-323.
4. Neubauer S, Gunz P, Hublin JJ (2009) The pattern of endocranial ontogenetic shape changes in humans. J Anat 215: 240-255.
5. Sturm RA (2009) Molecular genetics of human pigmentation diversity. Hum Mol Genet 18: R9-R17.
6. Shkoukani MA, Chen M, Vong A (2013) Cleft Lip - A Comprehensive Review. Front Pediatr 1: 53.
7. Kimonis V, Gold JA, Hoffman TL, Panchal J, Boyadjiev SA (2007) Genetics of Craniosynostosis. Semin Pediatr Neurol 14: 150-161.
8. Weinberg SM, Naidoo SD, Bardi KM, Brandon CA, Neiswanger K, et al. (2009) Face shape of unaffected parents with cleft affected offspring: combining three-dimensional surface imaging and geometric morphometrics. Orthod Craniofac Res 12: 271-281.
9. Coussens AK, Wilkinson CR, Hughes IP, Morris CP, van Daal A, et al. (2007) Unravelling the molecular control of calvarial suture fusion in children with craniosynostosis. BMC Genomics 8: 458.
10. Coussens AK, Hughes IP, Wilkinson CR, Morris CP, Anderson PJ, et al. (2008) Identification of genes differentially expressed by prematurely fused human sutures using a novel in vivo -in vitro approach. Differentiation 76: 531-545.
11. Boehringer S, van der Lijn F, Liu F, Gunther M, Sinigerova S, et al. (2011) Genetic determination of human facial morphology: links between cleft-lips and normal variation. Eur J Hum Genet 19: 1192-1197.
12. Nakamura A, Hattori M, Sakaki Y (1997) A novel gene isolated from human placenta located in Down syndrome critical region on chromosome 21. DNA Res 4: 321-324.
13. El Ghouzzi V, Legeai-Mallet L, Benoist-Lasselin C, Lajeunie E, Renier D, et al. (2001) Mutations in the basic domain and the loop-helix II junction of TWIST abolish DNA binding in Saethre-Chotzen syndrome. FEBS Lett 492: 112-118.
14. Kamath BM, Stolle C, Bason L, Colliton RP, Piccoli DA, et al. (2002) Craniosynostosis in Alagille syndrome. Am J Med Genet 112: 176-180.
15. Hood RL, Lines MA, Nikkel SM, Schwartzentruber J, Beaulieu C, et al. (2012) Mutations in SRCAP, encoding SNF2-related CREBBP activator protein, cause Floating-Harbor syndrome. Am J Hum Genet 90: 308-313.
16. Roper RJ, Baxter LL, Saran NG, Klinedinst DK, Beachy PA, et al. (2006) Defective cerebellar response to mitogenic Hedgehog signaling in Down's syndrome mice. Proc Natl Acad Sci U S A 103: 1452-1456.
17. Croonen EA, van der Burgt I, Kapusta L, Draaisma JM (2008) Electrocardiography in Noonan syndrome PTPN11 gene mutation--phenotype characterization. Am J Med Genet A 146: 350-353.
18. Fourie Z, Damstra J, Gerrits PO, Ren Y (2011) Evaluation of anthropometric accuracy and reliability using different three-dimensional scanning systems. Forensic Sci Int 207: 127-134.
19. Toma AM, Zhurov A, Playle R, Ong E, Richmond S (2009) Reproducibility of facial soft tissue landmarks on 3D laser-scanned facial images. Orthod Craniofac Res 12: 33-42.
20. Kovacs L, Zimmermann A, Brockmann G, Baurecht H, Schwenzer-Zimmerer K, et al. (2006) Accuracy and Precision of the Three-Dimensional Assessment of the Facial Surface Using a 3-D Laser Scanner. IEEE Trans Med Imaging 25: 742-754.
21. Coussens AK, van Daal A (2005) Linkage disequilibrium analysis identifies an FGFR1 haplotype-tag SNP associated with normal variation in craniofacial shape. Genomics 85: 563-573.
22. Paternoster L, Zhurov AI, Toma AM, Kemp JP, St Pourcain B, et al. (2012) Genome-wide association study of three-dimensional facial morphology identifies a variant in PAX3 associated with nasion position. Am J Hum Genet 90: 478-485.
23. Liu F, Van der Lijn F, Schurmann C, Zhu G, Chakravarty MM, et al. (2012) A genome-wide association study identifies five loci influencing facial morphology in Europeans. PLoS Genet 8: e1002932.
24. Michel S, Liang L, Depner M, Klopp N, Ruether A, et al. (2010) Unifying candidate gene and GWAS Approaches in Asthma. PLoS One 5: e13894.
25. Sun J, Jia P, Fanous AH, Webb BT, van den Oord EJ, et al. (2009) A multi-dimensional evidence-based candidate gene prioritization approach for complex diseases-schizophrenia as a case. Bioinformatics 25: 2595-6602.
26. Hrdy DB (1978) Analysis of hair samples of mummies from Semma South (Sudanese Nubia). Am J Phys Anthropol 49: 277-282.
27. Fitzpatrick TB (1988) The validity and practicality of sun-reactive skin types I through VI. Archives of Dermatology 124: 869.
28. Sturm RA, Larsson M (2009) Genetics of human iris colour and patterns. Pigment Cell Melanoma Res 22: 544-562
29. Farkas LG (1994) Anthropometry of the head and face. (2^nd edn). New York: Raven Press.
30. Website MMRW. The Jackson Laboratory, Bar Harbor, Maine.
31. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, et al. (2000) Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25: 25-29.
32. Online Mendelian Inheritance in Man O. McKusick-Nathans Institute of Genetic Medicine, Johns Hopkins University (Baltimore, MD).
33. Rebhan M, ChalifaCaspi V, Prilusky J, Lancet D (1997) GeneCards: Integrating information about genes, proteins and diseases. Trends in Genetics 13: 163.
34. Sabeti PC, Varilly P, Fry B, Lohmueller J, Hostetter E, et al. (2007) Genome-wide detection and characterization of positive selection in human populations. Nature 449: 913-918.
35. Zhou N, Wang L (2007) Effective selection of informative SNPs and classification on the HapMap genotype data. BMC Bioinformatics 8: 484.
36. Kosoy R, Nassir R, Tian C, White PA, Butler LM, et al. (2009) Ancestry informative marker sets for determining continental origin and admixture proportions in common populations in America. Hum Mutat 30: 69-78.
37. Paschou P, Lewis J, Javed A, Drineas P (2010) Ancestry informative markers for fine-scale individual assignment to worldwide populations. J Med Genet 47: 835-847.
38. Londin ER, Keller MA, Maista C, Smith G, Mamounas LA, et al. (2010) CoAIMs: a cost-effective panel of ancestry informative markers for determining continental origins. PLoS One 5: e13443.
39. Bulbul O, Filoglu G, Altuncul H, Aradas AF, Ruiz Y, et al. (2011) A SNP multiplex for the simultaneous prediction of biogeographic ancestry and pigmentation type. Forensic Sci Inter: Genet Suppl 3: e500-e501.
40. Tandon A, Patterson N, Reich D (2011) Ancestry informative marker panels for African Americans based on subsets of commercially available SNP arrays. Genet Epidemiol 35: 80-83.
41. Phillips C, Freire Aradas A, Kriegel AK, Fondevila M, Bulbul O, et al. (2013) Eurasiaplex: a forensic SNP assay for differentiating European and South Asian ancestries. Forensic Sci Int Genet 7: 359-366.
42. Gettings KB, Lai R, Johnson JL, Peck MA, Hart JA, et al. (2014) A 50-SNP assay for biogeographic ancestry and phenotype prediction in the U.S. population. Forensic Sci Int Genet 8: 101-108.
43. Kidd KK, Speed WC, Pakstis AJ, Furtado MR, Fang R, et al. (2014) Progress toward an efficient panel of SNPs for ancestry inference. Forensic Sci Int Genet 10: 23-32.
44. Amigo J, Salas A, Phillips C (2011) ENGINES: exploring single nucleotide variation in entire human genomes. BMC Bioinformatics 12: 105.
45. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, et al. (2010) A method and server for predicting damaging missense mutations. Nat Methods 7: 248-249.
46. Marinescu V, Kohane I, Riva A (2005) MAPPER: a search engine for the computational identification of putative transcription factor binding sites in multiple genomes. BMC Bioinformatics 6: 79.
47. Wang J, Ronaghi M, Chong SS, Lee CGL (2011) pfSNP: An integrated potentially functional SNP resource that facilitates hypotheses generation through knowledge syntheses. Human Mutation 32: 19-24.
48. Coassin S, Brandstätter A, Kronenberg F (2010) Lost in the space of bioinformatic tools: A constantly updated survival guide for genetic epidemiology. The GenEpi Toolbox. Atherosclerosis 209: 321-335.
49. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945-959.
50. Thorisson G, Smith A, Krishnan L, Stein L (2005) The International HapMap Project Web site. Genome Res 15: 1592.
51. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, et al. (2007) PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81: 559-575.
52. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, et al. (2006) Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet 38: 904-909.
53. Patterson N, Price AL, Reich D (2006) Population structure and eigenanalysis. PLoS Genet 2: e190.
54. UniProt Consortium (2014) Activities at the Universal Protein Resource (UniProt). Nucleic Acids Res 42: 191-198.
55. Rappaport N, Nativ N, Stelzer G, Twik M, Guan-Golan Y, et al. (2013) MalaCards: an integrated compendium for diseases and their annotation. Database (Oxford) 18.
56. Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, et al. (2010) The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res 38: 214-220.
57. Reference Genome Group of the Gene Ontology Consortium (2009) The Gene Ontology's Reference Genome Project: a unified framework for functional annotation across species. PLoS Comput Biol 5: e1000431.
58. Sherry ST, Ward MH, Kholodov M, Baker J, Phan L, et al. (2001) dbSNP: the NCBI database of genetic variation. Nucleic Acids Research 29: 308-311.
59. Durbin RM, Abecasis GR, Altshuler DL, Auton A, Brooks LD, et al. (2010) A map of human genome variation from population-scale sequencing. Nature 467: 1061-1073.
60. Chelala C, Khan A, Lemoine NR (2009) SNPnexus: a web database for functional annotation of newly discovered and public domain single nucleotide polymorphisms. Bioinformatics 25: 655-661.
61. Rajeevan H, Cheung KH, Gadagkar R, Stein S, Soundararajan U, et al. (2007) ALFRED: an allele frequency database for microevolutionary studies. Evol Bioinform Online 1: 1-10.
62. Johnson AD, Handsaker RE, Pulit SL, Nizzari MM, O'Donnell CJ, et al. (2008) SNAP: a web-based tool for identification and annotation of proxy SNPs using HapMap. Bioinformatics 24: 2938-2939.
63. Boyle AP, Hong EL, Hariharan M, Cheng Y, Schaub MA, et al. (2012) Annotation of functional variation in personal genomes using RegulomeDB. Genome Res 22: 1790-1797.
64. MacArthur J, Bowler E, Cerezo M, Gil L, Hall P, et al. (2017) The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog). Nucleic Acids Res 45: 896-901.
65. Visscher PM (2008) Sizing up human height variation. Nature Genetics 40: 489-490.
66. Frazer KA, Murray SS, Schork NJ, Topol EJ (2009) Human genetic variation and its contribution to complex traits. Nat Rev Genet 10: 241.
67. Wood AR, Esko T, Yang J, Vedantam S, Pers TH, et al. (2014) Defining the role of common variation in the genomic and biological architecture of adult human height. Nat Genet 46: 1173-1186.
68. Hallgrimsson B, Mio W, Marcucio RS, Spritz R (2014) Let's Face It - Complex Traits Are Just Not That Simple. PLoS Genet 10: e1004724.
69. Attanasio C, Nord AS, Zhu Y, Blow MJ, Li Z, et al. (2013) Fine tuning of craniofacial morphology by distant-acting enhancers. Science 342: 1241006.
70. Tak YG, Farnham PJ (2015) Making sense of GWAS: using epigenomics and genome engineering to understand the functional relevance of SNPs in non-coding regions of the human genome. Epigenetics Chromatin 8: 57.
71. Weston EM, Friday AE, Lio P (2007) Biometric evidence that sexual selection has shaped the hominin face. PLoS One 2: e710.
72. Panagiotou OA, Ioannidis JP, Genome-Wide Significance P (2012) What should the genome-wide significance threshold be? Empirical replication of borderline genetic associations. Int J Epidemiol 41: 273-286.
73. Dudbridge F, Gusnanto A (2008) Estimation of significance thresholds for genome wide association scans. Genet Epidemiol 32: 227-234.
74. Cooper RS, Tayo B, Zhu X (2008) Genome-wide association studies: implications for multiethnic samples. Hum Mol Genet 17: 151-155.
75. Barnholtz-Sloan JS, Chakraborty R, Sellers TA, Schwartz AG (2005) Examining population stratification via individual ancestry estimates versus self-reported race. Cancer Epidemiol Biomarkers Prev 14: 1545-1551.
76. Tang H, Quertermous T, Rodriguez B, Kardia SL, Zhu X, et al. (2005) Genetic structure, self-identified race/ethnicity, and confounding in case-control association studies. Am J Hum Genet 76: 268-275.
77. Price AL, Zaitlen NA, Reich D, Patterson N (2010) New approaches to population stratification in genome-wide association studies. Nat Rev Genet 11: 459-463.
78. Zuk O, Hechter E, Sunyaev SR, Lander ES (2012) The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci U S A 109: 1193-1198.
79. Makowsky R, Pajewski NM, Klimentidis YC, Vazquez AI, Duarte CW, et al. (2011) Beyond Missing Heritability: Prediction of Complex Traits. PLOS Genetics 7: e1002051.
80. Claes P, Liberton DK, Daniels K, Rosana KM, Quillen EE, et al. (2014) Modeling 3D facial shape from DNA. PLoS Genet 10: e1004224.
81. Lau E (2014) Non-coding RNA: Zooming in on lncRNA functions. Nat Rev Genet 15: 574-575.
82. Peschansky VJ, Wahlestedt C (2014) Non-coding RNAs as direct and indirect modulators of epigenetic regulation. Epigenetics 9: 3-12.
83. Wu H, Nord AS, Akiyama JA, Shoukry M, Afzal V, et al. (2014) Tissue-Specific RNA Expression Marks Distant-Acting Developmental Enhancers. PLoS Genet 10: e1004610.
84. Peng S, Tan J, Hu S, Zhou H, Guo J, et al. (2013) Detecting Genetic Association of Common Human Facial Morphological Variation Using High Density 3D Image Registration. PLoS Comput Biol 9: e1003375.
85. Claes P, Hill H, Shriver MD (2014) Toward DNA-based facial composites: preliminary results and validation. Forensic Sci Int Genet 13: 208-216.