Karlijne Indencleef1,2*
Jasmien Roosenboom3
Hanne Hoskens2,4
Julie D. White5
Mark D. Shriver5
Stephen Richmond6
Hilde Peeters4Eleanor Feingold3
Mary L. Marazita3,7
John R. Shaffer3,7
Seth M. Weinberg3,7
Greet Hens8
Peter Claes1,2*- 1Department of Electrical Engineering, ESAT/PSI, KU Leuven, Leuven, Belgium
- 2Medical Imaging Research Center, UZ Leuven, Leuven, Belgium
- 3Department of Oral Biology, Center for Craniofacial and Dental Genetics, School of Dental Medicine, University of Pittsburgh, Pittsburgh, PA, United States
- 4Department of Human Genetics, University Hospitals Leuven, Leuven, Belgium
- 5Department of Anthropology, The Pennsylvania State University, University Park, PA, United States
- 6Applied Clinical Research and Public Health, School of Dentistry, Cardiff University, College of Biomedical and Life Sciences, Heath Park, Cardiff, United Kingdom
- 7Department of Human Genetics, Graduate School of Public Health, University of Pittsburgh, Pittsburgh, PA, United States
- 8Department of Otorhinolaryngology, University Hospitals Leuven, Leuven, Belgium
Objectives: Orofacial clefting is one of the most prevalent craniofacial malformations. Previous research has demonstrated that unaffected relatives of patients with non-syndromic cleft lip with/without cleft palate (NSCL/P) show distinctive facial features, which can be an expression of underlying NSCL/P susceptibility genes. These results support the hypothesis that genes involved in the occurrence of a cleft also play a role in normal craniofacial development. In this study, we investigated the influence of genetic variants associated with NSCL/P on normal-range variation in facial shape.
Methods: A literature review of genome wide association studies (GWAS) investigating the genetic etiology of NSCL/P was performed, resulting in a list of 75 single nucleotide polymorphisms (SNPs) located in 38 genetic loci. Genotype data were available for 65 of these selected SNPs in three datasets with a combined sample size of 7,418 participants of European ancestry, whose 3D facial images were also available. The effect of each SNP was tested using a multivariate canonical correlation analysis (CCA) against 63 hierarchically-constructed facial segments in each of the three datasets and meta-analyzed. This allowed for the investigation of associations between SNPs known to be involved in NSCL/P and normal-range facial shape variations in a global-to-local perspective, without preselecting specific facial shape features or characteristics.
Results: Six NSCL/P SNPs showed significant associations with variation in normal-range facial morphology. rs6740960 showed significant effects in the chin area (p = 3.71 × 10−28). This SNP lies in a non-coding area. Another SNP, rs227731 near the NOG gene, showed a significant effect in the philtrum area (p = 1.96 × 10−16). Three SNPs showed significant effects on the shape of the nose. rs742071 (p = 8.71 × 10−14), rs34246903 (p = 6.87 × 10−12), and rs10512248 (p = 8.4 × 10−9). Respectively, these SNPs are annotated to PAX7, MSX1, and PTCH1. Finally, rs7590268, an intron variant of THADA, showed an effect in the shape of the supraorbital ridge (p = 3.84 × 10−7).
Conclusions: This study provides additional evidence NSCL/P-associated genetic variants influence normal-range craniofacial morphology, with significant effects observed for the chin, the nose, the supraorbital ridges and the philtrum area.
Introduction
Human facial features are highly variable and mostly genetically determined. Although craniofacial morphology is of interest to many scientists, its genetic architecture remains poorly understood (Roosenboom et al., 2016). Clinical dysmorphologists have hypothesized that genes responsible for syndromes with a distinctive facial phenotype are also involved in normal craniofacial development (Winter, 1996). Additionally, many loci found in genome wide association studies on human facial morphology were involved in syndromes affecting the face (Adhikari et al., 2016; Shaffer et al., 2016; Claes et al., 2018). Hence, investigating genes involved in craniofacial disorders constitutes a suitable approach to help unravel the genetic architecture of facial morphology.
With an incidence of 1.7 in 1,000 livebirths, non-syndromic cleft lip with or without cleft palate (NSCL/P) is one of the most common craniofacial anomalies, causing significant functional and psychological burden to the patient (Mossey et al., 2009). NSCL/P has a multifactorial etiology, which includes both genetic and environmental factors. Although the genetic background of NSCL/P is complex, significant progress has been made in the identification of NSCL/P susceptibility loci through genetic association studies (references listed in Table 1).
Table 1. Overview of lead-SNPs from the literature survey.
It has previously been hypothesized that NSCL/P genes can have an effect on normal facial morphology. Boehringer et al. (2011) found associations between genetic loci involved in NSCL/P and normal craniofacial traits, namely bizygomatic distance and nose width. Furthermore, several studies have found that unaffected relatives of patients with NSCL/P show distinctive facial characteristics in comparison to a control group, such as midfacial retrusion and broadening of the upper face, which could be defined as endophenotypes (Weinberg et al., 2008; Roosenboom et al., 2015).
In this study, we investigated the involvement of NSCL/P-associated genetic variants in normal-range facial variation. We used a candidate variant approach in combination with a new approach to study facial phenotypes based on spatially-dense, data-driven, global-to-local segmentations of facial 3D images. This segmentation approach was adopted from Claes et al. (2018) and offered two advantages. First, it allowed us to study the effects of candidate variants on facial shape in a hierarchical manner, providing complete coverage of the 3D facial surface at different levels of scale. Second, it allowed for an open-ended data-driven approach to establishing phenotypes, thereby avoiding preselection of phenotypic measurements. The latter aspect is in contrast to other cleft candidate gene studies (Boehringer et al., 2011; Liu et al., 2012) in which preselected phenotypic traits on a normal-range sample were studied.
Materials and Methods
Sample and Recruitment Details
Pittsburgh Dataset
Three datasets with 3D images and corresponding genetic data were used in this study. For the Pittsburgh sample, data were obtained from the 3D Facial Norms Database, which is a repository of 3D facial images and measurements (Weinberg et al., 2016). Participants were recruited in Pittsburgh, PA; Seattle, WA; Houston, TX; and Iowa City, IA. This dataset consists of 2,382 3D images with corresponding covariates of sex, age, weight, height, and genotype data. Participants ranged from 3 to 40 years old (mean age = 22; SD age = 9) and were of self-reported European ancestry. Individuals were excluded if they reported a personal or family history of any birth defect or syndrome affecting the head or face, a personal history of any significant trauma, surgery, or any medical condition that might alter the structure of the face. Of 2,382 participants with 3D images and covariate data, 42 were excluded based on having poor 3D image quality and 22 were excluded because of missing data on covariates. Based on the genotype data, relatives (n = 10) and genetic PCA outliers (n = 15) were identified and removed. The intersection of individuals with quality-controlled images, covariates and genotype data included 2,297 subjects.
Penn State Dataset
The data collected through Penn State consists of participants recruited in State College, PA; New York, NY; Urbana-Champaign, IL; Twinsburg, OH; Dublin, Ireland; Rome, Italy; Warsaw, Poland; and Porto, Portugal. The minimum age in this sample was 18, the maximum age 83 (mean age = 29; SD age = 14). Data on self-reported ancestry, body characteristics, age and sex as well as genotype data were obtained. Individuals were excluded if they reported a personal or family history of any birth defect or syndrome affecting the head or face, a personal history of any significant trauma, surgery, or any medical condition that might alter the structure of the face. From the entire Penn State dataset (n = 6,588), participants were excluded based on missing covariate data (n = 748) and quality control of the images resulted in the elimination of 52 participants. Based on the genotype data, European participants were selected (n = 1,614) and related individuals were removed (n = 59) (see Genotyping methods section), resulting in 1,555 participants for analysis. No genetic outliers were identified.
ALSPAC Dataset
The Avon Longitudinal Study of Parents and their Children (ALSPAC) is a longitudinal birth cohort in which pregnant women residing in Avon with an expected delivery date between 1 April 1991 and 31 December 1992 were recruited (Boyd et al., 2013; Fraser et al., 2013). At the time 14,541 pregnant women were recruited and DNA samples were collected for 11,343 children. Please note that the study website contains details of all the data that is available through a fully searchable data dictionary.1 Genome-wide data was available for 8,952 subjects of the B2261 study which is titled “Exploring distinctive facial features and their association with known candidate variants.” In addition to this, 4,731 3D-images were available along with corresponding data files containing information about sex, age, weight, height, ancestry, and other body characteristics. Participant ages ranged from 14 to 17 years old (mean age = 15; SD age = 0.5). Image quality control analysis resulted in the removal of 14 images of poor quality. 199 participants were removed due to self-reported non-European ancestry, 168 participants were removed because of missing covariate data and 726 individuals were removed because of relatedness. No genetic outliers were identified. The intersection of participants with quality-controlled images, covariates and genotype data included 3,566 individuals.
Genotyping
Genotype Quality Control and Population Structure
For both the Pittsburgh and Penn State sample, the genotype data were obtained as described in Claes et al. (2018). Pittsburgh participants were genotyped at the Center for Inherited Disease Research (CIDR) at Johns Hopkins University on the Illumina OmniExpress+ Exome v1.2 array plus 4,322 investigator-chosen SNPs included to capture variation in specific regions of interest based on previous studies of the genetics of facial variation (Shaffer et al., 2016). Genotypes were imputed to the 1000 Genomes Project Phase 3 reference panel (The 1000 Genomes Project Consortium et al., 2015), using SHAPEIT2 (Delaneau et al., 2013) for prephasing of haplotypes and IMPTUE2 for the imputation (Howie et al., 2009, 2011). Participants in the Penn State sample were either genotyped on the Illumina Human Hp200c1 BeadChip (IRB 32341) or the 23andMe v3 and v4 arrays (IRB 44929, 13103, 2503, 4320). In individuals with more than 500,000 variants, genotypes were prephased with SHAPEIT2 (Delaneau et al., 2013) and imputed to the 1000 Genomes Phase 3 reference using the Sanger Imputation Server2 with the Positional Burrows-Wheeler Transform (PBWT) pipeline (Durbin, 2014). For the ALSPAC sample, the participants were genotyped using the Illumina HumanHap550 quad genome-wide SNP genotyping platform by Sample Logistics and Genotyping Facilities at the Wellcome Trust Sanger Institute (Cambridge, UK) and the Laboratory Corporation of America (Burlington, NC, US), supported by 23andMe. Haplotypes were estimated using SHAPEIT2 (Delaneau et al., 2013) and imputed to the 1000 genomes reference panel (Phase 1, Version 3; Abecasis et al., 2012) using IMPUTE2 (Howie et al., 2009, 2011).
To select participants with primarily European ancestry in the Penn State Sample, an ADMIXTURE analysis was done with the 1000 Genomes Phase 3 dataset as the reference (Alexander et al., 2009). The estimated number of populations (k) was 5, determined by the cross-validation (CV) error for each k value. These results were then used to select samples with < 10% ancestry from all of the non-European admixture components. In the ASPAC sample, population stratification was assessed by multidimensional scaling analysis and compared with Hapmap II (release 22) European descent (CEU), Han Chinese, Japanese and Yoruba reference populations; all individuals with non-European ancestry were removed.
To assess population structure within the European subsets after removing non-European individuals, the same protocol was followed for the three datasets using PLINK 1.9 (Purcell et al., 2007). First, SNPs with a minor allele frequency (MAF) < 5% or more than 5% genotype data missing, were filtered out. Subsequently, SNPs were pruned for linkage disequilibrium (LD) with r2 set at 0.2 in a pairwise manner, with a moving window size of 50 variants shifting 5 variants each step. Subsequently, related individuals were identified and removed when the proportion of identity by descent (IBD) was higher than 0.125. Ancestry axes were determined with principal component analysis (PCA). Outliers (n = 15) were removed based on Z-scores calculated in the first 10 principal components. Z-scores higher than 6 indicated outliers, who were subsequently removed after which PCA was computed again.
Candidate Variant Selection
To select a set of NSCL/P candidate variants, 12 GWAS studies investigating the genetic etiology of NSCL/P were reviewed. Table 1 is a list of NSCL/P associated SNPs that have been selected based on a genome-wide significant association in at least one of the 12 GWAS studies. LD data for SNPs of the same locus were collected using the NCI NIH analysis tool LDmatrix in LDlink in the European populations (CEU, TSI, FIN, GBR, and IBS; Machiela and Chanock, 2015). With this data, SNP-pairs in perfect LD (r2 = 1 and D′ = 1) were detected and one SNP per pair was eliminated to avoid unnecessary multiple testing. The elimination of SNPs in LD resulted in a set of 75 lead-SNPs from 38 different loci (Table 1). Genotypes of 10 of these SNPs were absent in either of the three datasets (highlighted in gray). Therefore, 65 SNPs from 34 loci had been included in the analysis.
Phenotyping
Acquisition
Facial 3D surface images were acquired using two stereophotogrammetry systems and one laser scanning system. Facial surface data of the Pittsburgh sample were collected using the 3dMDface camera systems (3dMD, Atlanta, GA). For the Penn State sample, both 3dMDface camera systems and Vectra H1 (Canfield Scientific, Parsippany, NJ) camera systems were used. For the ALSPAC sample, a Konica Minolta Vivid 900 laser scanner (Konica Minolta Sensing Europe, Milton Keynes, UK) was used to take two high-resolution facial images which were subsequently processed, merged and registered using an algorithm implemented as a macro in Rapidform® software; INUS Technology Inc., Seoul, South Korea (Kau et al., 2004; Zhurov et al., 2005; Toma et al., 2008). Participants in all datasets were asked to have their mouth closed and to maintain a neutral facial expression during image capture.
Registration and Quality Control
3D surface images were imported into Matlab 2016b in.obj format to perform spatially dense registration (MeshMonk3). After importing the images, five positioning landmarks were indicated in the corners of the eye, the tip of the nose and the corners of the mouth to establish a crude alignment of the images. Subsequently, the images were cleaned by removing hair, ears, and any dissociated polygons. A symmetrical anthropometric mask (Claes et al., 2012) of 7,160 landmarks was then mapped onto the preprocessed images (Snyders et al., 2014). This resulted in homologous spatially dense configurations of quasi-landmarks.
After the registration, image quality control was performed to identify poorly remapped faces using two approaches. First, as described in Claes et al. (2018), outlier faces were identified by calculating Z-scores from the Mahalanobis distance between the mean face and each individual face. Faces with Z-scores higher than 2 were manually checked. Second, a score was calculated that reflects the missing data present in the image due to holes, spikes and other mesh artifacts that can be caused by facial hair or errors during the preprocessing steps, for example. Images with scores indicating a high amount of missing data, indicating large gaps in the mesh, were also manually checked. During the manual check, the images were either classified as images of poor quality or were preprocessed and mapped again.
The anthropometric mask (AM) is symmetrical relative to the sagittal plane, which allows reflected images to be created by changing the sign of the x-coordinate of the original mapped images. Both the original and the reflected remapped faces were then superimposed following a generalized Procrustes superimposition to eliminate differences in orientation, position and scale (Rohlf and Slice, 1990). Symmetrized images were created by averaging the original and the reflected images. All subsequent analyses were performed using these symmetrized images. Facial sizes were calculated by taking the centroid size of the spatially dense configurations.
Segmentation
To study the effects of NSCL/P candidate variants on facial shape in a detailed manner, a data-driven facial segmentation was performed (Claes et al., 2018). First, the images were corrected for the confounding factors age, age-squared, sex, weight, height, facial size and the first four ancestry axes using partial least-squares regression (PLSR; function plsregress from Matlab 2016b). Because of potential systematic differences in genotyping platforms or imputation efforts, the genomic ancestry axes were calculated on each dataset separately. Following this trend, the images were corrected in each dataset separately. After the correction, the segmentation was performed on all three datasets combined. Facial segments were defined by grouping vertices that are strongly correlated and connected using hierarchical spectral clustering, as described in Claes et al. (2018). The strength of covariation between quasi-landmarks was defined using Escoufier's RV coefficient (Escoufier, 1973), which is a scalar measure of strength of association between two groups of variables and is used in morphometric studies on biological shapes (Klingenberg, 2009). The RV coefficient allowed us to build a structural similarity matrix that defined the hierarchical construction of 63 facial segments, consisting of five levels (Figure 1). Subsequently, all segments independently were aligned using generalized Procrustes superimposition. To capture the major variance in the facial segments with fewer variables, a PCA was performed on each of the 63 segments in combination with parallel analysis (Hayton et al., 2004). Parallel analysis can be used to eliminate noisy or meaningless shape variations that result from sources of error, such as 3D image acquisition and/or quasi-landmark registration, as described in Claes et al. (2018). After combining the three datasets a large sample size (N = 7,418) is obtained which is beneficial to obtain a well-defined segmentation of landmark covariations and the determination of significant principal components (PCs) using parallel analysis. Moreover, facial shape data across all three datasets are now in the same shape space, enabling across-dataset analyses including a meta-analysis.
Figure 1. Global-to-local facial segmentation: 63 facial segments represented in blue. Starting from the global facial configuration in the center, the segment is split into two segments in level one, which are again split in two segments, which leads to four segments in level two. This is repeated up until level five which then contains 32 facial segments.
Meta-Analysis
Statistical Analysis
Each candidate SNP was tested against 63 facial segments each represented by multiple dimensions of variation (principal components). To this end we used canonical correlation analysis (CCA, canoncorr in Matlab 2016b) to test the effect of the 64 SNPs on facial shape under the reduced model. CCA is a multivariate testing framework which extracts the linear combination of PCs from a facial segment that has maximal correlation with the SNP being tested. The correlation is tested for significance based on Rao's F-test approximation (right tail, one-sided test; Olson, 1976). Using CCA, we tested each SNP individually under the additive genetic model in all three datasets. The reduced model was obtained after removing the linear effects of confounding factors age, age2, sex, weight, height, facial size and the first four ancestry axes using PLSR. Both the independent (SNP genotype) and the dependent (facial shape) variables were corrected for these covariates. Additionally, a partial least squares regression was performed (PLSR, plsregress in Matlab 2016b) in the same way as the CCA, to obtain the percentage of variance explained by the SNP in the facial segment that was tested.
Meta-Analysis
To maximize the statistical power from the combination of the three datasets we performed a meta-analysis in a round-robin fashion. Using CCA, each of the three datasets was in turn used to discover and define the linear combination of phenotypic variables that is maximally correlated with the SNP. This step of the meta-analysis will be referred to as the discovery. Subsequently, the other two datasets were projected onto the loadings obtained from CCA (the linear combination of phenotypic variables), creating univariate phenotypic variables which are then also tested for genotype-phenotype associations. This step of the meta-analysis will be referred to as the replication. The replication results in univariate phenotypic variables are measured in the replication datasets as a function of the discovered phenotypic trait in the CCA. In doing so, as proposed in Claes et al. (2018), the phenotypic trait discovered in CCA can be explicitly measured in the replication datasets enabling the combination of the statistical results across all datasets. The univariate phenotypic scores were statistically tested for association in a standard linear regression with the SNP genotypes as independent variables (function regstats in Matlab2016b). This function employs a t-statistic and a positive-sided p-value was obtained with the Student's T cumulative distribution function (function tcdf in Matlab2016b; Devroye, 1986).
After repeating the discovery plus two replications for each dataset in turn, we ended up with nine p-values (schematically represented in Figure 2). Each row contains a discovery p-value, determined by a multivariate CCA, and two replication p-values, each determined by a univariate linear regression. Row-wise each p-value was obtained from non-overlapping datasets and therefore independent and can be combined in a meta-analysis according to Stouffer et al.'s method 1949. This resulted in three meta-analyses p-values per segment, per SNP. Column-wise the p-values are not independent and therefore cannot be combined in the same manner.
Figure 2. Schematic representation of p-values used in the meta-analysis.
Multiple Testing Correction
Analyzing 65 SNPs in 63 separate facial segments introduces a multiple testing burden. The facial segments are overlapping and hierarchically constructed and are thus not completely independent. To determine the number of independent tests, we evaluated the eigenvalues of pairwise multivariate correlations between facial segments and of pairwise genotype correlations and determined 37 independent segments (Li and Ji, 2005). Some of the SNPs tested are in high LD, which resulted in 56 independent genetic variants. The round-robin meta-analysis resulted in three p-values, for which we also corrected. Therefore, a total Bonferroni correction for the effective numbers of independent segments, SNPs and round-robin sequences resulted in an adjusted significance threshold of 8.04 × 10−6 [i.e., 0.05/(37 × 56 × 3)].
Results
Meta-Analysis
The data-driven facial segmentation resulted in 63 facial segments that are hierarchically represented in Figure 1. At level two, four segments are identified, covering the nose, the mouth, the lower and the upper facial area, which coincides with the facial segmentation in Claes et al. (2018).
We identified six genetic loci involved in the etiology of NSCL/P that had significant effects on craniofacial morphology, namely rs742071 in 1p36, rs6740960, and rs7590268 in 2p21, rs34246903 in 4p16.2, rs10512248 in 9q22.32 and rs227731 in 17q22 (Table 2). These SNPs reached a p-value below the Bonferroni threshold of 8.04 × 10−6 in at least one facial segment in one of the three meta-analyses (Figure 3; Supplementary Figure 1). The percentage of variation explained by each SNP in a specific facial segment is reported in Table 2. Overall these percentages are low with the highest percentage of variation explained being 0.31% by rs6740960 in the chin area. Some variants, such as rs742071, rs6740960 and rs227731 show strong significance of association in the global segments as well as in the most local segments, with affected segments coming back in all three meta-analyses (Figure 3). Variants rs34246903 and rs10512248 both only show significant effects in the PITT and ALSPAC meta-analyses. Finally, rs7590268 is a variant that only shows effects in the local segments. Although there are no significant effects in the PENN STATE meta-analysis, the eyebrow region shows suggestive p-value signals.
Table 2. Discovery and meta-analysis results of significant SNPs in the most relevant segment.
Figure 3. Meta-analyses results plotting the –log10 p-value for the SNPs in each segment in each discovery dataset meta-analysis. A black circle surrounding the yellow circle shows significance (p-value < 8.04 × 10−6) in the segment encircled in the meta-analysis. From left to right, the discovery datasets are Pittsburgh, Penn State and ALSPAC.
The results of the meta-analyses of the remaining SNPs are provided in Supplementary Table 1. rs742071 in 1p36 is associated with an effect on the shape of the tip of the nose, with the angle between the nose and the lip being decreased in association with the risk allele (Figure 4). Two other SNPs (rs4920524 and rs9439713) showed very similar effects (Supplementary Table 1) and were in high LD with rs742071 (Supplementary Table 2). rs7590268 in 2p21, showed a significant effect in the supraorbital ridge, with a more protruded and pronounced shape associated with the minor allele. Another SNP in 2p21, rs6740960, showed an effect on the shape of the chin, where the major allele is associated with a more protruded but shorter in length chin. rs34246903 in 4p16.2 showed a significant effect in both the nose and philtrum region. With the major allele, the nasal ridge is narrowed, the width of the nose tip is decreased and the philtrum is protruded. rs10512248 in 9q22.32 showed an effect on the shape of the nose, where the angle between the nose and the lip decreases toward the major allele. The effect of rs227731 in 17q22 is located in the philtrum, with the minor allele associated with a more protruded philtrum. A short description of knowledge on the significant variants was provided in the Supplementary Text.
Figure 4. Visualization of the effects of the SNPs exaggerated in the direction of the minor allele (gray) and the major allele (orange), the risk allele is indicated with a *.
Discussion
This study combined the open-ended phenotyping approach proposed in Claes et al. (2018) with a candidate gene set-up to find possible effects of NSCL/P candidate genes on normal-range craniofacial morphology. To analyze the particular facial region in which the SNP is having the strongest evidence of effect, a global-to-local perspective was introduced by performing a data-driven segmentation (Claes et al., 2018). The segmentation in this study produced more compact segments in comparison to Claes et al. (2018), which is a consequence of the increase in sample size in this study. Another approach for a candidate gene association with facial morphology was described in Claes et al. (2014). Here, candidate variants were also tested on normal facial variation using spatially-dense configurations, but only on a global scale. Using this methodology three out of six SNP associations in this study would not have been discovered since they were only significant in local segments.
The replication strategy in Claes et al. (2018) was adapted toward a meta-analysis design. The meta-analysis was performed in a round-robin fashion, which is an elegant way to combine the power of three different datasets without having to select one as the sole discovery dataset. As expected, the behavior of the meta-analysis resulted in (1) a strengthening of the statistical evidence, if an effect was present in all three datasets, even if only suggestive evidence was found in each dataset separately and (2) a weakening of the statistical noise if a false positive was present in one of the datasets. The inclusion criterion for a significant signal in this study required a SNP to reach the Bonferroni adjusted threshold (8.04 × 10−6) in at least one segment for one of the three meta-analyses. To decrease the risk of identifying false positives, the Bonferroni threshold was determined by the number of independent segments, the number of independent genetic variants and, although not completely independent, the three meta-analyses.
In each meta-analysis a different dataset was used as the discovery dataset, this way our results indicate that the effects found were independent of the discovery dataset. Figure 3 shows that our results are influenced by sample size. In particular in the Penn State dataset, which has the lowest sample size (n = 1,555), rs7590268, rs34346903, and rs10512248 do not reach significance. Another possible explanation is that the true effects in the different datasets are not the same, for example due to small population differences between the datasets. This idea is analogous to between-study heterogeneity in genetic meta-analyses of quantitative traits (Magosi et al., 2017). In our study, we did not find heterogeneity in the effects found in the British (ALSPAC) and Euro-American (PITT and PENN) populations, this might be because we corrected for stratification within the European population, or because the effects found are in fact the same within the broader European and European-American population. To verify this, the directions of the SNP effects in shape space in each discovery analysis were visually inspected for qualitative coherence. Future analyses should include a quantitative analysis of these directions.
Our study is not the first to investigate the effects of genes involved in the etiology of NSCL/P on normal facial morphology. Two studies each analyzed the same set of 11 SNPs on different sets of phenotypic measurements (Boehringer et al., 2011; Liu et al., 2012). Boehringer et al. (2011) found a suggestive association between rs227731 (17q22) and nose width and the same variant was associated with normal facial morphology in Liu et al. (2012) (the associated phenotypic variables were not reported in this study). In our analysis, this SNP was associated with the shape of the philtrum (Figure 3). Liu et al. also found an association between rs7590268 (2p21) and facial morphology. This variant was associated with the shape of the brow ridge in our study (Figure 4). Both Boehringer et al. (2011) and Liu et al. (2012) found associations between 8q24 and craniofacial morphology. This locus is the strongest NSCL/P GWAS signal in European populations (Birnbaum et al., 2009; Beaty et al., 2010; Murray et al., 2012; Leslie et al., 2015). The association between 8q24 and facial morphology was not observed in this study. As for all variants tested that were not associated, true associations might not have been detected due to low power. Another reason for not finding the 8q24 association might be the difference in the methods applied. In both studies MRI-images and 2D photo images were used, in contrast to the use of 3D-surface images in our study. When correcting for multiple testing in Boehringer et al. (2011), the association between 8q24 and bizygomatic distance was no longer significant, which was discussed in the manuscript. Liu et al. corrected for the 48 phenotypes, but an additional correction for the number of SNPs tested was not considered in determining the threshold for statistical significance. Another difference is that both studies discussed used a limited number of landmarks from which distance measurements were derived. These distances cannot fully capture the complexity of human facial shape. In contrast, our use of spatially-dense quasi-landmarks strongly improves the description of facial morphology (Claes et al., 2012). In the case of Liu et al. (2012), principal components (PCs) were also derived from the sparse landmark configurations. These PCs represent a more complex configuration of facial shape than distances, yet they were examined separately instead of in a multivariate fashion. This preselection of single measurements causes information on the combination of measurements to be lost. In our study, a multivariate framework was applied in which the linear combination of PCs from a facial segment that are in maximal correlation with the SNP being tested are extracted. This way, the PCs can be combined without introducing another multiple testing problem.
Another study analyzing the genetic overlap between NSCL/P genetics and the genetics of normal craniofacial variation was conducted by Howe et al. (2018). Polygenic risk scores (PRS) were calculated from NSCL/P data and subsequently used as the variable of interest in an association with normal facial morphology, in which seven distances based on sparse landmarks were used as phenotypic features. The finding in this study was an association between an increased PRS and a decreased philtrum width. Although the use of PRS is distinct from a candidate gene approach, this study presents additional evidence to the hypothesis that genes involved in the development of NSCL/P have an effect on normal facial morphology as well. More specifically, in both Howe et al. and our study, variations in the philtrum were identified to be associated with NSCL/P susceptibility genes.
Additional evidence for the genetic overlap between cleft genetics and craniofacial morphology genetics lies in the hypothesis of facial endophenotypes for NSCL/P (Weinberg et al., 2008, 2009; Roosenboom et al., 2015). The NSCL/P facial endophenotype shows a retrusion of the midfacial region and a broadening of the upper facial area, which can be an expression of NSCL/P susceptibility genes in unaffected relatives but is also part of normal facial variation. In Weinberg et al. (2009), it was found that unaffected relatives show a decreased philtrum width in comparison to the controls. This result coincides with the association found between an increased cleft PRS and a decreased philtrum width in Howe et al. In our study, the philtrum is more protruded toward the risk allele in rs227731. Although distances are not easy to compare directly with a spatially-dense representation, when looking at the results of rs227731, one can imagine that philtrum protrusion leads to a smaller philtrum width and vice versa.
There is some evidence for an effect of NSCL/P susceptibility genes on the shape of the nose. Boehringer et al. (2011) and Liu et al. (2012) found one NSCL/P associated genetic variant (rs1258763 in 15q13) to be associated with nose width. This result was not replicated in this study. On the other hand, NSCL/P endophenotypic features have been found in the nose by Weinberg et al. (2008, 2009). In our study, most significant effects were found in the central facial area: three in the nose and one in the philtrum. This is not surprising for two reasons. First, when looking at human embryological development, in which the development of the nose and palate are tightly linked. The nose bridge and philtrum are formed by the fusion of the medial-nasal processes, which then fuse with the lateral nasal processes and the maxillary swellings to form the sides of the nose, alae, and the maxillae. The anterior palate is also formed with the fusion of the maxillary and medial-nasal process. The fusion of these prominences requires the coordinated growth of the oronasal prominences in a precise temporal-spatial sequence (Dixon et al., 2011). Thus, the genetic regulation mechanisms involved in lip and palate fusion likely also affects nasal and philtrum morphology. Second, the central facial features show a high genetic determination in genome wide studies (Claes et al., 2018).
The variants in this study were annotated to genes that have been shown to play a role in craniofacial development in previous studies (Table 2, Supplementary Text 1). PAX7, of which rs742071 is an intron variant, has been shown to play a role in neural crest development (Basch et al., 2006). Neural crest progenitors give rise to craniofacial cartilage, which interestingly is an important structural component in the tip of the nose, with which this variant was associated (Mansouri et al., 1996; Zalc et al., 2015). Another SNP, rs10512248, is an intron variant of PTCH1, which is an important factor in the Hedgehog signaling pathway (Aoto and Trainor, 2015). This pathway plays a fundamental role in craniofacial development in vertebrates (Xavier et al., 2016; Everson et al., 2018). Other variants found in this study were annotated to genes, such as MSX1 and NOG, which have been shown to play a role in craniofacial development (Satokata and Maas, 1994; Anderson et al., 2006). Many of these studies look at gene expression studies during embryonic development in animal models. Although these studies can complement ones like the current study, caution should be taken when comparing animal model studies with our study, since we are investigating different species at different stages (embryological vs. post-natal).
In this study we focused our analyses on participants with European ancestry. To analyze possible population similarities and differences, future analyses should include other populations. Furthermore, only lead SNPs from our literature survey were used, not all loci involved in cleft development have been discovered and this study only included variants discovered in NSCL/P GWAS studies. Thus, the candidate variant list used in this study is incomplete.
For future analyses, it would be interesting to calculate polygenic risk scores (PRS) for NSCL/P, and associate these against normal facial variation in an open-ended phenotyping approach. The PRS capture an individual's genetic propensity toward a trait and thus combines the effect sizes discovered in a GWAS into a univariate score capturing a risk for NSCL/P. Combining the information of multiple markers in one score analogous to this study would result in a facial phenotype associated with NSCL/P risk. Additionally, the risk-score could increase power in comparison to the use of individual SNPs. Another approach would be to calculate a quantitative measurement that indicates to what degree the endophenotype is present in each participant and to use this in a genetic mapping effort, using candidate variants, genome-wide SNPs or PRS. The methodology described in this study allows us to study the effects of genes with a known role in a certain condition on facial morphology. In the future, genes responsible for syndromes with a distinct facial phenotype can also be tested in this framework to see if variations in these genes are responsible for morphological changes in the face.
Conclusion
Many genetic syndromes and malformations are characterized by a distinctive facial phenotype (Jones et al., 2013). The underlying genes might also be involved in normal craniofacial development (Winter, 1996). Since NSCL/P is one of the most frequent congenital craniofacial malformations, its genetic background could also be involved in normal facial morphology. In this study, the open-ended data-driven phenotyping approach from Claes et al.'s (2018) GWAS on normal-range variation in facial shape was used to study the effects of NSCL/P candidate genes on normal-range craniofacial morphology in a global-to-local perspective. We identified six SNPs involved in NSCL/P with effects on the shape of the nose, chin and philtrum area in a non-clinical population. This study gives insight into the complex genetic architecture of normal-range craniofacial morphology. Furthermore, it provides evidence for the interplay between the genetic background of NSCL/P and craniofacial morphology.
Data Availability Statement
The datasets analyzed for this study are available through various sources. For the Pittsburgh dataset, the genotypic markers are available through the dbGaP controlled-access repository (http://www.ncbi.nlm.nih.gov/gap) at accession phs000949.v1.p1. The 3D facial images are available through the FaceBase Consortium (https://www.facebase.org/data/record/#1/isa:dataset/RID=14283) The participants making up the Penn State University dataset 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 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 M.D.S. The ALSPAC data will be made available to bona fide researchers on application to the ASLAPC Executive Committee.
Ethics Statement
Institutional Review Board (IRB) approval was obtained at each recruitment site, and all participants gave their written informed consent before participation; for children, written consent was obtained from a parent or legal guardian. For the Pittsburgh sample, the following local ethics approvals were obtained: University of Pittsburgh IRB PRO09060553 and RB0405013; UT Health Committee for the Protection of Human Subjects HSC-DB-09-0508; Seattle Children's IRB 12107; University of Iowa Human Subjects Office/IRB 200912764 and 200710721. For the Penn State sample, the following local ethics approvals were obtained: State College, PA (IRB 44929 and 4320); New York, NY (IRB 45727); Urbana-Champaign, IL (IRB 13103); Dublin, Ireland; Rome, Italy; Warsaw, Poland; and Porto, Portugal (IRB 32341); Twinsburg, OH (IRB 2503). For the ALSPAC sample, ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees.
Author Contributions
KI performed all analyses and wrote the first draft of the manuscript under supervision of PC, GH, and HP. PC, JR, EF, JS, SW, MM, JW, and MS conceptualized the design of the study. JR, EF, JS, SW, and MM organized the PITT cohort. JW and MS organized the PSU cohort and imputed the PSU genetic data. SR coordinated the collection of the ALSPAC images. JR, HH, and JW provided input throughout the analyses and the writing process. All authors contributed to manuscript revision, read and approved the submitted version.
Funding
This study was funded by the National Institutes of Health (1-RO1-DE027023), the Research Fund KU Leuven (BOF-C1, C14/15/081) and the Research Program of the The Research Foundation–Flanders (Belgium) (FWO, G078518N). HP is a Senior Clinical Investigator of The Research Foundation–Flanders (Belgium) (FWO). The collaborators at the University of Pittsburgh were supported by the National Institute for Dental and Craniofacial Research (see URLs) through the following grants: U01-DE020078, U01-DE020057, R01-DE016148, K99-DE02560, and 1-R01-DE027023. Funding for genotyping was provided by the National Human Genome Research Institute (see URLs): X01-HG007821 and X01-HG007485. Funding for initial genomic data cleaning by the University of Washington was provided by contract HHSN268201200008I from the National Institute for Dental and Craniofacial Research (http://www.nidcr.nih.gov/) awarded to the Center for Inherited Disease Research (CIDR). The collaborators at Penn State University were supported in part by grants from the Center for Human Evolution and Development at Penn State, the Science Foundation of Ireland Walton Fellowship (04.W4/B643), the US National Institute of Justice (see URLs; 2008-DN-BX-K125) and the US Department of Defense (see URLs). The UK Medical Research Council and Wellcome (Grant ref: 102215/2/13/2) and the University of Bristol provide core support for ALSPAC. This publication is the work of the authors and KI and PC will serve as guarantors for the contents of this paper. A comprehensive list of grants funding is available on the ALSPAC website (http://www.bristol.ac.uk/alspac/external/documents/grant-acknowledgements.pdf) This research was specifically funded by the University of Cardiff. For the ALSPAC sample, GWAS data was generated by Sample Logistics and Genotyping Facilities at Wellcome Sanger Institute and LabCorp (Laboratory Corporation of America) using support from 23andMe.
Conflict of Interest Statement
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.
Acknowledgments
We are extremely grateful to all the families who took part in this study, the midwives for their help in recruiting them, and the whole Pittsburgh, Penn State and ALSPAC team, which includes interviewers, computer and laboratory technicians, clerical workers, research scientists, volunteers, managers, receptionists and nurses.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2018.00502/full#supplementary-material
Footnotes
1. ^http://www.bris.ac.uk/alspac/researchers/data-access/data-dictionary/
2. ^https://imputation.sanger.ac.uk
3. ^Facial mapping software available on https://github.com/TheWebMonks/meshmonk
References
Abecasis, G. R., Auton, A., Brooks, L. D., DePristo, M. A., Durbin, R. M., Handsaker, R. E., et al. (2012). An integrated map of genetic variation from 1,092 human genomes. Nature 491, 56–65. doi: 10.1038/nature11632
Adhikari, K., Fuentes-Guajardo, M., Quinto-Sánchez, M., Mendoza-Revilla, J., Chacón-Duque, J. C., Acuña-Alonzo, V., et al. (2016). A genome-wide association scan implicates DCHS2, RUNX2, GLI3, PAX1 and EDAR in human facial variation. Nat. Commun. 7:11616. doi: 10.1038/ncomms11616
Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. doi: 10.1101/gr.094052.109
Anderson, R. M., Stottmann, R. W., Choi, M., and Klingensmith, J. (2006). Endogenous bone morphogenetic protein antagonists regulate mammalian neural crest generation and survival. Dev. Dyn. 235, 2507–2520. doi: 10.1002/dvdy.20891
Aoto, K., and Trainor, P. A. (2015). Co-Ordinated brain and craniofacial development depend upon patched1/XIAP regulation of cell survival. Hum. Mol. Genet. 24, 698–713. doi: 10.1093/hmg/ddu489
Basch, M. L., Bronner-Fraser, M., and García-Castro, M. I. (2006). Specification of the neural crest occurs during gastrulation and requires Pax7. Nature 441, 218–222. doi: 10.1038/nature04684
Beaty, T. H., Murray, J. C., Marazita, M. L., Munger, R. G., Ruczinski, I., Hetmanski, J. B., Liang, K. Y., et al. (2010). A genome-wide association study of cleft lip with and without cleft palate identifies risk variants near MAFB and ABCA4. Nat. Genet. 42, 525–529. doi: 10.1038/ng.580
Birnbaum, S., Ludwig, K. U., Reutter, H., Herms, S., Steffens, M., Rubini, M., et al. (2009). Key susceptibility locus for nonsyndromic cleft lip with or without cleft palate on chromosome 8q24. Nat. Genet. 41, 473–477. doi: 10.1038/ng.333
Boehringer, S., van der Lijn, F., Liu, F., Günther, M., Sinigerova, S., Nowak, S., et al. (2011). Genetic determination of human facial morphology: links between cleft-lips and normal variation. Eur.J. Hum. Genet. 19, 1192–1197. doi: 10.1038/ejhg.2011.110
Boyd, A., Golding, J., Macleod, J., Lawlor, D. A., Fraser, A., Henderson, J., et al. (2013). Cohort profile: the ‘children of the 90s’–the index offspring of the avon longitudinal study of parents and children. Int. J. Epidemiol. 42, 111–127. doi: 10.1093/ije/dys064
Claes, P., Liberton, D. K., Daniels, K., Rosana, K. M., Quillen, E. E., Pearson, L. N., et al. (2014). Modeling 3D facial shape from DNA. PLoS Genet. 10:e1004224. doi: 10.1371/journal.pgen.1004224
Claes, P., Roosenboom, J., White, J. D., Swigut, T., Sero, D., Li, J., et al. (2018). Genome-wide mapping of global-to-local genetic effects on human facial shape. Nat. Genet. 50, 414–423. doi: 10.1038/s41588-018-0057-4
Claes, P., Walters, M., and Clement, J. (2012). Improved facial outcome assessment using a 3D anthropometric mask. Int. J. Oral Maxillofac. Surg. 41, 324–330. doi: 10.1016/j.ijom.2011.10.019
Delaneau, O., Howie, B., Cox, A. J., Zagury, J.-F., and Marchini, J. (2013). Haplotype estimation using sequencing reads. Am. J. Hum. Genet. 93, 687–696. doi: 10.1016/j.ajhg.2013.09.002
Dixon, M. J., Marazita, M. L., Beaty, T. H., and Murray, J. C. (2011). Cleft Lip and Palate: understanding genetic and environmental influences. Nat. Rev. Genet. 12, 167–178. doi: 10.1038/nrg2933
Durbin, R. (2014). Efficient haplotype matching and storage using the positional burrows-wheeler transform (PBWT). Bioinformatics 30, 1266–1272. doi: 10.1093/bioinformatics/btu014
Escoufier, Y. (1973). Le traitement des variables vectorielles. Biometrics 29, 751–760. doi: 10.2307/2529140
Everson, J. L., Fink, D. M., Chung, H. M., Sun, M. R., and Lipinski, R. J. (2018). Identification of sonic hedgehog-regulated genes and biological processes in the cranial neural crest mesenchyme by comparative transcriptomics. BMC Genomics 19:497. doi: 10.1186/s12864-018-4885-5
Fraser, A., Macdonald-Wallis, C., Tilling, K., Boyd, A., Golding, J., Davey Smith, G., et al. (2013). Cohort Profile: the avon longitudinal study of parents and children: ALSPAC mothers cohort. Int. J. Epidemiol. 42, 97–110. doi: 10.1093/ije/dys066
Grant, S. F. A., Wang, K., Zhang, H., Glaberson, W., Annaiah, K., Kim, C. E., et al. (2009). A genome-wide association study identifies a locus for nonsyndromic cleft lip with or without cleft palate on 8q24. J. Pediatr. 155, 909–913. doi: 10.1016/j.jpeds.2009.06.020
Hayton, J. C., Allen, D. G., and Scarpello, V. (2004). Factor retention decisions in exploratory factor analysis: a tutorial on parallel analysis. Organizat. Res. Methods 7, 191–205. doi: 10.1177/1094428104263675
Howe, L. J., Lee, M. K., Sharp, G. C., Smith, G. D., St. Pourcain, B., Shaffer, J. R., et al. (2018). Investigating the shared genetics of non-syndromic cleft lip/palate and facial morphology. PLoS Genet. 14:e1007501. doi: 10.1371/journal.pgen.1007501
Howie, B., Marchini, J., and Stephens, M. (2011). Genotype Imputation with Thousands Of. G3 1, 457–470. doi: 10.1534/g3.111.001198
Howie, B. N., Donnelly, P., and Marchini, J. (2009). A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. 5:e1000529. doi: 10.1371/journal.pgen.1000529
Jones, K. L., Jones, M. C., and del Campo, M. (2013). Smith's Recognizable Patterns of Human Malformation E-Book. Philadelphia, PA: Elsevier Saunders.
Kau, C. H., Zhurov, A., Scheer, R., Bouwman, S., and Richmond, S. (2004). The feasibility of measuring three-dimensional facial morphology in children. Orthodont. Craniofac. Res. 7, 198–204. doi: 10.1111/j.1601-6343.2004.00289.x
Klingenberg, C. P. (2009). morphometric integration and modularity in configurations of landmarks: tools for evaluating a priori hypotheses. Evol. Dev. 11, 405–421. doi: 10.1111/j.1525-142X.2009.00347.x
Leslie, E. J., Carlson, J. C., Shaffer, J. R., Butali, A., Buxó, C. J., Castilla, E. E., et al. (2017). Genome-wide meta-analyses of nonsyndromic orofacial clefts identify novel associations between FOXE1 and all orofacial clefts, and TP63 and cleft lip with or without Cleft Palate. Hum. Genet. 136, 275–86. doi: 10.1007/s00439-016-1754-7
Leslie, E. J., Liu, H., Carlson, J. C., Shaffer, J. R., Feingold, E., Wehby, G., et al. (2016). A Genome-wide association study of nonsyndromic cleft palate identifies an etiologic missense variant in GRHL3. Am. J. Hum. Genet. 98, 744–754. doi: 10.1016/j.ajhg.2016.02.014
Leslie, E. J., Taub, M. A., Liu, H., Steinberg, K. M., Koboldt, D. C., Zhang, Q., et al. (2015). Identification of functional variants for cleft lip with or without cleft palate in or near PAX7, FGFR2, and NOG by targeted sequencing of GWAS loci. Am. J. Hum. Genet. 96, 397–411. doi: 10.1016/j.ajhg.2015.01.004
Li, J., and Ji, L. (2005). Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity 95, 221–227. doi: 10.1038/sj.hdy.6800717
Liu, F., van der Lijn, F., Schurmann, C., Zhu, G., Chakravarty, M. M., Hysi, P. G., et al. (2012). A Genome-wide association study identifies five loci influencing facial morphology in europeans. PLoS Genet. 8:e1002932. doi: 10.1371/journal.pgen.1002932
Ludwig, K. U., Ahmed, S. T., Böhmer, A. C., Sangani, N. B., Varghese, S., Klamt, J., et al. (2016). Meta-analysis reveals genome-wide significance at 15q13 for nonsyndromic clefting of both the lip and the palate, and functional analyses implicate GREM1 as a plausible causative gene. PLoS Genet. 12:e1005914. doi: 10.1371/journal.pgen.1005914
Ludwig, K. U., Böhmer, A. C., Bowes, J., Nikolić, M., Ishorst, N., Wyatt, N., et al. (2017). Imputation of orofacial clefting data identifies novel risk loci and sheds light on the genetic background of cleft lip ± cleft palate and cleft palate only. Hum. Mol. Genet. 26, 829–842. doi: 10.1093/hmg/ddx012
Ludwig, K. U., Mangold, E., Herms, S., Nowak, S., Reutter, H., Paul, A., et al. (2012). Genome-wide meta-analyses of nonsyndromic cleft lip with or without cleft palate identify six new risk loci. Nat. Genet. 44, 968–971. doi: 10.1038/ng.2360
Machiela, M. J., and Chanock, S. J. (2015). LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. Bioinformatics 31, 3555–3557. doi: 10.1093/bioinformatics/btv402
Magosi, L. E., Goel, A., Hopewell, J. C., Farrall, M. and on behalf of the CARDIoGRAMplusC4D Consortium (2017). Identifying systematic heterogeneity patterns in genetic association meta-analysis studies. PLoS Genet. 13:e1006755. doi: 10.1371/journal.pgen.1006755
Mangold, E., Ludwig, K. U., Birnbaum, S., Baluardo, C., Ferrian, M., Herms, S., et al. (2010). Genome-wide association study identifies two susceptibility loci for nonsyndromic cleft lip with or without cleft palate. Nat. Genet. 42, 24–26. doi: 10.1038/ng.506
Mansouri, A., Stoykova, A., Torres, M., and Gruss, P. (1996). Dysgenesis of cephalic neural crest derivatives in pax7−/− mutant mice. Development 122, 831–838.
Mossey, P. A., Little, J., Munger, R. G., Dixon, M. J., and Shaw, W. C. (2009). Cleft lip and palate. Lancet 374, 1773–1785. doi: 10.1016/S0140-6736(09)60695-4
Mostowska, A., Gaczkowska, A., Zukowski, K., Ludwig, K. U., Hozyasz, K. K., Wójcicki, P., et al. (2018). Common variants in DLG1 locus are associated with non-syndromic cleft lip with or without cleft palate. Clin. Genet. 93, 784–793. doi: 10.1111/cge.13141
Murray, T., Taub, M. A., Ruczinski, I., Scott, A. F., Hetmanski, J. B., Schwender, H., et al. (2012). Examining markers in 8q24 to explain differences in evidence for association with cleft lip with/without cleft palate between Asians and Europeans. Genet. Epidemiol. 36, 392–399. doi: 10.1002/gepi.21633
Olson, C. L. (1976). On choosing a test statistic in multivariate analysis of variance. Psychol. Bull. 83, 579–586. doi: 10.1037/0033-2909.83.4.579
Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795
Rohlf, F. J., and Slice, D. (1990). Extensions of the procrustes method for the optimal superimposition of landmarks. System. Zool. 39:40. doi: 10.2307/2992207
Roosenboom, J., Hens, G., Mattern, B. C., Shriver, M. D., and Claes, P. (2016). Exploring the underlying genetics of craniofacial morphology through various sources of knowledge. Res. Article 2016:3054578. doi: 10.1155/2016/3054578
Roosenboom, J., Saey, I., Peeters, H., Devriendt, K., Claes, P., and Hens, G. (2015). Facial characteristics and olfactory dysfunction: two endophenotypes related to nonsyndromic cleft lip and/or palate. BioMed. Res. Int. 2015:863429. doi: 10.1155/2015/863429
Satokata, I., and Maas, R. (1994). Msx1 deficient mice exhibit cleft palate and abnormalities of craniofacial and tooth development. Nat. Genet. 6, 348–356. doi: 10.1038/ng0494-348
Shaffer, J. R., Orlova, E., Lee, M. K., Leslie, E. J., Raffensperger, Z. D., Heike, C. L., et al. (2016). Genome-wide association study reveals multiple loci influencing normal human facial morphology. PLoS Genet. 12:e1006149. doi: 10.1371/journal.pgen.1006149
Snyders, J., Claes, P., Vandermeulen, D., and Suetens, P. (2014). Development and comparison of non-rigid surface registration algorithms and extensions. Technical report KUL/ESAT/PSI/1401, 1–55 (2014).
Stouffer, S. A., Suchman, E. A., Devinney, L. C., Star, S. A., and Williams, R. M. Jr. (1949). “The American Soldier: Adjustment during Army Life,” in Studies in Social Psychology in World War II, Vol. 1 (Oxford: Princeton Univ. Press).
Sun, Y., Huang, Y., Yin, A., Pan, Y., Wang, Y., Wang, C., et al. (2015). Genome-wide association study identifies a new susceptibility locus for cleft lip with or without a cleft palate. Nat. Commun. 6:6414. doi: 10.1038/ncomms7414
The 1000 Genomes Project Consortium, Richard Gibbs, A., Boerwinkle, E., Doddapaneni, H., Han, Y., Korchina, V., Kovar, C., et al. (2015). A global reference for human genetic variation. Nature 526, 68–74. doi: 10.1038/nature15393
Toma, A. M., Zhurov, A., Playle, R., and Richmond, S. (2008). A three-dimensional look for facial differences between males and females in a british-caucasian sample aged 151/2 years old. Orthodont. Craniofacial Res. 11, 180–185. doi: 10.1111/j.1601-6343.2008.00428.x
Weinberg, S. M., Naidoo, S. D., Bardi, K. M., Brandon, C. A., Neiswanger, K., Resick, J. M., et al. (2009). Face shape of unaffected parents with cleft affected offspring: combining three-dimensional surface imaging and geometric morphometrics. Orthodont. Craniofacial Res. 12, 271–281. doi: 10.1111/j.1601-6343.2009.01462.x
Weinberg, S. M., Neiswanger, K., Richtsmeier, J. T., Maher, B. S., Mooney, M. P., Siegel, M. I., et al. (2008). Three-dimensional morphometric analysis of craniofacial shape in the unaffected relatives of individuals with nonsyndromic orofacial clefts: a possible marker for genetic susceptibility. Am. J. Med. Genet. Part A. 146A, 409–420. doi: 10.1002/ajmg.a.32177
Weinberg, S. M., Raffensperger, Z. D., Kesterke, M. J., Heike, C. L., Cunningham, M. L., Hecht, J. T., et al. (2016). The 3D facial norms database: part 1. A web-based craniofacial anthropometric and image repository for the Clinical and Research Community. Cleft Palate-Craniofac. J. 53, e185–e197. doi: 10.1597/15-199
Xavier, G. M., Seppala, M., Barrell, W., Birjandi, A. A., Geoghegan, F., and Cobourne, M. T. (2016). Hedgehog receptor function during craniofacial development. Dev. Biol. 415, 198–215. doi: 10.1016/j.ydbio.2016.02.009
Yu, Y., Zuo, X., He, M., Gao, J., Fu, Y., Qin, C., et al. (2017). Genome-wide analyses of non-syndromic cleft lip with palate identify 14 novel loci and genetic heterogeneity. Nat. Commun. 8:14364. doi: 10.1038/ncomms14364
Keywords: craniofacial, morphology, NSCL/P, candidate gene, ALSPAC
Citation: Indencleef K, Roosenboom J, Hoskens H, White JD, Shriver MD, Richmond S, Peeters H, Feingold E, Marazita ML, Shaffer JR, Weinberg SM, Hens G and Claes P (2018) Six NSCL/P Loci Show Associations With Normal-Range Craniofacial Variation. Front. Genet. 9:502. doi: 10.3389/fgene.2018.00502
Received: 21 June 2018; Accepted: 05 October 2018;
Published: 25 October 2018.
Edited by:
Dana C. Crawford, Case Western Reserve University, United StatesReviewed by:
Kerstin U. Ludwig, Universitätsklinikum Bonn, GermanyJennifer L. Fish, University of Massachusetts Lowell, United States
Copyright © 2018 Indencleef, Roosenboom, Hoskens, White, Shriver, Richmond, Peeters, Feingold, Marazita, Shaffer, Weinberg, Hens and Claes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Karlijne Indencleef, karlijne.indencleef@kuleuven.be
Peter Claes, peter.claes@kuleuven.be