Orofacial clefting (OFC) is one of the most frequent congenital malformations with a substantial functional, financial, and mental health burden to the persons and families affected (Wehby and Cassell, 2010). Non-syndromic cleft lip with or without cleft palate (NSCL/P) has a multifactorial etiology, with both environmental and genetic factors playing roles (Mangold et al., 2011). To date, approximately 40 genetic risk loci have been identified for non-syndromic clefts through genome-wide association studies (GWAS; Thieme and Ludwig, 2017). Despite the substantial number of loci discovered, the sample sizes in these studies are relatively low when compared to GWA-studies of schizophrenia (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014) or height (Wood et al., 2014), which include over 100K individuals. Therefore, it is likely that additional NSCL/P risk loci are yet to be discovered. Furthermore, a thorough understanding of the roles that these loci play in development, and thus the etiology of OFC, is still lacking.

A comprehensive approach to studying the etiology of NSCL/P requires investigations from different angles. One approach is to study subclinical phenotypes that exist within the range of normal development and are associated with OFC. Such traits are referred to as endophenotypes, which can be present in the biological relatives of individuals with NSCL/P and are potentially an incomplete expression of susceptibility loci for OFC (Gottesman and Gould, 2003; Weinberg et al., 2006). A number of endophenotypes have been described for NSCL/P, for example, a mild disruption of the upper lip musculature (Neiswanger et al., 2007), the appearance of “whorls” on the surface of the lower lip (Neiswanger et al., 2009), dental anomalies (Howe et al., 2015), and altered facial shape. In particular, unaffected relatives have been shown to display midfacial retrusion and a broadening of the upper face versus comparison control groups (Weinberg et al., 2009; Roosenboom et al., 2015).

Extensive characterization of the NSCL/P facial endophenotype allows us to investigate whether the presence of facial endophenotypes indicates heightened genetic susceptibility for NSCL/P. The possibility that NSCL/P risk variants influence normal-range facial variation raises interesting questions about the intersection of the genetic architecture between normal and dysmorphic facial features. Some NSCL/P candidate genes have already been shown to impact normal facial morphology. Boehringer et al. (2011) found associations between two genetic loci, rs1258763 and rs987525, involved in NSCL/P and normal craniofacial traits, namely nose width and bizygomatic distance, respectively. In previous work, we found that six known NSCL/P risk variants affected nose, chin, eyebrow ridge, and philtrum 3D shape in a large unselected population (Indencleef et al., 2018). Moreover, in a recent GWAS of normal-range facial morphology in Europeans, there was strong evidence of enrichment for genes implicated in OFC and palate development (Claes et al., 2018).

In this study, we investigate the intersection of the genetic architectures of normal facial variation and NSCL/P, starting from explicitly defined NSCL/P facial endophenotypes. First, we quantify the NSCL/P facial endophenotype at multiple levels of scale, modeling both global and local facial shape variations in a sample of unaffected parents of individuals with NSCL/P (Claes et al., 2018). Next, we identify genetic loci affecting these discovered facial endophenotypic traits in a GWAS of individuals from the general population. Finally, we examine the extent to which the genetic architecture of the facial endophenotypic traits reflects a heightened susceptibility for OFC.

Facial shape in 59 out of 63 segments was found to be different between unaffected parents of individuals with NSCL/P and the unselected individuals at nominal significance (p < 0.05), as shown in Figure 2. In 52 segments, a significant difference was found below the Bonferroni-corrected threshold (p < 1.3 × 10^–3). Each shape regression returned a direction of effect, represented by a vector in multidimensional shape space, which we refer to as the endophenotypic trait. Endophenotypic traits that were identified in segments 1, 2, 6, and 7 are shown in Figure 3, respectively, covering the full face, the midface, the mandibular region, and the upper facial region. These endophenotypic traits illustrate a midfacial retrusion and a broadening of the upper face, a more protruded chin, and an upturned nose. All 59 global-to-local traits are visualized in Supplementary Figure 2.

We identified 29 independent genetic loci with a lead SNP associated with at least one of the 59 endophenotypic traits below a genome-wide significance threshold (p < 5 × 10^–8) (Figure 4). In total, 36 associations were significant at a genome-wide level as some of the 29 lead SNPs were associated with more than one trait. Of these 36 associations, 16 included traits in quadrant II (nose, Figure 1), 12 in quadrant I (philtrum), 5 in quadrant III (mandible), 2 in quadrant IV (upper face), and lastly 1 association in segment 2 (midface, Figure 5). Nine out of 29 lead SNPs also passed the study-wide threshold (p < 8.47 × 10^–10). SNP information and meta-analysis statistics of the 29 best associations (lowest p-value of all associated traits) are given in Table 2. Supplementary Table 1 contains additional information on the regression statistics and information for all 36 significant associations. Twenty-two of the lead SNPs lie in proximity (±500 kb) to lead SNPs identified in a prior normal facial variation GWAS (White et al., 2020a). Sixteen lead SNPs showed strong evidence for a role in risk for OFC, another 12 showed weaker evidence (Supplementary Table 2).

Associations of the 59 endophenotypic traits with polygenic risk scores calculated across different p-value inclusion thresholds are visualized in Figure 6. In the United States cohort, the strongest associations occurred in polygenic risk scores calculated with inclusion thresholds around 0.1. The endophenotypic traits showing the strongest associations are present in the forehead, eye, and nose tip regions (Figure 6). The strongest associations in the United Kingdom cohort were found between the endophenotypic traits and polygenic risk scores calculated using a p-value inclusion threshold around 0.001. The traits responsible for these associations are localized in the philtrum area (Figure 6).

Orofacial clefting has a wide range of expression and severity (Marazita, 2012). CL/P can occur bilaterally or unilaterally and is generally considered distinct from cleft palate only (CPO), based on embryological and epidemiological considerations (Roosenboom et al., 2017). The range of severity also includes minor defects and subclinical phenotypes, which are on the mild side of the severity spectrum. One of those subclinical phenotypes includes subtle alterations in facial shape—a facial endophenotype indicative of NSCL/P risk. Unaffected members of families with a history of OFC have been found to show a broadening of the upper face and midfacial retrusion in comparison to control groups with no family history of clefting (Weinberg et al., 2009; Roosenboom et al., 2015, 2017). These endophenotypic features were confirmed in this work, where unaffected parents also showed an upturned nose and a protrusion of the chin.

Although several studies have confirmed the presence of a facial endophenotype for NSCL/P in unaffected relatives, a biological cause during embryonic development remains speculative. Lip formation happens around week 6 of gestation, and CL/P occurs when the nasal and maxillary prominences fail to fuse (Marazita, 2012). This could be the result of a failure of prominence growth and/or alignment due to the shape of the face and the direction of growth of the prominences. Mouse strains at higher risk for spontaneous clefting have been shown to exhibit unique patterns in the orientation of the facial prominences forming the lip and palate (Trasler, 1968; Juriloff and Trasler, 1976; Millicovsky et al., 1982; Young et al., 2007; Parsons et al., 2008). Mouse studies including mutants in TFAP2A [MIM: 107580], a gene implicated in CL/P, have suggested small gene expression changes altering the growth and shape of the facial prominences as a possible cause for the development of a cleft (Green et al., 2015). Accordingly, cleft risk variants might influence the shape of embryological structures in a way that complicates proper alignment for fusion. If fusion succeeds however, the altered shape could result in the facial endophenotype. This idea is supported by our endophenotype GWAS, in which the lead SNPs were predominantly associated with traits in the nasal and philtrum area (80% of the genome-wide associations were with traits in this area), the regions developed from the nasal and maxillary prominences. Thus, there is a clear intersection between normal facial development and shape alterations potentially leading to malformations such as CL/P.

In our GWAS of endophenotypes, 22 of the 29 lead SNPs among associated loci were found in vicinity (±500 kb) of the lead SNPs from a previous GWAS of normal facial variation (White et al., 2020a). The overlap of the endophenotypes with normal facial variation is substantial and expected. The endophenotypic traits were established by analyzing unaffected parents’ facial shape, yet tested in a sample of normal facial variation where a family history of cleft was absent or unknown. An intersection between the genetic pathways of NSCL/P and normal facial variation was already implicated by the influence of NSCL/P risk variants on normal facial variation. Boehringer et al. found associations between two NSCL/P loci and bizygomatic distance and nose width (Boehringer et al., 2011). In previous work, we found six NSCL/P risk variants affecting nose, chin, eyebrow ridge, and philtrum 3D shape in the larger unselected population (Indencleef et al., 2018). One of those six variants, rs6740960, has been shown to play a role in both NSCL/P and normal facial variation in multiple studies. In a normal facial variation GWAS conducted by our group, rs6740960 was associated with shape in the mandibular area (Claes et al., 2018). The same variant was associated with NSCL/P in a multi-ethnic meta-analysis (Ludwig et al., 2017). In this work, rs4952552 (17 kb from rs6740960) was associated with segment 43, which covers small regions on the outer sides of the nostrils. Our results indicate a strong overlap in the underlying genetics of the NSCL/P endophenotype and normal facial variation, with rs6740960 as an example of this intersection.

The facial endophenotype is hypothesized to be an expression of NSCL/P susceptibility genes. In this work, 16 of the 29 lead SNPs also have substantial evidence for a role in OFC and 12 other lead SNPs show weaker evidence (Supplementary Table 2). Of special interest are the seven lead SNPs implicated with OFC that did not show an overlap with normal facial variation, which will be discussed here. One of these variants is rs56399962, an SNP located in 1p31.3. A clinical report covering a patient with a deletion in 1p31 described a unique form of OFC manifesting as a facial midline defect in both the lower and upper lip (Y1ld1r1m et al., 2014). In a family based study, rs787499, an SNP 87 kb from rs56399962, was associated with NSCL/P (Camargo et al., 2012). This variant is also 404 kb from WLS [MIM: 611514], of which expressed mRNA and protein have been found in the embryonic development of the lip and palate in mice (Yu et al., 2010). WLS deletion in mice has resulted in abnormal palate development (Fu et al., 2011). In zebrafish, this gene was shown to be involved in palate morphogenesis (Rochard et al., 2016). These findings show that rs56399962 is located in a possible risk locus for NSCL/P, with WLS as a candidate gene. Another SNP associated with NSCL/P endophenotypes in our study, rs2578135, is 169 kb from a transcript variant of GRID2 [MIM: 602368], a gene that showed suggestive evidence (p-value = 7.9 × 10^–7) with NSCL/P in a GWAS with an Asian cohort (Beaty et al., 2010). GRID2 has also shown evidence for a gene-smoking interaction effect for CL/P (Beaty et al., 2013). Another SNP, rs57106444, which is 250 kb from rs831227, showed a suggestive association (p-value = 5.34 × 10^–7) in a POFC multi-ethnic GWAS for NSCL/P (Leslie et al., 2016). PIK3R1 [MIM: 171833], a gene 168 kb from rs57106444, was differentially expressed in association with NSCL/P (Wang et al., 2016). F13A1 [MIM: 134570] is a gene located 72 kb from rs2326766 in 6p25.1, which was associated with an endophenotypic trait in the philtrum. This gene was suggested as a candidate gene in a NSCL/P linkage study, but this result has not been replicated (Eiberg et al., 1987). Although 6p23 has been a candidate locus for NSCL/P in various studies, five individuals with isolated 6p25 deletions have been reported with a high or cleft palate (Carinci et al., n.d.; Scapoli et al., 1997; Linhares et al., 2015). The variant identified in 15q22.2 (rs2326766) lies in the same chromosomal band but still 5MB from a SNP that was suggestively associated with NSCL/P (de Aquino et al., 2014). Another variant, rs113689916, located in 17p13.3 and associated with an endophenotypic trait in the philtrum, lies within a microduplication site for which the syndromic phenotype shows CL/P (Curry et al., 2013; Tucker and Escobar, 2014; Ibitoye et al., 2015). Lastly, rs3091552 is associated with an endophenotypic trait in the nose and showed a suggestive association (p = 2.5 × 10^–5) in a previous NSCL/P GWAS (Leslie et al., 2017). The SNP is located in 20q13.12, the duplication of which has caused cleft lip/palate in two clinical cases (Blanc et al., 2008). One familial case of Galloway-Mowat syndrome [MIM: 617730] due to a novel TP53RK [MIM: 608679] mutation, a gene within 500 kb of rs3091552, also showed a cleft palate (Hyun et al., 2018). Although these seven variants have not yet been significantly associated with either normal facial variation or NSCL/P in a GWAS, these results imply that they are new candidate risk variants for NSCL/P in a European population, with a potential additional role in normal facial variation.

The strong overlap of the endophenotypic lead SNPs and genetics of OFC indicates that the endophenotype is an expression of NSCL/P susceptibility genes and that our quantification method has captured this heightened susceptibility successfully. In further support of this observation, we have provided additional evidence by calculating polygenic risk scores for NSCL/P and subjecting them to an association with the endophenotypic traits. In both the United States and United Kingdom cohorts, endophenotypic traits showed associations with the polygenic risk scores, indicating an overall heightened susceptibility for NSCL/P. Although associations were found in both United States and United Kingdom cohorts, the specific traits associated were different. The United Kingdom cohort showed associations in the philtrum, an important area in lip formation during development. In the same cohort, Howe et al. (2018) found an association between philtrum width and an increased PRS for NSCL/P. These results indicate that variants affecting the event of lip fusion can also influence philtrum shape. The United States cohort showed associations in the orbital region and in the nasal tip. During embryological development, the nasal prominence is involved in lip formation. Even though the forehead region does not seem to be directly involved in lip formation, it has previously been identified as a region involved in the facial endophenotype (Weinberg et al., 2009; Roosenboom et al., 2015). The discordances of traits associated with the PRS scores between the United States and United Kingdom cohorts could be present due to a couple of reasons. First, in contrast to the United States, the United Kingdom cohort only includes individuals of 15 years old, which do not have fully developed facial features. Second, the summary statistics used to calculate the PRS are based on United States individuals. Although all are of European ancestry, genetic differences between populations of United States and United Kingdom are likely to exist and influence the polygenic risk scores predicted. Lastly, the sample sizes for calculating PRS are relatively small, and the summary statistics were based on a GWAS of N = 1,411 individuals.

In this work, we tested a range of inclusion thresholds for the polygenic risk scores. However, the common practice is to determine an optimal inclusion threshold in another independent dataset (Dudbridge, 2013), which was not available in this work. For most polygenic traits, the more the SNPs included in the polygenic risk model, the higher the predictive power (Dudbridge, 2013). This often leads to relatively high optimal inclusion thresholds such as p < 0.5. In a study by Ludwig et al. (2017), the optimum inclusion threshold for NSCL/P polygenic risk scores was p < 0.2. In contrast, Howe et al. (2018) found P < 10^–5 to be the optimum inclusion p-value and thus included a lower number of SNPs. Consequently, the optimum p-value inclusion threshold seems to be strongly dependent on the trait as well as the datasets used. Linkage disequilibrium score regression is a method in which genetic correlations between different traits are estimated, yet unfortunately our sample sizes did not meet the criteria for this analysis.

In this study, we focused our analyses on individuals of European ancestry. To analyze possible population similarities and differences in both OFC genetics as well as normal facial genetics, future analyses should focus on individuals of other ancestries. For example in the context of NSCL/P GWA-studies, evidence for the gene desert 8q24 is consistently strong in European ancestry groups, but less so in Asian ancestry groups (Beaty et al., 2016). Inclusion of more diverse populations in human genetics studies is of great importance, because the current bias for European ancestry can translate into a poor disease prediction and treatment for individuals of under-represented ancestries (Sirugo et al., 2019). Another limitation of this study is the lack of additional datasets to replicate the genetic associations found, as well as the polygenic risk scores. More specifically in the PRS analysis, additional datasets would be most beneficial to identify the optimum p-value inclusion threshold. A third limitation is that the unselected individuals are extracted from different databases in which extensive filtering was necessary. Additionally, only the unselected individuals included from our collaboration with Pittsburgh were checked for a family history of OFC, and excluded if it was indeed so. Lastly, no functional analyses were performed to further investigate the role of the loci discovered in OFC. This would be of great interest for future analyses.

This study is the first to identify facial endophenotypic traits for NSCL/P based on a comprehensive global-to-local facial phenotyping approach. We have used these traits to examine their underlying genetic architecture and observed a substantial link between the endophenotypic traits and NSCL/P genetics, with 16 of our lead SNPs possessing strong literature evidence for an association with NSCL/P and an additional 12 SNPs having weaker literature evidence. We identified seven variants as possible new candidate risk loci for NSCL/P. Furthermore, endophenotypic traits also indicated an overall heightened NSCL/P susceptibility based on polygenic risk scores. This study delivers additional evidence for the existence of the facial endophenotype as an expression of NSCL/P susceptibility genes. Our study shows that 76% of our lead SNPs are also associated with normal facial variation. Therefore, and in conclusion, this study contributes to the understanding of the intersection of NSCL/P genetics and the genetic architecture of normal facial variation.

Publicly available datasets were analyzed in this study. The data can be found here: All of the genotypic markers for the 3DFN dataset are available to the research community through the dbGaP controlled-access repository (http://www.ncbi.nlm.nih.gov/gap) at accession #phs000929.v1.p1. The raw source data for the phenotypes—the 3D facial surface models in.obj format—are available through the FaceBase Consortium (https://www.facebase.org) at accession #FB00000491.01. Access to these 3D facial surface models requires proper institutional ethics approval and approval from the FaceBase data access committee. Additional details can be requested from SMW. The participants making up the PSU and IUPUI cohorts were not collected with broad data sharing consent. Given the highly identifiable nature of both facial and genomic information and unresolved issues regarding risk to participants, we opted for a more conservative approach to participant recruitment. Broad data sharing of the raw data from these collections would thus be in legal and ethical violation of the informed consent obtained from the participants. This restriction is not because of any personal or commercial interests. Additional details can be requested from MS and SW for the PSU and IUPUI datasets, respectively. The ALSPAC (United Kingdom) data will be made available to bona fide researchers on application to the ALSPAC Executive Committee (http://www.bris.ac.uk/alspac/researchers/data-access). Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees. The ALSPAC study website contains details of all the data that are available through a fully searchable data dictionary and variable search tool: http://www.bris.ac.uk/alspac/researchers/our-data/.

The studies involving human participants were reviewed and approved by ALSPAC Ethics and Law Committee. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin. Written informed consent was obtained from the individual(s), and minor(s)’ legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

KI performed all analyses and wrote the first draft of the manuscript under the supervision of PC. PC, JS, SMW, and MM conceptualized the design of the study. JS, SMW, MM, RE, SW, JW, and MS organized and imputed the EURO data. ML, GW, LM, JH, RL, KC, FD, JS, MM, and SMW collected and curated the PARENT data. SR coordinated the collection of the ALSPAC images. SN, CL, ML, HH, and JW provided input throughout the analyses and the writing process. All authors contributed to manuscript revision, and read and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer AR-M declared a past co-authorship with one of the authors, GH, to the handling editor